AliPhysics  1811c8f (1811c8f)
trendingCluster.C
Go to the documentation of this file.
1 
27 
28 /* $Id: $ */
29 //--------------------------------------------------
30 
31 //
32 // Author: Yaxian Mao
33 // Modified M. Germain
34 //
35 
36 #include <TFile.h>
37 #include <TH1F.h>
38 #include <TH2F.h>
39 #include <TH3D.h>
40 #include <Riostream.h>
41 #include <TCanvas.h>
42 #include <TGraphErrors.h>
43 #include <TGrid.h>
44 #include <TFileMerger.h>
45 #include <TMultiGraph.h>
46 #include <TROOT.h>
47 #include <TString.h>
48 #include <TStyle.h>
49 #include <TLegend.h>
50 #include <TGridCollection.h>
51 #include <TROOT.h>
52 #include <TGridResult.h>
53 #include <TClonesArray.h>
54 #include <TObjString.h>
55 #include <stdio.h>
56 
57 #include <string>
58 #include <fstream>
59 #include <iostream>
60 #include <sstream> // Required for stringstreams
61 
62 using namespace std;
63 
64 void trendingCluster(TString fCalorimeter = "EMCAL", TString period = "LHC11h", TString pass = "pass1_HLT", const Int_t n = 10, TString fTrigger = "MB")
65 {
66 
67  FILE * pFile;
68 
69  TString file = "";
70  if (fTrigger=="EMC") file = "/scratch/alicehp2/germain/QA/"+period+"/"+ pass + "/runlist.txt" ;
71  else file = "/scratch/alicehp2/germain/QA/"+period+"/"+ pass + "/runlistMB.txt" ;
72 
73  //cout<<file<<endl;
74  pFile = fopen(file.Data(), "r"); //open the text file where include the run list and correct run index
75 
76  cout<<file<<endl;
77  cout << " fcalo: " << fCalorimeter << "; period: " << period << "; pass: " << pass << " trigger "<<fTrigger<< endl;
78  char outfilename [100] ;
79 
80  Int_t index[500];
81  Int_t p;
82  Int_t q;
83  Int_t ncols;
84  Int_t nlines = 0 ;
85  Int_t RunId[500] ;
86  Double_t x[500] ;
87  Double_t xrun[500] ;
88  //TString label;
89 
90  while (1){
91  ncols = fscanf(pFile,"%d %d ",&p,&q);
92  if (ncols< 0) break;
93  x[nlines]=p;
94  index[nlines]=p;
95  RunId[nlines]=q;
96  xrun[nlines]=1.*q;
97  nlines++;
98  }
99  fclose(pFile);
100  const Int_t nRun = nlines ;
101 
102  Double_t xe[nRun] ;
103  Double_t nEvent[nRun] ;
104  // Double_t TimeMean[nRun] ;
105  // Double_t TimeRMS[nRun] ;
106  Double_t CellMean[nRun] ;
107  Double_t CellRMS[nRun] ;
108  Double_t ClusterMean[nRun] ;
109  Double_t ClusterRMS[nRun] ;
110  Double_t EtotalMean[nRun] ;//total energy deposited per event
111  Double_t EtotalRMS[nRun] ;
112 
113  Double_t CellPerClusterMean[nRun] ;
114  Double_t CellPerClusterRMS[nRun] ;
115  Double_t ECell1Mean[nRun] ;//total energy deposited per event without 1 cell clusters
116  Double_t ECell1RMS[nRun] ;
117 
118  TFile * f ;
119  TDirectoryFile * dir;
120  TList * outputList;
121 
122  TH1F * fhNEvents;
123  TH1F * fhE;
124  TH1F * fhNClusters;
125  TH1F * fhNCells;
126  TH2F * fhIM ;
127  TH3D * fhNCellsPerCluster ;
128  TH2F * fhTimeAmp ;
129 
130  TH1F * NCells[n];
131  TH1F * NClusters[n];
132  TH2F * NCellsPerCluster[n];
133  TH1F * E[n];
134 
135  TGraphErrors * AverNcellsSM[n];
136  TGraphErrors * AverNclustersSM[n];
137  TGraphErrors * AverNcellsPerClusterSM[n];
138  TGraphErrors * AverESM[n];
139  TGraphErrors * AverTimeSM[n];
140  TGraphErrors * AverMggSM[n];
141  TGraphErrors * AverEcell1SM[n];
142 
143  Double_t CellMeanSM[n][nRun] ;
144  Double_t CellRMSSM[n][nRun] ;
145  Double_t ClusterMeanSM[n][nRun] ;
146  Double_t ClusterRMSSM[n][nRun] ;
147  Double_t EtotalMeanSM[n][nRun] ;//total energy deposited per event
148  Double_t EtotalRMSSM[n][nRun] ;
149  Double_t CellPerClusterMeanSM[n][nRun] ;
150  Double_t CellPerClusterRMSSM[n][nRun] ;
151  Double_t ECell1MeanSM[n][nRun] ;//total energy deposited per event without 1 cell clusters
152  Double_t ECell1RMSSM[n][nRun] ;
153 
154 
155  TString namefile = "/scratch/alicehp2/germain/QA/"+period+"/"+pass+"/"+ fCalorimeter + period + pass + fTrigger+"data.txt";
156 
157  fstream QAData(namefile, ios::out); //write the QA check values at the end
158 
159  cout << " namefile " << namefile << endl;
160  cout << " nRun " << nRun << " index(nRun) " << index[nRun-1] << endl;
161 
162  for(Int_t i = 0 ; i < nRun ; i++){
163 
164  xe[i] = 0.5 ;
165 
166  TString name = "/scratch/alicehp2/germain/QA/"+period +"/"+ pass + "/";
167 
168  name += RunId[i] ;
169  name += ".root";
170  f = TFile::Open(name.Data(),"read") ;
171  if (!f) continue;
172 
173  // cout << " i = " << i << " file opend Output" << RunId[i]<< endl;
174 
175  if(fTrigger=="EMC"){ dir = (TDirectoryFile *)f->Get("CaloQA_EMC7");
176  outputList = (TList*)dir->Get("CaloQA_EMC7");
177  }
178  else{
179  dir = (TDirectoryFile *)f->Get("CaloQA_default");
180  outputList = (TList*)dir->Get("CaloQA_default");
181  }
182 
183  //define the averages for checking Histograms
184  //
185 
186  fhNEvents =(TH1F *)outputList->FindObject("hNEvents");
187  nEvent[i]=fhNEvents->GetEntries();
188  cout << " Run " << RunId[i]<< " nevent " << nEvent[i]<< endl;
189  if( nEvent[i] == 0) continue ;
190  if( nEvent[i] < 2) continue ;
191  fhE = (TH1F *)outputList->FindObject(fCalorimeter+"_hE");
192 
193  Double_t energy = 0. ;
194 
195  for(Int_t ibin = fhE->FindBin(0.3) ; ibin <fhE->FindBin(50.) ; ibin++){ //Starting from 0.3eV
196  energy+=fhE->GetBinCenter(ibin)*fhE->GetBinContent(ibin);
197  }
198  EtotalMean[i]=energy/fhE->Integral(fhE->FindBin(0.3), fhE->FindBin(50.)) ;
199  EtotalRMS[i]=fhE->GetMeanError();
200 
201 
202 
203  //for single module check
204  TH1F*fhNCells =0;
205  TH1F*fhNClusters =0;
206  for(Int_t ism = 0 ; ism < n ; ism++)
207  {
208  TString nameNCell = Form("%s_hNCells_Mod%d",fCalorimeter.Data(),ism);
209 
210  TString nameNCluster = Form("%s_hNClusters_Mod%d",fCalorimeter.Data(),ism);
211  TString nameNCellPerCluster = Form("%s_hNCellsPerCluster_Mod%d",fCalorimeter.Data(),ism);
212  TString nameE = Form("%s_hE_Mod%d",fCalorimeter.Data(),ism);
213 
214  NCells[ism] = (TH1F*)outputList->FindObject(nameNCell.Data());
215  NClusters[ism] = (TH1F*)outputList->FindObject(nameNCluster.Data());
216  NCellsPerCluster[ism] = (TH2F*)outputList->FindObject(nameNCellPerCluster.Data());
217  E[ism] = (TH1F*)outputList->FindObject(nameE.Data());
218  CellMeanSM[ism][i]=NCells[ism]->GetMean();
219  CellRMSSM[ism][i]=NCells[ism]->GetMeanError();
220  ClusterMeanSM[ism][i]=NClusters[ism]->GetMean();
221  ClusterRMSSM[ism][i]=NClusters[ism]->GetMeanError();
222  CellPerClusterMeanSM[ism][i]=NCellsPerCluster[ism]->GetMean(2);
223  CellPerClusterRMSSM[ism][i]=NCellsPerCluster[ism]->GetMeanError(2);
224  // cout<<"SM = "<<ism<<" Mean : = "<<CellPerClusterMeanSM[ism][i]<<endl ;
225  ECell1MeanSM[ism][i] =NCellsPerCluster[ism]->ProjectionX("",2,300,"")->Integral(5,100)/(nEvent[i]);
226  ECell1RMSSM[ism][i] =NCellsPerCluster[ism]->ProjectionX("",2,300,"")->GetMeanError();
227  Double_t energySM = 0. ;
228  for(Int_t ibin = E[ism]->FindBin(0.3) ; ibin <E[ism]->FindBin(50.) ; ibin++){ //starting from 0.3GeV
229  energySM+=E[ism]->GetBinCenter(ibin)*(E[ism]->GetBinContent(ibin));
230  }
231  EtotalMeanSM[ism][i]=energySM/(E[ism]->Integral(E[ism]->FindBin(0.3),E[ism]->FindBin(50.)));
232 
233  EtotalRMSSM[ism][i]=E[ism]->GetMeanError();
234 
235  if(ism==0) {
236  fhNCells = (TH1F*)NCells[ism]->Clone("NCells");
237  fhNClusters = (TH1F*)NClusters[ism]->Clone("NClusters");
238  }
239  else {
240  fhNCells->Add(NCells[ism],1);
241  fhNClusters->Add(NClusters[ism],1);
242  }
243  } //per SM loop
244  ClusterMean[i]=fhNClusters->GetMean();
245  ClusterRMS[i]=fhNClusters->GetMeanError();
246  CellMean[i]=fhNCells->GetMean();
247  CellRMS[i]=fhNCells->GetMeanError();
248  outputList->Clear() ;
249  dir->Close();
250  f->Close();
251  f=NULL;
252  dir=NULL;
253  outputList=NULL;
254 
255  //if you want to write all the QA check values at the end, otherwise just comment out below
256  // becareful with different detectors as the check output are different since different modules/SM for different detetcor
257 
258  QAData <<i+1<<" "<< RunId[i] <<" "<< nEvent[i]
259  // <<" "<< RunId[i]<<" "<<CellMean[i]
260  // <<" "<< CellMeanSM[0][i] <<" "<< CellMeanSM[1][i]
261  // <<" "<<CellMeanSM[2][i]
262  // <<" "<< RunId[i] <<" "<<ClusterMean[i]
263  // <<" "<< ClusterMeanSM[0][i]<<" "<< ClusterMeanSM[1][i]
264  // <<" "<< ClusterMeanSM[2][i]
265  // <<" "<< RunId[i]<<" "<< CellPerClusterMean[i]
266  // <<" "<< CellPerClusterMeanSM[0][i]<<" "<< CellPerClusterMeanSM[1][i]
267  // <<" "<< CellPerClusterMeanSM[2][i]
268  // <<" "<< RunId[i]<<" "<< EtotalMean[i]
269  // <<" "<< EtotalMeanSM[0][i] <<" "<< EtotalMeanSM[1][i]
270  // <<" "<< EtotalMeanSM[2][i]
271 
272  <<"\n" ;
273 
274  } // end loop on nrun
275 
276 
277  QAData.close();
278 
279  TString base = "/scratch/alicehp2/germain/QA/";
280  base += period ;
281  base += "/";
282  base += pass ;
283  base += "/";
284 
285  base += fTrigger;
286  TString ClusterAverages ; ClusterAverages = base + "ClusterAverages.gif";
287  TString Entries ; Entries = base + "Nentries.gif";
288  TString ClusterAveragesEnergy ; ClusterAveragesEnergy = base + "ClusterAveragesEnergy.gif";
289  TString ClusterAveragesEntries ; ClusterAveragesEntries = base + "ClusterAveragesEntries.gif";
290  TString ClusterAveragesCells ; ClusterAveragesCells = base + "ClusterAveragescells.gif";
291 
292 
293  cout << "c11 nEvents" << endl;
294  //just for the canvas defination
295 
296  cout << " index(0)" << index[nRun-1] << " index(20) " << index[20] <<endl;
297  TH1F * dummy = new TH1F("dummy", "dummy", nRun, 0., nRun+0.5);
298  // TH1F * dummy = new TH1F("dummy", "dummy", index[nRun-1], 0., index[nRun-1]+0.5);
299  dummy->SetTitle("") ;
300  dummy->SetStats(kFALSE) ;
301  dummy->SetAxisRange(0., nRun, "X") ;
302 
303  for(Int_t i = 0 ; i < nRun ; i++){
304  TString label = " ";
305  label+=RunId[i];
306  cout <<" run "<< RunId[i] << " label " <<label << endl;
307 
308  dummy->GetXaxis()->SetBinLabel(i+1,label.Data());
309  dummy->GetXaxis()->LabelsOption("v");
310  }
311 
312 
313 
314  //number of events passes physics selection for each run
315  TCanvas * c11 = new TCanvas("nEvents", "nEvents", 1000, 500);
316  c11->SetFillColor(0);
317  c11->SetBorderSize(0);
318  c11->SetFrameBorderMode(0);
319  gStyle->SetOptStat(0);
320  c11->SetLogy();
321  c11->SetGrid();
322  // dummy->GetXaxis()->SetTitle("RUN");
323  dummy->GetXaxis()->SetTitleOffset(0.05);
324  dummy->GetYaxis()->SetTitle("N_{events}");
325  dummy->SetMinimum(1.) ; //should addjust based on the statistics
326  dummy->SetMaximum(1.e6) ; //should addjust based on the statistics
327  dummy->Draw();
328  TGraph * nEvents = new TGraph(nRun, x, nEvent);
329  nEvents->SetMarkerStyle(20);
330  nEvents->SetMarkerColor(1);
331  nEvents->SetLineColor(2);
332  nEvents->Draw("same PL") ;
333  c11->Update();
334 
335 
336  if (fTrigger=="MB")sprintf(outfilename,"nEventMB.gif");
337  if (fTrigger=="EMC")sprintf(outfilename,"nEventEMC.gif");
338 
339 
340  c11->SaveAs(Entries);
341 
342 
343  cout << "c1 Aver NCell" << endl;
344 
345  TCanvas * c1 = new TCanvas("AverNCell", "AverNCell", 600, 600);
346  // c1->SetLogy();
347  c1->SetFillColor(0);
348  c1->SetBorderSize(0);
349  c1->SetFrameBorderMode(0);
350  gStyle->SetOptStat(0);
351  TH1F * h1 = (TH1F*)dummy->Clone("");
352  h1->GetXaxis()->SetTitle("RUN Index");
353  h1->GetYaxis()->SetTitle("<N_{cells}>");
354  h1->SetMinimum(0.) ;
355  h1->SetMaximum(10.) ;
356  if(fCalorimeter=="EMCAL") h1->SetMaximum(5.) ;
357  h1->Draw();
358 
359  TGraphErrors * AverNcells = new TGraphErrors(nRun, x, CellMean, xe, CellRMS);
360 
361  AverNcells->SetMarkerColor(1);
362  AverNcells->SetMarkerStyle(20);
363  AverNcells->Draw("same P") ;
364  for(Int_t ism = 0 ; ism < n ; ism++){
365  AverNcellsSM[ism] = new TGraphErrors(nRun, x, CellMeanSM[ism], xe, CellRMSSM[ism]);
366  AverNcellsSM[ism]->SetMarkerColor(ism+2);
367  AverNcellsSM[ism]->SetMarkerStyle(21+ism);
368  AverNcellsSM[ism]->Draw("same P");
369  }
370 
371 
372 
373  TLegend * l1 = new TLegend(0.4, 0.6, 0.75, 0.85);
374  l1->SetFillColor(0);
375  l1->SetBorderSize(0);
376  l1->SetTextSize(0.02);
377  l1->AddEntry(AverNcells, "<# of cells>", "") ;
378  l1->AddEntry(AverNcells, Form("det = %s",fCalorimeter.Data()), "") ;
379  l1->AddEntry(AverNcells,"average", "p");
380  for(Int_t ism = 0 ; ism < n ; ism++){
381  TString projname = Form("SM_ %d",ism);
382  l1->AddEntry(AverNcellsSM[ism],projname.Data(), "p");
383  }
384  l1->Draw("same");
385  c1->Update();
386 
387 
388  TCanvas * c200 = new TCanvas("ClusterAverages", "ClusterAverages", 1000, 500);
389  c200->SetFillColor(0);
390  c200->SetBorderSize(0);
391  c200->SetFrameBorderMode(0);
392  c200->SetGrid();
393 
394  gPad->SetLeftMargin(0.08);
395  gPad->SetRightMargin(0.02);
396  gPad->SetGrid();
397 
398  TH1F * h2 = (TH1F*)dummy->Clone("");
399 
400  h2->GetYaxis()->SetTitle("<N_{clusters}>/event");
401 
402  if (fTrigger=="EMC") h2->SetMaximum(70) ;
403  else h2->SetMaximum(50) ;
404  h2->SetMinimum(10) ; // for Pb Pb
405  h2->Draw();
406 
407  TGraphErrors * AverNclusters = new TGraphErrors(nRun, x, ClusterMean, xe, ClusterRMS);
408  AverNclusters->SetMarkerStyle(20);
409  AverNclusters->SetMarkerColor(1);
410  AverNclusters->Draw("same P") ;
411 
412  for(Int_t ism = 0 ; ism < n ; ism++){
413  AverNclustersSM[ism] = new TGraphErrors(nRun, x, ClusterMeanSM[ism], xe, ClusterRMSSM[ism]);
414  AverNclustersSM[ism]->SetMarkerColor(ism+2);
415  AverNclustersSM[ism]->SetMarkerStyle(21+ism);
416 
417  AverNclustersSM[ism]->Draw("same P");
418  }
419 
420 
421  TLegend * l200 = new TLegend(0.4, 0.6, 0.75, 0.85);
422  l200->SetFillColor(0);
423  l200->SetBorderSize(0);
424  l200->SetTextSize(0.02);
425  l200->AddEntry(AverNclusters, "<# of clusters>", "") ;
426  l200->AddEntry(AverNclusters, Form("det = %s",fCalorimeter.Data()), "") ;
427  l200->AddEntry(AverNclusters,"average", "p");
428 
429  for(Int_t ism = 0 ; ism < n ; ism++){
430  TString projname = Form("SM_ %d",ism);
431  l200->AddEntry(AverNclustersSM[ism],projname.Data(), "p");
432  }
433 
434  l200->Draw("same");
435  c200->Update();
436  c200->SaveAs(ClusterAveragesEntries);
437 
438  TCanvas * c201 = new TCanvas("ClusterAveragesEnergy", "ClusterAveragesEnergy", 1000, 500);
439  c201->SetFillColor(0);
440  c201->SetBorderSize(0);
441  c201->SetFrameBorderMode(0);
442  c201->SetGrid();
443 
444 
445  gPad->SetLeftMargin(0.08);
446  gPad->SetRightMargin(0.02);
447  gPad->SetGrid();
448 
449 
450  TH1F * h3 = (TH1F*)dummy->Clone("");
451 
452  h3->GetYaxis()->SetTitle("<E> (GeV)");
453  h3->SetMinimum(0.2) ;
454  h3->SetMaximum(1.) ;
455  h3->Draw();
456 
457  TGraphErrors * AverE = new TGraphErrors(nRun, x, EtotalMean, xe, EtotalRMS);
458  AverE->SetMarkerStyle(20);
459  AverE->SetMarkerColor(1);
460  AverE->Draw("same P");
461 
462  for(Int_t ism = 0 ; ism < n ; ism++){
463  AverESM[ism] = new TGraphErrors(nRun, x, EtotalMeanSM[ism], xe, EtotalRMSSM[ism]);
464  AverESM[ism]->SetMarkerColor(ism+2);
465  AverESM[ism]->SetMarkerStyle(21+ism);
466  AverESM[ism]->Draw("same P");
467 
468  }
469 
470  TLegend * l3 = new TLegend(0.4, 0.6, 0.75, 0.85);
471  l3->SetFillColor(0);
472  l3->SetBorderSize(0);
473  l3->SetTextSize(0.02);
474  l3->AddEntry(AverE, "<E>", "") ;
475  l3->AddEntry(AverE, Form("det = %s",fCalorimeter.Data()), "") ;
476  l3->AddEntry(AverE,"average", "p");
477 
478  for(Int_t ism = 0 ; ism < n ; ism++){
479  TString projname = Form("SM_ %d",ism);
480  l3->AddEntry(AverESM[ism],projname.Data(), "p");
481  }
482 
483  l3->Draw("same");
484 
485  c201->SaveAs(ClusterAveragesEnergy);
486 
487 
488  TCanvas * c202 = new TCanvas("ClusterAveragesCells", "ClusterAveragesCells", 1000, 500);
489  c202->SetFillColor(0);
490  c202->SetBorderSize(0);
491  c202->SetFrameBorderMode(0);
492  c202->SetGrid();
493 
494  gPad->SetLeftMargin(0.08);
495  gPad->SetRightMargin(0.02);
496 
497  gPad->SetGrid();
498 
499 
500  TH1F * h4 = (TH1F*)dummy->Clone("");
501 
502  h4->GetYaxis()->SetTitle("<N_{CellsPerCluster}>");
503  h4->SetMinimum(1.5) ;
504  h4->SetMaximum(5.5) ;
505  h4->Draw();
506 
507 
508  TGraphErrors * AverCellPerCluster = new TGraphErrors(nRun, x, CellPerClusterMean, xe, CellPerClusterRMS);
509  AverCellPerCluster->SetMarkerStyle(20);
510  AverCellPerCluster->SetMarkerColor(1);
511 
512  for(Int_t ism = 0 ; ism < n ; ism++){
513  AverNcellsPerClusterSM[ism] = new TGraphErrors(nRun, x, CellPerClusterMeanSM[ism], xe, CellPerClusterRMSSM[ism]);
514  AverNcellsPerClusterSM[ism]->SetMarkerColor(ism+2);
515  AverNcellsPerClusterSM[ism]->SetMarkerStyle(21+ism);
516  AverNcellsPerClusterSM[ism]->Draw("same P");
517 
518  }
519 
520  TLegend * l4 = new TLegend(0.4, 0.6, 0.75, 0.85);
521  l4->SetFillColor(0);
522  l4->SetBorderSize(0);
523  l4->SetTextSize(0.02);
524  l4->AddEntry(AverCellPerCluster, "<# of cells per cluster>", "") ;
525  l4->AddEntry(AverCellPerCluster, Form("det = %s",fCalorimeter.Data()), "") ;
526  l4->AddEntry(AverCellPerCluster,"average", "p");
527 
528  for(Int_t ism = 0 ; ism < n ; ism++){
529  TString projname = Form("SM_ %d",ism);
530  l4->AddEntry(AverNcellsPerClusterSM[ism],projname.Data(), "p");
531  }
532 
533  l4->Draw("same");
534 
535  c202->SaveAs(ClusterAveragesCells);
536 
537 }
double Double_t
Definition: External.C:58
Definition: External.C:236
energy
Definition: HFPtSpectrum.C:44
int Int_t
Definition: External.C:63
Definition: External.C:252
Double_t nEvents
plot quality messages
TFile * file
TList with histograms for a given trigger.
TDirectoryFile * dir