AliPhysics  66e96a0 (66e96a0)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 fNEnergybins(1000),
61 fMinEnergyBin(0.), fMaxEnergyBin(100.),
62 fNTimeBins(1000), fMinTimeBin(0.), fMaxTimeBin(1000.),
63 // Temporal
64 fMomentum1(), fMomentum2(), fMomentum12(),
65 // Histograms
66 fHmgg(0x0), fHmggDifferentSM(0x0),
67 fHmggMaskFrame(0x0), fHmggDifferentSMMaskFrame(0x0),
68 fHOpeningAngle(0x0), fHOpeningAngleDifferentSM(0x0),
69 fHAsymmetry(0x0), fHAsymmetryDifferentSM(0x0),
70 fhNEvents(0x0),
71 fhClusterTime(0x0), fhClusterPairDiffTime(0x0)
72 {
73  for(Int_t iMod=0; iMod < AliEMCALGeoParams::fgkEMCALModules; iMod++)
74  {
75  for(Int_t iX=0; iX<24; iX++)
76  {
77  for(Int_t iZ=0; iZ<48; iZ++)
78  {
79  fHmpi0[iMod][iZ][iX] = 0 ;
80  fhEnergy[iMod][iZ][iX] = 0 ;
81  }
82  }
83  }
84 
85  fVertex[0]=fVertex[1]=fVertex[2]=-1000;
86 
87  fHTpi0[0]= 0 ;
88  fHTpi0[1]= 0 ;
89  fHTpi0[2]= 0 ;
90  fHTpi0[3]= 0 ;
91 
93  fMaskCellColumns[0] = 6 ; fMaskCellColumns[1] = 7 ; fMaskCellColumns[2] = 8 ;
94  fMaskCellColumns[3] = 35; fMaskCellColumns[4] = 36; fMaskCellColumns[5] = 37;
95  fMaskCellColumns[6] = 12+AliEMCALGeoParams::fgkEMCALCols; fMaskCellColumns[7] = 13+AliEMCALGeoParams::fgkEMCALCols;
96  fMaskCellColumns[8] = 40+AliEMCALGeoParams::fgkEMCALCols; fMaskCellColumns[9] = 41+AliEMCALGeoParams::fgkEMCALCols;
97  fMaskCellColumns[10]= 42+AliEMCALGeoParams::fgkEMCALCols;
98 
99  for(Int_t iSMPair = 0; iSMPair < AliEMCALGeoParams::fgkEMCALModules/2; iSMPair++)
100  {
101  fHmggPairSameSectorSM[iSMPair] = 0;
102  fHmggPairSameSectorSMMaskFrame[iSMPair] = 0;
104  }
105 
106  for(Int_t iSMPair = 0; iSMPair < AliEMCALGeoParams::fgkEMCALModules-2; iSMPair++)
107  {
108  fHmggPairSameSideSM[iSMPair] = 0;
109  fHmggPairSameSideSMMaskFrame[iSMPair] = 0;
110  fhClusterPairDiffTimeSameSide[iSMPair] = 0;
111  }
112 
113  for(Int_t iSM = 0; iSM < AliEMCALGeoParams::fgkEMCALModules; iSM++)
114  {
115  fHmggSM[iSM] = 0;
116  fHmggSM_Zone1[iSM] = 0;
117  fHmggSM_Zone2[iSM] = 0;
118  fHmggSM_Zone3[iSM] = 0;
119  fHmggSM_Zone4[iSM] = 0;
120  fHmggSM_Zone5[iSM] = 0;
121  fHmggSM_Zone6[iSM] = 0;
122  fHmggSM_Zone7[iSM] = 0;
123  fHmggSMMaskFrame[iSM] = 0;
124  fHOpeningAngleSM[iSM] = 0;
125  fHOpeningAnglePairSM[iSM] = 0;
126  fHAsymmetrySM[iSM] = 0;
127  fHAsymmetryPairSM[iSM] = 0;
128  fhTowerDecayPhotonHit[iSM] = 0;
129  fhTowerDecayPhotonEnergy[iSM] = 0;
132  fMatrix[iSM] = 0x0;
133  fhClusterTimeSM[iSM] = 0;
134  fhTopoClusterAmpCase0[iSM] =0;
135  fhTopoClusterAmpCase1[iSM] =0;
136  fhTopoClusterAmpCase2[iSM] =0;
137  fhTopoClusterAmpCase3[iSM] =0;
143  }
144 }
145 
151 //______________________________________________________________________________________________
153 AliAnalysisTaskSE(name),
154 fEMCALGeo(0x0), fLoadMatrices(0),
155 fEMCALGeoName("EMCAL_COMPLETE12SMV1_DCAL_8SM"),
156 fTriggerName("EMC"),
157 fRecoUtils(new AliEMCALRecoUtils),
158 fOADBFilePath(""), fCalibFilePath(""),
159 fCorrectClusters(kFALSE), fRecalPosition(kTRUE),
160 fCaloClustersArr(0x0), fEMCALCells(0x0),
161 //fCuts(0x0),
162 fOutputContainer(0x0),
163 fVertex(), fFilteredInput(kFALSE),
164 fImportGeometryFromFile(1), fImportGeometryFilePath(""),
165 fEmin(0.5), fEmax(15.),
166 fL0min(0.01), fL0max(0.5),
167 fDTimeCut(100.), fTimeMax(1000000), fTimeMin(-1000000),
168 fAsyCut(1.), fMinNCells(2), fGroupNCells(0),
169 fLogWeight(4.5), fSameSM(kFALSE),
170 fNMaskCellColumns(11), fMaskCellColumns(0x0),
171 fInvMassCutMin(110.), fInvMassCutMax(160.),
172 // Histograms binning
173 fNbins(300),
174 fMinBin(0.), fMaxBin(300.),
175 fNEnergybins(1000),
176 fMinEnergyBin(0.), fMaxEnergyBin(100.),
177 fNTimeBins(1000), fMinTimeBin(0.), fMaxTimeBin(1000.),
178 // Temporal
179 fMomentum1(), fMomentum2(), fMomentum12(),
180 // Histograms
181 fHmgg(0x0), fHmggDifferentSM(0x0),
182 fHmggMaskFrame(0x0), fHmggDifferentSMMaskFrame(0x0),
183 fHOpeningAngle(0x0), fHOpeningAngleDifferentSM(0x0),
184 fHAsymmetry(0x0), fHAsymmetryDifferentSM(0x0),
185 fhNEvents(0x0),
186 fhClusterTime(0x0), fhClusterPairDiffTime(0x0)
187 {
188  for(Int_t iMod=0; iMod < AliEMCALGeoParams::fgkEMCALModules; iMod++)
189  {
190  for(Int_t iX=0; iX<24; iX++)
191  {
192  for(Int_t iZ=0; iZ<48; iZ++)
193  {
194  fHmpi0[iMod][iZ][iX] = 0 ;
195  fhEnergy[iMod][iZ][iX] = 0 ;
196  }
197  }
198  }
199 
200  fVertex[0]=fVertex[1]=fVertex[2]=-1000;
201 
202  fHTpi0[0]= 0 ;
203  fHTpi0[1]= 0 ;
204  fHTpi0[2]= 0 ;
205  fHTpi0[3]= 0 ;
206 
207  fMaskCellColumns = new Int_t[fNMaskCellColumns];
208  fMaskCellColumns[0] = 6 ; fMaskCellColumns[1] = 7 ; fMaskCellColumns[2] = 8 ;
209  fMaskCellColumns[3] = 35; fMaskCellColumns[4] = 36; fMaskCellColumns[5] = 37;
210  fMaskCellColumns[6] = 12+AliEMCALGeoParams::fgkEMCALCols; fMaskCellColumns[7] = 13+AliEMCALGeoParams::fgkEMCALCols;
211  fMaskCellColumns[8] = 40+AliEMCALGeoParams::fgkEMCALCols; fMaskCellColumns[9] = 41+AliEMCALGeoParams::fgkEMCALCols;
212  fMaskCellColumns[10]= 42+AliEMCALGeoParams::fgkEMCALCols;
213 
214  for(Int_t iSMPair = 0; iSMPair < AliEMCALGeoParams::fgkEMCALModules/2; iSMPair++)
215  {
216  fHmggPairSameSectorSM[iSMPair] = 0;
217  fHmggPairSameSectorSMMaskFrame[iSMPair] = 0;
219  }
220 
221  for(Int_t iSMPair = 0; iSMPair < AliEMCALGeoParams::fgkEMCALModules-2; iSMPair++)
222  {
223  fHmggPairSameSideSM[iSMPair] = 0;
224  fHmggPairSameSideSMMaskFrame[iSMPair] = 0;
225  fhClusterPairDiffTimeSameSide[iSMPair] = 0;
226  }
227 
228  for(Int_t iSM = 0; iSM < AliEMCALGeoParams::fgkEMCALModules; iSM++)
229  {
230  fHmggSM[iSM] = 0;
231  fHmggSM_Zone1[iSM] = 0;
232  fHmggSM_Zone2[iSM] = 0;
233  fHmggSM_Zone3[iSM] = 0;
234  fHmggSM_Zone4[iSM] = 0;
235  fHmggSM_Zone5[iSM] = 0;
236  fHmggSM_Zone6[iSM] = 0;
237  fHmggSM_Zone7[iSM] = 0;
238  fHmggSMMaskFrame[iSM] = 0;
239  fHOpeningAngleSM[iSM] = 0;
240  fHOpeningAnglePairSM[iSM] = 0;
241  fHAsymmetrySM[iSM] = 0;
242  fHAsymmetryPairSM[iSM] = 0;
243  fhTowerDecayPhotonHit[iSM] = 0;
244  fhTowerDecayPhotonEnergy[iSM] = 0;
247  fMatrix[iSM] = 0x0;
248  fhClusterTimeSM[iSM] = 0;
249  fhTopoClusterAmpCase0[iSM] =0;
250  fhTopoClusterAmpCase1[iSM] =0;
251  fhTopoClusterAmpCase2[iSM] =0;
252  fhTopoClusterAmpCase3[iSM] =0;
258  }
259 
260  DefineOutput(1, TList::Class());
261 //DefineOutput(2, TList::Class()); // will contain cuts or local params
262 }
263 
266 //_____________________________________________________________________________
268 {
269  if(fOutputContainer)
270  {
271  fOutputContainer->Delete() ;
272  delete fOutputContainer ;
273  }
274 
275  if(fEMCALGeo) delete fEMCALGeo ;
276  if(fRecoUtils) delete fRecoUtils ;
277  if(fNMaskCellColumns) delete [] fMaskCellColumns;
278 }
279 
283 //____________________________________________________________
285 {
286  if(fRecoUtils->GetParticleType()!=AliEMCALRecoUtils::kPhoton)
287  AliFatal(Form("Wrong particle type for cluster position recalculation! = %d\n", fRecoUtils->GetParticleType()));
288 
289  AliDebug(1,Form("It will use fLogWeight %.3f",fLogWeight));
290 
291  Float_t pos[]={0,0,0};
292 
293  for(Int_t iClu=0; iClu < fCaloClustersArr->GetEntriesFast(); iClu++)
294  {
295  AliVCluster *c1 = (AliVCluster *) fCaloClustersArr->At(iClu);
296 
297  Float_t e1i = c1->E(); // cluster energy before correction
298  if (e1i < fEmin) continue;
299  else if (e1i > fEmax) continue;
300 
301  else if (c1->GetNCells() < fMinNCells) continue;
302 
303  else if (c1->GetM02() < fL0min || c1->GetM02() > fL0max) continue;
304 
305  if(fRecoUtils->ClusterContainsBadChannel(fEMCALGeo, c1->GetCellsAbsId(), c1->GetNCells())) continue;
306 
307  if(DebugLevel() > 2)
308  {
309  AliInfo(Form("Std : i %d, E %f, dispersion %f, m02 %f, m20 %f\n",c1->GetID(),c1->E(),c1->GetDispersion(),c1->GetM02(),c1->GetM20()));
310  c1->GetPosition(pos);
311  AliInfo(Form("Std : i %d, x %f, y %f, z %f\n",c1->GetID(), pos[0], pos[1], pos[2]));
312  }
313 
314  // Correct cluster energy and position if requested, and not corrected previously, by default Off
315  if(fRecoUtils->IsRecalibrationOn())
316  {
317  fRecoUtils->RecalibrateClusterEnergy(fEMCALGeo, c1, fEMCALCells);
318  fRecoUtils->RecalculateClusterShowerShapeParameters(fEMCALGeo, fEMCALCells,c1);
319  fRecoUtils->RecalculateClusterPID(c1);
320  }
321 
322  AliDebug(2,Form("Energy: after recalibration %f",c1->E()));
323 
324  // Recalculate cluster position
325  if ( fRecalPosition ) fRecoUtils->RecalculateClusterPosition(fEMCALGeo, fEMCALCells,c1);
326 
327  // Correct Non-Linearity
328  c1->SetE(fRecoUtils->CorrectClusterEnergyLinearity(c1));
329 
330  AliDebug(2,Form("after linearity correction %f",c1->E()));
331 
332  // In case of MC analysis, to match resolution/calibration in real data
333  //c1->SetE(fRecoUtils->SmearClusterEnergy(c1)); // Not needed anymore
334 
335  AliDebug(2,Form("after smearing %f\n",c1->E()));
336 
337  if(DebugLevel() > 2)
338  {
339  AliInfo(Form("Cor : i %d, E %f, dispersion %f, m02 %f, m20 %f\n",c1->GetID(),c1->E(),c1->GetDispersion(),c1->GetM02(),c1->GetM20()));
340  c1->GetPosition(pos);
341  AliInfo(Form("Cor : i %d, x %f, y %f, z %f\n",c1->GetID(), pos[0], pos[1], pos[2]));
342  }
343  } // cluster loop
344 }
345 
349 //__________________________________________________________
351 {
352  Int_t absId1 = -1;
353  Int_t iSupMod1 = -1;
354  Int_t iphi1 = -1;
355  Int_t ieta1 = -1;
356  Int_t absId2 = -1;
357  Int_t iSupMod2 = -1;
358  Int_t iphi2 = -1;
359  Int_t ieta2 = -1;
360  Bool_t shared = kFALSE;
361 
362  Float_t pos[] = {0,0,0};
363 
364  Int_t bc = InputEvent()->GetBunchCrossNumber();
365  Int_t nSM = (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules();
366 
367  for(Int_t iClu=0; iClu<fCaloClustersArr->GetEntriesFast()-1; iClu++)
368  {
369  AliVCluster *c1 = (AliVCluster *) fCaloClustersArr->At(iClu);
370 
371  // Exclude bad channels
372  if(fRecoUtils->ClusterContainsBadChannel(fEMCALGeo, c1->GetCellsAbsId(), c1->GetNCells())) continue;
373 
374  Float_t e1i = c1->E(); // cluster energy before correction
375 
376  if (e1i < fEmin) continue;
377  else if (e1i > fEmax) continue;
378 
379  else if (!fRecoUtils->IsGoodCluster(c1,fEMCALGeo,fEMCALCells,bc)) continue;
380 
381  else if (c1->GetNCells() < fMinNCells) continue;
382 
383  else if (c1->GetM02() < fL0min || c1->GetM02() > fL0max) continue;
384 
385  if(DebugLevel() > 2)
386  {
387  AliInfo(Form("IMA : i %d, E %f, dispersion %f, m02 %f, m20 %f",c1->GetID(),e1i,c1->GetDispersion(),c1->GetM02(),c1->GetM20()));
388  c1->GetPosition(pos);
389  AliInfo(Form("IMA : i %d, x %f, y %f, z %f",c1->GetID(), pos[0], pos[1], pos[2]));
390  }
391 
392  fRecoUtils->GetMaxEnergyCell(fEMCALGeo, fEMCALCells,c1,absId1,iSupMod1,ieta1,iphi1,shared);
393 
394  c1->GetMomentum(fMomentum1,fVertex);
395 
396  // Check if cluster is in fidutial region, not too close to borders
397  Bool_t in1 = fRecoUtils->CheckCellFiducialRegion(fEMCALGeo, c1, fEMCALCells);
398 
399  // Clusters not facing frame structures
400  Bool_t mask1 = MaskFrameCluster(iSupMod1, ieta1);
401  //if(mask1) printf("Reject eta %d SM %d\n",ieta1, iSupMod1);
402 
403  Double_t time1 = c1->GetTOF()*1.e9;
404 
405  if(fSelectOnlyCellSignalOutOfCollision && (time1 < fTimeMax || time1 > fTimeMin)) continue;
406  else if(!fSelectOnlyCellSignalOutOfCollision && (time1 > fTimeMax || time1 < fTimeMin)) continue;
407 
408  fhClusterTime ->Fill(c1->E(),time1);
409  fhClusterTimeSM[iSupMod1]->Fill(c1->E(),time1);
410 
411  if(fClusterTopology)
412  {
413  Int_t iPosInNoisyQuartet = FindPositionInNoisyQuartet(iphi1,ieta1,iSupMod1);
414  AliEMCALGeometry* geom = AliEMCALGeometry::GetInstance();
415 
416  for(Int_t iCell = 0; iCell < c1->GetNCells(); iCell++)
417  {
418  Int_t iSupMod = -1, iIeta =-1, iIphi =-1, iTower =-1;
419 
420  Int_t CellID = c1->GetCellsAbsId()[iCell];
421  geom->GetCellIndex(CellID,iSupMod,iTower,iIphi,iIeta);
422  Float_t AmpFraction = c1->GetCellAmplitudeFraction(CellID);
423  Float_t amp = fEMCALCells->GetCellAmplitude(CellID);
424 
425  switch (iPosInNoisyQuartet) {
426  case 0:
427  fhTopoClusterAmpCase0[iSupMod1]->Fill(iIeta-ieta1,iIphi-iphi1,amp);
428  fhTopoClusterAmpFractionCase0[iSupMod1]->Fill(iIeta-ieta1,iIphi-iphi1,AmpFraction);
429  break;
430  case 1:
431  fhTopoClusterAmpCase1[iSupMod1]->Fill(iIeta-ieta1,iIphi-iphi1,amp);
432  fhTopoClusterAmpFractionCase1[iSupMod1]->Fill(iIeta-ieta1,iIphi-iphi1,AmpFraction);
433  break;
434  case 2:
435  fhTopoClusterAmpCase2[iSupMod1]->Fill(iIeta-ieta1,iIphi-iphi1,amp);
436  fhTopoClusterAmpFractionCase2[iSupMod1]->Fill(iIeta-ieta1,iIphi-iphi1,AmpFraction);
437  break;
438  case 3:
439  fhTopoClusterAmpCase3[iSupMod1]->Fill(iIeta-ieta1,iIphi-iphi1,amp);
440  fhTopoClusterAmpFractionCase3[iSupMod1]->Fill(iIeta-ieta1,iIphi-iphi1,AmpFraction);
441  break;
442  default:
443  break;
444  }
445  }
446  }
447 
448  // Combine cluster with other clusters and get the invariant mass
449  for (Int_t jClu=iClu+1; jClu < fCaloClustersArr->GetEntriesFast(); jClu++)
450  {
451  AliAODCaloCluster *c2 = (AliAODCaloCluster *) fCaloClustersArr->At(jClu);
452 
453  Float_t e2i = c2->E();
454  if (e2i < fEmin) continue;
455  else if (e2i > fEmax) continue;
456 
457  else if (!fRecoUtils->IsGoodCluster(c2,fEMCALGeo,fEMCALCells,bc))continue;
458 
459  else if (c2->GetNCells() < fMinNCells) continue;
460 
461  else if (c2->GetM02() < fL0min || c2->GetM02() > fL0max) continue;
462 
463  fRecoUtils->GetMaxEnergyCell(fEMCALGeo, fEMCALCells,c2,absId2,iSupMod2,ieta2,iphi2,shared);
464 
465  c2->GetMomentum(fMomentum2,fVertex);
466 
468  Float_t invmass = fMomentum12.M()*1000;
469 
470  //Asimetry cut
471  Float_t asym = TMath::Abs(fMomentum1.E()-fMomentum2.E())/(fMomentum1.E()+fMomentum2.E());
472 
473  if(asym > fAsyCut) continue;
474 
475  //Time cut
476  Double_t time2 = c2->GetTOF()*1.e9;
477 
478  if(fSelectOnlyCellSignalOutOfCollision && (time2 < fTimeMax || time2 > fTimeMin)) continue;
479  else if(!fSelectOnlyCellSignalOutOfCollision && (time2 > fTimeMax || time2 < fTimeMin)) continue;
480 
481  fhClusterPairDiffTime->Fill(fMomentum12.E(),time1-time2);
482  if(TMath::Abs(time1-time2) > fDTimeCut) continue;
483 
484  if(invmass < fMaxBin && invmass > fMinBin )
485  {
486  //Check if cluster is in fidutial region, not too close to borders
487  Bool_t in2 = fRecoUtils->CheckCellFiducialRegion(fEMCALGeo, c2, fEMCALCells);
488 
489  // Clusters not facing frame structures
490  Bool_t mask2 = MaskFrameCluster(iSupMod2, ieta2);
491  //if(mask2) printf("Reject eta %d SM %d\n",ieta2, iSupMod2);
492 
493  if(fSelectOnlyPhotonsInDifferentSM && (iSupMod1 == iSupMod2)) continue;
494 
495  if(in1 && in2)
496  {
497  fHmgg->Fill(invmass,fMomentum12.Pt());
498 
499  if(iSupMod1==iSupMod2)
500  {
501  fHmggSM[iSupMod1]->Fill(invmass,fMomentum12.Pt());
502  fhClusterPairDiffTimeSameSM[iSupMod1]->Fill(fMomentum12.E(),time1-time2);
503 
504  //Is in zone number i
505  Bool_t zone1 = IsInZone1(iSupMod1,ieta1,iphi1);
506  Bool_t zone2 = IsInZone2(iSupMod1,ieta1,iphi1);
507  Bool_t zone3 = IsInZone3(iSupMod1,ieta1,iphi1);
508  Bool_t zone4 = IsInZone4(iSupMod1,ieta1,iphi1);
509  Bool_t zone5 = IsInZone5(iSupMod1,ieta1,iphi1);
510  Bool_t zone6 = IsInZone6(iSupMod1,ieta1,iphi1);
511  Bool_t zone7 = IsInZone7(iSupMod1,ieta1,iphi1);
512 
513 
514  if(zone1) fHmggSM_Zone1[iSupMod1]->Fill(invmass,fMomentum12.Pt());
515  if(zone2) fHmggSM_Zone2[iSupMod1]->Fill(invmass,fMomentum12.Pt());
516  if(zone3) fHmggSM_Zone3[iSupMod1]->Fill(invmass,fMomentum12.Pt());
517  if(zone4) fHmggSM_Zone4[iSupMod1]->Fill(invmass,fMomentum12.Pt());
518  if(zone5) fHmggSM_Zone5[iSupMod1]->Fill(invmass,fMomentum12.Pt());
519  if(zone6) fHmggSM_Zone6[iSupMod1]->Fill(invmass,fMomentum12.Pt());
520  if(zone7) fHmggSM_Zone7[iSupMod1]->Fill(invmass,fMomentum12.Pt());
521 
522  }
523  else
524  fHmggDifferentSM ->Fill(invmass,fMomentum12.Pt());
525 
526  // Same sector
527  Int_t j=0;
528  for(Int_t i = 0; i < nSM/2; i++)
529  {
530  j=2*i;
531  if((iSupMod1==j && iSupMod2==j+1) || (iSupMod1==j+1 && iSupMod2==j))
532  {
533  fHmggPairSameSectorSM[i]->Fill(invmass,fMomentum12.Pt());
534  fhClusterPairDiffTimeSameSector[i]->Fill(fMomentum12.E(),time1-time2);
535  }
536  }
537 
538  // Same side
539  for(Int_t i = 0; i < nSM-2; i++)
540  {
541  if((iSupMod1==i && iSupMod2==i+2) || (iSupMod1==i+2 && iSupMod2==i))
542  {
543  fHmggPairSameSideSM[i]->Fill(invmass,fMomentum12.Pt());
544  fhClusterPairDiffTimeSameSide[i]->Fill(fMomentum12.E(),time1-time2);
545  }
546  }
547 
548 
549  if(!mask1 && !mask2)
550  {
551  fHmggMaskFrame->Fill(invmass,fMomentum12.Pt());
552 
553  if(iSupMod1==iSupMod2) fHmggSMMaskFrame[iSupMod1]->Fill(invmass,fMomentum12.Pt());
554  else fHmggDifferentSMMaskFrame ->Fill(invmass,fMomentum12.Pt());
555 
556  // Same sector
557  j=0;
558  for(Int_t i = 0; i < nSM/2; i++)
559  {
560  j=2*i;
561  if((iSupMod1==j && iSupMod2==j+1) || (iSupMod1==j+1 && iSupMod2==j)) fHmggPairSameSectorSMMaskFrame[i]->Fill(invmass,fMomentum12.Pt());
562  }
563 
564  // Same side
565  for(Int_t i = 0; i < nSM-2; i++)
566  {
567  if((iSupMod1==i && iSupMod2==i+2) || (iSupMod1==i+2 && iSupMod2==i)) fHmggPairSameSideSMMaskFrame[i]->Fill(invmass,fMomentum12.Pt());
568  }
569 
570  }// Pair not facing frame
571 
572  if(invmass > fInvMassCutMin && invmass < fInvMassCutMax) //restrict to clusters really close to pi0 peak
573  {
574 
575  // Check time of cells in both clusters, and fill time histogram
576  for(Int_t icell = 0; icell < c1->GetNCells(); icell++)
577  {
578  Int_t absID = c1->GetCellAbsId(icell);
579  fHTpi0[bc%4]->Fill(absID, fEMCALCells->GetCellTime(absID)*1.e9);
580  }
581 
582  for(Int_t icell = 0; icell < c2->GetNCells(); icell++)
583  {
584  Int_t absID = c2->GetCellAbsId(icell);
585  fHTpi0[bc%4]->Fill(absID, fEMCALCells->GetCellTime(absID)*1.e9);
586  }
587 
588  //Opening angle of 2 photons
589  Float_t opangle = fMomentum1.Angle(fMomentum2.Vect())*TMath::RadToDeg();
590  //printf("*******>>>>>>>> In PEAK pt %f, angle %f \n",fMomentum12.Pt(),opangle);
591 
592 
593  fHOpeningAngle ->Fill(opangle,fMomentum12.Pt());
594  fHAsymmetry ->Fill(asym,fMomentum12.Pt());
595 
596  if(iSupMod1==iSupMod2)
597  {
598  fHOpeningAngleSM[iSupMod1] ->Fill(opangle,fMomentum12.Pt());
599  fHAsymmetrySM[iSupMod1] ->Fill(asym,fMomentum12.Pt());
600  }
601  else
602  {
603  fHOpeningAngleDifferentSM ->Fill(opangle,fMomentum12.Pt());
604  fHAsymmetryDifferentSM ->Fill(asym,fMomentum12.Pt());
605  }
606 
607  if((iSupMod1==0 && iSupMod2==2) || (iSupMod1==2 && iSupMod2==0))
608  {
609  fHOpeningAnglePairSM[0] ->Fill(opangle,fMomentum12.Pt());
610  fHAsymmetryPairSM[0] ->Fill(asym,fMomentum12.Pt());
611 
612  }
613  if((iSupMod1==1 && iSupMod2==3) || (iSupMod1==3 && iSupMod2==1))
614  {
615  fHOpeningAnglePairSM[1] ->Fill(opangle,fMomentum12.Pt());
616  fHAsymmetryPairSM[1] ->Fill(asym,fMomentum12.Pt());
617  }
618 
619  if((iSupMod1==0 && iSupMod2==1) || (iSupMod1==1 && iSupMod2==0))
620  {
621  fHOpeningAnglePairSM[2] ->Fill(opangle,fMomentum12.Pt());
622  fHAsymmetryPairSM[2] ->Fill(asym,fMomentum12.Pt());
623  }
624  if((iSupMod1==2 && iSupMod2==3) || (iSupMod1==3 && iSupMod2==2))
625  {
626  fHOpeningAnglePairSM[3] ->Fill(opangle,fMomentum12.Pt());
627  fHAsymmetryPairSM[3] ->Fill(asym,fMomentum12.Pt());
628  }
629 
630  }// pair in 100 < mass < 160
631 
632  }//in acceptance cuts
633 
634  //In case of filling only channels with second cluster in same SM
635  if(fSameSM && iSupMod1!=iSupMod2) continue;
636 
637  if (fGroupNCells == 0)
638  {
639  fHmpi0[iSupMod1][ieta1][iphi1]->Fill(invmass);
640  fHmpi0[iSupMod2][ieta2][iphi2]->Fill(invmass);
641 
642  if (fCellEnergyHiso) fhEnergy[iSupMod1][ieta1][iphi1]->Fill(fMomentum1.E());
643  if (fCellEnergyHiso) fhEnergy[iSupMod2][ieta2][iphi2]->Fill(fMomentum2.E());
644 
645  if(invmass > fInvMassCutMin && invmass < fInvMassCutMax)//restrict to clusters really close to pi0 peak
646  {
647  fhTowerDecayPhotonHit [iSupMod1]->Fill(ieta1,iphi1);
648  fhTowerDecayPhotonEnergy [iSupMod1]->Fill(ieta1,iphi1,fMomentum1.E());
649  fhTowerDecayPhotonAsymmetry[iSupMod1]->Fill(ieta1,iphi1,asym);
650 
651  fhTowerDecayPhotonHit [iSupMod2]->Fill(ieta2,iphi2);
652  fhTowerDecayPhotonEnergy [iSupMod2]->Fill(ieta2,iphi2,fMomentum2.E());
653  fhTowerDecayPhotonAsymmetry[iSupMod2]->Fill(ieta2,iphi2,asym);
654 
655  if(!mask1)fhTowerDecayPhotonHitMaskFrame[iSupMod1]->Fill(ieta1,iphi1);
656  if(!mask2)fhTowerDecayPhotonHitMaskFrame[iSupMod2]->Fill(ieta2,iphi2);
657 
658  }// pair in mass of pi0
659  }
660  else
661  {
662  //printf("Regroup N %d, eta1 %d, phi1 %d, eta2 %d, phi2 %d \n",fGroupNCells, ieta1, iphi1, ieta2, iphi2);
663  for (Int_t i = -fGroupNCells; i < fGroupNCells+1; i++)
664  {
665  for (Int_t j = -fGroupNCells; j < fGroupNCells+1; j++)
666  {
667  Int_t absId11 = fEMCALGeo->GetAbsCellIdFromCellIndexes(iSupMod1, iphi1+j, ieta1+i);
668  Int_t absId22 = fEMCALGeo->GetAbsCellIdFromCellIndexes(iSupMod2, iphi2+j, ieta2+i);
669 
670  Bool_t ok1 = kFALSE;
671  Bool_t ok2 = kFALSE;
672 
673  for(Int_t icell = 0; icell < c1->GetNCells(); icell++)
674  {
675  if(c1->GetCellsAbsId()[icell] == absId11) ok1=kTRUE;
676  }
677 
678  for(Int_t icell = 0; icell < c2->GetNCells(); icell++)
679  {
680  if(c2->GetCellsAbsId()[icell] == absId22) ok2=kTRUE;
681  }
682 
683  if(ok1 && (ieta1+i >= 0) && (iphi1+j >= 0) && (ieta1+i < 48) && (iphi1+j < 24))
684  {
685  fHmpi0[iSupMod1][ieta1+i][iphi1+j]->Fill(invmass);
686  if(fCellEnergyHiso) fhEnergy[iSupMod1][ieta1+i][iphi1+j]->Fill(fMomentum1.E());
687  }
688 
689  if(ok2 && (ieta2+i >= 0) && (iphi2+j >= 0) && (ieta2+i < 48) && (iphi2+j < 24))
690  {
691  fHmpi0[iSupMod2][ieta2+i][iphi2+j]->Fill(invmass);
692  if(fCellEnergyHiso) fhEnergy[iSupMod2][ieta2+i][iphi2+j]->Fill(fMomentum2.E());
693  }
694  }// j loop
695  }//i loop
696  }//group cells
697 
698  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",
699  iSupMod1,iphi1,ieta1,iSupMod2,iphi2,ieta2,fMomentum12.M(),e1i,c1->E(),e2i,c2->E()));
700  }
701  }
702  } // end of loop over EMCAL clusters
703 }
704 
705 
710 //______________________________________________________________________
712 {
713  if ( !fRecoUtils->IsRecalibrationOn() || fCalibFilePath == "" ) return ;
714 
715  if(!fEMCALGeo)fEMCALGeo = AliEMCALGeometry::GetInstance(fEMCALGeoName) ;
716 
717  TFile * calibFactorsFile = TFile::Open(fCalibFilePath.Data());
718 
719  if ( !calibFactorsFile ) AliFatal("Cannot recover the calibration factors");
720 
721  for(Int_t ism = 0; ism < fEMCALGeo->GetNumberOfSuperModules(); ism++)
722  {
723  TH2F * histo = (TH2F*) calibFactorsFile->Get(Form("EMCALRecalFactors_SM%d",ism));
724  printf("In AliAnalysisTaskEMCALPi0CalibSelection::InitEnergyCalibrationFactors \n ---Recover calibration factor for : EMCALRecalFactors_SM%d %p\n",ism,histo);
725 
726  if ( histo )
727  fRecoUtils->SetEMCALChannelRecalibrationFactors(ism,histo);
728  else
729  AliWarning(Form("Null histogram with calibration factors for SM%d, 1 will be used for the full SM!",ism));
730  }
731 }
732 
736 //________________________________________________________________
738 {
739  Int_t runnumber = InputEvent()->GetRunNumber() ;
740 
741  //
742  // Load default geo matrices if requested
743  if(fImportGeometryFromFile && !gGeoManager)
744  {
745  if(fImportGeometryFilePath=="") // If not specified, set location depending on run number
746  {
747  // "$ALICE_ROOT/EVE/alice-data/default_geo.root"
748  if (runnumber < 140000) fImportGeometryFilePath = "$ALICE_PHYSICS/OADB/EMCAL/geometry_2010.root";
749  else if(runnumber < 171000) fImportGeometryFilePath = "$ALICE_PHYSICS/OADB/EMCAL/geometry_2011.root";
750  else if(runnumber < 198000) fImportGeometryFilePath = "$ALICE_PHYSICS/OADB/EMCAL/geometry_2012.root"; // 2012-2013
751  else fImportGeometryFilePath = "$ALICE_PHYSICS/OADB/EMCAL/geometry_2015.root"; // >= 2015
752  }
753 
754  AliInfo(Form("Import %s",fImportGeometryFilePath.Data()));
755 
756  TGeoManager::Import(fImportGeometryFilePath) ; // default need file "geometry.root" in local dir!!!!
757  }
758 
759  //
760  if(fLoadMatrices)
761  {
762  AliInfo("Load user defined EMCAL geometry matrices");
763  // OADB if available
764  AliOADBContainer emcGeoMat("AliEMCALgeo");
765 
766  if(fOADBFilePath=="") fOADBFilePath = "$ALICE_PHYSICS/OADB/EMCAL" ;
767 
768  emcGeoMat.InitFromFile(Form("%s/EMCALlocal2master.root",fOADBFilePath.Data()),"AliEMCALgeo");
769 
770  TObjArray *matEMCAL=(TObjArray*)emcGeoMat.GetObject(runnumber,"EmcalMatrices");
771 
772  for(Int_t mod = 0; mod < (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules(); mod++)
773  {
774  if (!fMatrix[mod]) // Get it from OADB
775  {
776  AliDebug(1,Form("EMCAL matrices SM %d, %p",mod,((TGeoHMatrix*) matEMCAL->At(mod))));
777  //((TGeoHMatrix*) matEMCAL->At(mod))->Print();
778 
779  fMatrix[mod] = (TGeoHMatrix*) matEMCAL->At(mod) ;
780  }
781 
782  if(fMatrix[mod])
783  {
784  if(DebugLevel() > 1)
785  fMatrix[mod]->Print();
786 
787  fEMCALGeo->SetMisalMatrix(fMatrix[mod],mod) ;
788  }
789  else if(gGeoManager)
790  {
791  AliWarning(Form("Set matrix for SM %d from gGeoManager",mod));
792  fEMCALGeo->SetMisalMatrix(fEMCALGeo->GetMatrixForSuperModuleFromGeoManager(mod),mod) ;
793  }
794  else
795  {
796  AliError(Form("Alignment matrix for SM %d is not available",mod));
797  }
798  }//SM loop
799  }//Load matrices
800  else if(!gGeoManager)
801  {
802  AliInfo("Get geo matrices from data");
803  //Still not implemented in AOD, just a workaround to be able to work at least with ESDs
804  if(!strcmp(InputEvent()->GetName(),"AliAODEvent"))
805  {
806  AliWarning("Use ideal geometry, values geometry matrix not kept in AODs");
807  }//AOD
808  else
809  {
810  AliDebug(1,"AliAnalysisTaskEMCALClusterize Load Misaligned matrices");
811 
812  for(Int_t mod=0; mod < (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules(); mod++)
813  {
814  if(InputEvent()->GetEMCALMatrix(mod))
815  {
816  if(DebugLevel() > 1)
817  InputEvent()->GetEMCALMatrix(mod)->Print();
818 
819  fEMCALGeo->SetMisalMatrix(InputEvent()->GetEMCALMatrix(mod),mod) ;
820  }
821 
822  }
823  }// ESD
824  }// Load matrices from Data
825  else if(gGeoManager) // Load default matrices
826  {
827  for(Int_t mod = 0; mod < (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules(); mod++)
828  {
829  AliWarning(Form("Set matrix for SM %d from gGeoManager",mod));
830  fEMCALGeo->SetMisalMatrix(fEMCALGeo->GetMatrixForSuperModuleFromGeoManager(mod),mod) ;
831  }
832  } // gGeoManager matrices
833 
834 }
835 
839 //______________________________________________________________________
841 {
842  if(!fRecoUtils->IsRunDepRecalibrationOn()) return;
843 
844  AliOADBContainer *contRFTD=new AliOADBContainer("");
845 
846  contRFTD->InitFromFile(Form("%s/EMCALTemperatureCorrCalib.root",fOADBFilePath.Data()),"AliEMCALRunDepTempCalibCorrections");
847 
848  Int_t runnumber = InputEvent()->GetRunNumber() ;
849 
850  TH1S *htd=(TH1S*)contRFTD->GetObject(runnumber);
851 
852  //If it did not exist for this run, get closes one
853  if (!htd)
854  {
855  AliWarning(Form("No TemperatureCorrCalib Objects for run: %d",runnumber));
856 
857  // let's get the closest runnumber instead then..
858  Int_t lower = 0;
859  Int_t ic = 0;
860  Int_t maxEntry = contRFTD->GetNumberOfEntries();
861 
862  while ( (ic < maxEntry) && (contRFTD->UpperLimit(ic) < runnumber) )
863  {
864  lower = ic;
865  ic++;
866  }
867 
868  Int_t closest = lower;
869  if ( (ic<maxEntry) &&
870  (contRFTD->LowerLimit(ic)-runnumber) < (runnumber - contRFTD->UpperLimit(lower)) )
871  {
872  closest = ic;
873  }
874 
875  AliWarning(Form("TemperatureCorrCalib Objects found closest id %d from run: %d",
876  closest, contRFTD->LowerLimit(closest)));
877 
878  htd = (TH1S*) contRFTD->GetObjectByIndex(closest);
879  }
880 
881  // Fill parameters
882  if(htd)
883  {
884  AliInfo("Recalibrate (Temperature) EMCAL");
885 
886  Int_t nSM = fEMCALGeo->GetNumberOfSuperModules();
887 
888  for (Int_t ism = 0; ism < nSM; ++ism)
889  {
890  for (Int_t icol = 0; icol < 48; ++icol)
891  {
892  for (Int_t irow = 0; irow < 24; ++irow)
893  {
894  Float_t factor = fRecoUtils->GetEMCALChannelRecalibrationFactor(ism,icol,irow);
895 
896  Int_t absID = fEMCALGeo->GetAbsCellIdFromCellIndexes(ism, irow, icol); // original calibration factor
897 
898  AliDebug(3,Form(" ism %d, icol %d, irow %d,absID %d - Calib factor %1.5f - ",ism, icol, irow, absID, factor));
899 
900  factor *= htd->GetBinContent(absID) / 10000. ; // correction dependent on T
901 
902  fRecoUtils->SetEMCALChannelRecalibrationFactor(ism,icol,irow,factor);
903 
904  AliDebug(3,Form("T factor %1.5f - final factor %1.5f",
905  htd->GetBinContent(absID) / 10000.,
906  fRecoUtils->GetEMCALChannelRecalibrationFactor(ism,icol,irow)));
907  } // columns
908  } // rows
909  } // SM loop
910  }
911  else AliInfo("Do NOT recalibrate EMCAL with T variations, no params TH1");
912 
913  delete contRFTD;
914 }
915 
916 
919 //___________________________________________________________________
921 {
922  if(!fEMCALGeo)fEMCALGeo = AliEMCALGeometry::GetInstance(fEMCALGeoName) ;
923  Int_t nSM = (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules();
924 
925  fOutputContainer = new TList();
926  const Int_t buffersize = 255;
927  char hname[buffersize], htitl[buffersize], htitlEnergy[buffersize];
928 
929  fhNEvents = new TH1I("hNEvents", "Number of analyzed events" , 1 , 0 , 1 ) ;
931 
932  fHmgg = new TH2F("hmgg","2-cluster invariant mass",fNbins,fMinBin,fMaxBin,100,0,10);
933  fHmgg->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
934  fHmgg->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
935  fOutputContainer->Add(fHmgg);
936 
937  fHmggDifferentSM = new TH2F("hmggDifferentSM","2-cluster invariant mass, different SM",fNbins,fMinBin,fMaxBin,100,0,10);
938  fHmggDifferentSM->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
939  fHmggDifferentSM->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
941 
942  fHOpeningAngle = new TH2F("hopang","2-cluster opening angle",100,0.,50.,100,0,10);
943  fHOpeningAngle->SetXTitle("#alpha_{#gamma #gamma}");
944  fHOpeningAngle->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
946 
947  fHOpeningAngleDifferentSM = new TH2F("hopangDifferentSM","2-cluster opening angle, different SM",100,0,50.,100,0,10);
948  fHOpeningAngleDifferentSM->SetXTitle("#alpha_{#gamma #gamma}");
949  fHOpeningAngleDifferentSM->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
951 
952  fHAsymmetry = new TH2F("hasym","2-cluster opening angle",100,0.,1.,100,0,10);
953  fHAsymmetry->SetXTitle("a");
954  fHAsymmetry->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
956 
957  fHAsymmetryDifferentSM = new TH2F("hasymDifferentSM","2-cluster opening angle, different SM",100,0,1.,100,0,10);
958  fHAsymmetryDifferentSM->SetXTitle("a");
959  fHAsymmetryDifferentSM->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
961 
962  //TString pairname[] = {"A side (0-2)", "C side (1-3)","Row 0 (0-1)", "Row 1 (2-3)"};
963 
964  fHmggMaskFrame = new TH2F("hmggMaskFrame","2-cluster invariant mass, frame masked",fNbins,fMinBin,fMaxBin,100,0,10);
965  fHmggMaskFrame->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
966  fHmggMaskFrame->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
968 
969  fHmggDifferentSMMaskFrame = new TH2F("hmggDifferentSMMaskFrame","2-cluster invariant mass, different SM, frame masked",
970  fNbins,fMinBin,fMaxBin,100,0,10);
971  fHmggDifferentSMMaskFrame->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
972  fHmggDifferentSMMaskFrame->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
974 
975  for(Int_t iSM = 0; iSM < nSM; iSM++)
976  {
977  snprintf(hname, buffersize, "hmgg_SM%d",iSM);
978  snprintf(htitl, buffersize, "Two-gamma inv. mass for super mod %d",iSM);
979  fHmggSM[iSM] = new TH2F(hname,htitl,fNbins,fMinBin,fMaxBin,100,0,10);
980  fHmggSM[iSM]->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
981  fHmggSM[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
982  fOutputContainer->Add(fHmggSM[iSM]);
983 
984  snprintf(hname, buffersize, "hmgg_SM%d_MaskFrame",iSM);
985  snprintf(htitl, buffersize, "Two-gamma inv. mass for super mod %d",iSM);
986  fHmggSMMaskFrame[iSM] = new TH2F(hname,htitl,fNbins,fMinBin,fMaxBin,100,0,10);
987  fHmggSMMaskFrame[iSM]->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
988  fHmggSMMaskFrame[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
990 
991  snprintf(hname, buffersize, "hmgg_SM%d_Zone1",iSM);
992  snprintf(htitl, buffersize, "Two-gamma inv. mass for super mod %d in zone 1",iSM);
993  fHmggSM_Zone1[iSM] = new TH2F(hname,htitl,fNbins,fMinBin,fMaxBin,100,0,10);
994  fHmggSM_Zone1[iSM]->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
995  fHmggSM_Zone1[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
996  fOutputContainer->Add(fHmggSM_Zone1[iSM]);
997 
998  snprintf(hname, buffersize, "hmgg_SM%d_Zone2",iSM);
999  snprintf(htitl, buffersize, "Two-gamma inv. mass for super mod %d in zone 2",iSM);
1000  fHmggSM_Zone2[iSM] = new TH2F(hname,htitl,fNbins,fMinBin,fMaxBin,100,0,10);
1001  fHmggSM_Zone2[iSM]->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
1002  fHmggSM_Zone2[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
1003  fOutputContainer->Add(fHmggSM_Zone2[iSM]);
1004 
1005  snprintf(hname, buffersize, "hmgg_SM%d_Zone3",iSM);
1006  snprintf(htitl, buffersize, "Two-gamma inv. mass for super mod %d in zone 3",iSM);
1007  fHmggSM_Zone3[iSM] = new TH2F(hname,htitl,fNbins,fMinBin,fMaxBin,100,0,10);
1008  fHmggSM_Zone3[iSM]->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
1009  fHmggSM_Zone3[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
1010  fOutputContainer->Add(fHmggSM_Zone3[iSM]);
1011 
1012  snprintf(hname, buffersize, "hmgg_SM%d_Zone4",iSM);
1013  snprintf(htitl, buffersize, "Two-gamma inv. mass for super mod %d in zone 4",iSM);
1014  fHmggSM_Zone4[iSM] = new TH2F(hname,htitl,fNbins,fMinBin,fMaxBin,100,0,10);
1015  fHmggSM_Zone4[iSM]->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
1016  fHmggSM_Zone4[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
1017  fOutputContainer->Add(fHmggSM_Zone4[iSM]);
1018 
1019  snprintf(hname, buffersize, "hmgg_SM%d_Zone5",iSM);
1020  snprintf(htitl, buffersize, "Two-gamma inv. mass for super mod %d in zone 5",iSM);
1021  fHmggSM_Zone5[iSM] = new TH2F(hname,htitl,fNbins,fMinBin,fMaxBin,100,0,10);
1022  fHmggSM_Zone5[iSM]->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
1023  fHmggSM_Zone5[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
1024  fOutputContainer->Add(fHmggSM_Zone5[iSM]);
1025 
1026  snprintf(hname, buffersize, "hmgg_SM%d_Zone6",iSM);
1027  snprintf(htitl, buffersize, "Two-gamma inv. mass for super mod %d in zone 6",iSM);
1028  fHmggSM_Zone6[iSM] = new TH2F(hname,htitl,fNbins,fMinBin,fMaxBin,100,0,10);
1029  fHmggSM_Zone6[iSM]->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
1030  fHmggSM_Zone6[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
1031  fOutputContainer->Add(fHmggSM_Zone6[iSM]);
1032 
1033  snprintf(hname, buffersize, "hmgg_SM%d_Zone7",iSM);
1034  snprintf(htitl, buffersize, "Two-gamma inv. mass for super mod %d in zone 7",iSM);
1035  fHmggSM_Zone7[iSM] = new TH2F(hname,htitl,fNbins,fMinBin,fMaxBin,100,0,10);
1036  fHmggSM_Zone7[iSM]->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
1037  fHmggSM_Zone7[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
1038  fOutputContainer->Add(fHmggSM_Zone7[iSM]);
1039 
1040  if(iSM < nSM/2)
1041  {
1042  snprintf(hname,buffersize, "hmgg_PairSameSectorSM%d",iSM);
1043  snprintf(htitl,buffersize, "Two-gamma inv. mass for SM pair Sector: %d",iSM);
1044  fHmggPairSameSectorSM[iSM] = new TH2F(hname,htitl,fNbins,fMinBin,fMaxBin,100,0,10);
1045  fHmggPairSameSectorSM[iSM]->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
1046  fHmggPairSameSectorSM[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
1048 
1049  snprintf(hname,buffersize, "hmgg_PairSameSectorSM%d_MaskFrame",iSM);
1050  snprintf(htitl,buffersize, "Two-gamma inv. mass for SM pair Sector: %d",iSM);
1051  fHmggPairSameSectorSMMaskFrame[iSM] = new TH2F(hname,htitl,fNbins,fMinBin,fMaxBin,100,0,10);
1052  fHmggPairSameSectorSMMaskFrame[iSM]->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
1053  fHmggPairSameSectorSMMaskFrame[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
1055 
1056  fhClusterPairDiffTimeSameSector[iSM] = new TH2F(Form("hClusterPairDiffTimeSameSector%d",iSM),
1057  Form("cluster pair time difference vs E, Sector %d",iSM),
1058  100,0,10, 200,-100,100);
1059  fhClusterPairDiffTimeSameSector[iSM]->SetXTitle("E_{pair} (GeV)");
1060  fhClusterPairDiffTimeSameSector[iSM]->SetYTitle("#Delta t (ns)");
1062  }
1063 
1064  if(iSM < nSM-2)
1065  {
1066  snprintf(hname,buffersize, "hmgg_PairSameSideSM%d",iSM);
1067  snprintf(htitl,buffersize, "Two-gamma inv. mass for SM pair Sector: %d",iSM);
1068  fHmggPairSameSideSM[iSM] = new TH2F(hname,htitl,fNbins,fMinBin,fMaxBin,100,0,10);
1069  fHmggPairSameSideSM[iSM]->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
1070  fHmggPairSameSideSM[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
1072 
1073  snprintf(hname,buffersize, "hmgg_PairSameSideSM%d_MaskFrame",iSM);
1074  snprintf(htitl,buffersize, "Two-gamma inv. mass for SM pair Sector: %d",iSM);
1075  fHmggPairSameSideSMMaskFrame[iSM] = new TH2F(hname,htitl,fNbins,fMinBin,fMaxBin,100,0,10);
1076  fHmggPairSameSideSMMaskFrame[iSM]->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
1077  fHmggPairSameSideSMMaskFrame[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
1079 
1080  fhClusterPairDiffTimeSameSide[iSM] = new TH2F(Form("hClusterPairDiffTimeSameSide%d",iSM),
1081  Form("cluster pair time difference vs E, Side %d",iSM),
1082  100,0,10, 200,-100,100);
1083  fhClusterPairDiffTimeSameSide[iSM]->SetXTitle("E_{pair} (GeV)");
1084  fhClusterPairDiffTimeSameSide[iSM]->SetYTitle("#Delta t (ns)");
1086  }
1087 
1088  snprintf(hname, buffersize, "hopang_SM%d",iSM);
1089  snprintf(htitl, buffersize, "Opening angle for super mod %d",iSM);
1090  fHOpeningAngleSM[iSM] = new TH2F(hname,htitl,100,0.,50.,100,0,10);
1091  fHOpeningAngleSM[iSM]->SetXTitle("#alpha_{#gamma #gamma} (deg)");
1092  fHOpeningAngleSM[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
1094 
1095  snprintf(hname,buffersize, "hopang_PairSM%d",iSM);
1096  snprintf(htitl,buffersize, "Opening angle for SM pair: %d",iSM);
1097  fHOpeningAnglePairSM[iSM] = new TH2F(hname,htitl,100,0.,50.,100,0,10);
1098  fHOpeningAnglePairSM[iSM]->SetXTitle("#alpha_{#gamma #gamma} (deg)");
1099  fHOpeningAnglePairSM[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
1101 
1102  snprintf(hname, buffersize, "hasym_SM%d",iSM);
1103  snprintf(htitl, buffersize, "Asymmetry for super mod %d",iSM);
1104  fHAsymmetrySM[iSM] = new TH2F(hname,htitl,100,0.,1.,100,0,10);
1105  fHAsymmetrySM[iSM]->SetXTitle("a");
1106  fHAsymmetrySM[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
1107  fOutputContainer->Add(fHAsymmetrySM[iSM]);
1108 
1109  snprintf(hname,buffersize, "hasym_PairSM%d",iSM);
1110  snprintf(htitl,buffersize, "Asymmetry for SM pair: %d",iSM);
1111  fHAsymmetryPairSM[iSM] = new TH2F(hname,htitl,100,0.,1.,100,0,10);
1112  fHAsymmetryPairSM[iSM]->SetXTitle("a");
1113  fHAsymmetryPairSM[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
1114  fOutputContainer->Add(fHAsymmetryPairSM[iSM]);
1115 
1116  Int_t colmax = 48;
1117  Int_t rowmax = 24;
1118 
1119  fhTowerDecayPhotonHit[iSM] = new TH2F (Form("hTowerDecPhotonHit_Mod%d",iSM),
1120  Form("Entries in grid of cells in Module %d",iSM),
1121  colmax+2,-1.5,colmax+0.5, rowmax+2,-1.5,rowmax+0.5);
1122  fhTowerDecayPhotonHit[iSM]->SetYTitle("row (phi direction)");
1123  fhTowerDecayPhotonHit[iSM]->SetXTitle("column (eta direction)");
1125 
1126  fhTowerDecayPhotonEnergy[iSM] = new TH2F (Form("hTowerDecPhotonEnergy_Mod%d",iSM),
1127  Form("Accumulated energy in grid of cells in Module %d",iSM),
1128  colmax+2,-1.5,colmax+0.5, rowmax+2,-1.5,rowmax+0.5);
1129  fhTowerDecayPhotonEnergy[iSM]->SetYTitle("row (phi direction)");
1130  fhTowerDecayPhotonEnergy[iSM]->SetXTitle("column (eta direction)");
1132 
1133  fhTowerDecayPhotonAsymmetry[iSM] = new TH2F (Form("hTowerDecPhotonAsymmetry_Mod%d",iSM),
1134  Form("Accumulated asymmetry in grid of cells in Module %d",iSM),
1135  colmax+2,-1.5,colmax+0.5, rowmax+2,-1.5,rowmax+0.5);
1136  fhTowerDecayPhotonAsymmetry[iSM]->SetYTitle("row (phi direction)");
1137  fhTowerDecayPhotonAsymmetry[iSM]->SetXTitle("column (eta direction)");
1139 
1140  fhTowerDecayPhotonHitMaskFrame[iSM] = new TH2F (Form("hTowerDecPhotonHit_Mod%d_MaskFrame",iSM),Form("Entries in grid of cells in Module %d",iSM),
1141  colmax+2,-1.5,colmax+0.5, rowmax+2,-1.5,rowmax+0.5);
1142  fhTowerDecayPhotonHitMaskFrame[iSM]->SetYTitle("row (phi direction)");
1143  fhTowerDecayPhotonHitMaskFrame[iSM]->SetXTitle("column (eta direction)");
1145 
1146  fhClusterTimeSM[iSM] = new TH2F(Form("hClusterTime_SM%d",iSM),"cluster time vs E",100,0,10, 200,-1000,1000);
1147  fhClusterTimeSM[iSM]->SetXTitle("E (GeV)");
1148  fhClusterTimeSM[iSM]->SetYTitle("t (ns)");
1149  fOutputContainer->Add(fhClusterTimeSM[iSM]);
1150 
1151  fhClusterPairDiffTimeSameSM[iSM] = new TH2F(Form("hClusterPairDiffTimeSameSM%d",iSM),
1152  Form("cluster pair time difference vs E, SM %d",iSM),
1153  100,0,10, 200,-100,100);
1154  fhClusterPairDiffTimeSameSM[iSM]->SetXTitle("E (GeV)");
1155  fhClusterPairDiffTimeSameSM[iSM]->SetYTitle("#Delta t (ns)");
1157 
1158 
1159  if(fClusterTopology)
1160  {
1161 
1162  fhTopoClusterAmpCase0[iSM] = new TH2F(Form("hTopoClusterAmpCase0SM%d",iSM),
1163  Form("cluster topology for cluster in position 0 in noisy quartet, SM %d",iSM),
1164  20,-10,10, 20,-10,10);
1165  fhTopoClusterAmpCase0[iSM]->SetXTitle("column");
1166  fhTopoClusterAmpCase0[iSM]->SetYTitle("row");
1168 
1169  fhTopoClusterAmpCase1[iSM] = new TH2F(Form("hTopoClusterAmpCase1SM%d",iSM),
1170  Form("cluster topology for cluster in position 1 in noisy quartet, SM %d",iSM),
1171  20,-10,10, 20,-10,10);
1172  fhTopoClusterAmpCase1[iSM]->SetXTitle("column");
1173  fhTopoClusterAmpCase1[iSM]->SetYTitle("row");
1175 
1176  fhTopoClusterAmpCase2[iSM] = new TH2F(Form("hTopoClusterAmpCase2SM%d",iSM),
1177  Form("cluster topology for cluster in position 2 in noisy quartet, SM %d",iSM),
1178  20,-10,10, 20,-10,10);
1179  fhTopoClusterAmpCase2[iSM]->SetXTitle("column");
1180  fhTopoClusterAmpCase2[iSM]->SetYTitle("row");
1182 
1183  fhTopoClusterAmpCase3[iSM] = new TH2F(Form("hTopoClusterAmpCase3SM%d",iSM),
1184  Form("cluster topology for cluster in position 3 in noisy quartet, SM %d",iSM),
1185  20,-10,10, 20,-10,10);
1186  fhTopoClusterAmpCase3[iSM]->SetXTitle("column");
1187  fhTopoClusterAmpCase3[iSM]->SetYTitle("row");
1189 
1190 
1191  fhTopoClusterAmpFractionCase0[iSM] = new TH2F(Form("hTopoClusterAmpFractionCase0SM%d",iSM),
1192  Form("cluster topology for cluster in position 0 in noisy quartet, SM %d",iSM),
1193  20,-10,10, 20,-10,10);
1194  fhTopoClusterAmpFractionCase0[iSM]->SetXTitle("column");
1195  fhTopoClusterAmpFractionCase0[iSM]->SetYTitle("row");
1197 
1198  fhTopoClusterAmpFractionCase1[iSM] = new TH2F(Form("hTopoClusterAmpFractionCase1SM%d",iSM),
1199  Form("cluster topology for cluster in position 1 in noisy quartet, SM %d",iSM),
1200  20,-10,10, 20,-10,10);
1201  fhTopoClusterAmpFractionCase1[iSM]->SetXTitle("column");
1202  fhTopoClusterAmpFractionCase1[iSM]->SetYTitle("row");
1204 
1205  fhTopoClusterAmpFractionCase2[iSM] = new TH2F(Form("hTopoClusterAmpFractionCase2SM%d",iSM),
1206  Form("cluster topology for cluster in position 2 in noisy quartet, SM %d",iSM),
1207  20,-10,10, 20,-10,10);
1208  fhTopoClusterAmpFractionCase2[iSM]->SetXTitle("column");
1209  fhTopoClusterAmpFractionCase2[iSM]->SetYTitle("row");
1211 
1212  fhTopoClusterAmpFractionCase3[iSM] = new TH2F(Form("hTopoClusterAmpFractionCase3SM%d",iSM),
1213  Form("cluster topology for cluster in position 3 in noisy quartet, SM %d",iSM),
1214  20,-10,10, 20,-10,10);
1215  fhTopoClusterAmpFractionCase3[iSM]->SetXTitle("column");
1216  fhTopoClusterAmpFractionCase3[iSM]->SetYTitle("row");
1218  }
1219  }
1220 
1221  Int_t nchannels = nSM*AliEMCALGeoParams::fgkEMCALRows*AliEMCALGeoParams::fgkEMCALCols;
1222  for(Int_t ibc = 0; ibc < 4; ibc++)
1223  {
1224  fHTpi0[ibc] = new TH2F(Form("hTime_BC%d",ibc),Form("Time of cell clusters under pi0 peak, bunch crossing %d",ibc),
1225  nchannels,0,nchannels, fNTimeBins,fMinTimeBin,fMaxTimeBin);
1226  fOutputContainer->Add(fHTpi0[ibc]);
1227  fHTpi0[ibc]->SetYTitle("time (ns)");
1228  fHTpi0[ibc]->SetXTitle("abs. Id. ");
1229  }
1230 
1231  fhClusterTime = new TH2F("hClusterTime","cluster time vs E",100,0,10, 100,0,1000);
1232  fhClusterTime->SetXTitle("E (GeV)");
1233  fhClusterTime->SetYTitle("t (ns)");
1235 
1236  fhClusterPairDiffTime = new TH2F("hClusterPairDiffTime","cluster pair time difference vs E",100,0,10, 800,-400,400);
1237  fhClusterPairDiffTime->SetXTitle("E_{pair} (GeV)");
1238  fhClusterPairDiffTime->SetYTitle("#Delta t (ns)");
1240 
1241  for(Int_t iMod=0; iMod < nSM; iMod++)
1242  {
1243  for(Int_t iRow=0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++)
1244  {
1245  for(Int_t iCol=0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++)
1246  {
1247  snprintf(hname,buffersize, "%d_%d_%d",iMod,iCol,iRow);
1248  snprintf(htitl,buffersize, "Two-gamma inv. mass for super mod %d, cell(col,row)=(%d,%d)",iMod,iCol,iRow);
1249  fHmpi0[iMod][iCol][iRow] = new TH1F(hname,htitl,fNbins,fMinBin,fMaxBin);
1250  fHmpi0[iMod][iCol][iRow]->SetXTitle("mass (MeV/c^{2})");
1251  fOutputContainer->Add(fHmpi0[iMod][iCol][iRow]);
1252 
1253  if(fCellEnergyHiso)
1254  {
1255  snprintf(htitlEnergy,buffersize, "Energy for super mod %d, cell(col,row)=(%d,%d)",iMod,iCol,iRow);
1256  fhEnergy[iMod][iCol][iRow] = new TH1F(hname,htitlEnergy,fNEnergybins,fMinEnergyBin,fMaxEnergyBin);
1257  fhEnergy[iMod][iCol][iRow]->SetXTitle("E (GeV)");
1258  fOutputContainer->Add(fhEnergy[iMod][iCol][iRow]);
1259  }
1260  }
1261  }
1262  }
1263 
1264  fOutputContainer->SetOwner(kTRUE);
1265 
1266  PostData(1,fOutputContainer);
1267 
1268 // // cuts container, set in terminate but init and post here
1269 //
1270 // fCuts = new TList();
1271 //
1272 // fCuts ->SetOwner(kTRUE);
1273 //
1274 // PostData(2, fCuts);
1275 }
1276 
1282 //______________________________________________________________________________________________________
1284 {
1285  Int_t icol = ieta;
1286  if(iSM%2) icol+=48; // Impair SM, shift index [0-47] to [48-96]
1287 
1289  {
1290  for (Int_t imask = 0; imask < fNMaskCellColumns; imask++)
1291  {
1292  if(icol==fMaskCellColumns[imask]) return kTRUE;
1293  }
1294  }
1295 
1296  return kFALSE;
1297 }
1298 
1299 
1308 //______________________________________________________________________________________________________
1309 Bool_t AliAnalysisTaskEMCALPi0CalibSelection::IsInZone1(Int_t iSupMod, Int_t ieta, Int_t iphi)
1310 {
1311  Int_t icol = ieta;
1312  Int_t irow = iphi;
1313 
1314 // printf("SM = %i\n",iSupMod);
1315 // printf("icol = %i\n",icol);
1316 // printf("irow = %i\n",irow);
1317 
1318  if(iSupMod%2) //Odd SM
1319  {
1320  if((irow >= 2 && irow <= 21) && ((icol >= 42 && icol <= 46) || (icol >= 13 && icol <= 37) || (icol >= 1 && icol <= 9)))
1321  {
1322  if((irow >= 2 && irow <= 3) || (irow >= 20 && irow <= 21))
1323  {
1324 // printf("zone 1\n");
1325  return kTRUE;
1326  }
1327  }
1328  }
1329  else //Even SM
1330  {
1331  if((irow >= 2 && irow <= 21) && ((icol >= 1 && icol <= 5) || (icol >= 10 && icol <= 34) || (icol >= 38 && icol <= 46)))
1332  {
1333  if((irow >= 2 && irow <= 3) || (irow >= 20 && irow <= 21))
1334  {
1335 // printf("zone 1\n");
1336  return kTRUE;
1337  }
1338  }
1339  }
1340 
1341  return kFALSE;
1342 }
1343 
1344 
1353 //______________________________________________________________________________________________________
1354 Bool_t AliAnalysisTaskEMCALPi0CalibSelection::IsInZone2(Int_t iSupMod, Int_t ieta, Int_t iphi)
1355 {
1356  Int_t icol = ieta;
1357  Int_t irow = iphi;
1358 
1359 // printf("SM = %i\n",iSupMod);
1360 // printf("icol = %i\n",icol);
1361 // printf("irow = %i\n",irow);
1362 
1363  if(iSupMod%2) //Odd SM
1364  {
1365  if((irow >= 2 && irow <= 21) && ((icol >= 42 && icol <= 46) || (icol >= 13 && icol <= 37) || (icol >= 1 && icol <= 9)))
1366  {
1367  if((irow >= 2 && irow <= 3) || (irow >= 20 && irow <= 21))
1368  {
1369  return kFALSE;
1370  }
1371  else
1372  {
1373 // printf("zone 2\n");
1374  return kTRUE;
1375  }
1376  }
1377  }
1378  else //Even SM
1379  {
1380  if((irow >= 2 && irow <= 21) && ((icol >= 1 && icol <= 5) || (icol >= 10 && icol <= 34) || (icol >= 38 && icol <= 46)))
1381  {
1382  if((irow >= 2 && irow <= 3) || (irow >= 20 && irow <= 21))
1383  {
1384  return kFALSE;
1385  }
1386  else
1387  {
1388 // printf("zone 2\n");
1389  return kTRUE;
1390  }
1391  }
1392  }
1393 
1394  return kFALSE;
1395 }
1396 
1397 
1406 //______________________________________________________________________________________________________
1407 Bool_t AliAnalysisTaskEMCALPi0CalibSelection::IsInZone3(Int_t iSupMod, Int_t ieta, Int_t iphi)
1408 {
1409  Int_t icol = ieta;
1410  Int_t irow = iphi;
1411 
1412 // printf("SM = %i\n",iSupMod);
1413 // printf("icol = %i\n",icol);
1414 // printf("irow = %i\n",irow);
1415 
1416  if(iSupMod%2) //Odd SM
1417  {
1418  if((irow >= 2 && irow <= 21) && ((icol >= 42 && icol <= 46) || (icol >= 13 && icol <= 37) || (icol >= 1 && icol <= 9)))
1419  {
1420  if((icol >= 1 && icol <= 3) || (icol >= 44 && icol <= 46))
1421  {
1422 // printf("zone 3\n");
1423  return kTRUE;
1424  }
1425  }
1426  }
1427  else //Even SM
1428  {
1429  if((irow >= 2 && irow <= 21) && ((icol >= 1 && icol <= 5) || (icol >= 10 && icol <= 34) || (icol >= 38 && icol <= 46)))
1430  {
1431  if((icol >= 1 && icol <= 3) || (icol >= 44 && icol <= 46))
1432  {
1433 // printf("zone 3\n");
1434  return kTRUE;
1435  }
1436  }
1437  }
1438 
1439  return kFALSE;
1440 }
1441 
1442 
1451 //______________________________________________________________________________________________________
1452 Bool_t AliAnalysisTaskEMCALPi0CalibSelection::IsInZone4(Int_t iSupMod, Int_t ieta, Int_t iphi)
1453 {
1454  Int_t icol = ieta;
1455  Int_t irow = iphi;
1456 
1457 // printf("SM = %i\n",iSupMod);
1458 // printf("icol = %i\n",icol);
1459 // printf("irow = %i\n",irow);
1460 
1461  if(iSupMod%2) //Odd SM
1462  {
1463  if((irow >= 2 && irow <= 21) && ((icol >= 42 && icol <= 46) || (icol >= 13 && icol <= 37) || (icol >= 1 && icol <= 9)))
1464  {
1465  if((icol >= 1 && icol <= 3) || (icol >= 44 && icol <= 46))
1466  {
1467  return kFALSE;
1468  }
1469  else
1470  {
1471 // printf("zone 4\n");
1472  return kTRUE;
1473  }
1474  }
1475  }
1476  else //Even SM
1477  {
1478  if((irow >= 2 && irow <= 21) && ((icol >= 1 && icol <= 5) || (icol >= 10 && icol <= 34) || (icol >= 38 && icol <= 46)))
1479  {
1480  if((icol >= 1 && icol <= 3) || (icol >= 44 && icol <= 46))
1481  {
1482  return kFALSE;
1483  }
1484  else
1485  {
1486 // printf("zone 4\n");
1487  return kTRUE;
1488  }
1489  }
1490  }
1491 
1492  return kFALSE;
1493 }
1494 
1495 
1505 //______________________________________________________________________________________________________
1506 Bool_t AliAnalysisTaskEMCALPi0CalibSelection::IsInZone5(Int_t iSupMod, Int_t ieta, Int_t iphi)
1507 {
1508  Int_t icol = ieta;
1509  Int_t irow = iphi;
1510 
1511  //Center of ellipse
1512  Float_t col0 = 47/2;
1513  Float_t row0 = 23/2;
1514 
1515  //Parameters
1516  Float_t a = 3-col0;
1517  Float_t b = 2-row0;
1518 
1519  if(((icol-col0)*(icol-col0)) / (a*a) + ((irow-row0)*(irow-row0) / (b*b)) > 1)
1520  {
1521  return kTRUE;
1522  }
1523  else
1524  {
1525  return kFALSE;
1526  }
1527 }
1528 
1529 
1539 //______________________________________________________________________________________________________
1540 Bool_t AliAnalysisTaskEMCALPi0CalibSelection::IsInZone6(Int_t iSupMod, Int_t ieta, Int_t iphi)
1541 {
1542  Int_t icol = ieta;
1543  Int_t irow = iphi;
1544 
1545  //Center of ellipse
1546  Float_t col0 = 47/2;
1547  Float_t row0 = 23/2;
1548 
1549  //Paramters
1550  Float_t aLarge = 3-col0;
1551  Float_t bLarge = 2-row0;
1552  Float_t aSmall = 16-col0;
1553  Float_t bSmall = 7-row0;
1554 
1555  if((((icol-col0)*(icol-col0)) / (aLarge*aLarge) + ((irow-row0)*(irow-row0) / (bLarge*bLarge)) < 1) && (((icol-col0)*(icol-col0)) / (aSmall*aSmall) + ((irow-row0)*(irow-row0) / (bSmall*bSmall)) > 1))
1556  {
1557  return kTRUE;
1558  }
1559  else
1560  {
1561  return kFALSE;
1562  }
1563 }
1564 
1565 
1575 //______________________________________________________________________________________________________
1576 Bool_t AliAnalysisTaskEMCALPi0CalibSelection::IsInZone7(Int_t iSupMod, Int_t ieta, Int_t iphi)
1577 {
1578  Int_t icol = ieta;
1579  Int_t irow = iphi;
1580 
1581  //Center of ellipse
1582  Float_t col0 = 47/2;
1583  Float_t row0 = 23/2;
1584 
1585  //Paramters
1586  Float_t a = 16-col0;
1587  Float_t b = 7-row0;
1588 
1589  if(((icol-col0)*(icol-col0)) / (a*a) + ((irow-row0)*(irow-row0) / (b*b)) < 1)
1590  {
1591  return kTRUE;
1592  }
1593  else
1594  {
1595  return kFALSE;
1596  }
1597 
1598 }
1599 
1600 
1606 //__________________________________________________________________________
1608 {
1609  // Event selection
1610 
1611  if(fTriggerName!="")
1612  {
1613  AliESDEvent* esdevent = dynamic_cast<AliESDEvent*> (InputEvent());
1614  AliAODEvent* aodevent = dynamic_cast<AliAODEvent*> (InputEvent());
1615 
1616  TString triggerClass = "";
1617  if (esdevent) triggerClass = esdevent->GetFiredTriggerClasses();
1618  else if(aodevent) triggerClass = aodevent->GetFiredTriggerClasses();
1619 
1620  AliDebug(1,Form("Event %d, FiredClass %s",
1621  (Int_t)Entry(),(((AliESDEvent*)InputEvent())->GetFiredTriggerClasses()).Data()));
1622 
1623  if(!triggerClass.Contains(fTriggerName))
1624  {
1625  AliDebug(1,"Reject event!");
1626  return;
1627  }
1628  else
1629  AliDebug(1,"Accept event!");
1630  }
1631 
1632  // Get the input event
1633 
1634  AliVEvent* event = 0;
1635  if(fFilteredInput) event = AODEvent();
1636  else event = InputEvent();
1637 
1638  if(!event)
1639  {
1640  AliWarning("Input event not available!");
1641  return;
1642  }
1643 
1644  AliDebug(1,Form("<<< %s: Event %d >>>",event->GetName(), (Int_t)Entry()));
1645 
1646  // Get the primary vertex
1647 
1648  event->GetPrimaryVertex()->GetXYZ(fVertex) ;
1649 
1650  AliDebug(1,Form("Vertex: (%.3f,%.3f,%.3f)",fVertex[0],fVertex[1],fVertex[2]));
1651 
1652  //Int_t runNum = aod->GetRunNumber();
1653  //if(DebugLevel() > 1) printf("Run number: %d\n",runNum);
1654 
1655  fhNEvents->Fill(0); //Count the events to be analyzed
1656 
1657  // Acccess once the geometry matrix and temperature corrections and calibration coefficients
1658  if(fhNEvents->GetEntries() == 1)
1659  {
1661 
1663 
1664 // InitEnergyCalibrationFactors();
1665  }
1666 
1667  //Get the list of clusters and cells
1668  fEMCALCells = event->GetEMCALCells();
1669 
1670  fCaloClustersArr = new TRefArray();
1671  event->GetEMCALClusters(fCaloClustersArr);
1672 
1673  AliDebug(1,Form("N CaloClusters: %d - N CaloCells %d",fCaloClustersArr->GetEntriesFast(), fEMCALCells->GetNumberOfCells()));
1674 
1675  // Apply non linearity, new calibration, T calibration to the clusters
1676  if( fCorrectClusters )
1677  CorrectClusters();
1678 
1679  FillHistograms();
1680 
1681  delete fCaloClustersArr;
1682 
1683  PostData(1,fOutputContainer);
1684 }
1685 
1688 //_____________________________________________________
1690 {
1691  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",
1693 
1694  printf("Group %d cells\n", fGroupNCells) ;
1695 
1696  printf("Cluster maximal cell away from border at least %d cells\n", fRecoUtils->GetNumberOfCellsFromEMCALBorder()) ;
1697 
1698  printf("Histograms: bins %d; energy range: %2.2f < E < %2.2f MeV\n",fNbins,fMinBin,fMaxBin) ;
1699 
1700  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",
1701  fRecoUtils->IsBadChannelsRemovalSwitchedOn(),fFilteredInput,fCorrectClusters, fRecalPosition, fSameSM) ;
1702 
1703  printf("OADB path : %s\n",fOADBFilePath .Data());
1704  printf("Calibration path : %s\n",fCalibFilePath.Data());
1705 
1706  printf("EMCAL Geometry name: < %s >, Load Matrices %d\n",fEMCALGeoName.Data(), fLoadMatrices) ;
1707 
1708  if(fLoadMatrices) { for(Int_t ism = 0; ism < AliEMCALGeoParams::fgkEMCALModules; ism++) if(fMatrix[ism]) fMatrix[ism]->Print() ; }
1709 }
1710 
1714 //_____________________________________________________________________
1716 {
1717  if(n > fNMaskCellColumns)
1718  {
1719  delete [] fMaskCellColumns ;
1720 
1721  fMaskCellColumns = new Int_t[n] ;
1722  }
1723 
1724  fNMaskCellColumns = n ;
1725 }
1726 
1731 //___________________________________________________________________________________
1733 {
1734  if(ipos < fNMaskCellColumns) fMaskCellColumns[ipos] = icol ;
1735  else AliWarning("Mask column not set, position larger than allocated set size first") ;
1736 }
1737 
1743 //___________________________________________________________________________________
1745 {
1746  Int_t iPos;
1747 
1748  if(icol%2 == 0)
1749  {
1750  if(irow%8 < 4)
1751  {
1752  iPos = 0;
1753  }
1754  else if(irow%8 < 8)
1755  {
1756  iPos = 2;
1757  }
1758  else iPos = -1;
1759 
1760  }
1761  else
1762  {
1763  if(irow%8 < 4)
1764  {
1765  iPos = 1;
1766  }
1767  else if(irow%8 < 8)
1768  {
1769  iPos = 3;
1770  }
1771  else iPos = -1;
1772  }
1773 
1774  return iPos;
1775 }
1776 
1779 //______________________________________________________________
1781 {
1782  AliDebug(1,"Not implemented");
1783 // const Int_t buffersize = 255;
1784 // char onePar[buffersize] ;
1785 
1786 // 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",
1787 // fEmin,fEmax, fL0min, fL0max, fMinNCells, fAsyCut, fDTimeCut, fTimeMin, fTimeMax, fInvMassCutMin, fInvMassCutMax) ;
1788 // fCuts->Add(new TObjString(onePar));
1789 // snprintf(onePar,buffersize, "Group %d cells;", fGroupNCells) ;
1790 // fCuts->Add(new TObjString(onePar));
1791 // snprintf(onePar,buffersize, "Cluster maximal cell away from border at least %d cells;", fRecoUtils->GetNumberOfCellsFromEMCALBorder()) ;
1792 // fCuts->Add(new TObjString(onePar));
1793 // snprintf(onePar,buffersize, "Histograms, Mass bins %d; energy range: %2.2f < E < %2.2f GeV;",fNbins,fMinBin,fMaxBin) ;
1794 // fCuts->Add(new TObjString(onePar));
1795 // snprintf(onePar,buffersize, "Histograms, Time bins %d; energy range: %2.2f < E < %2.2f GeV;",fNTimeBins,fMinTimeBin,fMaxTimeBin) ;
1796 // fCuts->Add(new TObjString(onePar));
1797 // 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 ",
1798 // fRecoUtils->IsBadChannelsRemovalSwitchedOn(),fFilteredInput,fCorrectClusters, fRecalPosition, fSameSM) ;
1799 // fCuts->Add(new TObjString(onePar));
1800 // snprintf(onePar,buffersize, "EMCAL Geometry name: < %s >, Load Matrices? %d",fEMCALGeoName.Data(),fLoadMatrices) ;
1801 // fCuts->Add(new TObjString(onePar));
1802 //
1803 // // Post Data
1804 // PostData(2, fCuts);
1805 }
1806 
TH2F * fHmggSM_Zone4[AliEMCALGeoParams::fgkEMCALModules]
! Two-cluster invariant mass per SM in zone 4.
Float_t fDTimeCut
Maximum difference between time of cluster pairs (ns).
TH2F * fHmggPairSameSectorSMMaskFrame[AliEMCALGeoParams::fgkEMCALModules/2]
! Two-cluster invariant mass per Pair, mask clusters facing frames.
TGeoHMatrix * fMatrix[AliEMCALGeoParams::fgkEMCALModules]
Bool_t IsInZone3(Int_t iSupMod, Int_t ieta, Int_t iphi)
TH2F * fhClusterPairDiffTimeSameSide[AliEMCALGeoParams::fgkEMCALModules-2]
! Diference in time of clusters same side.
Bool_t IsInZone6(Int_t iSupMod, Int_t ieta, Int_t iphi)
TH2F * fHmggSM_Zone3[AliEMCALGeoParams::fgkEMCALModules]
! Two-cluster invariant mass per SM in zone 3.
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.
Bool_t IsInZone4(Int_t iSupMod, Int_t ieta, Int_t iphi)
TString fImportGeometryFilePath
Path fo geometry.root file.
TH2F * fhTopoClusterAmpCase0[AliEMCALGeoParams::fgkEMCALModules]
! Cell amplitude map for type 0 cluster in noisy quartet
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.
Bool_t fSelectOnlyPhotonsInDifferentSM
Select only pairs of photons that are not in the same SM.
Float_t fLogWeight
Logarithmic weight used in cluster recalibration.
Bool_t IsInZone7(Int_t iSupMod, Int_t ieta, Int_t iphi)
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.
Bool_t IsInZone5(Int_t iSupMod, Int_t ieta, Int_t iphi)
Int_t FindPositionInNoisyQuartet(Int_t irow, Int_t icol, Int_t iSM)
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 * fHmggSM_Zone5[AliEMCALGeoParams::fgkEMCALModules]
! Two-cluster invariant mass per SM in zone 5.
TH2F * fHOpeningAngleSM[AliEMCALGeoParams::fgkEMCALModules]
! Two-cluster opening angle vs pt per SM,with mass close to pi0.
TLorentzVector fMomentum2
Cluster kinematics, temporal.
TH2F * fhTopoClusterAmpCase2[AliEMCALGeoParams::fgkEMCALModules]
! Cell amplitude map for type 2 cluster in noisy quartet
Float_t fInvMassCutMin
Minimum mass cut for clusters to fill time or other histograms.
TH2F * fhTopoClusterAmpCase3[AliEMCALGeoParams::fgkEMCALModules]
! Cell amplitude map for type 3 cluster in noisy quartet
Float_t fMinEnergyBin
Minimum energy bins of cell energy histograms.
Bool_t IsInZone2(Int_t iSupMod, Int_t ieta, Int_t iphi)
Int_t fNEnergybins
N energy bins of cell energy 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 * fhTopoClusterAmpCase1[AliEMCALGeoParams::fgkEMCALModules]
! Cell amplitude map for type 1 cluster in noisy quartet
Bool_t fSelectOnlyCellSignalOutOfCollision
Select cells / clusters that are due to noise, i.e. signal in EMCal that happens not during collision...
TH2F * fHmggSM_Zone6[AliEMCALGeoParams::fgkEMCALModules]
! Two-cluster invariant mass per SM in zone 6.
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.
TH1F * fhEnergy[AliEMCALGeoParams::fgkEMCALModules][AliEMCALGeoParams::fgkEMCALCols][AliEMCALGeoParams::fgkEMCALRows]
! Energy distribution for each cell.
Float_t fMinBin
Minimum mass bins of invariant mass histograms.
TH2F * fhTopoClusterAmpFractionCase2[AliEMCALGeoParams::fgkEMCALModules]
! Cell amplitude fraction map for type 2 cluster in noisy quartet
Bool_t fRecalPosition
Switch on/off cluster position calculation, in case alignment matrices are not available.
Float_t fMaxEnergyBin
Maximum energy bins of cell energy histograms.
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 * fhTopoClusterAmpFractionCase3[AliEMCALGeoParams::fgkEMCALModules]
! Cell amplitude fraction map for type 3 cluster in noisy quartet
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 * fHmggSM_Zone2[AliEMCALGeoParams::fgkEMCALModules]
! Two-cluster invariant mass per SM in zone 2.
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...
ClassImp(AliAnalysisTaskCRC) AliAnalysisTaskCRC
TH2F * fHmggPairSameSideSM[AliEMCALGeoParams::fgkEMCALModules-2]
! Two-cluster invariant mass per Pair.
TH2F * fhTopoClusterAmpFractionCase0[AliEMCALGeoParams::fgkEMCALModules]
! Cell amplitude fraction map for type 0 cluster in noisy quartet
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.
Bool_t IsInZone1(Int_t iSupMod, Int_t ieta, Int_t iphi)
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 * fHmggSM_Zone1[AliEMCALGeoParams::fgkEMCALModules]
! Two-cluster invariant mass per SM in zone 1.
TH2F * fhTopoClusterAmpFractionCase1[AliEMCALGeoParams::fgkEMCALModules]
! Cell amplitude fraction map for type 1 cluster in noisy quartet
TH2F * fHmggSM_Zone7[AliEMCALGeoParams::fgkEMCALModules]
! Two-cluster invariant mass per SM in zone 7.
TH2F * fHmggPairSameSectorSM[AliEMCALGeoParams::fgkEMCALModules/2]
! Two-cluster invariant mass per Pair.