20 #if !defined(__CINT__) || defined(__MAKECINT__)
26 #include <Riostream.h>
28 #include <TGraphErrors.h>
31 #include <TFileMerger.h>
32 #include <TMultiGraph.h>
36 #include <TGridCollection.h>
37 #include <TGridResult.h>
38 #include <TClonesArray.h>
39 #include <TObjString.h>
53 void Draw(Int_t cell[], Int_t iBC, Int_t nBC,
const Int_t cellref=151)
55 gROOT ->SetStyle(
"Plain");
56 gStyle->SetOptStat(0);
57 gStyle->SetFillColor(kWhite);
58 gStyle->SetTitleFillColor(kWhite);
59 gStyle->SetPalette(1);
61 char out[120];
char title[100];
char name[100];
62 TString slide =
"Cells ";
65 slide+=cell[iBC+nBC-1];
66 sprintf(out,
"%d-%d.gif",cell[iBC],cell[iBC+nBC-1]);
68 TH2 *hCellAmplitude = (TH2*) gFile->Get(
"hCellAmplitude");
69 TH1 *hCellref = hCellAmplitude->ProjectionY(
"badcells",cellref+1,cellref+1);
72 TCanvas *c1 =
new TCanvas(
"badcells",
"badcells",1000,750) ;
74 if (nBC > 6) c1->Divide(3,3) ;
75 else if (nBC > 3) c1->Divide(3,2) ;
79 TLegend *leg =
new TLegend(0.7, 0.7, 0.9, 0.9);
80 for(i=0; i<nBC ; i++){
81 sprintf(name,
"Cell %d",cell[iBC+i]) ;
82 TH1 *hCell = hCellAmplitude->ProjectionY(name,cell[iBC+i]+1,cell[iBC+i]+1);
84 c1->cd(i%9 + 1)->SetLogy();
85 sprintf(title,
"Cell %d Entries : %d",cell[iBC+i], (Int_t)hCell->GetEntries()) ;
86 hCell->SetLineColor(2) ;
88 hCell->SetMaximum(1e5);
90 hCell->SetAxisRange(0.,4.);
91 hCell->GetXaxis()->SetTitle(
"E (GeV)");
92 hCell->GetYaxis()->SetTitle(
"N Entries");
93 hCellref->SetAxisRange(0.,4.);
94 hCell->SetLineWidth(1) ;
95 hCellref->SetLineWidth(1) ;
96 hCell->SetTitle(title);
97 hCellref->SetLineColor(1) ;
99 leg->AddEntry(hCellref,
"reference",
"l");
100 leg->AddEntry(hCell,
"current",
"l");;
103 hCellref->Draw(
"same") ;
108 TString PdfName =
"BadCellsCandidate.pdf";
111 c1->Print(PdfName.Data());}
114 c1->Print(PdfName.Data());}
115 else c1->Print(PdfName.Data());
132 void Convert(TString fCalorimeter =
"EMCAL", TString period =
"LHC11h", TString pass =
"pass1_HLT", TString trigger=
"default")
134 TH2F * hCellAmplitude =
new TH2F(
"hCellAmplitude",
"Cell Amplitude",11520,0,11520,200,0,10);
135 TH1D * hNEventsProcessedPerRun =
new TH1D(
"hNEventsProcessedPerRun",
"Number of processed events vs run number",200000,100000,300000);
137 TString
file =
"/scratch/alicehp2/mas/analyse/QA/"+period+
"/"+ pass +
"/runlistMB.txt" ;
139 pFile = fopen(file.Data(),
"r");
142 cout <<
" fcalo: " << fCalorimeter <<
"; period: " << period <<
"; pass: " << pass <<
" trigger "<<trigger<< endl;
156 ncols = fscanf(pFile,
"%d %d ",&p,&q);
166 const Int_t nRun = nlines ;
170 TString direct =
"CaloQA_";
173 for(Int_t i = 0 ; i < nRun ; i++) {
174 base =
"/scratch/alicehp2/germain/QA/";
175 BCfile = base + period ;
184 infile = base +
".root" ;
185 TFile *f = TFile::Open(infile);
188 TDirectoryFile *dir = (TDirectoryFile *)f->Get(direct);
189 TList *outputList = (TList*)dir->Get(direct);
192 hAmpId =(TH2F *)outputList->FindObject(
"EMCAL_hAmpId");
193 hNEvents =(TH2F *)outputList->FindObject(
"hNEvents");
194 Nentr = (Int_t)hNEvents->GetEntries();
195 if (Nentr<100) continue ;
196 hNEventsProcessedPerRun->SetBinContent(RunId[i]-100000,(Double_t)Nentr);
198 for(ix=1;ix<=200;ix++){
199 for(iy=1; iy<=11520; iy++){
201 content += hAmpId->GetBinContent(ix,iy);
202 content += hCellAmplitude->GetBinContent(iy,ix);
204 hCellAmplitude->SetBinContent(iy,ix,content);
209 if(i==0){ cout<<
"Merging/Converting procedure ..." ; cout.flush();}
210 else { cout<<
"..." ; cout.flush();}
211 outputList->Delete();
218 TFile *BCF = TFile::Open(BCfile,
"recreate");
219 hNEventsProcessedPerRun->Write();
220 hCellAmplitude->Write();
222 cout<<
"DONE !"<<endl;
237 void Process(Int_t *pflag[11520][7], TH1* inhisto, Double_t Nsigma = 4.,
238 Int_t dnbins = 200, Double_t dmaxval = -1.)
240 Int_t crit = *pflag[0][0] ;
241 Double_t goodmax= 0. ;
242 Double_t goodmin= 0. ;
247 dmaxval = inhisto->GetMaximum()*1.01;
248 if(crit==2 && dmaxval > 1) dmaxval =1. ;
251 TH1 *distrib =
new TH1F(Form(
"%sDistr",inhisto->GetName()),
"", dnbins, inhisto->GetMinimum(), dmaxval);
252 distrib->SetXTitle(inhisto->GetYaxis()->GetTitle());
253 distrib->SetYTitle(
"Entries");
256 for (Int_t c = 1; c <= inhisto->GetNbinsX(); c++)
257 distrib->Fill(inhisto->GetBinContent(c));
260 TCanvas *c1 =
new TCanvas(inhisto->GetName(),inhisto->GetName(), 800,400);
264 gPad->SetLeftMargin(0.14);
265 gPad->SetRightMargin(0.06);
267 inhisto->SetTitleOffset(1.7,
"Y");
271 gPad->SetLeftMargin(0.14);
272 gPad->SetRightMargin(0.10);
277 for (i = 2; i <= distrib->GetNbinsX(); i++){
278 if(distrib->GetBinContent(higherbin) < distrib->GetBinContent(i))
282 for(i=higherbin ; i<=dnbins ; i++)
283 if(distrib->GetBinContent(i)<2)
break ;
284 goodmax = distrib->GetBinCenter(i);
286 for(i=higherbin ; i>0 ; i--)
287 if(distrib->GetBinContent(i)<2)
break ;
288 goodmin = distrib->GetBinLowEdge(i);
290 TF1 *fit2 =
new TF1(
"fit2",
"gaus");
292 distrib->Fit(fit2,
"0LQEM",
"", goodmin, goodmax);
293 Double_t sig, mean, chi2ndf;
294 mean = fit2->GetParameter(1);
295 sig = fit2->GetParameter(2);
296 chi2ndf = fit2->GetChisquare()/fit2->GetNDF();
297 goodmin = mean - Nsigma*sig ;
298 goodmax = mean + Nsigma*sig ;
301 TLine *lline =
new TLine(goodmin, 0, goodmin, distrib->GetMaximum());
302 lline->SetLineColor(kOrange);
303 lline->SetLineStyle(7);
306 TLine *rline =
new TLine(goodmax, 0, goodmax, distrib->GetMaximum());
307 rline->SetLineColor(kOrange);
308 rline->SetLineStyle(7);
312 TLegend *leg =
new TLegend(0.60,0.82,0.9,0.88);
313 leg->AddEntry(lline,
"Good region boundary",
"l");
318 TString name =
"criteria-" ;
326 for (Int_t c = 1; c <= inhisto->GetNbinsX(); c++) {
327 cel=(Int_t)(inhisto->GetBinLowEdge(c)+0.1);
328 if (inhisto->GetBinContent(c) < goodmin) {
332 else if (inhisto->GetBinContent(c) > goodmax) {
350 void TestCellEandN(Int_t *pflag[11520][7], Double_t Emin = 0.1, Double_t Nsigma = 4.,
char* hname =
"hCellAmplitude", Int_t dnbins = 200)
352 TH2 *hCellAmplitude = (TH2*) gFile->Get(hname);
355 Int_t ncells = hCellAmplitude->GetNbinsX();
356 Double_t amin = hCellAmplitude->GetXaxis()->GetXmin();
357 Double_t amax = hCellAmplitude->GetXaxis()->GetXmax();
359 TH1* hCellEtotal =
new TH1F(Form(
"%s_hCellEtotal_E%.2f",hname,Emin),
360 Form(
"Total deposited energy, E > %.2f GeV",Emin), ncells,amin,amax);
361 hCellEtotal->SetXTitle(
"AbsId");
362 hCellEtotal->SetYTitle(
"Energy, GeV");
364 TH1F *hCellNtotal =
new TH1F(Form(
"%s_hCellNtotal_E%.2f",hname,Emin),
365 Form(
"Number of entries per events, E > %.2f GeV",Emin), ncells,amin,amax);
366 hCellNtotal->SetXTitle(
"AbsId");
367 hCellNtotal->SetYTitle(
"Entries");
369 TH1F *hCellEtoNtotal =
new TH1F(Form(
"%s_hCellEtoNtotal_E%.2f",hname,Emin),
370 Form(
"Average energy per hit, E > %.2f GeV",Emin), ncells,amin,amax);
371 hCellEtoNtotal->SetXTitle(
"AbsId");
372 hCellEtoNtotal->SetYTitle(
"Energy, GeV");
374 TH1* hNEventsProcessedPerRun = (TH1*) gFile->Get(
"hNEventsProcessedPerRun");
375 Double_t totalevents = hNEventsProcessedPerRun->Integral(1, hNEventsProcessedPerRun->GetNbinsX());
378 for (Int_t c = 1; c <= ncells; c++) {
383 for (Int_t j = 1; j <= hCellAmplitude->GetNbinsY(); j++) {
384 Double_t E = hCellAmplitude->GetYaxis()->GetBinCenter(j);
385 Double_t N = hCellAmplitude->GetBinContent(c, j);
386 if (E < Emin)
continue;
391 hCellEtotal->SetBinContent(c, Esum);
392 hCellNtotal->SetBinContent(c, Nsum/totalevents);
395 hCellEtoNtotal->SetBinContent(c, Esum/Nsum);
399 delete hCellAmplitude;
403 Process(pflag, hCellEtoNtotal, Nsigma, dnbins );
405 Process(pflag, hCellNtotal, Nsigma, dnbins );
415 void TestCellShapes(Int_t *pflag[11520][7], Double_t fitEmin, Double_t fitEmax, Double_t Nsigma =4.,
char* hname =
"hCellAmplitude", Int_t dnbins = 1000)
417 TH2 *hCellAmplitude = (TH2*) gFile->Get(Form(
"%s",hname));
420 Int_t ncells = hCellAmplitude->GetNbinsX();
421 Double_t amin = hCellAmplitude->GetXaxis()->GetXmin();
422 Double_t amax = hCellAmplitude->GetXaxis()->GetXmax();
425 TH1 *hFitA =
new TH1F(Form(
"hFitA_%s",hname),
"Fit A value", ncells,amin,amax);
426 hFitA->SetXTitle(
"AbsId");
427 hFitA->SetYTitle(
"A");
429 TH1 *hFitB =
new TH1F(Form(
"hFitB_%s",hname),
"Fit B value", ncells,amin,amax);
430 hFitB->SetXTitle(
"AbsId");
431 hFitB->SetYTitle(
"B");
433 TH1 *hFitChi2Ndf =
new TH1F(Form(
"hFitChi2Ndf_%s",hname),
"Fit #chi^{2}/ndf value", ncells,amin,amax);
434 hFitChi2Ndf->SetXTitle(
"AbsId");
435 hFitChi2Ndf->SetYTitle(
"#chi^{2}/ndf");
437 Double_t maxval1=0., maxval2=0., maxval3=0.;
438 Double_t prev=0., MSA=0., AvA = 0. ;
439 Double_t prev2=0., MSB=0., AvB = 0. ;
440 Double_t prev3=0., MSki2=0., Avki2 = 0. ;
443 for (Int_t k = 1; k <= ncells; k++) {
444 TF1 *fit =
new TF1(
"fit",
"[0]*exp(-[1]*x)/x^2");
445 TH1 *hCell = hCellAmplitude->ProjectionY(
"",k,k);
446 if (hCell->GetEntries() == 0)
continue;
448 hCell->Fit(fit,
"0QEM",
"", fitEmin, fitEmax);
451 if(fit->GetParameter(0) < 5000.){
452 hFitA->SetBinContent(k, fit->GetParameter(0));
454 AvA += fit->GetParameter(0);
455 if(k==2999) maxval1 = AvA/3000. ;
456 if (prev < fit->GetParameter(0)) MSA += fit->GetParameter(0) - prev;
457 else MSA -= (fit->GetParameter(0) - prev) ;
458 prev = fit->GetParameter(0);
464 if((fit->GetParameter(0) - maxval1) > 0. && (fit->GetParameter(0) - maxval1) < (MSA/1000.))
466 maxval1 = fit->GetParameter(0);
470 else hFitA->SetBinContent(k, 5000.);
474 if(fit->GetParameter(1) < 5000.){
475 hFitB->SetBinContent(k, fit->GetParameter(1));
477 AvB += fit->GetParameter(1);
478 if(k==2999) maxval2 = AvB/3000. ;
479 if (prev2 < fit->GetParameter(1)) MSB += fit->GetParameter(1) - prev2;
480 else MSB -= (fit->GetParameter(1) - prev2) ;
481 prev2 = fit->GetParameter(1);
487 if((fit->GetParameter(1) - maxval2) > 0. && (fit->GetParameter(1) - maxval2) < (MSB/1000.))
489 maxval2 = fit->GetParameter(1);
493 else hFitB->SetBinContent(k, 5000.);
496 if (fit->GetNDF() != 0 ) ki2 = fit->GetChisquare()/fit->GetNDF();
500 hFitChi2Ndf->SetBinContent(k, ki2);
503 if(k==2999) maxval3 = Avki2/3000. ;
504 if (prev3 < ki2) MSki2 += ki2 - prev3;
505 else MSki2 -= (ki2 - prev3) ;
512 if((ki2 - maxval3) > 0. && (ki2 - maxval3) < (MSki2/1000.))
518 else hFitChi2Ndf->SetBinContent(k, 1000.);
524 delete hCellAmplitude;
533 Process(pflag, hFitChi2Ndf, Nsigma, dnbins, maxval3);
537 Process(pflag, hFitA, Nsigma, dnbins, maxval1);
541 Process(pflag, hFitB, Nsigma, dnbins, maxval2);
552 TH2 *hCellAmplitude = (TH2*) gFile->Get(
"hCellAmplitude");
555 for (Int_t c = 1; c <= 11520; c++) {
558 for (Int_t l = 1; l <= hCellAmplitude->GetNbinsY(); l++) {
559 Double_t N = hCellAmplitude->GetBinContent(c, l);
562 if(Nsum < 0.5 && *pexclu[c-1]!=1) *pexclu[c-1]=1;
567 delete hCellAmplitude;
578 TH2 *hCellAmplitude = (TH2*) gFile->Get(
"hCellAmplitude");
580 for(Int_t i =0; i<nbc; i++){
581 for(Int_t j=0; j<= hCellAmplitude->GetNbinsY() ;j++){
582 hCellAmplitude->SetBinContent(filter[i]+1,j,0) ; }}
584 TH1* hNEventsProcessedPerRun = (TH1*) gFile->Get(
"hNEventsProcessedPerRun");
586 TFile *tf =
new TFile(
"filter.root",
"recreate");
587 hCellAmplitude->Write();
588 hNEventsProcessedPerRun->Write();
591 delete hCellAmplitude;
delete hNEventsProcessedPerRun;
609 Double_t Emax=1.0, TString
file =
"none")
612 Int_t *pexclu[11520] ;
614 Int_t *pflag[11520][7] ;
615 Int_t flag[11520][7];
620 TString output, bilan;
621 if(criterum == 7) bilan =
"Results.txt" ;
622 output.Form(
"Criterum-%d_Emin-%.2f.txt",criterum,Emin);
623 for(i=0;i<11520;i++) { exclu[i]=0; pexclu[i] =&exclu[i];
624 for(j=0;j<7;j++) { flag[i][j] =1 ; pflag[i][j] = &flag[i][j];}}
625 flag[0][0]=criterum ;
631 cout<<
"Excluded/dead cells : "<<endl;
632 for(i=0;i<11520;i++) {
if(exclu[i]!=0) {cout<<i<<
", " ; nb++;}}
633 cout<<
"("<<nb<<
")"<<endl; nb=0;}
638 cout<<
"FINAL RESULTS"<<endl;
639 ofstream fichier(bilan, ios::out | ios::trunc);
641 fichier<<
"Dead cells : "<<endl;
642 cout<<
"Dead cells : "<<endl;
643 for(i=0;i<11520;i++) {
644 if(exclu[i]!=0) {fichier<<i<<
", " ; cout<<i<<
", " ; nb++;}}
645 fichier<<
"("<<nb<<
")"<<endl; cout<<
"("<<nb<<
")"<<endl; nb=0;
647 TFile::Open(
"filter.root");
649 fichier<<
"Bad cells candidates : "<<endl; cout<<
"Bad cells candidates : "<<endl;
650 for(i=0;i<11520;i++) {
651 if(exclu[i]!=0) {bad[nb]=i; fichier<<i<<
", " ; cout<<i<<
", " ;
652 nb++;
if(nb==999){ cout<<
"TO MUCH BAD CELLS"<<endl ;
break;}}}
653 fichier<<
"("<<nb<<
")"<<endl; cout<<
"("<<nb<<
")"<<endl;}
660 for(w=0; (w*9)<=nb; w++) {
661 if(9<=(nb-w*9)) c = 9 ;
672 else if (criterum < 6)
677 if(criterum < 6) { nb=0;
678 cout<<
"bad by lower value : "<<endl;
679 for(i=0;i<11520;i++) {
680 if(flag[i][criterum]==0 && exclu[i]==0){nb++;
681 cout<<i<<
", " ;}} cout<<
"("<<nb<<
")"<<endl; nb=0;
683 cout<<
"bad by higher value : "<<endl;
684 for(i=0;i<11520;i++) {
685 if(flag[i][criterum]==2 && exclu[i]==0) {nb++;
686 cout<<i<<
", " ;}} cout<<
"("<<nb<<
")"<<endl; nb=0;
688 cout<<
"total bad "<<endl;
689 for(i=0;i<11520;i++) {
690 if(flag[i][criterum]!=1 && exclu[i]==0) {
693 cout<<i<<
", " ; }} cout<<
"("<<nb<<
")"<<endl;
700 ofstream fichier(output, ios::out | ios::trunc);
703 fichier <<
"criterum : "<<criterum<<
", Emin = "<<Emin<<
" GeV"<<
", Emax = "<<Emax<<
" GeV"<<endl;
704 fichier<<
"bad by lower value : "<<endl;
705 for(i=0;i<11520;i++) {
706 if(flag[i][criterum]==0 && exclu[i]==0){nb++;
707 fichier<<i<<
", " ;}} fichier<<
"("<<nb<<
")"<<endl; nb=0;
709 fichier<<
"bad by higher value : "<<endl;
710 for(i=0;i<11520;i++) {
711 if(flag[i][criterum]==2 && exclu[i]==0) {nb++;
712 fichier<<i<<
", " ;}} fichier<<
"("<<nb<<
")"<<endl; nb=0;
714 fichier<<
"total bad "<<endl;
715 for(i=0;i<11520;i++) {
716 if(flag[i][criterum]!=1 && exclu[i]==0) {
719 fichier<<i<<
", " ; }} fichier<<
"("<<nb<<
")"<<endl;
723 cerr <<
"opening error" << endl;
739 if(trigger==
"default")
743 TFile::Open(
"filter.root");
745 TFile::Open(
"filter.root");
747 TFile::Open(
"filter.root");
749 TFile::Open(
"filter.root");
758 TFile::Open(
"filter.root");
760 TFile::Open(
"filter.root");
762 TFile::Open(
"filter.root");
764 TFile::Open(
"filter.root");
779 TString pass =
"pass1_HLT", TString trigger=
"default")
781 Convert(fCalorimeter, period, pass, trigger);
782 TString inputfile =
"/scratch/alicehp2/germain/QA/" + period;
783 inputfile += trigger;
784 inputfile +=
".root";
void TestCellEandN(Int_t *pflag[11520][7], Double_t Emin=0.1, Double_t Nsigma=4., char *hname="hCellAmplitude", Int_t dnbins=200)
void Process(Int_t *pflag[11520][7], TH1 *inhisto, Double_t Nsigma=4., Int_t dnbins=200, Double_t dmaxval=-1.)
void TestCellShapes(Int_t *pflag[11520][7], Double_t fitEmin, Double_t fitEmax, Double_t Nsigma=4., char *hname="hCellAmplitude", Int_t dnbins=1000)
void PeriodAnalysis(Int_t criterum=7, Double_t Nsigma=4.0, Double_t Emin=0.1, Double_t Emax=1.0, TString file="none")
void BCAnalysis(TString file, TString trigger="default")
void BadChannelAnalysis(TString fCalorimeter="EMCAL", TString period="LHC11h", TString pass="pass1_HLT", TString trigger="default")
void Draw(Int_t cell[], Int_t iBC, Int_t nBC, const Int_t cellref=151)
void KillCells(Int_t filter[], Int_t nbc)
void Convert(TString fCalorimeter="EMCAL", TString period="LHC11h", TString pass="pass1_HLT", TString trigger="default")
void ExcludeCells(Int_t *pexclu[11520])