AliPhysics  vAN-20151012 (2287573)
 All Classes Namespaces Files Functions Variables Enumerations Enumerator Macros
AliAnalysisTaskEMCALPi0CalibSelection.cxx
Go to the documentation of this file.
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 
4  * Permission to use, copy, modify and distribute this software and its *
5  * documentation strictly for non-commercial purposes is hereby granted *
6  * without fee, provided that the above copyright notice appears in all *
7  * copies and that both the copyright notice and this permission notice *
8  * appear in the supporting documentation. The authors make no claims *
9  * about the suitability of this software for any purpose. It is *
10  * provided "as is" without express or implied warranty. *
11  **************************************************************************/
12 
13 // Root
14 #include <TRefArray.h>
15 #include <TList.h>
16 #include <TH1F.h>
17 #include <TGeoManager.h>
18 #include <TFile.h>
19 
20 // AliRoot
22 #include "AliAODEvent.h"
23 #include "AliESDEvent.h"
24 #include "AliEMCALGeometry.h"
25 #include "AliVCluster.h"
26 #include "AliVCaloCells.h"
27 #include "AliEMCALRecoUtils.h"
28 #include "AliOADBContainer.h"
29 
33 
36 //______________________________________________________________________________________________
38 AliAnalysisTaskSE(),
39 fEMCALGeo(0x0), fLoadMatrices(0),
40 fEMCALGeoName("EMCAL_COMPLETE12SMV1_DCAL_8SM"),
41 fTriggerName("EMC"),
42 fRecoUtils(new AliEMCALRecoUtils),
43 fOADBFilePath(""), fCalibFilePath(""),
44 fCorrectClusters(kFALSE), fRecalPosition(kTRUE),
45 fCaloClustersArr(0x0), fEMCALCells(0x0),
46 //fCuts(0x0),
47 fOutputContainer(0x0),
48 fVertex(), fFilteredInput(kFALSE),
49 fImportGeometryFromFile(1), fImportGeometryFilePath(""),
50 fEmin(0.5), fEmax(15.),
51 fL0min(0.01), fL0max(0.5),
52 fDTimeCut(100.), fTimeMax(1000000), fTimeMin(-1000000),
53 fAsyCut(1.), fMinNCells(2), fGroupNCells(0),
54 fLogWeight(4.5), fSameSM(kFALSE),
55 fNMaskCellColumns(11), fMaskCellColumns(0x0),
56 fInvMassCutMin(110.), fInvMassCutMax(160.),
57 // Histograms binning
58 fNbins(300),
59 fMinBin(0.), fMaxBin(300.),
60 fNTimeBins(1000), fMinTimeBin(0.), fMaxTimeBin(1000.),
61 // Temporal
62 fMomentum1(), fMomentum2(), fMomentum12(),
63 // Histograms
64 fHmgg(0x0), fHmggDifferentSM(0x0),
65 fHmggMaskFrame(0x0), fHmggDifferentSMMaskFrame(0x0),
66 fHOpeningAngle(0x0), fHOpeningAngleDifferentSM(0x0),
67 fHAsymmetry(0x0), fHAsymmetryDifferentSM(0x0),
68 fhNEvents(0x0),
69 fhClusterTime(0x0), fhClusterPairDiffTime(0x0)
70 {
71  for(Int_t iMod=0; iMod < AliEMCALGeoParams::fgkEMCALModules; iMod++)
72  {
73  for(Int_t iX=0; iX<24; iX++)
74  {
75  for(Int_t iZ=0; iZ<48; iZ++)
76  {
77  fHmpi0[iMod][iZ][iX] = 0 ;
78  }
79  }
80  }
81 
82  fVertex[0]=fVertex[1]=fVertex[2]=-1000;
83 
84  fHTpi0[0]= 0 ;
85  fHTpi0[1]= 0 ;
86  fHTpi0[2]= 0 ;
87  fHTpi0[3]= 0 ;
88 
90  fMaskCellColumns[0] = 6 ; fMaskCellColumns[1] = 7 ; fMaskCellColumns[2] = 8 ;
91  fMaskCellColumns[3] = 35; fMaskCellColumns[4] = 36; fMaskCellColumns[5] = 37;
92  fMaskCellColumns[6] = 12+AliEMCALGeoParams::fgkEMCALCols; fMaskCellColumns[7] = 13+AliEMCALGeoParams::fgkEMCALCols;
93  fMaskCellColumns[8] = 40+AliEMCALGeoParams::fgkEMCALCols; fMaskCellColumns[9] = 41+AliEMCALGeoParams::fgkEMCALCols;
94  fMaskCellColumns[10]= 42+AliEMCALGeoParams::fgkEMCALCols;
95 
96  for(Int_t iSMPair = 0; iSMPair < AliEMCALGeoParams::fgkEMCALModules/2; iSMPair++)
97  {
98  fHmggPairSameSectorSM[iSMPair] = 0;
99  fHmggPairSameSectorSMMaskFrame[iSMPair] = 0;
101  }
102 
103  for(Int_t iSMPair = 0; iSMPair < AliEMCALGeoParams::fgkEMCALModules-2; iSMPair++)
104  {
105  fHmggPairSameSideSM[iSMPair] = 0;
106  fHmggPairSameSideSMMaskFrame[iSMPair] = 0;
107  fhClusterPairDiffTimeSameSide[iSMPair] = 0;
108  }
109 
110  for(Int_t iSM = 0; iSM < AliEMCALGeoParams::fgkEMCALModules; iSM++)
111  {
112  fHmggSM[iSM] = 0;
113  fHmggSMMaskFrame[iSM] = 0;
114  fHOpeningAngleSM[iSM] = 0;
115  fHOpeningAnglePairSM[iSM] = 0;
116  fHAsymmetrySM[iSM] = 0;
117  fHAsymmetryPairSM[iSM] = 0;
118  fhTowerDecayPhotonHit[iSM] = 0;
119  fhTowerDecayPhotonEnergy[iSM] = 0;
122  fMatrix[iSM] = 0x0;
123  fhClusterTimeSM[iSM] = 0;
125  }
126 }
127 
133 //______________________________________________________________________________________________
135 AliAnalysisTaskSE(name),
136 fEMCALGeo(0x0), fLoadMatrices(0),
137 fEMCALGeoName("EMCAL_COMPLETE12SMV1_DCAL_8SM"),
138 fTriggerName("EMC"),
139 fRecoUtils(new AliEMCALRecoUtils),
140 fOADBFilePath(""), fCalibFilePath(""),
141 fCorrectClusters(kFALSE), fRecalPosition(kTRUE),
142 fCaloClustersArr(0x0), fEMCALCells(0x0),
143 //fCuts(0x0),
144 fOutputContainer(0x0),
145 fVertex(), fFilteredInput(kFALSE),
146 fImportGeometryFromFile(1), fImportGeometryFilePath(""),
147 fEmin(0.5), fEmax(15.),
148 fL0min(0.01), fL0max(0.5),
149 fDTimeCut(100.), fTimeMax(1000000), fTimeMin(-1000000),
150 fAsyCut(1.), fMinNCells(2), fGroupNCells(0),
151 fLogWeight(4.5), fSameSM(kFALSE),
152 fNMaskCellColumns(11), fMaskCellColumns(0x0),
153 fInvMassCutMin(110.), fInvMassCutMax(160.),
154 // Histograms binning
155 fNbins(300),
156 fMinBin(0.), fMaxBin(300.),
157 fNTimeBins(1000), fMinTimeBin(0.), fMaxTimeBin(1000.),
158 // Temporal
159 fMomentum1(), fMomentum2(), fMomentum12(),
160 // Histograms
161 fHmgg(0x0), fHmggDifferentSM(0x0),
162 fHmggMaskFrame(0x0), fHmggDifferentSMMaskFrame(0x0),
163 fHOpeningAngle(0x0), fHOpeningAngleDifferentSM(0x0),
164 fHAsymmetry(0x0), fHAsymmetryDifferentSM(0x0),
165 fhNEvents(0x0),
166 fhClusterTime(0x0), fhClusterPairDiffTime(0x0)
167 {
168  for(Int_t iMod=0; iMod < AliEMCALGeoParams::fgkEMCALModules; iMod++)
169  {
170  for(Int_t iX=0; iX<24; iX++)
171  {
172  for(Int_t iZ=0; iZ<48; iZ++)
173  {
174  fHmpi0[iMod][iZ][iX] = 0 ;
175  }
176  }
177  }
178 
179  fVertex[0]=fVertex[1]=fVertex[2]=-1000;
180 
181  fHTpi0[0]= 0 ;
182  fHTpi0[1]= 0 ;
183  fHTpi0[2]= 0 ;
184  fHTpi0[3]= 0 ;
185 
186  fMaskCellColumns = new Int_t[fNMaskCellColumns];
187  fMaskCellColumns[0] = 6 ; fMaskCellColumns[1] = 7 ; fMaskCellColumns[2] = 8 ;
188  fMaskCellColumns[3] = 35; fMaskCellColumns[4] = 36; fMaskCellColumns[5] = 37;
189  fMaskCellColumns[6] = 12+AliEMCALGeoParams::fgkEMCALCols; fMaskCellColumns[7] = 13+AliEMCALGeoParams::fgkEMCALCols;
190  fMaskCellColumns[8] = 40+AliEMCALGeoParams::fgkEMCALCols; fMaskCellColumns[9] = 41+AliEMCALGeoParams::fgkEMCALCols;
191  fMaskCellColumns[10]= 42+AliEMCALGeoParams::fgkEMCALCols;
192 
193  for(Int_t iSMPair = 0; iSMPair < AliEMCALGeoParams::fgkEMCALModules/2; iSMPair++)
194  {
195  fHmggPairSameSectorSM[iSMPair] = 0;
196  fHmggPairSameSectorSMMaskFrame[iSMPair] = 0;
198  }
199 
200  for(Int_t iSMPair = 0; iSMPair < AliEMCALGeoParams::fgkEMCALModules-2; iSMPair++)
201  {
202  fHmggPairSameSideSM[iSMPair] = 0;
203  fHmggPairSameSideSMMaskFrame[iSMPair] = 0;
204  fhClusterPairDiffTimeSameSide[iSMPair] = 0;
205  }
206 
207  for(Int_t iSM = 0; iSM < AliEMCALGeoParams::fgkEMCALModules; iSM++)
208  {
209  fHmggSM[iSM] = 0;
210  fHmggSMMaskFrame[iSM] = 0;
211  fHOpeningAngleSM[iSM] = 0;
212  fHOpeningAnglePairSM[iSM] = 0;
213  fHAsymmetrySM[iSM] = 0;
214  fHAsymmetryPairSM[iSM] = 0;
215  fhTowerDecayPhotonHit[iSM] = 0;
216  fhTowerDecayPhotonEnergy[iSM] = 0;
219  fMatrix[iSM] = 0x0;
220  fhClusterTimeSM[iSM] = 0;
222  }
223 
224  DefineOutput(1, TList::Class());
225 //DefineOutput(2, TList::Class()); // will contain cuts or local params
226 }
227 
230 //_____________________________________________________________________________
232 {
233  if(fOutputContainer)
234  {
235  fOutputContainer->Delete() ;
236  delete fOutputContainer ;
237  }
238 
239  if(fEMCALGeo) delete fEMCALGeo ;
240  if(fRecoUtils) delete fRecoUtils ;
241  if(fNMaskCellColumns) delete [] fMaskCellColumns;
242 }
243 
247 //____________________________________________________________
249 {
250  if(fRecoUtils->GetParticleType()!=AliEMCALRecoUtils::kPhoton)
251  AliFatal(Form("Wrong particle type for cluster position recalculation! = %d\n", fRecoUtils->GetParticleType()));
252 
253  AliDebug(1,Form("It will use fLogWeight %.3f",fLogWeight));
254 
255  Float_t pos[]={0,0,0};
256 
257  for(Int_t iClu=0; iClu < fCaloClustersArr->GetEntriesFast(); iClu++)
258  {
259  AliVCluster *c1 = (AliVCluster *) fCaloClustersArr->At(iClu);
260 
261  Float_t e1i = c1->E(); // cluster energy before correction
262  if (e1i < fEmin) continue;
263  else if (e1i > fEmax) continue;
264 
265  else if (c1->GetNCells() < fMinNCells) continue;
266 
267  else if (c1->GetM02() < fL0min || c1->GetM02() > fL0max) continue;
268 
269  if(fRecoUtils->ClusterContainsBadChannel(fEMCALGeo, c1->GetCellsAbsId(), c1->GetNCells())) continue;
270 
271  if(DebugLevel() > 2)
272  {
273  AliInfo(Form("Std : i %d, E %f, dispersion %f, m02 %f, m20 %f\n",c1->GetID(),c1->E(),c1->GetDispersion(),c1->GetM02(),c1->GetM20()));
274  c1->GetPosition(pos);
275  AliInfo(Form("Std : i %d, x %f, y %f, z %f\n",c1->GetID(), pos[0], pos[1], pos[2]));
276  }
277 
278  // Correct cluster energy and position if requested, and not corrected previously, by default Off
279  if(fRecoUtils->IsRecalibrationOn())
280  {
281  fRecoUtils->RecalibrateClusterEnergy(fEMCALGeo, c1, fEMCALCells);
282  fRecoUtils->RecalculateClusterShowerShapeParameters(fEMCALGeo, fEMCALCells,c1);
283  fRecoUtils->RecalculateClusterPID(c1);
284  }
285 
286  AliDebug(2,Form("Energy: after recalibration %f",c1->E()));
287 
288  // Recalculate cluster position
289  if ( fRecalPosition ) fRecoUtils->RecalculateClusterPosition(fEMCALGeo, fEMCALCells,c1);
290 
291  // Correct Non-Linearity
292  c1->SetE(fRecoUtils->CorrectClusterEnergyLinearity(c1));
293 
294  AliDebug(2,Form("after linearity correction %f",c1->E()));
295 
296  // In case of MC analysis, to match resolution/calibration in real data
297  //c1->SetE(fRecoUtils->SmearClusterEnergy(c1)); // Not needed anymore
298 
299  AliDebug(2,Form("after smearing %f\n",c1->E()));
300 
301  if(DebugLevel() > 2)
302  {
303  AliInfo(Form("Cor : i %d, E %f, dispersion %f, m02 %f, m20 %f\n",c1->GetID(),c1->E(),c1->GetDispersion(),c1->GetM02(),c1->GetM20()));
304  c1->GetPosition(pos);
305  AliInfo(Form("Cor : i %d, x %f, y %f, z %f\n",c1->GetID(), pos[0], pos[1], pos[2]));
306  }
307  } // cluster loop
308 }
309 
313 //__________________________________________________________
315 {
316  Int_t absId1 = -1;
317  Int_t iSupMod1 = -1;
318  Int_t iphi1 = -1;
319  Int_t ieta1 = -1;
320  Int_t absId2 = -1;
321  Int_t iSupMod2 = -1;
322  Int_t iphi2 = -1;
323  Int_t ieta2 = -1;
324  Bool_t shared = kFALSE;
325 
326  Float_t pos[] = {0,0,0};
327 
328  Int_t bc = InputEvent()->GetBunchCrossNumber();
329  Int_t nSM = (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules();
330 
331  for(Int_t iClu=0; iClu<fCaloClustersArr->GetEntriesFast()-1; iClu++)
332  {
333  AliVCluster *c1 = (AliVCluster *) fCaloClustersArr->At(iClu);
334 
335  // Exclude bad channels
336  if(fRecoUtils->ClusterContainsBadChannel(fEMCALGeo, c1->GetCellsAbsId(), c1->GetNCells())) continue;
337 
338  Float_t e1i = c1->E(); // cluster energy before correction
339 
340  if (e1i < fEmin) continue;
341  else if (e1i > fEmax) continue;
342 
343  else if (!fRecoUtils->IsGoodCluster(c1,fEMCALGeo,fEMCALCells,bc)) continue;
344 
345  else if (c1->GetNCells() < fMinNCells) continue;
346 
347  else if (c1->GetM02() < fL0min || c1->GetM02() > fL0max) continue;
348 
349  if(DebugLevel() > 2)
350  {
351  AliInfo(Form("IMA : i %d, E %f, dispersion %f, m02 %f, m20 %f",c1->GetID(),e1i,c1->GetDispersion(),c1->GetM02(),c1->GetM20()));
352  c1->GetPosition(pos);
353  AliInfo(Form("IMA : i %d, x %f, y %f, z %f",c1->GetID(), pos[0], pos[1], pos[2]));
354  }
355 
356  fRecoUtils->GetMaxEnergyCell(fEMCALGeo, fEMCALCells,c1,absId1,iSupMod1,ieta1,iphi1,shared);
357 
358  c1->GetMomentum(fMomentum1,fVertex);
359 
360  // Check if cluster is in fidutial region, not too close to borders
361  Bool_t in1 = fRecoUtils->CheckCellFiducialRegion(fEMCALGeo, c1, fEMCALCells);
362 
363  // Clusters not facing frame structures
364  Bool_t mask1 = MaskFrameCluster(iSupMod1, ieta1);
365  //if(mask1) printf("Reject eta %d SM %d\n",ieta1, iSupMod1);
366 
367  Double_t time1 = c1->GetTOF()*1.e9;
368 
369  if(time1 > fTimeMax || time1 < fTimeMin) continue;
370 
371  fhClusterTime ->Fill(c1->E(),time1);
372  fhClusterTimeSM[iSupMod1]->Fill(c1->E(),time1);
373 
374  // Combine cluster with other clusters and get the invariant mass
375  for (Int_t jClu=iClu+1; jClu < fCaloClustersArr->GetEntriesFast(); jClu++)
376  {
377  AliAODCaloCluster *c2 = (AliAODCaloCluster *) fCaloClustersArr->At(jClu);
378 
379  Float_t e2i = c2->E();
380  if (e2i < fEmin) continue;
381  else if (e2i > fEmax) continue;
382 
383  else if (!fRecoUtils->IsGoodCluster(c2,fEMCALGeo,fEMCALCells,bc))continue;
384 
385  else if (c2->GetNCells() < fMinNCells) continue;
386 
387  else if (c2->GetM02() < fL0min || c2->GetM02() > fL0max) continue;
388 
389  fRecoUtils->GetMaxEnergyCell(fEMCALGeo, fEMCALCells,c2,absId2,iSupMod2,ieta2,iphi2,shared);
390 
391  c2->GetMomentum(fMomentum2,fVertex);
392 
394  Float_t invmass = fMomentum12.M()*1000;
395 
396  //Asimetry cut
397  Float_t asym = TMath::Abs(fMomentum1.E()-fMomentum2.E())/(fMomentum1.E()+fMomentum2.E());
398 
399  if(asym > fAsyCut) continue;
400 
401  //Time cut
402  Double_t time2 = c2->GetTOF()*1.e9;
403 
404  if(time2 > fTimeMax || time2 < fTimeMin) continue;
405 
406  fhClusterPairDiffTime->Fill(fMomentum12.E(),time1-time2);
407  if(TMath::Abs(time1-time2) > fDTimeCut) continue;
408 
409  if(invmass < fMaxBin && invmass > fMinBin )
410  {
411  //Check if cluster is in fidutial region, not too close to borders
412  Bool_t in2 = fRecoUtils->CheckCellFiducialRegion(fEMCALGeo, c2, fEMCALCells);
413 
414  // Clusters not facing frame structures
415  Bool_t mask2 = MaskFrameCluster(iSupMod2, ieta2);
416  //if(mask2) printf("Reject eta %d SM %d\n",ieta2, iSupMod2);
417 
418  if(in1 && in2)
419  {
420  fHmgg->Fill(invmass,fMomentum12.Pt());
421 
422  if(iSupMod1==iSupMod2)
423  {
424  fHmggSM[iSupMod1]->Fill(invmass,fMomentum12.Pt());
425  fhClusterPairDiffTimeSameSM[iSupMod1]->Fill(fMomentum12.E(),time1-time2);
426  }
427  else
428  fHmggDifferentSM ->Fill(invmass,fMomentum12.Pt());
429 
430  // Same sector
431  Int_t j=0;
432  for(Int_t i = 0; i < nSM/2; i++)
433  {
434  j=2*i;
435  if((iSupMod1==j && iSupMod2==j+1) || (iSupMod1==j+1 && iSupMod2==j))
436  {
437  fHmggPairSameSectorSM[i]->Fill(invmass,fMomentum12.Pt());
438  fhClusterPairDiffTimeSameSector[i]->Fill(fMomentum12.E(),time1-time2);
439  }
440  }
441 
442  // Same side
443  for(Int_t i = 0; i < nSM-2; i++)
444  {
445  if((iSupMod1==i && iSupMod2==i+2) || (iSupMod1==i+2 && iSupMod2==i))
446  {
447  fHmggPairSameSideSM[i]->Fill(invmass,fMomentum12.Pt());
448  fhClusterPairDiffTimeSameSide[i]->Fill(fMomentum12.E(),time1-time2);
449  }
450  }
451 
452 
453  if(!mask1 && !mask2)
454  {
455  fHmggMaskFrame->Fill(invmass,fMomentum12.Pt());
456 
457  if(iSupMod1==iSupMod2) fHmggSMMaskFrame[iSupMod1]->Fill(invmass,fMomentum12.Pt());
458  else fHmggDifferentSMMaskFrame ->Fill(invmass,fMomentum12.Pt());
459 
460  // Same sector
461  j=0;
462  for(Int_t i = 0; i < nSM/2; i++)
463  {
464  j=2*i;
465  if((iSupMod1==j && iSupMod2==j+1) || (iSupMod1==j+1 && iSupMod2==j)) fHmggPairSameSectorSMMaskFrame[i]->Fill(invmass,fMomentum12.Pt());
466  }
467 
468  // Same side
469  for(Int_t i = 0; i < nSM-2; i++)
470  {
471  if((iSupMod1==i && iSupMod2==i+2) || (iSupMod1==i+2 && iSupMod2==i)) fHmggPairSameSideSMMaskFrame[i]->Fill(invmass,fMomentum12.Pt());
472  }
473 
474  }// Pair not facing frame
475 
476  if(invmass > fInvMassCutMin && invmass < fInvMassCutMax) //restrict to clusters really close to pi0 peak
477  {
478 
479  // Check time of cells in both clusters, and fill time histogram
480  for(Int_t icell = 0; icell < c1->GetNCells(); icell++)
481  {
482  Int_t absID = c1->GetCellAbsId(icell);
483  fHTpi0[bc%4]->Fill(absID, fEMCALCells->GetCellTime(absID)*1.e9);
484  }
485 
486  for(Int_t icell = 0; icell < c2->GetNCells(); icell++)
487  {
488  Int_t absID = c2->GetCellAbsId(icell);
489  fHTpi0[bc%4]->Fill(absID, fEMCALCells->GetCellTime(absID)*1.e9);
490  }
491 
492  //Opening angle of 2 photons
493  Float_t opangle = fMomentum1.Angle(fMomentum2.Vect())*TMath::RadToDeg();
494  //printf("*******>>>>>>>> In PEAK pt %f, angle %f \n",fMomentum12.Pt(),opangle);
495 
496 
497  fHOpeningAngle ->Fill(opangle,fMomentum12.Pt());
498  fHAsymmetry ->Fill(asym,fMomentum12.Pt());
499 
500  if(iSupMod1==iSupMod2)
501  {
502  fHOpeningAngleSM[iSupMod1] ->Fill(opangle,fMomentum12.Pt());
503  fHAsymmetrySM[iSupMod1] ->Fill(asym,fMomentum12.Pt());
504  }
505  else
506  {
507  fHOpeningAngleDifferentSM ->Fill(opangle,fMomentum12.Pt());
508  fHAsymmetryDifferentSM ->Fill(asym,fMomentum12.Pt());
509  }
510 
511  if((iSupMod1==0 && iSupMod2==2) || (iSupMod1==2 && iSupMod2==0))
512  {
513  fHOpeningAnglePairSM[0] ->Fill(opangle,fMomentum12.Pt());
514  fHAsymmetryPairSM[0] ->Fill(asym,fMomentum12.Pt());
515 
516  }
517  if((iSupMod1==1 && iSupMod2==3) || (iSupMod1==3 && iSupMod2==1))
518  {
519  fHOpeningAnglePairSM[1] ->Fill(opangle,fMomentum12.Pt());
520  fHAsymmetryPairSM[1] ->Fill(asym,fMomentum12.Pt());
521  }
522 
523  if((iSupMod1==0 && iSupMod2==1) || (iSupMod1==1 && iSupMod2==0))
524  {
525  fHOpeningAnglePairSM[2] ->Fill(opangle,fMomentum12.Pt());
526  fHAsymmetryPairSM[2] ->Fill(asym,fMomentum12.Pt());
527  }
528  if((iSupMod1==2 && iSupMod2==3) || (iSupMod1==3 && iSupMod2==2))
529  {
530  fHOpeningAnglePairSM[3] ->Fill(opangle,fMomentum12.Pt());
531  fHAsymmetryPairSM[3] ->Fill(asym,fMomentum12.Pt());
532  }
533 
534  }// pair in 100 < mass < 160
535 
536  }//in acceptance cuts
537 
538  //In case of filling only channels with second cluster in same SM
539  if(fSameSM && iSupMod1!=iSupMod2) continue;
540 
541  if (fGroupNCells == 0)
542  {
543  fHmpi0[iSupMod1][ieta1][iphi1]->Fill(invmass);
544  fHmpi0[iSupMod2][ieta2][iphi2]->Fill(invmass);
545 
546  if(invmass > fInvMassCutMin && invmass < fInvMassCutMax)//restrict to clusters really close to pi0 peak
547  {
548  fhTowerDecayPhotonHit [iSupMod1]->Fill(ieta1,iphi1);
549  fhTowerDecayPhotonEnergy [iSupMod1]->Fill(ieta1,iphi1,fMomentum1.E());
550  fhTowerDecayPhotonAsymmetry[iSupMod1]->Fill(ieta1,iphi1,asym);
551 
552  fhTowerDecayPhotonHit [iSupMod2]->Fill(ieta2,iphi2);
553  fhTowerDecayPhotonEnergy [iSupMod2]->Fill(ieta2,iphi2,fMomentum2.E());
554  fhTowerDecayPhotonAsymmetry[iSupMod2]->Fill(ieta2,iphi2,asym);
555 
556  if(!mask1)fhTowerDecayPhotonHitMaskFrame[iSupMod1]->Fill(ieta1,iphi1);
557  if(!mask2)fhTowerDecayPhotonHitMaskFrame[iSupMod2]->Fill(ieta2,iphi2);
558 
559  }// pair in mass of pi0
560  }
561  else
562  {
563  //printf("Regroup N %d, eta1 %d, phi1 %d, eta2 %d, phi2 %d \n",fGroupNCells, ieta1, iphi1, ieta2, iphi2);
564  for (Int_t i = -fGroupNCells; i < fGroupNCells+1; i++)
565  {
566  for (Int_t j = -fGroupNCells; j < fGroupNCells+1; j++)
567  {
568  Int_t absId11 = fEMCALGeo->GetAbsCellIdFromCellIndexes(iSupMod1, iphi1+j, ieta1+i);
569  Int_t absId22 = fEMCALGeo->GetAbsCellIdFromCellIndexes(iSupMod2, iphi2+j, ieta2+i);
570 
571  Bool_t ok1 = kFALSE;
572  Bool_t ok2 = kFALSE;
573 
574  for(Int_t icell = 0; icell < c1->GetNCells(); icell++)
575  {
576  if(c1->GetCellsAbsId()[icell] == absId11) ok1=kTRUE;
577  }
578 
579  for(Int_t icell = 0; icell < c2->GetNCells(); icell++)
580  {
581  if(c2->GetCellsAbsId()[icell] == absId22) ok2=kTRUE;
582  }
583 
584  if(ok1 && (ieta1+i >= 0) && (iphi1+j >= 0) && (ieta1+i < 48) && (iphi1+j < 24))
585  {
586  fHmpi0[iSupMod1][ieta1+i][iphi1+j]->Fill(invmass);
587  }
588 
589  if(ok2 && (ieta2+i >= 0) && (iphi2+j >= 0) && (ieta2+i < 48) && (iphi2+j < 24))
590  {
591  fHmpi0[iSupMod2][ieta2+i][iphi2+j]->Fill(invmass);
592  }
593  }// j loop
594  }//i loop
595  }//group cells
596 
597  AliDebug(1,Form("Mass in (SM%d,%d,%d) and (SM%d,%d,%d): %.3f GeV E1_i=%f E1_ii=%f E2_i=%f E2_ii=%f\n",
598  iSupMod1,iphi1,ieta1,iSupMod2,iphi2,ieta2,fMomentum12.M(),e1i,c1->E(),e2i,c2->E()));
599  }
600  }
601  } // end of loop over EMCAL clusters
602 }
603 
604 
609 //______________________________________________________________________
611 {
612  if ( !fRecoUtils->IsRecalibrationOn() || fCalibFilePath == "" ) return ;
613 
614  TFile * calibFactorsFile = TFile::Open(fCalibFilePath);
615 
616  if ( !calibFactorsFile ) AliFatal("Cannot recover the calibration factors");
617 
618  for(Int_t ism = 0; ism < fEMCALGeo->GetNumberOfSuperModules(); ism++)
619  {
620  TH2F * histo = (TH2F*) calibFactorsFile->Get(Form("EMCALRecalFactors_SM%d",ism));
621 
622  if ( histo )
623  fRecoUtils->SetEMCALChannelRecalibrationFactors(ism,histo);
624  else
625  AliWarning(Form("Null histogram with calibration factors for SM%d, 1 will be used for the full SM!",ism));
626  }
627 }
628 
632 //________________________________________________________________
634 {
635  Int_t runnumber = InputEvent()->GetRunNumber() ;
636 
637  //
638  // Load default geo matrices if requested
639  if(fImportGeometryFromFile && !gGeoManager)
640  {
641  if(fImportGeometryFilePath=="") // If not specified, set location depending on run number
642  {
643  // "$ALICE_ROOT/EVE/alice-data/default_geo.root"
644  if (runnumber < 140000) fImportGeometryFilePath = "$ALICE_PHYSICS/OADB/EMCAL/geometry_2010.root";
645  else if(runnumber < 171000) fImportGeometryFilePath = "$ALICE_PHYSICS/OADB/EMCAL/geometry_2011.root";
646  else if(runnumber < 198000) fImportGeometryFilePath = "$ALICE_PHYSICS/OADB/EMCAL/geometry_2012.root"; // 2012-2013
647  else fImportGeometryFilePath = "$ALICE_PHYSICS/OADB/EMCAL/geometry_2015.root"; // >= 2015
648  }
649 
650  AliInfo(Form("Import %s",fImportGeometryFilePath.Data()));
651 
652  TGeoManager::Import(fImportGeometryFilePath) ; // default need file "geometry.root" in local dir!!!!
653  }
654 
655  //
656  if(fLoadMatrices)
657  {
658  AliInfo("Load user defined EMCAL geometry matrices");
659  // OADB if available
660  AliOADBContainer emcGeoMat("AliEMCALgeo");
661 
662  if(fOADBFilePath=="") fOADBFilePath = "$ALICE_PHYSICS/OADB/EMCAL" ;
663 
664  emcGeoMat.InitFromFile(Form("%s/EMCALlocal2master.root",fOADBFilePath.Data()),"AliEMCALgeo");
665 
666  TObjArray *matEMCAL=(TObjArray*)emcGeoMat.GetObject(runnumber,"EmcalMatrices");
667 
668  for(Int_t mod = 0; mod < (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules(); mod++)
669  {
670  if (!fMatrix[mod]) // Get it from OADB
671  {
672  AliDebug(1,Form("EMCAL matrices SM %d, %p",mod,((TGeoHMatrix*) matEMCAL->At(mod))));
673  //((TGeoHMatrix*) matEMCAL->At(mod))->Print();
674 
675  fMatrix[mod] = (TGeoHMatrix*) matEMCAL->At(mod) ;
676  }
677 
678  if(fMatrix[mod])
679  {
680  if(DebugLevel() > 1)
681  fMatrix[mod]->Print();
682 
683  fEMCALGeo->SetMisalMatrix(fMatrix[mod],mod) ;
684  }
685  else if(gGeoManager)
686  {
687  AliWarning(Form("Set matrix for SM %d from gGeoManager",mod));
688  fEMCALGeo->SetMisalMatrix(fEMCALGeo->GetMatrixForSuperModuleFromGeoManager(mod),mod) ;
689  }
690  else
691  {
692  AliError(Form("Alignment matrix for SM %d is not available",mod));
693  }
694  }//SM loop
695  }//Load matrices
696  else if(!gGeoManager)
697  {
698  AliInfo("Get geo matrices from data");
699  //Still not implemented in AOD, just a workaround to be able to work at least with ESDs
700  if(!strcmp(InputEvent()->GetName(),"AliAODEvent"))
701  {
702  AliWarning("Use ideal geometry, values geometry matrix not kept in AODs");
703  }//AOD
704  else
705  {
706  AliDebug(1,"AliAnalysisTaskEMCALClusterize Load Misaligned matrices");
707 
708  for(Int_t mod=0; mod < (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules(); mod++)
709  {
710  if(InputEvent()->GetEMCALMatrix(mod))
711  {
712  if(DebugLevel() > 1)
713  InputEvent()->GetEMCALMatrix(mod)->Print();
714 
715  fEMCALGeo->SetMisalMatrix(InputEvent()->GetEMCALMatrix(mod),mod) ;
716  }
717 
718  }
719  }// ESD
720  }// Load matrices from Data
721  else if(gGeoManager) // Load default matrices
722  {
723  for(Int_t mod = 0; mod < (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules(); mod++)
724  {
725  AliWarning(Form("Set matrix for SM %d from gGeoManager",mod));
726  fEMCALGeo->SetMisalMatrix(fEMCALGeo->GetMatrixForSuperModuleFromGeoManager(mod),mod) ;
727  }
728  } // gGeoManager matrices
729 
730 }
731 
735 //______________________________________________________________________
737 {
738  if(!fRecoUtils->IsRunDepRecalibrationOn()) return;
739 
740  AliOADBContainer *contRFTD=new AliOADBContainer("");
741 
742  contRFTD->InitFromFile(Form("%s/EMCALTemperatureCorrCalib.root",fOADBFilePath.Data()),"AliEMCALRunDepTempCalibCorrections");
743 
744  Int_t runnumber = InputEvent()->GetRunNumber() ;
745 
746  TH1S *htd=(TH1S*)contRFTD->GetObject(runnumber);
747 
748  //If it did not exist for this run, get closes one
749  if (!htd)
750  {
751  AliWarning(Form("No TemperatureCorrCalib Objects for run: %d",runnumber));
752 
753  // let's get the closest runnumber instead then..
754  Int_t lower = 0;
755  Int_t ic = 0;
756  Int_t maxEntry = contRFTD->GetNumberOfEntries();
757 
758  while ( (ic < maxEntry) && (contRFTD->UpperLimit(ic) < runnumber) )
759  {
760  lower = ic;
761  ic++;
762  }
763 
764  Int_t closest = lower;
765  if ( (ic<maxEntry) &&
766  (contRFTD->LowerLimit(ic)-runnumber) < (runnumber - contRFTD->UpperLimit(lower)) )
767  {
768  closest = ic;
769  }
770 
771  AliWarning(Form("TemperatureCorrCalib Objects found closest id %d from run: %d",
772  closest, contRFTD->LowerLimit(closest)));
773 
774  htd = (TH1S*) contRFTD->GetObjectByIndex(closest);
775  }
776 
777  // Fill parameters
778  if(htd)
779  {
780  AliInfo("Recalibrate (Temperature) EMCAL");
781 
782  Int_t nSM = fEMCALGeo->GetNumberOfSuperModules();
783 
784  for (Int_t ism = 0; ism < nSM; ++ism)
785  {
786  for (Int_t icol = 0; icol < 48; ++icol)
787  {
788  for (Int_t irow = 0; irow < 24; ++irow)
789  {
790  Float_t factor = fRecoUtils->GetEMCALChannelRecalibrationFactor(ism,icol,irow);
791 
792  Int_t absID = fEMCALGeo->GetAbsCellIdFromCellIndexes(ism, irow, icol); // original calibration factor
793 
794  AliDebug(3,Form(" ism %d, icol %d, irow %d,absID %d - Calib factor %1.5f - ",ism, icol, irow, absID, factor));
795 
796  factor *= htd->GetBinContent(absID) / 10000. ; // correction dependent on T
797 
798  fRecoUtils->SetEMCALChannelRecalibrationFactor(ism,icol,irow,factor);
799 
800  AliDebug(3,Form("T factor %1.5f - final factor %1.5f",
801  htd->GetBinContent(absID) / 10000.,
802  fRecoUtils->GetEMCALChannelRecalibrationFactor(ism,icol,irow)));
803  } // columns
804  } // rows
805  } // SM loop
806  }
807  else AliInfo("Do NOT recalibrate EMCAL with T variations, no params TH1");
808 
809  delete contRFTD;
810 }
811 
812 
815 //___________________________________________________________________
817 {
818  fEMCALGeo = AliEMCALGeometry::GetInstance(fEMCALGeoName) ;
819  Int_t nSM = (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules();
820 
821  fOutputContainer = new TList();
822  const Int_t buffersize = 255;
823  char hname[buffersize], htitl[buffersize];
824 
825  fhNEvents = new TH1I("hNEvents", "Number of analyzed events" , 1 , 0 , 1 ) ;
827 
828  fHmgg = new TH2F("hmgg","2-cluster invariant mass",fNbins,fMinBin,fMaxBin,100,0,10);
829  fHmgg->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
830  fHmgg->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
831  fOutputContainer->Add(fHmgg);
832 
833  fHmggDifferentSM = new TH2F("hmggDifferentSM","2-cluster invariant mass, different SM",fNbins,fMinBin,fMaxBin,100,0,10);
834  fHmggDifferentSM->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
835  fHmggDifferentSM->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
837 
838  fHOpeningAngle = new TH2F("hopang","2-cluster opening angle",100,0.,50.,100,0,10);
839  fHOpeningAngle->SetXTitle("#alpha_{#gamma #gamma}");
840  fHOpeningAngle->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
842 
843  fHOpeningAngleDifferentSM = new TH2F("hopangDifferentSM","2-cluster opening angle, different SM",100,0,50.,100,0,10);
844  fHOpeningAngleDifferentSM->SetXTitle("#alpha_{#gamma #gamma}");
845  fHOpeningAngleDifferentSM->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
847 
848  fHAsymmetry = new TH2F("hasym","2-cluster opening angle",100,0.,1.,100,0,10);
849  fHAsymmetry->SetXTitle("a");
850  fHAsymmetry->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
852 
853  fHAsymmetryDifferentSM = new TH2F("hasymDifferentSM","2-cluster opening angle, different SM",100,0,1.,100,0,10);
854  fHAsymmetryDifferentSM->SetXTitle("a");
855  fHAsymmetryDifferentSM->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
857 
858  //TString pairname[] = {"A side (0-2)", "C side (1-3)","Row 0 (0-1)", "Row 1 (2-3)"};
859 
860  fHmggMaskFrame = new TH2F("hmggMaskFrame","2-cluster invariant mass, frame masked",fNbins,fMinBin,fMaxBin,100,0,10);
861  fHmggMaskFrame->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
862  fHmggMaskFrame->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
864 
865  fHmggDifferentSMMaskFrame = new TH2F("hmggDifferentSMMaskFrame","2-cluster invariant mass, different SM, frame masked",
866  fNbins,fMinBin,fMaxBin,100,0,10);
867  fHmggDifferentSMMaskFrame->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
868  fHmggDifferentSMMaskFrame->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
870 
871  for(Int_t iSM = 0; iSM < nSM; iSM++)
872  {
873  snprintf(hname, buffersize, "hmgg_SM%d",iSM);
874  snprintf(htitl, buffersize, "Two-gamma inv. mass for super mod %d",iSM);
875  fHmggSM[iSM] = new TH2F(hname,htitl,fNbins,fMinBin,fMaxBin,100,0,10);
876  fHmggSM[iSM]->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
877  fHmggSM[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
878  fOutputContainer->Add(fHmggSM[iSM]);
879 
880  snprintf(hname, buffersize, "hmgg_SM%d_MaskFrame",iSM);
881  snprintf(htitl, buffersize, "Two-gamma inv. mass for super mod %d",iSM);
882  fHmggSMMaskFrame[iSM] = new TH2F(hname,htitl,fNbins,fMinBin,fMaxBin,100,0,10);
883  fHmggSMMaskFrame[iSM]->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
884  fHmggSMMaskFrame[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
886 
887  if(iSM < nSM/2)
888  {
889  snprintf(hname,buffersize, "hmgg_PairSameSectorSM%d",iSM);
890  snprintf(htitl,buffersize, "Two-gamma inv. mass for SM pair Sector: %d",iSM);
891  fHmggPairSameSectorSM[iSM] = new TH2F(hname,htitl,fNbins,fMinBin,fMaxBin,100,0,10);
892  fHmggPairSameSectorSM[iSM]->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
893  fHmggPairSameSectorSM[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
895 
896  snprintf(hname,buffersize, "hmgg_PairSameSectorSM%d_MaskFrame",iSM);
897  snprintf(htitl,buffersize, "Two-gamma inv. mass for SM pair Sector: %d",iSM);
898  fHmggPairSameSectorSMMaskFrame[iSM] = new TH2F(hname,htitl,fNbins,fMinBin,fMaxBin,100,0,10);
899  fHmggPairSameSectorSMMaskFrame[iSM]->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
900  fHmggPairSameSectorSMMaskFrame[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
902 
903  fhClusterPairDiffTimeSameSector[iSM] = new TH2F(Form("hClusterPairDiffTimeSameSector%d",iSM),
904  Form("cluster pair time difference vs E, Sector %d",iSM),
905  100,0,10, 200,-100,100);
906  fhClusterPairDiffTimeSameSector[iSM]->SetXTitle("E_{pair} (GeV)");
907  fhClusterPairDiffTimeSameSector[iSM]->SetYTitle("#Delta t (ns)");
909  }
910 
911  if(iSM < nSM-2)
912  {
913  snprintf(hname,buffersize, "hmgg_PairSameSideSM%d",iSM);
914  snprintf(htitl,buffersize, "Two-gamma inv. mass for SM pair Sector: %d",iSM);
915  fHmggPairSameSideSM[iSM] = new TH2F(hname,htitl,fNbins,fMinBin,fMaxBin,100,0,10);
916  fHmggPairSameSideSM[iSM]->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
917  fHmggPairSameSideSM[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
919 
920  snprintf(hname,buffersize, "hmgg_PairSameSideSM%d_MaskFrame",iSM);
921  snprintf(htitl,buffersize, "Two-gamma inv. mass for SM pair Sector: %d",iSM);
922  fHmggPairSameSideSMMaskFrame[iSM] = new TH2F(hname,htitl,fNbins,fMinBin,fMaxBin,100,0,10);
923  fHmggPairSameSideSMMaskFrame[iSM]->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
924  fHmggPairSameSideSMMaskFrame[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
926 
927  fhClusterPairDiffTimeSameSide[iSM] = new TH2F(Form("hClusterPairDiffTimeSameSide%d",iSM),
928  Form("cluster pair time difference vs E, Side %d",iSM),
929  100,0,10, 200,-100,100);
930  fhClusterPairDiffTimeSameSide[iSM]->SetXTitle("E_{pair} (GeV)");
931  fhClusterPairDiffTimeSameSide[iSM]->SetYTitle("#Delta t (ns)");
933  }
934 
935  snprintf(hname, buffersize, "hopang_SM%d",iSM);
936  snprintf(htitl, buffersize, "Opening angle for super mod %d",iSM);
937  fHOpeningAngleSM[iSM] = new TH2F(hname,htitl,100,0.,50.,100,0,10);
938  fHOpeningAngleSM[iSM]->SetXTitle("#alpha_{#gamma #gamma} (deg)");
939  fHOpeningAngleSM[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
941 
942  snprintf(hname,buffersize, "hopang_PairSM%d",iSM);
943  snprintf(htitl,buffersize, "Opening angle for SM pair: %d",iSM);
944  fHOpeningAnglePairSM[iSM] = new TH2F(hname,htitl,100,0.,50.,100,0,10);
945  fHOpeningAnglePairSM[iSM]->SetXTitle("#alpha_{#gamma #gamma} (deg)");
946  fHOpeningAnglePairSM[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
948 
949  snprintf(hname, buffersize, "hasym_SM%d",iSM);
950  snprintf(htitl, buffersize, "Asymmetry for super mod %d",iSM);
951  fHAsymmetrySM[iSM] = new TH2F(hname,htitl,100,0.,1.,100,0,10);
952  fHAsymmetrySM[iSM]->SetXTitle("a");
953  fHAsymmetrySM[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
954  fOutputContainer->Add(fHAsymmetrySM[iSM]);
955 
956  snprintf(hname,buffersize, "hasym_PairSM%d",iSM);
957  snprintf(htitl,buffersize, "Asymmetry for SM pair: %d",iSM);
958  fHAsymmetryPairSM[iSM] = new TH2F(hname,htitl,100,0.,1.,100,0,10);
959  fHAsymmetryPairSM[iSM]->SetXTitle("a");
960  fHAsymmetryPairSM[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
962 
963  Int_t colmax = 48;
964  Int_t rowmax = 24;
965 
966  fhTowerDecayPhotonHit[iSM] = new TH2F (Form("hTowerDecPhotonHit_Mod%d",iSM),
967  Form("Entries in grid of cells in Module %d",iSM),
968  colmax+2,-1.5,colmax+0.5, rowmax+2,-1.5,rowmax+0.5);
969  fhTowerDecayPhotonHit[iSM]->SetYTitle("row (phi direction)");
970  fhTowerDecayPhotonHit[iSM]->SetXTitle("column (eta direction)");
972 
973  fhTowerDecayPhotonEnergy[iSM] = new TH2F (Form("hTowerDecPhotonEnergy_Mod%d",iSM),
974  Form("Accumulated energy in grid of cells in Module %d",iSM),
975  colmax+2,-1.5,colmax+0.5, rowmax+2,-1.5,rowmax+0.5);
976  fhTowerDecayPhotonEnergy[iSM]->SetYTitle("row (phi direction)");
977  fhTowerDecayPhotonEnergy[iSM]->SetXTitle("column (eta direction)");
979 
980  fhTowerDecayPhotonAsymmetry[iSM] = new TH2F (Form("hTowerDecPhotonAsymmetry_Mod%d",iSM),
981  Form("Accumulated asymmetry in grid of cells in Module %d",iSM),
982  colmax+2,-1.5,colmax+0.5, rowmax+2,-1.5,rowmax+0.5);
983  fhTowerDecayPhotonAsymmetry[iSM]->SetYTitle("row (phi direction)");
984  fhTowerDecayPhotonAsymmetry[iSM]->SetXTitle("column (eta direction)");
986 
987  fhTowerDecayPhotonHitMaskFrame[iSM] = new TH2F (Form("hTowerDecPhotonHit_Mod%d_MaskFrame",iSM),Form("Entries in grid of cells in Module %d",iSM),
988  colmax+2,-1.5,colmax+0.5, rowmax+2,-1.5,rowmax+0.5);
989  fhTowerDecayPhotonHitMaskFrame[iSM]->SetYTitle("row (phi direction)");
990  fhTowerDecayPhotonHitMaskFrame[iSM]->SetXTitle("column (eta direction)");
992 
993  fhClusterTimeSM[iSM] = new TH2F(Form("hClusterTime_SM%d",iSM),"cluster time vs E",100,0,10, 100,0,1000);
994  fhClusterTimeSM[iSM]->SetXTitle("E (GeV)");
995  fhClusterTimeSM[iSM]->SetYTitle("t (ns)");
997 
998  fhClusterPairDiffTimeSameSM[iSM] = new TH2F(Form("hClusterPairDiffTimeSameSM%d",iSM),
999  Form("cluster pair time difference vs E, SM %d",iSM),
1000  100,0,10, 200,-100,100);
1001  fhClusterPairDiffTimeSameSM[iSM]->SetXTitle("E (GeV)");
1002  fhClusterPairDiffTimeSameSM[iSM]->SetYTitle("#Delta t (ns)");
1004  }
1005 
1006  Int_t nchannels = nSM*AliEMCALGeoParams::fgkEMCALRows*AliEMCALGeoParams::fgkEMCALCols;
1007  for(Int_t ibc = 0; ibc < 4; ibc++)
1008  {
1009  fHTpi0[ibc] = new TH2F(Form("hTime_BC%d",ibc),Form("Time of cell clusters under pi0 peak, bunch crossing %d",ibc),
1010  nchannels,0,nchannels, fNTimeBins,fMinTimeBin,fMaxTimeBin);
1011  fOutputContainer->Add(fHTpi0[ibc]);
1012  fHTpi0[ibc]->SetYTitle("time (ns)");
1013  fHTpi0[ibc]->SetXTitle("abs. Id. ");
1014  }
1015 
1016  fhClusterTime = new TH2F("hClusterTime","cluster time vs E",100,0,10, 100,0,1000);
1017  fhClusterTime->SetXTitle("E (GeV)");
1018  fhClusterTime->SetYTitle("t (ns)");
1020 
1021  fhClusterPairDiffTime = new TH2F("hClusterPairDiffTime","cluster pair time difference vs E",100,0,10, 800,-400,400);
1022  fhClusterPairDiffTime->SetXTitle("E_{pair} (GeV)");
1023  fhClusterPairDiffTime->SetYTitle("#Delta t (ns)");
1025 
1026  for(Int_t iMod=0; iMod < nSM; iMod++)
1027  {
1028  for(Int_t iRow=0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++)
1029  {
1030  for(Int_t iCol=0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++)
1031  {
1032  snprintf(hname,buffersize, "%d_%d_%d",iMod,iCol,iRow);
1033  snprintf(htitl,buffersize, "Two-gamma inv. mass for super mod %d, cell(col,row)=(%d,%d)",iMod,iCol,iRow);
1034  fHmpi0[iMod][iCol][iRow] = new TH1F(hname,htitl,fNbins,fMinBin,fMaxBin);
1035  fHmpi0[iMod][iCol][iRow]->SetXTitle("mass (MeV/c^{2})");
1036  fOutputContainer->Add(fHmpi0[iMod][iCol][iRow]);
1037  }
1038  }
1039  }
1040 
1041  fOutputContainer->SetOwner(kTRUE);
1042 
1043  PostData(1,fOutputContainer);
1044 
1045 // // cuts container, set in terminate but init and post here
1046 //
1047 // fCuts = new TList();
1048 //
1049 // fCuts ->SetOwner(kTRUE);
1050 //
1051 // PostData(2, fCuts);
1052 }
1053 
1059 //______________________________________________________________________________________________________
1061 {
1062  Int_t icol = ieta;
1063  if(iSM%2) icol+=48; // Impair SM, shift index [0-47] to [48-96]
1064 
1066  {
1067  for (Int_t imask = 0; imask < fNMaskCellColumns; imask++)
1068  {
1069  if(icol==fMaskCellColumns[imask]) return kTRUE;
1070  }
1071  }
1072 
1073  return kFALSE;
1074 }
1075 
1081 //__________________________________________________________________________
1083 {
1084  // Event selection
1085 
1086  if(fTriggerName!="")
1087  {
1088  AliESDEvent* esdevent = dynamic_cast<AliESDEvent*> (InputEvent());
1089  AliAODEvent* aodevent = dynamic_cast<AliAODEvent*> (InputEvent());
1090 
1091  TString triggerClass = "";
1092  if (esdevent) triggerClass = esdevent->GetFiredTriggerClasses();
1093  else if(aodevent) triggerClass = aodevent->GetFiredTriggerClasses();
1094 
1095  AliDebug(1,Form("Event %d, FiredClass %s",
1096  (Int_t)Entry(),(((AliESDEvent*)InputEvent())->GetFiredTriggerClasses()).Data()));
1097 
1098  if(!triggerClass.Contains(fTriggerName))
1099  {
1100  AliDebug(1,"Reject event!");
1101  return;
1102  }
1103  else
1104  AliDebug(1,"Accept event!");
1105  }
1106 
1107  // Get the input event
1108 
1109  AliVEvent* event = 0;
1110  if(fFilteredInput) event = AODEvent();
1111  else event = InputEvent();
1112 
1113  if(!event)
1114  {
1115  AliWarning("Input event not available!");
1116  return;
1117  }
1118 
1119  AliDebug(1,Form("<<< %s: Event %d >>>",event->GetName(), (Int_t)Entry()));
1120 
1121  // Get the primary vertex
1122 
1123  event->GetPrimaryVertex()->GetXYZ(fVertex) ;
1124 
1125  AliDebug(1,Form("Vertex: (%.3f,%.3f,%.3f)",fVertex[0],fVertex[1],fVertex[2]));
1126 
1127  //Int_t runNum = aod->GetRunNumber();
1128  //if(DebugLevel() > 1) printf("Run number: %d\n",runNum);
1129 
1130  fhNEvents->Fill(0); //Count the events to be analyzed
1131 
1132  // Acccess once the geometry matrix and temperature corrections and calibration coefficients
1133  if(fhNEvents->GetEntries() == 1)
1134  {
1136 
1138 
1140  }
1141 
1142  //Get the list of clusters and cells
1143  fEMCALCells = event->GetEMCALCells();
1144 
1145  fCaloClustersArr = new TRefArray();
1146  event->GetEMCALClusters(fCaloClustersArr);
1147 
1148  AliDebug(1,Form("N CaloClusters: %d - N CaloCells %d",fCaloClustersArr->GetEntriesFast(), fEMCALCells->GetNumberOfCells()));
1149 
1150  // Apply non linearity, new calibration, T calibration to the clusters
1151  if( fCorrectClusters )
1152  CorrectClusters();
1153 
1154  FillHistograms();
1155 
1156  delete fCaloClustersArr;
1157 
1158  PostData(1,fOutputContainer);
1159 }
1160 
1163 //_____________________________________________________
1165 {
1166  printf("Cluster cuts: %2.2f < E < %2.2f GeV; number of cells > %d; Assymetry < %1.2f, pair time diff < %2.2f, %2.2f < t < %2.2f ns\n",
1168 
1169  printf("Group %d cells\n", fGroupNCells) ;
1170 
1171  printf("Cluster maximal cell away from border at least %d cells\n", fRecoUtils->GetNumberOfCellsFromEMCALBorder()) ;
1172 
1173  printf("Histograms: bins %d; energy range: %2.2f < E < %2.2f MeV\n",fNbins,fMinBin,fMaxBin) ;
1174 
1175  printf("Switchs:\n \t Remove Bad Channels? %d; Use filtered input? %d; Correct Clusters? %d, and their position? %d \n \t Mass per channel same SM clusters? %d\n",
1176  fRecoUtils->IsBadChannelsRemovalSwitchedOn(),fFilteredInput,fCorrectClusters, fRecalPosition, fSameSM) ;
1177 
1178  printf("OADB path : %s\n",fOADBFilePath .Data());
1179  printf("Calibration path : %s\n",fCalibFilePath.Data());
1180 
1181  printf("EMCAL Geometry name: < %s >, Load Matrices %d\n",fEMCALGeoName.Data(), fLoadMatrices) ;
1182 
1183  if(fLoadMatrices) { for(Int_t ism = 0; ism < AliEMCALGeoParams::fgkEMCALModules; ism++) if(fMatrix[ism]) fMatrix[ism]->Print() ; }
1184 }
1185 
1189 //_____________________________________________________________________
1191 {
1192  if(n > fNMaskCellColumns)
1193  {
1194  delete [] fMaskCellColumns ;
1195 
1196  fMaskCellColumns = new Int_t[n] ;
1197  }
1198 
1199  fNMaskCellColumns = n ;
1200 }
1201 
1206 //___________________________________________________________________________________
1208 {
1209  if(ipos < fNMaskCellColumns) fMaskCellColumns[ipos] = icol ;
1210  else AliWarning("Mask column not set, position larger than allocated set size first") ;
1211 }
1212 
1215 //______________________________________________________________
1217 {
1218  AliDebug(1,"Not implemented");
1219 // const Int_t buffersize = 255;
1220 // char onePar[buffersize] ;
1221 
1222 // snprintf(onePar,buffersize, "Custer cuts: %2.2f < E < %2.2f GeV; %2.2f < Lambda0_2 < %2.2f GeV; min number of cells %d; Assymetry cut %1.2f, time1-time2 < %2.2f; %2.2f < T < %2.2f ns; %3.1f < Mass < %3.1f",
1223 // fEmin,fEmax, fL0min, fL0max, fMinNCells, fAsyCut, fDTimeCut, fTimeMin, fTimeMax, fInvMassCutMin, fInvMassCutMax) ;
1224 // fCuts->Add(new TObjString(onePar));
1225 // snprintf(onePar,buffersize, "Group %d cells;", fGroupNCells) ;
1226 // fCuts->Add(new TObjString(onePar));
1227 // snprintf(onePar,buffersize, "Cluster maximal cell away from border at least %d cells;", fRecoUtils->GetNumberOfCellsFromEMCALBorder()) ;
1228 // fCuts->Add(new TObjString(onePar));
1229 // snprintf(onePar,buffersize, "Histograms, Mass bins %d; energy range: %2.2f < E < %2.2f GeV;",fNbins,fMinBin,fMaxBin) ;
1230 // fCuts->Add(new TObjString(onePar));
1231 // snprintf(onePar,buffersize, "Histograms, Time bins %d; energy range: %2.2f < E < %2.2f GeV;",fNTimeBins,fMinTimeBin,fMaxTimeBin) ;
1232 // fCuts->Add(new TObjString(onePar));
1233 // snprintf(onePar,buffersize, "Switchs: Remove Bad Channels? %d; Use filtered input? %d; Correct Clusters? %d and their position? %d, Mass per channel same SM clusters? %d ",
1234 // fRecoUtils->IsBadChannelsRemovalSwitchedOn(),fFilteredInput,fCorrectClusters, fRecalPosition, fSameSM) ;
1235 // fCuts->Add(new TObjString(onePar));
1236 // snprintf(onePar,buffersize, "EMCAL Geometry name: < %s >, Load Matrices? %d",fEMCALGeoName.Data(),fLoadMatrices) ;
1237 // fCuts->Add(new TObjString(onePar));
1238 //
1239 // // Post Data
1240 // PostData(2, fCuts);
1241 }
1242 
Float_t fDTimeCut
Maximum difference between time of cluster pairs (ns).
ClassImp(AliAnalysisTaskTriggerRates) AliAnalysisTaskTriggerRates
TH2F * fHmggPairSameSectorSMMaskFrame[AliEMCALGeoParams::fgkEMCALModules/2]
! Two-cluster invariant mass per Pair, mask clusters facing frames.
TGeoHMatrix * fMatrix[AliEMCALGeoParams::fgkEMCALModules]
TH2F * fhClusterPairDiffTimeSameSide[AliEMCALGeoParams::fgkEMCALModules-2]
! Diference in time of clusters same side.
Int_t fNTimeBins
N time bins of invariant mass histograms.
TH2F * fhTowerDecayPhotonHit[AliEMCALGeoParams::fgkEMCALModules]
! Cells ordered in column/row for different module, number of times a decay photon hits...
Float_t fInvMassCutMax
Maximum mass cut for clusters to fill time or other histograms.
TH2F * fHAsymmetryPairSM[AliEMCALGeoParams::fgkEMCALModules]
! Two-cluster asymmetry vs pt per Pair,with mass close to pi0.
TH1F * fHmpi0[AliEMCALGeoParams::fgkEMCALModules][AliEMCALGeoParams::fgkEMCALCols][AliEMCALGeoParams::fgkEMCALRows]
< Two-cluster invariant mass assigned to each cell.
TString fImportGeometryFilePath
Path fo geometry.root file.
Bool_t fSameSM
Combine clusters in channels on same SM.
Bool_t fLoadMatrices
Matrices set from configuration, not get from geometry.root or from ESDs/AODs.
AliEMCALRecoUtils * fRecoUtils
Access to reconstruction utilities.
AliEMCALGeometry * fEMCALGeo
! EMCAL geometry pointer.
Float_t fMaxTimeBin
Maximum time bins of invariant mass histograms.
Float_t fLogWeight
Logarithmic weight used in cluster recalibration.
TH2F * fHAsymmetry
! Two-cluster asymmetry vs pt of pair, with mass close to pi0.
This task provides the input for the EMCal energy calibration with pi0 invariant mass analysis per ch...
Float_t fMinTimeBin
Minimum time bins of invariant mass histograms.
TH2F * fHOpeningAngleDifferentSM
! Two-cluster opening angle vs pt of pair, each cluster in different SM, with mass close to pi0...
TH2F * fhClusterPairDiffTimeSameSM[AliEMCALGeoParams::fgkEMCALModules]
! Diference in time of clusters same SM.
TH2F * fHOpeningAngle
! Two-cluster opening angle vs pt of pair, with mass close to pi0.
TH2F * fhClusterTimeSM[AliEMCALGeoParams::fgkEMCALModules]
! Timing of clusters vs energy per SM.
TH2F * fHOpeningAngleSM[AliEMCALGeoParams::fgkEMCALModules]
! Two-cluster opening angle vs pt per SM,with mass close to pi0.
TLorentzVector fMomentum2
Cluster kinematics, temporal.
Float_t fInvMassCutMin
Minimum mass cut for clusters to fill time or other histograms.
TH2F * fHAsymmetrySM[AliEMCALGeoParams::fgkEMCALModules]
! Two-cluster asymmetry vs pt per SM,with mass close to pi0.
TString fCalibFilePath
Full path with file with energy calibration factors per channel from previous iteration.
TH2F * fHmggSM[AliEMCALGeoParams::fgkEMCALModules]
! Two-cluster invariant mass per SM.
TH2F * fhTowerDecayPhotonHitMaskFrame[AliEMCALGeoParams::fgkEMCALModules]
! Cells ordered in column/row for different module, number of times a decay photon hits...
TH2F * fHmgg
! Two-cluster invariant mass vs pt of pair.
TH2F * fHTpi0[4]
! Time of cell under pi0 mass, for 4 bunch crossings.
AliAnalysisTaskEMCALPi0CalibSelection()
Default constructor. Arrays initialization is done here.
Float_t fMinBin
Minimum mass bins of invariant mass histograms.
Bool_t fRecalPosition
Switch on/off cluster position calculation, in case alignment matrices are not available.
TString fTriggerName
Trigger name must contain this name.
TH2F * fHmggDifferentSMMaskFrame
! Two-cluster invariant mass vs pt of pair, each cluster in different SM,mask clusters facing frames...
TH2F * fhClusterPairDiffTime
! Diference in time of clusters.
Bool_t Data(TH1F *h, Double_t *rangefit, Bool_t writefit, Double_t &sgn, Double_t &errsgn, Double_t &bkg, Double_t &errbkg, Double_t &sgnf, Double_t &errsgnf, Double_t &sigmafit, Int_t &status)
Int_t fNbins
N mass bins of invariant mass histograms.
TH2F * fhClusterTime
! Timing of clusters vs energy.
TLorentzVector fMomentum12
Cluster pair kinematics, temporal.
TH2F * fhTowerDecayPhotonAsymmetry[AliEMCALGeoParams::fgkEMCALModules]
! Cells ordered in column/row for different module, accumulated asymmetry in the tower by decay photo...
TH2F * fHmggPairSameSideSMMaskFrame[AliEMCALGeoParams::fgkEMCALModules-2]
! Two-cluster invariant mass per Pair, mask clusters facing frames.
Bool_t fCorrectClusters
Correct clusters energy, position etc.
void UserCreateOutputObjects()
Create output container, init geometry.
TH2F * fHmggSMMaskFrame[AliEMCALGeoParams::fgkEMCALModules]
! Two-cluster invariant mass per SM, mask clusters facing frames.
void Terminate(Option_t *opt)
Create cuts/param objects and publish to slot. Comment out for the moment.
TH2F * fhTowerDecayPhotonEnergy[AliEMCALGeoParams::fgkEMCALModules]
! Cells ordered in column/row for different module, accumulated energy in the tower by decay photons...
TH2F * fHmggPairSameSideSM[AliEMCALGeoParams::fgkEMCALModules-2]
! Two-cluster invariant mass per Pair.
Float_t fMaxBin
Maximum mass bins of invariant mass histograms.
TH2F * fHmggDifferentSM
! Two-cluster invariant mass vs pt of pair, each cluster in different SM.
TH2F * fhClusterPairDiffTimeSameSector[AliEMCALGeoParams::fgkEMCALModules/2]
! Diference in time of clusters same sector.
TH1I * fhNEvents
! Number of events counter histogram.
TLorentzVector fMomentum1
Cluster kinematics, temporal.
TH2F * fHmggMaskFrame
! Two-cluster invariant mass vs pt of pair, mask clusters facing frames.
Bool_t fImportGeometryFromFile
Import geometry settings in geometry.root file.
TString fOADBFilePath
Default path $ALICE_PHYSICS/OADB/EMCAL, if needed change.
TH2F * fHOpeningAnglePairSM[AliEMCALGeoParams::fgkEMCALModules]
! Two-cluster opening angle vs pt per Pair,with mass close to pi0.
Bool_t fFilteredInput
Read input produced with filter.
TH2F * fHAsymmetryDifferentSM
! Two-cluster asymmetry vs pt of pair, each cluster in different SM, with mass close to pi0...
TH2F * fHmggPairSameSectorSM[AliEMCALGeoParams::fgkEMCALModules/2]
! Two-cluster invariant mass per Pair.