AliPhysics  648edd6 (648edd6)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
helperMacrosBC.C
Go to the documentation of this file.
1 
21 // --- ROOT system ---
22 #include <Riostream.h>
23 #include <TFile.h>
24 #include <TH1.h>
25 #include <TStyle.h>
26 #include <TROOT.h>
27 #include <TCanvas.h>
28 
29 
30 // --- ANALYSIS system ---
31 #include "AliEMCALGeometry.h" //include when compile
32 #include "AliCalorimeterUtils.h" //include when compile
33 #include "AliAODEvent.h" //include when compile
34 
35 
36 //definition of methods
37 void SetHisto(TH1 *Histo,TString Xtitel,TString Ytitel,Bool_t longhisto);
38 Bool_t IsCellMaskedByHand(Int_t cell, std::vector<Int_t> cellVector);
39 void CreateCellCompPDF(TH2F* hAmpIDMasked, std::vector<Int_t> cellVector, TH1* goodCellsMerged, TH1* goodCellsRbR, TString pdfName);
40 void Plot2DCells(TString Block, Int_t runNo, std::vector<Int_t> cellVectorRbR, std::vector<Int_t> cellVectorMerge);
41 
46 void Get_RowCollumnID(Int_t runNum= 244411,Int_t inputCellID=-1,Int_t inputRow=-1,Int_t inputCollumn=-1,Int_t inputSM=-1)
47 {
48  //......................................................
49  //..Initialize EMCal/DCal geometry
50  AliCalorimeterUtils* fCaloUtils = new AliCalorimeterUtils();
51  //..Create a dummy event for the CaloUtils
52  AliAODEvent* aod = new AliAODEvent();
53  fCaloUtils->SetRunNumber(runNum);
54  fCaloUtils->AccessGeometry(aod);
55  AliEMCALGeometry * geom = fCaloUtils->GetEMCALGeometry();
56 
57  //..get row collumn from cell ID
58  Int_t cellColumn=0,cellRow=0;
59  Int_t cellColumnAbs=0,cellRowAbs=0;
60  Int_t cellID=0;
61  Int_t trash;
62 
63  cout<<"...................................................."<<endl;
64  cout<<""<<endl;
65  if(inputCellID!=-1)
66  {
67  fCaloUtils->GetModuleNumberCellIndexesAbsCaloMap(inputCellID,0,cellColumn,cellRow,trash,cellColumnAbs,cellRowAbs);
68  cout<<"Cell Id provided: "<<inputCellID<<endl;
69  cout<<"This corresponds to absolute row: "<<cellRowAbs<<" and absolute collumn: "<<cellColumnAbs<<endl;
70  }
71  else if(inputRow!=-1 && inputCollumn!=-1 && inputSM!=-1)
72  {
73  cout<<"Supermodule provided: "<<inputSM<<endl;
74  cout<<"Absolute row provided: "<< inputRow<<" and absolute collumn provided : "<<inputCollumn<<endl;
75  cellID=geom->GetAbsCellIdFromCellIndexes(inputSM,inputRow,inputCollumn);
76  cout<<"This corresponds to Cell Id: "<<cellID<<endl;
77  }
78  else cout<<"need more information!"<<endl;
79  cout<<""<<endl;
80  cout<<"...................................................."<<endl;
81 }
82 
88 //________________________________________________________________________
89 void Compare2Blocks(TString period="LHC15n",Int_t trainNo=603,Int_t versionA=0, Int_t versionB=1)
90 {
91  gROOT->ProcessLine("gErrorIgnoreLevel = kWarning;"); //..to supress a lot of standard output
92  gStyle->SetOptStat(0); //..Do not plot stat boxes
93  gStyle->SetPadLeftMargin(0.13);
94  gStyle->SetPadRightMargin(0.1);
95  gStyle->SetPadBottomMargin(0.13);
96  gStyle->SetPadTopMargin(0.02);
97 
98  Int_t noOfCells=17674;
99  Int_t nBadCellMerged =noOfCells;
100  Int_t nBadCellRbR =noOfCells;
101 
102  //..............................................
103  //..manually disable cells
104  std::vector<Int_t> badcellsBlock1;
105 
106  //badcellsBlock1.insert(badcellsBlock1.end(),{6644,6655,10140,12036,12037,12038,12039,12040,12041,12926,13067,13066,13125});
107  //badcellsBlock1.insert(badcellsBlock1.end(),{13133,13483,13971,13978,14116,14118,14122,14411,14593,14599,14600,14606,14699});
108 
109  std::vector<Int_t> vOnlyMaskedInMergedA;
110  std::vector<Int_t> vOnlyMaskedInMergedB;
111 
112  //......................................................
113  //..Get the .root file from analyzing runs as 1 block together
114  //......................................................
115  cout<<"** Open file A with merged runlist analysis: "<<endl;
116  TString pathA = Form("./AnalysisOutput/%s/Train_%i/Version%i",period.Data(),trainNo,versionA);
117  TString rootFileNameA= Form("%s_INT7_Histograms_V%i.root",period.Data(),versionA);
118  TFile* outputRootA = TFile::Open(Form("%s/%s",pathA.Data(),rootFileNameA.Data()));
119  if(!outputRootA)cout<<"File "<<outputRootA->GetName()<<" does not exist"<<endl;
120  else cout<<"file A: "<<outputRootA->GetName()<<endl;
121 
122  //..get the necessary histograms
123  TH2F* hCellAmplitudeA =(TH2F*)outputRootA->Get("hCellAmplitude");
124  TH1F* hnEventsA =(TH1F*)outputRootA->Get("hNEvents");
125  hCellAmplitudeA->Scale(hnEventsA->GetBinContent(1));
126  TH2F* hCellAmplitudeBlockMaskA=(TH2F*)hCellAmplitudeA->Clone("hCellAmplitudeMask");
127  TH1F* hCellFlagA =(TH1F*)outputRootA->Get("fhCellFlag");
128 
129  //......................................................
130  //..Get the .root file with run-by-run bad channel mask
131  //......................................................
132  cout<<endl;
133  cout<<"**Open file B with merged runlist analysis: "<<endl;
134  TString pathB = Form("./AnalysisOutput/%s/Train_%i/Version%i",period.Data(),trainNo,versionB);
135  TString rootFileNameB= Form("%s_INT7_Histograms_V%i.root",period.Data(),versionB);
136  TFile* outputRootB = TFile::Open(Form("%s/%s",pathB.Data(),rootFileNameB.Data()));
137  if(!outputRootB)cout<<"File "<<outputRootB->GetName()<<" does not exist"<<endl;
138  else cout<<"file B: "<<outputRootB->GetName()<<endl;
139 
140  //..get the necessary histograms
141  TH2F* hCellAmplitudeB =(TH2F*)outputRootB->Get("hCellAmplitude");
142  TH1F* hnEventsB =(TH1F*)outputRootB->Get("hNEvents");
143  hCellAmplitudeB->Scale(hnEventsB->GetBinContent(1));
144  TH2F* hCellAmplitudeBlockMaskB=(TH2F*)hCellAmplitudeB->Clone("hCellAmplitudeMask");
145  TH1F* hCellFlagB =(TH1F*)outputRootB->Get("fhCellFlag");
146 
147  //......................................................
148  //..mask the bad cells according to both versions
149  //......................................................
150  Int_t maskA=0;
151  Int_t maskB=0;
152 
153  for(Int_t ic = 0; ic < noOfCells; ic++)
154  {
155  //......................................................
156  //..Part A - analyzed as one merged runblock
157  //......................................................
158  maskA=0;
159  maskA=IsCellMaskedByHand(ic,badcellsBlock1);
160 
161  //..mask the bad cells
162  for (Int_t amp = 1; amp <= hCellAmplitudeBlockMaskA->GetNbinsX(); amp++)
163  {
164  if(hCellFlagA->GetBinContent(ic+1)>0 || maskA==1)
165  {
166  hCellAmplitudeBlockMaskA->SetBinContent(amp,ic+1,0);
167  if(amp==1)nBadCellMerged--;
168  }
169  }
170 
171  //......................................................
172  //..Part B - analyzed run-by-run, corrected and masked in block
173  //......................................................
174  maskB=0;
175  maskB=IsCellMaskedByHand(ic,badcellsBlock1);
176 
177  //..mask the bad cells
178  for (Int_t amp = 1; amp <= hCellAmplitudeBlockMaskB->GetNbinsX(); amp++)
179  {
180  if(hCellFlagB->GetBinContent(ic+1)>0 || maskB==1)
181  {
182  hCellAmplitudeBlockMaskB->SetBinContent(amp,ic+1,0);
183  if(amp==1)nBadCellRbR--;
184  }
185  }
186 
187  //......................................................
188  //..Compare the different channels that are marked in the two versions
189  //......................................................
190  if(!IsCellMaskedByHand(ic,badcellsBlock1))
191  {
192  if(hCellFlagA->GetBinContent(ic+1)>0 && hCellFlagB->GetBinContent(ic+1)==0)vOnlyMaskedInMergedA.push_back(ic);
193  if(hCellFlagA->GetBinContent(ic+1)==0 && hCellFlagB->GetBinContent(ic+1)>0) vOnlyMaskedInMergedB.push_back(ic);
194  }
195  }
196  //..merged runblock
197  TH1* projMaskedCellsA = hCellAmplitudeBlockMaskA->ProjectionX("MaskedCellsMergedBlockA");
198  TH1* projMaskedCellsB = hCellAmplitudeBlockMaskB->ProjectionX("MaskedCellsMergedBlockB");
199 
200  //......................................................
201  //..Plot results
202  //......................................................
203  TCanvas* C1 = new TCanvas("-1-","Projections of A and B",900,900);
204  C1->Divide(2,2);
205  C1->cd(1)->SetLogy();
206  SetHisto(projMaskedCellsA,"","hits/event",0);
207  projMaskedCellsA->DrawCopy("hist");
208  projMaskedCellsB->SetLineColor(8);
209  projMaskedCellsB->DrawCopy("same hist");
210  C1->cd(2);
211  projMaskedCellsA->Divide(projMaskedCellsB);
212  SetHisto(projMaskedCellsA,"","Merged Version A/Merged Version B",0);
213  projMaskedCellsA->DrawCopy("hist");
214  projMaskedCellsA->Multiply(projMaskedCellsB);
215 
216  TCanvas* C3 = new TCanvas("-3-","2D of A and B",900,900);
217  C3->Divide(2,2);
218  C3->cd(1)->SetLogz();
219  SetHisto(hCellAmplitudeBlockMaskA,"","cell ID",0);
220  hCellAmplitudeBlockMaskA->DrawCopy("colz");
221  C3->cd(2)->SetLogz();
222  SetHisto(hCellAmplitudeBlockMaskB,"","cell ID",0);
223  hCellAmplitudeBlockMaskB->DrawCopy("colz");
224 
225 
226  //......................................................
227  //..Print out compared cells and plot the spectra
228  //......................................................
229  projMaskedCellsA->Scale(1.0/nBadCellMerged);
230  projMaskedCellsB->Scale(1.0/nBadCellRbR);
231 
232 
233  cout<<" Cells masked in version A and not in version B ("<<vOnlyMaskedInMergedA.size()<<"):"<<endl;
234  for(Int_t i=0; i<(Int_t)vOnlyMaskedInMergedA.size();i++)
235  {
236  cout<<vOnlyMaskedInMergedA.at(i)<<","<<flush;
237  }
238  cout<<endl;
239  CreateCellCompPDF(hCellAmplitudeBlockMaskB,vOnlyMaskedInMergedA,projMaskedCellsA,projMaskedCellsB,"./cOnlyMergedBlockA.pdf");
240  cout<<" Cells masked in version B and not in version A ("<<vOnlyMaskedInMergedB.size()<<"):"<<endl;
241  for(Int_t i=0; i<(Int_t)vOnlyMaskedInMergedB.size();i++)
242  {
243  cout<<vOnlyMaskedInMergedB.at(i)<<","<<flush;
244  }
245  cout<<endl;
246  CreateCellCompPDF(hCellAmplitudeBlockMaskA,vOnlyMaskedInMergedB,projMaskedCellsA,projMaskedCellsB,"./cOnlyMergedBlockB.pdf");
247 
248  //......................................................
249  //..build two dimensional histogram with cells rejected from
250  //..the one or the other method
251  //......................................................
252  Plot2DCells("A",244917,vOnlyMaskedInMergedB,vOnlyMaskedInMergedA);
253 }
254 //-----------------------------------------------------------
255 // All functions below this point are helping the
256 // Compare2Blocks function
257 //-----------------------------------------------------------
258 
262 void Plot2DCells(TString Block, Int_t runNo, std::vector<Int_t> cellVectorRbR, std::vector<Int_t> cellVectorMerge)
263 {
264  //......................................................
265  //..Initialize EMCal/DCal geometry
266  AliCalorimeterUtils* fCaloUtils = new AliCalorimeterUtils();
267  //..Create a dummy event for the CaloUtils
268  AliAODEvent* aod = new AliAODEvent();
269  fCaloUtils->SetRunNumber(runNo);
270  fCaloUtils->AccessGeometry(aod);
271  //......................................................
272  //..setings for the 2D histogram
273  Int_t nModules=fCaloUtils->GetEMCALGeometry()->GetNumberOfSuperModules();
274  Int_t fNMaxColsAbs = 2*48;
275  Int_t fNMaxRowsAbs = Int_t (nModules/2)*24; //multiply by number of supermodules
276 
277  TH2F *plot2D_RbR = new TH2F(Form("Block%s_MaskedRbR",Block.Data()),Form("Block%s_MaskedRbR",Block.Data()),fNMaxColsAbs+1,-0.5,fNMaxColsAbs+0.5, fNMaxRowsAbs+1,-0.5,fNMaxRowsAbs+0.5);
278  plot2D_RbR->GetXaxis()->SetTitle("cell column (#eta direction)");
279  plot2D_RbR->GetYaxis()->SetTitle("cell row (#phi direction)");
280  TH2F *plot2D_Merge = new TH2F(Form("Block%s_MaskedMerge",Block.Data()),Form("Block%s_MaskedMerge",Block.Data()),fNMaxColsAbs+1,-0.5,fNMaxColsAbs+0.5, fNMaxRowsAbs+1,-0.5,fNMaxRowsAbs+0.5);
281  plot2D_Merge->GetXaxis()->SetTitle("cell column (#eta direction)");
282  plot2D_Merge->GetYaxis()->SetTitle("cell row (#phi direction)");
283 
284  Int_t cellColumn=0,cellRow=0;
285  Int_t cellColumnAbs=0,cellRowAbs=0;
286  Int_t trash;
287 
288  for(Int_t i = 0; i < (Int_t)cellVectorRbR.size(); i++)
289  {
290  Int_t cell=cellVectorRbR.at(i);
291  //..Do that only for cell ids also accepted by the code
292  if(!fCaloUtils->GetEMCALGeometry()->CheckAbsCellId(cell))continue;
293 
294  //..Get Row and Collumn for cell ID c
295  fCaloUtils->GetModuleNumberCellIndexesAbsCaloMap(cell,0,cellColumn,cellRow,trash,cellColumnAbs,cellRowAbs);
296  if(cellColumnAbs> fNMaxColsAbs || cellRowAbs>fNMaxRowsAbs)
297  {
298  cout<<"Problem! wrong calculated number of max col and max rows"<<endl;
299  cout<<"current col: "<<cellColumnAbs<<", max col"<<fNMaxColsAbs<<endl;
300  cout<<"current row: "<<cellRowAbs<<", max row"<<fNMaxRowsAbs<<endl;
301  }
302  plot2D_RbR->Fill(cellColumnAbs,cellRowAbs);
303  }
304  for(Int_t i = 0; i < (Int_t)cellVectorMerge.size(); i++)
305  {
306  Int_t cell=cellVectorMerge.at(i);
307  //..Do that only for cell ids also accepted by the code
308  if(!fCaloUtils->GetEMCALGeometry()->CheckAbsCellId(cell))continue;
309 
310  //..Get Row and Collumn for cell ID c
311  fCaloUtils->GetModuleNumberCellIndexesAbsCaloMap(cell,0,cellColumn,cellRow,trash,cellColumnAbs,cellRowAbs);
312  if(cellColumnAbs> fNMaxColsAbs || cellRowAbs>fNMaxRowsAbs)
313  {
314  cout<<"Problem! wrong calculated number of max col and max rows"<<endl;
315  cout<<"current col: "<<cellColumnAbs<<", max col"<<fNMaxColsAbs<<endl;
316  cout<<"current row: "<<cellRowAbs<<", max row"<<fNMaxRowsAbs<<endl;
317  }
318  plot2D_Merge->Fill(cellColumnAbs,cellRowAbs,1);
319  }
320  //. . . . . . . . . . . . . . . . . . . .
321  TCanvas *c1 = new TCanvas(Form("2DMapForBlock%s",Block.Data()),Form("2DMapForBlock%s",Block.Data()),900,500);
322  c1->ToggleEventStatus();
323  c1->Divide(2);
324  c1->cd(1);
325  plot2D_RbR->Draw("colz");
326  c1->cd(2);
327  plot2D_Merge->Draw("colz");
328 }
332 Bool_t IsCellMaskedByHand(Int_t cell, std::vector<Int_t> cellVector)
333 {
334  Bool_t bad=0;
335  for(Int_t i=0; i<(Int_t)cellVector.size();i++)
336  {
337  if(cell==cellVector.at(i))bad=1;
338  }
339 
340  return bad;
341 }
345 void CreateCellCompPDF(TH2F* hAmpIDMasked, std::vector<Int_t> cellVector,TH1* goodCellsMerged, TH1* goodCellsRbR, TString pdfName)
346 {
347  Int_t NoOfCells=cellVector.size();
348  Bool_t firstCanvas=0;
349  TString name;
350  /*TLatex* textA = new TLatex(0.2,0.8,"*test*");
351  textA->SetTextSize(0.08);
352  textA->SetTextColor(1);
353  textA->SetNDC();
354  */
355  for(Int_t cell=0;cell<NoOfCells;cell++)
356  {
357  TString internal_pdfName=pdfName;
358  TCanvas *c1;
359  if(cell%9==0)
360  {
361  c1 = new TCanvas(Form("badcells%i",cell),"badcells",1000,750);
362  if(cellVector.size() > 6) c1->Divide(3,3);
363  else if (cellVector.size() > 3) c1->Divide(3,2);
364  else c1->Divide(3,1);
365  }
366  TH1 *hCell = hAmpIDMasked->ProjectionX(Form("Cell %d",cellVector.at(cell)),cellVector.at(cell)+1,cellVector.at(cell)+1);
367  TH1 *hCell2 = (TH1*)hCell->Clone("hCell2");
368 
369  c1->cd(cell%9 + 1);
370  hCell->Divide(goodCellsRbR);
371  hCell2->Divide(goodCellsMerged);
372 
373  hCell->SetLineColor(kBlue-8);
374  hCell2->SetLineColor(kRed-9);
375  hCell->GetXaxis()->SetTitle("E (GeV)");
376  hCell->GetYaxis()->SetTitle("cell/mean of good");
377  hCell->GetXaxis()->SetRangeUser(0.,10.);
378  hCell->SetLineWidth(1) ;
379  hCell->SetTitle(Form("Cell No. %d",cellVector.at(cell)));
380  hCell->Draw("hist");
381  hCell2->DrawCopy("same hist");
382 
383  //textA->SetTitle(Form("Cell No. %d",cellVector.at(cell)));
384  //textA->Draw();
385 
386  //..page is full or end of loop
387  if(cell%9==8 ||cell == NoOfCells-1)
388  {
389  if(cell == NoOfCells-1)
390  {
391  //internal_pdfName +=")";
392  c1->Print(Form("%s)",pdfName.Data()));
393  }
394  else if(firstCanvas==0)
395  {
396  internal_pdfName +="(";
397  c1->Print(internal_pdfName.Data());
398  firstCanvas=1;
399  }
400  else
401  {
402  c1->Print(internal_pdfName.Data());
403  }
404  delete c1;
405  }
406  }
407 
408 }
413 void SetHisto(TH1 *Histo,TString Xtitel,TString Ytitel,Bool_t longhisto)
414 {
415  Histo->SetStats(0);
416  Histo->SetTitle("");
417  if(longhisto==0)
418  {
419  Histo->GetYaxis()->SetTitleOffset(1.4);
420  Histo->GetXaxis()->SetTitleOffset(1.4);
421  Histo->GetXaxis()->SetLabelSize(0.05);
422  Histo->GetYaxis()->SetLabelSize(0.05);
423  Histo->GetXaxis()->SetTitleSize(0.045);
424  Histo->GetYaxis()->SetTitleSize(0.045);
425  Histo->GetXaxis()->SetNdivisions(505);
426  Histo->GetYaxis()->SetNdivisions(505);
427  }
428 
429  if(longhisto==1)
430  {
431  Histo->GetYaxis()->SetTitleOffset(0.2);
432  Histo->GetXaxis()->SetTitleOffset(1.0);
433  //if(big==1) Histo->GetYaxis()->SetLabelOffset(0.015);
434  //if(big==1) Histo->GetXaxis()->SetLabelOffset(0.015);
435  Histo->GetXaxis()->SetLabelSize(0.07);
436  Histo->GetYaxis()->SetLabelSize(0.07);
437  Histo->GetXaxis()->SetTitleSize(0.08);
438  Histo->GetYaxis()->SetTitleSize(0.08);
439  //Histo->GetXaxis()->CenterTitle();
440  //Histo->GetYaxis()->CenterTitle();
441  Histo->GetXaxis()->SetNdivisions(520);
442  Histo->GetYaxis()->SetNdivisions(10);
443  }
444 
445  Histo->GetXaxis()->SetLabelFont(42);
446  Histo->GetYaxis()->SetLabelFont(42);
447  Histo->GetXaxis()->SetTitleFont(42);
448  Histo->GetYaxis()->SetTitleFont(42);
449  if(Xtitel!="")Histo->GetXaxis()->SetTitle(Xtitel);
450  if(Ytitel!="")Histo->GetYaxis()->SetTitle(Ytitel);
451 
452  Histo->SetLineColor(1);
453  Histo->SetMarkerColor(1);
454  Histo->SetMarkerStyle(20);
455  Histo->SetMarkerSize(0.5);
456 }
Plot2DCells
void Plot2DCells(TString Block, Int_t runNo, std::vector< Int_t > cellVectorRbR, std::vector< Int_t > cellVectorMerge)
Definition: helperMacrosBC.C:262
AliCalorimeterUtils::SetRunNumber
void SetRunNumber(Int_t run)
Definition: AliCalorimeterUtils.h:391
TH2F
Definition: External.C:236
IsCellMaskedByHand
Bool_t IsCellMaskedByHand(Int_t cell, std::vector< Int_t > cellVector)
Definition: helperMacrosBC.C:332
SetHisto
void SetHisto(TH1 *Histo, TString Xtitel, TString Ytitel, Bool_t longhisto)
Definition: helperMacrosBC.C:413
Compare2Blocks
void Compare2Blocks(TString period="LHC15n", Int_t trainNo=603, Int_t versionA=0, Int_t versionB=1)
Definition: helperMacrosBC.C:89
AliCalorimeterUtils::GetEMCALGeometry
AliEMCALGeometry * GetEMCALGeometry() const
Definition: AliCalorimeterUtils.h:104
Int_t
int Int_t
Definition: External.C:63
pdfName
const char * pdfName
Definition: DrawAnaELoss.C:30
Get_RowCollumnID
void Get_RowCollumnID(Int_t runNum=244411, Int_t inputCellID=-1, Int_t inputRow=-1, Int_t inputCollumn=-1, Int_t inputSM=-1)
Definition: helperMacrosBC.C:46
Bool_t
bool Bool_t
Definition: External.C:53
AliCalorimeterUtils
Class with utils specific to calorimeter clusters/cells.
Definition: AliCalorimeterUtils.h:44
AliCalorimeterUtils::AccessGeometry
void AccessGeometry(AliVEvent *inputEvent)
Definition: AliCalorimeterUtils.cxx:439
AliCalorimeterUtils.h
TH1
Definition: External.C:196
AliCalorimeterUtils::GetModuleNumberCellIndexesAbsCaloMap
Int_t GetModuleNumberCellIndexesAbsCaloMap(Int_t absId, Int_t calo, Int_t &icol, Int_t &irow, Int_t &iRCU, Int_t &icolAbs, Int_t &irowAbs) const
Definition: AliCalorimeterUtils.cxx:1343
CreateCellCompPDF
void CreateCellCompPDF(TH2F *hAmpIDMasked, std::vector< Int_t > cellVector, TH1 *goodCellsMerged, TH1 *goodCellsRbR, TString pdfName)
Definition: helperMacrosBC.C:345