add new time cut, time histograms
[u/mrichter/AliRoot.git] / PWG4 / CaloCalib / AliAnalysisTaskEMCALPi0CalibSelection.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: Boris Polishchuk                                               *
5  * Adapted to AOD reading by Gustavo Conesa                               *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 //---------------------------------------------------------------------------// 
17 //                                                                           //
18 // Fill histograms (one per cell) with two-cluster invariant mass            //
19 // using calibration coefficients of the previous iteration.                 //
20 // Histogram for a given cell is filled if the most energy of one cluster    //
21 // is deposited in this cell and the other cluster could be anywherein EMCAL.//
22 //                                                                           //
23 //---------------------------------------------------------------------------//
24
25 //#include <cstdlib>
26 //#include <Riostream.h>
27 // Root 
28 #include "TLorentzVector.h"
29 #include "TRefArray.h"
30 #include "TList.h"
31 #include "TH1F.h"
32 #include <TGeoManager.h>
33
34 // AliRoot
35 #include "AliAnalysisTaskEMCALPi0CalibSelection.h"
36 #include "AliAODEvent.h"
37 #include "AliESDEvent.h"
38 #include "AliEMCALGeometry.h"
39 #include "AliVCluster.h"
40 #include "AliVCaloCells.h"
41 #include "AliEMCALRecoUtils.h"
42 //#include "AliEMCALAodCluster.h"
43 //#include "AliEMCALCalibData.h"
44
45 ClassImp(AliAnalysisTaskEMCALPi0CalibSelection)
46
47
48 //__________________________________________________
49 AliAnalysisTaskEMCALPi0CalibSelection::AliAnalysisTaskEMCALPi0CalibSelection(const char* name) :
50   AliAnalysisTaskSE(name),fEMCALGeo(0x0),//fCalibData(0x0), 
51   fEmin(0.5), fEmax(15.), fDTimeCut(20.), fAsyCut(1.),fMinNCells(2), fGroupNCells(0),
52   fLogWeight(4.5), fSameSM(kFALSE), fFilteredInput(kFALSE),
53   fCorrectClusters(kFALSE), fEMCALGeoName("EMCAL_COMPLETEV1"), 
54   fRecoUtils(new AliEMCALRecoUtils),
55   fNbins(300), fMinBin(0.), fMaxBin(300.),fOutputContainer(0x0),
56   fHmgg(0x0),           fHmggDifferentSM(0x0), 
57   fHmggMaskFrame(0x0),  fHmggDifferentSMMaskFrame(0x0), 
58   fHOpeningAngle(0x0),  fHOpeningAngleDifferentSM(0x0),  
59   fHIncidentAngle(0x0), fHIncidentAngleDifferentSM(0x0),
60   fHAsymmetry(0x0),     fHAsymmetryDifferentSM(0x0),  
61   fhNEvents(0x0),fCuts(0x0),fLoadMatrices(0),
62   fNMaskCellColumns(11), fMaskCellColumns(0x0),
63   fhClusterTime(0x0), fhClusterPairDiffTime(0x0)
64 {
65   //Named constructor which should be used.
66   
67   for(Int_t iMod=0; iMod < AliEMCALGeoParams::fgkEMCALModules; iMod++) {
68     for(Int_t iX=0; iX<24; iX++) {
69       for(Int_t iZ=0; iZ<48; iZ++) {
70         fHmpi0[iMod][iZ][iX]=0;
71       }
72     } 
73   }
74   
75   fMaskCellColumns = new Int_t[fNMaskCellColumns];
76   fMaskCellColumns[0] = 6 ;  fMaskCellColumns[1] = 7 ;  fMaskCellColumns[2] = 8 ; 
77   fMaskCellColumns[3] = 35;  fMaskCellColumns[4] = 36;  fMaskCellColumns[5] = 37; 
78   fMaskCellColumns[6] = 12+AliEMCALGeoParams::fgkEMCALCols; fMaskCellColumns[7] = 13+AliEMCALGeoParams::fgkEMCALCols;
79   fMaskCellColumns[8] = 40+AliEMCALGeoParams::fgkEMCALCols; fMaskCellColumns[9] = 41+AliEMCALGeoParams::fgkEMCALCols; 
80   fMaskCellColumns[10]= 42+AliEMCALGeoParams::fgkEMCALCols; 
81   
82   for(Int_t iSMPair = 0; iSMPair < AliEMCALGeoParams::fgkEMCALModules/2; iSMPair++) {
83     fHmggPairSameSectorSM[iSMPair]          = 0;
84     fHmggPairSameSectorSMMaskFrame[iSMPair] = 0;
85     fhClusterPairDiffTimeSameSector[iSMPair]= 0;
86   }
87   for(Int_t iSMPair = 0; iSMPair < AliEMCALGeoParams::fgkEMCALModules-2; iSMPair++){ 
88     fHmggPairSameSideSM[iSMPair]            = 0;
89     fHmggPairSameSideSMMaskFrame[iSMPair]   = 0;
90     fhClusterPairDiffTimeSameSide[iSMPair]  = 0;
91   }
92   
93   for(Int_t iSM = 0; iSM < AliEMCALGeoParams::fgkEMCALModules; iSM++) {
94     fHmggSM[iSM]                     = 0;
95     fHmggSMMaskFrame[iSM]            = 0;
96     fHOpeningAngleSM[iSM]            = 0;
97     fHOpeningAnglePairSM[iSM]        = 0;
98     fHAsymmetrySM[iSM]               = 0;
99     fHAsymmetryPairSM[iSM]           = 0;
100     fHIncidentAngleSM[iSM]           = 0;
101     fHIncidentAnglePairSM[iSM]       = 0;
102     fhTowerDecayPhotonHit[iSM]       = 0;
103     fhTowerDecayPhotonEnergy[iSM]    = 0;
104     fhTowerDecayPhotonAsymmetry[iSM] = 0;
105     fhTowerDecayPhotonHitMaskFrame[iSM]= 0;
106     fMatrix[iSM]                     = 0x0;
107     fhClusterTimeSM[iSM]             = 0;
108     fhClusterPairDiffTimeSameSM[iSM] = 0;
109   }
110   
111   DefineOutput(1, TList::Class());
112   DefineOutput(2, TList::Class());  // will contain cuts or local params
113   
114 }
115
116 //__________________________________________________
117 AliAnalysisTaskEMCALPi0CalibSelection::~AliAnalysisTaskEMCALPi0CalibSelection()
118 {
119   //Destructor.
120   
121   if(fOutputContainer){
122     fOutputContainer->Delete() ; 
123     delete fOutputContainer ;
124   }
125   
126   if(fEMCALGeo)         delete fEMCALGeo  ;
127   if(fRecoUtils)        delete fRecoUtils ;
128   if(fNMaskCellColumns) delete [] fMaskCellColumns;
129   
130 }
131
132 //_____________________________________________________
133 void AliAnalysisTaskEMCALPi0CalibSelection::LocalInit()
134 {
135   // Local Initialization
136   
137   // Create cuts/param objects and publish to slot
138   const Int_t buffersize = 255;
139   char onePar[buffersize] ;
140   fCuts = new TList();
141   
142   snprintf(onePar,buffersize, "Custer cuts: %2.2f < E < %2.2f GeV; min number of cells %d; Assymetry cut %1.2f", fEmin,fEmax, fMinNCells, fAsyCut) ;
143   fCuts->Add(new TObjString(onePar));
144   snprintf(onePar,buffersize, "Group %d cells;", fGroupNCells) ;
145   fCuts->Add(new TObjString(onePar));
146   snprintf(onePar,buffersize, "Cluster maximal cell away from border at least %d cells;", fRecoUtils->GetNumberOfCellsFromEMCALBorder()) ;
147   fCuts->Add(new TObjString(onePar));
148   snprintf(onePar,buffersize, "Histograms: bins %d; energy range: %2.2f < E < %2.2f GeV;",fNbins,fMinBin,fMaxBin) ;
149   fCuts->Add(new TObjString(onePar));
150   snprintf(onePar,buffersize, "Switchs: Remove Bad Channels? %d; Use filtered input? %d;  Correct Clusters? %d, Mass per channel same SM clusters? %d ",
151            fRecoUtils->IsBadChannelsRemovalSwitchedOn(),fFilteredInput,fCorrectClusters, fSameSM) ;
152   fCuts->Add(new TObjString(onePar));
153   snprintf(onePar,buffersize, "EMCAL Geometry name: < %s >, Load Matrices? %d",fEMCALGeoName.Data(),fLoadMatrices) ;
154   fCuts->Add(new TObjString(onePar));
155   
156   fCuts ->SetOwner(kTRUE);
157   
158   // Post Data
159   PostData(2, fCuts);
160   
161 }
162
163 //__________________________________________________
164 void AliAnalysisTaskEMCALPi0CalibSelection::UserCreateOutputObjects()
165 {
166   //Create output container, init geometry 
167   
168   fEMCALGeo =  AliEMCALGeometry::GetInstance(fEMCALGeoName) ;   
169   Int_t nSM = (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules();
170   
171   fOutputContainer = new TList();
172   const Int_t buffersize = 255;
173   char hname[buffersize], htitl[buffersize];
174   
175   for(Int_t iMod=0; iMod < nSM; iMod++) {
176     for(Int_t iRow=0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
177       for(Int_t iCol=0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++) {
178         snprintf(hname,buffersize, "%d_%d_%d",iMod,iCol,iRow);
179         snprintf(htitl,buffersize, "Two-gamma inv. mass for super mod %d, cell(col,row)=(%d,%d)",iMod,iCol,iRow);
180         fHmpi0[iMod][iCol][iRow] = new TH1F(hname,htitl,fNbins,fMinBin,fMaxBin);
181         fOutputContainer->Add(fHmpi0[iMod][iCol][iRow]);
182       }
183     }
184   }
185   
186   fHmgg = new TH2F("hmgg","2-cluster invariant mass",fNbins,fMinBin,fMaxBin,100,0,10);
187   fHmgg->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
188   fHmgg->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
189   fOutputContainer->Add(fHmgg);
190
191   fHmggDifferentSM = new TH2F("hmggDifferentSM","2-cluster invariant mass, different SM",fNbins,fMinBin,fMaxBin,100,0,10);
192   fHmggDifferentSM->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
193   fHmggDifferentSM->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
194   fOutputContainer->Add(fHmggDifferentSM);
195
196   fHOpeningAngle = new TH2F("hopang","2-cluster opening angle",100,0.,50.,100,0,10);
197   fHOpeningAngle->SetXTitle("#alpha_{#gamma #gamma}");
198   fHOpeningAngle->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
199   fOutputContainer->Add(fHOpeningAngle);
200   
201   fHOpeningAngleDifferentSM = new TH2F("hopangDifferentSM","2-cluster opening angle, different SM",100,0,50.,100,0,10);
202   fHOpeningAngleDifferentSM->SetXTitle("#alpha_{#gamma #gamma}");
203   fHOpeningAngleDifferentSM->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
204   fOutputContainer->Add(fHOpeningAngleDifferentSM);
205   
206   fHIncidentAngle = new TH2F("hinang","#gamma incident angle in SM",100,0.,20.,100,0,10);
207   fHIncidentAngle->SetXTitle("#alpha_{#gamma SM center}");
208   fHIncidentAngle->SetYTitle("p_{T #gamma} (GeV/c)");
209   fOutputContainer->Add(fHIncidentAngle);
210   
211   fHIncidentAngleDifferentSM = new TH2F("hinangDifferentSM","#gamma incident angle in SM, different SM pair",100,0,20.,100,0,10);
212   fHIncidentAngleDifferentSM->SetXTitle("#alpha_{#gamma - SM center}");
213   fHIncidentAngleDifferentSM->SetYTitle("p_{T #gamma} (GeV/c)");
214   fOutputContainer->Add(fHIncidentAngleDifferentSM);
215   
216   fHAsymmetry = new TH2F("hasym","2-cluster opening angle",100,0.,1.,100,0,10);
217   fHAsymmetry->SetXTitle("a");
218   fHAsymmetry->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
219   fOutputContainer->Add(fHAsymmetry);
220   
221   fHAsymmetryDifferentSM = new TH2F("hasymDifferentSM","2-cluster opening angle, different SM",100,0,1.,100,0,10);
222   fHAsymmetryDifferentSM->SetXTitle("a");
223   fHAsymmetryDifferentSM->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
224   fOutputContainer->Add(fHAsymmetryDifferentSM);
225   
226   
227   //TString pairname[] = {"A side (0-2)", "C side (1-3)","Row 0 (0-1)", "Row 1 (2-3)"};
228   
229   fHmggMaskFrame = new TH2F("hmggMaskFrame","2-cluster invariant mass, frame masked",fNbins,fMinBin,fMaxBin,100,0,10);
230   fHmggMaskFrame->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
231   fHmggMaskFrame->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
232   fOutputContainer->Add(fHmggMaskFrame);
233   
234   fHmggDifferentSMMaskFrame = new TH2F("hmggDifferentSMMaskFrame","2-cluster invariant mass, different SM, frame masked",
235                                        fNbins,fMinBin,fMaxBin,100,0,10);
236   fHmggDifferentSMMaskFrame->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
237   fHmggDifferentSMMaskFrame->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
238   fOutputContainer->Add(fHmggDifferentSMMaskFrame);
239   
240   
241   for(Int_t iSM = 0; iSM < nSM; iSM++) {
242     
243     snprintf(hname, buffersize, "hmgg_SM%d",iSM);
244     snprintf(htitl, buffersize, "Two-gamma inv. mass for super mod %d",iSM);
245     fHmggSM[iSM] = new TH2F(hname,htitl,fNbins,fMinBin,fMaxBin,100,0,10);
246     fHmggSM[iSM]->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
247     fHmggSM[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
248     fOutputContainer->Add(fHmggSM[iSM]);
249     
250     snprintf(hname, buffersize, "hmgg_SM%d_MaskFrame",iSM);
251     snprintf(htitl, buffersize, "Two-gamma inv. mass for super mod %d",iSM);
252     fHmggSMMaskFrame[iSM] = new TH2F(hname,htitl,fNbins,fMinBin,fMaxBin,100,0,10);
253     fHmggSMMaskFrame[iSM]->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
254     fHmggSMMaskFrame[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
255     fOutputContainer->Add(fHmggSMMaskFrame[iSM]);
256     
257     
258     if(iSM < nSM/2){
259       snprintf(hname,buffersize, "hmgg_PairSameSectorSM%d",iSM);
260       snprintf(htitl,buffersize, "Two-gamma inv. mass for SM pair Sector: %d",iSM);
261       fHmggPairSameSectorSM[iSM] = new TH2F(hname,htitl,fNbins,fMinBin,fMaxBin,100,0,10);
262       fHmggPairSameSectorSM[iSM]->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
263       fHmggPairSameSectorSM[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
264       fOutputContainer->Add(fHmggPairSameSectorSM[iSM]);
265       
266       snprintf(hname,buffersize, "hmgg_PairSameSectorSM%d_MaskFrame",iSM);
267       snprintf(htitl,buffersize, "Two-gamma inv. mass for SM pair Sector: %d",iSM);
268       fHmggPairSameSectorSMMaskFrame[iSM] = new TH2F(hname,htitl,fNbins,fMinBin,fMaxBin,100,0,10);
269       fHmggPairSameSectorSMMaskFrame[iSM]->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
270       fHmggPairSameSectorSMMaskFrame[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
271       fOutputContainer->Add(fHmggPairSameSectorSMMaskFrame[iSM]);
272
273       fhClusterPairDiffTimeSameSector[iSM] = new TH2F(Form("hClusterPairDiffTimeSameSector%d",iSM),
274                                                       Form("cluster pair time difference vs E, Sector %d",iSM),
275                                                       100,0,10, 200,-100,100);
276       fhClusterPairDiffTimeSameSector[iSM]->SetXTitle("E_{pair} (GeV)");
277       fhClusterPairDiffTimeSameSector[iSM]->SetYTitle("#Delta t (ns)");
278       fOutputContainer->Add(fhClusterPairDiffTimeSameSector[iSM]);
279
280
281     }
282     
283     if(iSM < nSM-2){
284       snprintf(hname,buffersize, "hmgg_PairSameSideSM%d",iSM);
285       snprintf(htitl,buffersize, "Two-gamma inv. mass for SM pair Sector: %d",iSM);
286       fHmggPairSameSideSM[iSM] = new TH2F(hname,htitl,fNbins,fMinBin,fMaxBin,100,0,10);
287       fHmggPairSameSideSM[iSM]->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
288       fHmggPairSameSideSM[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
289       fOutputContainer->Add(fHmggPairSameSideSM[iSM]);
290       
291       snprintf(hname,buffersize, "hmgg_PairSameSideSM%d_MaskFrame",iSM);
292       snprintf(htitl,buffersize, "Two-gamma inv. mass for SM pair Sector: %d",iSM);
293       fHmggPairSameSideSMMaskFrame[iSM] = new TH2F(hname,htitl,fNbins,fMinBin,fMaxBin,100,0,10);
294       fHmggPairSameSideSMMaskFrame[iSM]->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
295       fHmggPairSameSideSMMaskFrame[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
296       fOutputContainer->Add(fHmggPairSameSideSMMaskFrame[iSM]);   
297
298       fhClusterPairDiffTimeSameSide[iSM] = new TH2F(Form("hClusterPairDiffTimeSameSide%d",iSM),
299                                                     Form("cluster pair time difference vs E,  Side %d",iSM),
300                                                     100,0,10, 200,-100,100);
301       fhClusterPairDiffTimeSameSide[iSM]->SetXTitle("E_{pair} (GeV)");
302       fhClusterPairDiffTimeSameSide[iSM]->SetYTitle("#Delta t (ns)");
303       fOutputContainer->Add(fhClusterPairDiffTimeSameSide[iSM]);
304
305     }
306     
307     snprintf(hname, buffersize, "hopang_SM%d",iSM);
308     snprintf(htitl, buffersize, "Opening angle for super mod %d",iSM);
309     fHOpeningAngleSM[iSM] = new TH2F(hname,htitl,100,0.,50.,100,0,10);
310     fHOpeningAngleSM[iSM]->SetXTitle("#alpha_{#gamma #gamma} (deg)");
311     fHOpeningAngleSM[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
312     fOutputContainer->Add(fHOpeningAngleSM[iSM]);
313     
314     snprintf(hname,buffersize, "hopang_PairSM%d",iSM);
315     snprintf(htitl,buffersize, "Opening angle for SM pair: %d",iSM);
316     fHOpeningAnglePairSM[iSM] = new TH2F(hname,htitl,100,0.,50.,100,0,10);
317     fHOpeningAnglePairSM[iSM]->SetXTitle("#alpha_{#gamma #gamma} (deg)");
318     fHOpeningAnglePairSM[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
319     fOutputContainer->Add(fHOpeningAnglePairSM[iSM]);    
320     
321     snprintf(hname, buffersize, "hinang_SM%d",iSM);
322     snprintf(htitl, buffersize, "Incident angle for super mod %d",iSM);
323     fHIncidentAngleSM[iSM] = new TH2F(hname,htitl,100,0.,20.,100,0,10);
324     fHIncidentAngleSM[iSM]->SetXTitle("#alpha_{#gamma - SM center} (deg)");
325     fHIncidentAngleSM[iSM]->SetYTitle("p_{T #gamma} (GeV/c)");
326     fOutputContainer->Add(fHIncidentAngleSM[iSM]);
327     
328     snprintf(hname,buffersize, "hinang_PairSM%d",iSM);
329     snprintf(htitl,buffersize, "Incident angle for SM pair: %d",iSM);
330     fHIncidentAnglePairSM[iSM] = new TH2F(hname,htitl,100,0.,20.,100,0,10);
331     fHIncidentAnglePairSM[iSM]->SetXTitle("#alpha_{#gamma - SM center} (deg)");
332     fHIncidentAnglePairSM[iSM]->SetYTitle("p_{T #gamma} (GeV/c)");
333     fOutputContainer->Add(fHIncidentAnglePairSM[iSM]);   
334     
335     snprintf(hname, buffersize, "hasym_SM%d",iSM);
336     snprintf(htitl, buffersize, "Asymmetry for super mod %d",iSM);
337     fHAsymmetrySM[iSM] = new TH2F(hname,htitl,100,0.,1.,100,0,10);
338     fHAsymmetrySM[iSM]->SetXTitle("a");
339     fHAsymmetrySM[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
340     fOutputContainer->Add(fHAsymmetrySM[iSM]);
341     
342     snprintf(hname,buffersize, "hasym_PairSM%d",iSM);
343     snprintf(htitl,buffersize, "Asymmetry for SM pair: %d",iSM);
344     fHAsymmetryPairSM[iSM] = new TH2F(hname,htitl,100,0.,1.,100,0,10);
345     fHAsymmetryPairSM[iSM]->SetXTitle("a");
346     fHAsymmetryPairSM[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
347     fOutputContainer->Add(fHAsymmetryPairSM[iSM]);    
348     
349     Int_t colmax = 48;
350     Int_t rowmax = 24;
351     
352     fhTowerDecayPhotonHit[iSM]  = new TH2F (Form("hTowerDecPhotonHit_Mod%d",iSM),
353                                             Form("Entries in grid of cells in Module %d",iSM), 
354                                             colmax+2,-1.5,colmax+0.5, rowmax+2,-1.5,rowmax+0.5); 
355     fhTowerDecayPhotonHit[iSM]->SetYTitle("row (phi direction)");
356     fhTowerDecayPhotonHit[iSM]->SetXTitle("column (eta direction)");
357     fOutputContainer->Add(fhTowerDecayPhotonHit[iSM]);
358     
359     fhTowerDecayPhotonEnergy[iSM]  = new TH2F (Form("hTowerDecPhotonEnergy_Mod%d",iSM),
360                                                Form("Accumulated energy in grid of cells in Module %d",iSM), 
361                                                colmax+2,-1.5,colmax+0.5, rowmax+2,-1.5,rowmax+0.5); 
362     fhTowerDecayPhotonEnergy[iSM]->SetYTitle("row (phi direction)");
363     fhTowerDecayPhotonEnergy[iSM]->SetXTitle("column (eta direction)");
364     fOutputContainer->Add(fhTowerDecayPhotonEnergy[iSM]);
365     
366     fhTowerDecayPhotonAsymmetry[iSM]  = new TH2F (Form("hTowerDecPhotonAsymmetry_Mod%d",iSM),
367                                                   Form("Accumulated asymmetry in grid of cells in Module %d",iSM), 
368                                                   colmax+2,-1.5,colmax+0.5, rowmax+2,-1.5,rowmax+0.5); 
369     fhTowerDecayPhotonAsymmetry[iSM]->SetYTitle("row (phi direction)");
370     fhTowerDecayPhotonAsymmetry[iSM]->SetXTitle("column (eta direction)");
371     fOutputContainer->Add(fhTowerDecayPhotonAsymmetry[iSM]);
372     
373     fhTowerDecayPhotonHitMaskFrame[iSM]  = new TH2F (Form("hTowerDecPhotonHit_Mod%d_MaskFrame",iSM),Form("Entries in grid of cells in Module %d",iSM), 
374                                             colmax+2,-1.5,colmax+0.5, rowmax+2,-1.5,rowmax+0.5); 
375     fhTowerDecayPhotonHitMaskFrame[iSM]->SetYTitle("row (phi direction)");
376     fhTowerDecayPhotonHitMaskFrame[iSM]->SetXTitle("column (eta direction)");
377     fOutputContainer->Add(fhTowerDecayPhotonHitMaskFrame[iSM]);
378
379     fhClusterTimeSM[iSM] = new TH2F(Form("hClusterTime_SM%d",iSM),"cluster time vs E",100,0,10, 100,0,1000);
380     fhClusterTimeSM[iSM]->SetXTitle("E (GeV)");
381     fhClusterTimeSM[iSM]->SetYTitle("t (ns)");
382     fOutputContainer->Add(fhClusterTimeSM[iSM]);
383     
384     fhClusterPairDiffTimeSameSM[iSM] = new TH2F(Form("hClusterPairDiffTimeSameSM%d",iSM),
385                                                 Form("cluster pair time difference vs E, SM %d",iSM),
386                                                 100,0,10, 200,-100,100);
387     fhClusterPairDiffTimeSameSM[iSM]->SetXTitle("E (GeV)");
388     fhClusterPairDiffTimeSameSM[iSM]->SetYTitle("#Delta t (ns)");
389     fOutputContainer->Add(fhClusterPairDiffTimeSameSM[iSM]);
390
391   }  
392   
393   fhClusterTime = new TH2F("hClusterTime","cluster time vs E",100,0,10, 100,0,1000);
394   fhClusterTime->SetXTitle("E (GeV)");
395   fhClusterTime->SetYTitle("t (ns)");
396   fOutputContainer->Add(fhClusterTime);
397
398   fhClusterPairDiffTime = new TH2F("hClusterPairDiffTime","cluster pair time difference vs E",100,0,10, 200,-100,100);
399   fhClusterPairDiffTime->SetXTitle("E_{pair} (GeV)");
400   fhClusterPairDiffTime->SetYTitle("#Delta t (ns)");
401   fOutputContainer->Add(fhClusterPairDiffTime);
402
403   
404   fhNEvents        = new TH1I("hNEvents", "Number of analyzed events"   , 1 , 0 , 1  ) ;
405   fOutputContainer->Add(fhNEvents);
406   
407   fOutputContainer->SetOwner(kTRUE);
408   
409   //  fCalibData = new AliEMCALCalibData();
410   
411   PostData(1,fOutputContainer);
412
413 }
414
415 //__________________________________________________
416 Bool_t AliAnalysisTaskEMCALPi0CalibSelection::MaskFrameCluster(const Int_t iSM,  const Int_t ieta) const {
417   //Check if cell is in one of the regions where we have significant amount of material in front of EMCAL
418   
419   Int_t icol = ieta;
420   if(iSM%2) icol+=48; // Impair SM, shift index [0-47] to [48-96]
421   
422   if (fNMaskCellColumns && fMaskCellColumns) {
423     for (Int_t imask = 0; imask < fNMaskCellColumns; imask++) {
424       if(icol==fMaskCellColumns[imask]) return kTRUE;
425     }
426   }
427   
428   return kFALSE;
429   
430 }
431
432 //__________________________________________________
433 void AliAnalysisTaskEMCALPi0CalibSelection::UserExec(Option_t* /* option */)
434 {
435   //Analysis per event.
436   
437   if(fRecoUtils->GetParticleType()!=AliEMCALRecoUtils::kPhoton){
438     printf("Wrong particle type for cluster position recalculation! = %d\n", fRecoUtils->GetParticleType());
439     abort();
440   }
441   
442   if(!(((AliESDEvent*)InputEvent())->GetFiredTriggerClasses()).Contains("EMC")) {
443     //printf("Reject Event %d, FiredClass %s\n",(Int_t)Entry(),(((AliESDEvent*)InputEvent())->GetFiredTriggerClasses()).Data());
444     return;
445   }
446
447   fhNEvents->Fill(0); //Event analyzed
448   
449   //Get the input event
450   AliVEvent* event = 0;
451   if(fFilteredInput) event = AODEvent();
452   else               event = InputEvent();
453   
454   if(!event) {
455     printf("Input event not available!\n");
456     return;
457   }
458   
459   if(DebugLevel() > 1) 
460     printf("AliAnalysisTaskEMCALPi0CalibSelection <<< %s: Event %d >>>\n",event->GetName(), (Int_t)Entry());
461   
462   
463   //Get the primary vertex
464   Double_t v[3];
465   event->GetPrimaryVertex()->GetXYZ(v) ;
466   
467   if(DebugLevel() > 1) printf("AliAnalysisTaskEMCALPi0CalibSelection Vertex: (%.3f,%.3f,%.3f)\n",v[0],v[1],v[2]);
468   
469   //Int_t runNum = aod->GetRunNumber();
470   //if(DebugLevel() > 1) printf("Run number: %d\n",runNum);
471   
472   Int_t nSM = (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules();
473   //Get the matrix with geometry information
474   if(fhNEvents->GetEntries()==1){
475     if(fLoadMatrices){
476       printf("AliAnalysisTaskEMCALPi0CalibSelection::UserExec() - Load user defined geometry matrices\n");
477       for(Int_t mod=0; mod < nSM ; mod++){
478         if(fMatrix[mod]){
479           if(DebugLevel() > 1) 
480             fMatrix[mod]->Print();
481           fEMCALGeo->SetMisalMatrix(fMatrix[mod],mod) ;   
482         }
483       }//SM loop
484     }//Load matrices
485     else if(!gGeoManager){
486       printf("AliAnalysisTaskEMCALPi0CalibSelection::UserExec() - Get geo matrices from data\n");
487       //Still not implemented in AOD, just a workaround to be able to work at least with ESDs   
488       if(!strcmp(event->GetName(),"AliAODEvent")) {
489         if(DebugLevel() > 1) 
490           printf("AliAnalysisTaskEMCALPi0CalibSelection Use ideal geometry, values geometry matrix not kept in AODs.\n");
491       }//AOD
492       else {    
493         if(DebugLevel() > 1) printf("AliAnalysisTaskEMCALPi0CalibSelection Load Misaligned matrices. \n");
494         AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event) ;
495         if(!esd) {
496           printf("AliAnalysisTaskEMCALPi0CalibSelection::UserExec() - This event does not contain ESDs?");
497           return;
498         }
499         for(Int_t mod=0; mod < nSM; mod++){
500           if(DebugLevel() > 1) 
501             esd->GetEMCALMatrix(mod)->Print();
502           if(esd->GetEMCALMatrix(mod)) fEMCALGeo->SetMisalMatrix(esd->GetEMCALMatrix(mod),mod) ;
503         } 
504       }//ESD
505     }//Load matrices from Data 
506   }//first event
507   
508   if(DebugLevel() > 1) printf("AliAnalysisTaskEMCALPi0CalibSelection Will use fLogWeight %.3f .\n",fLogWeight);
509   Int_t absId1   = -1;
510   Int_t iSupMod1 = -1;
511   Int_t iphi1    = -1;
512   Int_t ieta1    = -1;
513   Int_t absId2   = -1;
514   Int_t iSupMod2 = -1;
515   Int_t iphi2    = -1;
516   Int_t ieta2    = -1;
517   Bool_t shared  = kFALSE;
518   
519   TLorentzVector p1;
520   TLorentzVector p2;
521   TLorentzVector p12;
522   
523   //Get the list of clusters
524   TRefArray * caloClustersArr  = new TRefArray();
525   event->GetEMCALClusters(caloClustersArr);
526   const Int_t kNumberOfEMCALClusters   = caloClustersArr->GetEntries() ;
527   if(DebugLevel() > 1) printf("AliAnalysisTaskEMCALPi0CalibSelection - N CaloClusters: %d \n", kNumberOfEMCALClusters);
528   
529   // Get EMCAL cells
530   AliVCaloCells *emCells = event->GetEMCALCells();
531   
532   // loop over EMCAL clusters
533   //----------------------------------------------------------
534   // First recalibrate and recalculate energy and position
535   Float_t pos[]={0,0,0};
536   if(fCorrectClusters){
537     for(Int_t iClu=0; iClu<kNumberOfEMCALClusters; iClu++) {
538       AliVCluster *c1 = (AliVCluster *) caloClustersArr->At(iClu);
539       
540       Float_t e1i = c1->E();   // cluster energy before correction   
541       if      (e1i < fEmin) continue;
542       else if (e1i > fEmax) continue;
543       else if (c1->GetNCells() < fMinNCells) continue; 
544
545       if(fRecoUtils->ClusterContainsBadChannel(fEMCALGeo, c1->GetCellsAbsId(), c1->GetNCells())) continue;      
546       
547       if(DebugLevel() > 2)
548       { 
549         printf("Std  : i %d, E %f, dispersion %f, m02 %f, m20 %f\n",c1->GetID(),c1->E(),c1->GetDispersion(),c1->GetM02(),c1->GetM20());
550         c1->GetPosition(pos);
551         printf("Std  : i %d, x %f, y %f, z %f\n",c1->GetID(), pos[0], pos[1], pos[2]);
552       }
553       
554       //Correct cluster energy and position if requested, and not corrected previously, by default Off
555       if(fRecoUtils->IsRecalibrationOn())       {
556         fRecoUtils->RecalibrateClusterEnergy(fEMCALGeo, c1, emCells);
557         fRecoUtils->RecalculateClusterShowerShapeParameters(fEMCALGeo, emCells,c1);
558         fRecoUtils->RecalculateClusterPID(c1);
559       }
560       if(DebugLevel() > 2) 
561         printf("Energy: after recalibration %f; ",c1->E());
562       
563       // Recalculate cluster position
564       fRecoUtils->RecalculateClusterPosition(fEMCALGeo, emCells,c1);
565       
566       // Correct Non-Linearity
567      
568       c1->SetE(fRecoUtils->CorrectClusterEnergyLinearity(c1));
569       if(DebugLevel() > 2) 
570         printf("after linearity correction %f\n",c1->E());
571       
572       if(DebugLevel() > 2)
573       { 
574         printf("Cor  : i %d, E %f, dispersion %f, m02 %f, m20 %f\n",c1->GetID(),c1->E(),c1->GetDispersion(),c1->GetM02(),c1->GetM20());
575         c1->GetPosition(pos);
576         printf("Cor  : i %d, x %f, y %f, z %f\n",c1->GetID(), pos[0], pos[1], pos[2]);
577       }    
578     }    
579   }
580   
581   //----------------------------------------------------------
582   //Now the invariant mass analysis with the corrected clusters  
583   for(Int_t iClu=0; iClu<kNumberOfEMCALClusters-1; iClu++) {
584     
585     AliVCluster *c1 = (AliVCluster *) caloClustersArr->At(iClu);
586     if(fRecoUtils->ClusterContainsBadChannel(fEMCALGeo, c1->GetCellsAbsId(), c1->GetNCells())) continue;        
587     
588     Float_t e1i = c1->E();   // cluster energy before correction   
589     if      (e1i < fEmin) continue;
590     else if (e1i > fEmax) continue;
591     else if (c1->GetNCells() < fMinNCells) continue; 
592     
593     if(DebugLevel() > 2)
594     { 
595       printf("IMA  : i %d, E %f, dispersion %f, m02 %f, m20 %f\n",c1->GetID(),e1i,c1->GetDispersion(),c1->GetM02(),c1->GetM20());
596       c1->GetPosition(pos);
597       printf("IMA  : i %d, x %f, y %f, z %f\n",c1->GetID(), pos[0], pos[1], pos[2]);
598     }
599     
600     fRecoUtils->GetMaxEnergyCell(fEMCALGeo, emCells,c1,absId1,iSupMod1,ieta1,iphi1,shared);
601     c1->GetMomentum(p1,v);
602     
603     //Check if cluster is in fidutial region, not too close to borders
604     Bool_t in1 = fRecoUtils->CheckCellFiducialRegion(fEMCALGeo, c1, emCells);
605     // Clusters not facing frame structures
606     Bool_t mask1 = MaskFrameCluster(iSupMod1, ieta1);
607     //if(mask1) printf("Reject eta %d SM %d\n",ieta1, iSupMod1);
608
609     Double_t time1 = c1->GetTOF()*1.e9;
610     fhClusterTime            ->Fill(c1->E(),time1);
611     fhClusterTimeSM[iSupMod1]->Fill(c1->E(),time1);
612     
613     // Combine cluster with other clusters and get the invariant mass
614     for (Int_t jClu=iClu+1; jClu<kNumberOfEMCALClusters; jClu++) {
615       AliAODCaloCluster *c2 = (AliAODCaloCluster *) caloClustersArr->At(jClu);
616       
617       if(fRecoUtils->ClusterContainsBadChannel(fEMCALGeo, c2->GetCellsAbsId(), c2->GetNCells())) continue;      
618       
619       Float_t e2i = c2->E();
620       if      (e2i < fEmin) continue;
621       else if (e2i > fEmax) continue;
622       else if (c2->GetNCells() < fMinNCells) continue; 
623      
624       
625       fRecoUtils->GetMaxEnergyCell(fEMCALGeo, emCells,c2,absId2,iSupMod2,ieta2,iphi2,shared);
626       c2->GetMomentum(p2,v);
627
628       p12 = p1+p2;
629       Float_t invmass = p12.M()*1000; 
630       //printf("*** mass %f\n",invmass);
631
632       //Asimetry cut      
633       Float_t asym = TMath::Abs(p1.E()-p2.E())/(p1.E()+p2.E());
634       //printf("asymmetry %f\n",asym);
635       if(asym > fAsyCut) continue;
636
637       //Time cut
638       Double_t time2 = c2->GetTOF()*1.e9;
639       fhClusterPairDiffTime->Fill(p12.E(),time1-time2);
640       if(TMath::Abs(time1-time2) > fDTimeCut) continue;
641
642       if(invmass < fMaxBin && invmass > fMinBin ){
643         
644         //Check if cluster is in fidutial region, not too close to borders
645         Bool_t in2 = fRecoUtils->CheckCellFiducialRegion(fEMCALGeo, c2, emCells);
646         
647         // Clusters not facing frame structures
648         Bool_t mask2 = MaskFrameCluster(iSupMod2, ieta2);         
649         //if(mask2) printf("Reject eta %d SM %d\n",ieta2, iSupMod2);
650
651         if(in1 && in2){
652           
653           fHmgg->Fill(invmass,p12.Pt()); 
654           
655           if(iSupMod1==iSupMod2) {
656             fHmggSM[iSupMod1]->Fill(invmass,p12.Pt()); 
657             fhClusterPairDiffTimeSameSM[iSupMod1]->Fill(p12.E(),time1-time2);
658           }
659           else                   
660             fHmggDifferentSM ->Fill(invmass,p12.Pt());
661           
662           // Same sector
663           Int_t j=0;
664           for(Int_t i = 0; i < nSM/2; i++){
665             j=2*i;
666             if((iSupMod1==j && iSupMod2==j+1) || (iSupMod1==j+1 && iSupMod2==j)) {
667               fHmggPairSameSectorSM[i]->Fill(invmass,p12.Pt());
668               fhClusterPairDiffTimeSameSector[i]->Fill(p12.E(),time1-time2);
669             } 
670           }
671           
672           // Same side
673           for(Int_t i = 0; i < nSM-2; i++){
674             if((iSupMod1==i && iSupMod2==i+2) || (iSupMod1==i+2 && iSupMod2==i)) {
675               fHmggPairSameSideSM[i]->Fill(invmass,p12.Pt()); 
676               fhClusterPairDiffTimeSameSide[i]->Fill(p12.E(),time1-time2);
677             }
678           }
679
680           
681           if(!mask1 && !mask2){
682           
683             fHmggMaskFrame->Fill(invmass,p12.Pt()); 
684             
685             if(iSupMod1==iSupMod2) fHmggSMMaskFrame[iSupMod1]->Fill(invmass,p12.Pt()); 
686             else                   fHmggDifferentSMMaskFrame ->Fill(invmass,p12.Pt());
687             
688             // Same sector
689             j=0;
690             for(Int_t i = 0; i < nSM/2; i++){
691               j=2*i;
692               if((iSupMod1==j && iSupMod2==j+1) || (iSupMod1==j+1 && iSupMod2==j)) fHmggPairSameSectorSMMaskFrame[i]->Fill(invmass,p12.Pt()); 
693             }
694             
695             // Same side
696             for(Int_t i = 0; i < nSM-2; i++){
697               if((iSupMod1==i && iSupMod2==i+2) || (iSupMod1==i+2 && iSupMod2==i)) fHmggPairSameSideSMMaskFrame[i]->Fill(invmass,p12.Pt()); 
698             }
699             
700           }// Pair not facing frame
701           
702           
703           if(invmass > 100. && invmass < 160.){//restrict to clusters really close to pi0 peak
704             
705             //Opening angle of 2 photons
706             Float_t opangle = p1.Angle(p2.Vect())*TMath::RadToDeg();
707             //printf("*******>>>>>>>> In PEAK pt %f, angle %f \n",p12.Pt(),opangle);
708             
709             //Incident angle of each photon
710             Float_t inangle1 =0., inangle2=0.;
711             if(gGeoManager){
712               Float_t posSM1cen[3]={0.,0.,0.};
713               Float_t depth = fRecoUtils->GetDepth(p1.Energy(),fRecoUtils->GetParticleType(),iSupMod1);
714               fEMCALGeo->RecalculateTowerPosition(11.5, 23.5, iSupMod1, depth, fRecoUtils->GetMisalTransShiftArray(),fRecoUtils->GetMisalRotShiftArray(),posSM1cen); 
715               Float_t posSM2cen[3]={0.,0.,0.}; 
716               depth = fRecoUtils->GetDepth(p2.Energy(),fRecoUtils->GetParticleType(),iSupMod2);
717               fEMCALGeo->RecalculateTowerPosition(11.5, 23.5, iSupMod2, depth, fRecoUtils->GetMisalTransShiftArray(),fRecoUtils->GetMisalRotShiftArray(),posSM2cen); 
718               //printf("SM1 %d pos (%2.3f,%2.3f,%2.3f) \n",iSupMod1,posSM1cen[0],posSM1cen[1],posSM1cen[2]);
719               //printf("SM2 %d pos (%2.3f,%2.3f,%2.3f) \n",iSupMod2,posSM2cen[0],posSM2cen[1],posSM2cen[2]);
720               
721               TVector3 vecSM1cen(posSM1cen[0]-v[0],posSM1cen[1]-v[1],posSM1cen[2]-v[2]); 
722               TVector3 vecSM2cen(posSM2cen[0]-v[0],posSM2cen[1]-v[1],posSM2cen[2]-v[2]); 
723               inangle1 = p1.Angle(vecSM1cen)*TMath::RadToDeg();
724               inangle2 = p2.Angle(vecSM2cen)*TMath::RadToDeg();
725               //printf("Incident angle: cluster 1 %2.3f; cluster 2 %2.3f\n",inangle1,inangle2);
726             }
727             
728             fHOpeningAngle ->Fill(opangle,p12.Pt()); 
729             fHIncidentAngle->Fill(inangle1,p1.Pt()); 
730             fHIncidentAngle->Fill(inangle2,p2.Pt()); 
731             fHAsymmetry    ->Fill(asym,p12.Pt()); 
732             
733             if(iSupMod1==iSupMod2) {
734               fHOpeningAngleSM[iSupMod1] ->Fill(opangle,p12.Pt());
735               fHIncidentAngleSM[iSupMod1]->Fill(inangle1,p1.Pt());
736               fHIncidentAngleSM[iSupMod1]->Fill(inangle2,p2.Pt());
737               fHAsymmetrySM[iSupMod1]    ->Fill(asym,p12.Pt());
738             }
739             else{      
740               fHOpeningAngleDifferentSM  ->Fill(opangle,p12.Pt());
741               fHIncidentAngleDifferentSM ->Fill(inangle1,p1.Pt());
742               fHIncidentAngleDifferentSM ->Fill(inangle2,p2.Pt());
743               fHAsymmetryDifferentSM     ->Fill(asym,p12.Pt());
744             }
745             
746             if((iSupMod1==0 && iSupMod2==2) || (iSupMod1==2 && iSupMod2==0)) {
747               fHOpeningAnglePairSM[0] ->Fill(opangle,p12.Pt()); 
748               fHIncidentAnglePairSM[0]->Fill(inangle1,p1.Pt());
749               fHIncidentAnglePairSM[0]->Fill(inangle2,p2.Pt());
750               fHAsymmetryPairSM[0]    ->Fill(asym,p12.Pt());
751               
752             } 
753             if((iSupMod1==1 && iSupMod2==3) || (iSupMod1==3 && iSupMod2==1)) {
754               fHOpeningAnglePairSM[1] ->Fill(opangle,p12.Pt()); 
755               fHIncidentAnglePairSM[1]->Fill(inangle1,p1.Pt());
756               fHIncidentAnglePairSM[1]->Fill(inangle2,p2.Pt());
757               fHAsymmetryPairSM[1]    ->Fill(asym,p12.Pt());
758               
759             }
760             
761             if((iSupMod1==0 && iSupMod2==1) || (iSupMod1==1 && iSupMod2==0)) {
762               fHOpeningAnglePairSM[2] ->Fill(opangle,p12.Pt()); 
763               fHIncidentAnglePairSM[2]->Fill(inangle1,p1.Pt());
764               fHIncidentAnglePairSM[2]->Fill(inangle2,p2.Pt());
765               fHAsymmetryPairSM[2]    ->Fill(asym,p12.Pt());
766               
767               
768             }
769             if((iSupMod1==2 && iSupMod2==3) || (iSupMod1==3 && iSupMod2==2)) {
770               fHOpeningAnglePairSM[3] ->Fill(opangle,p12.Pt()); 
771               fHIncidentAnglePairSM[3]->Fill(inangle1,p1.Pt());
772               fHIncidentAnglePairSM[3]->Fill(inangle2,p2.Pt());
773               fHAsymmetryPairSM[3]    ->Fill(asym,p12.Pt());
774             }
775             
776           }// pair in 100 < mass < 160
777           
778         }//in acceptance cuts
779         
780         //In case of filling only channels with second cluster in same SM
781         if(fSameSM && iSupMod1!=iSupMod2) continue;
782         
783         if (fGroupNCells == 0){
784           fHmpi0[iSupMod1][ieta1][iphi1]->Fill(invmass);
785           fHmpi0[iSupMod2][ieta2][iphi2]->Fill(invmass);
786           
787           if(invmass > 100. && invmass < 160.){//restrict to clusters really close to pi0 peak
788             fhTowerDecayPhotonHit      [iSupMod1]->Fill(ieta1,iphi1);
789             fhTowerDecayPhotonEnergy   [iSupMod1]->Fill(ieta1,iphi1,p1.E());
790             fhTowerDecayPhotonAsymmetry[iSupMod1]->Fill(ieta1,iphi1,asym);
791             
792             fhTowerDecayPhotonHit      [iSupMod2]->Fill(ieta2,iphi2);
793             fhTowerDecayPhotonEnergy   [iSupMod2]->Fill(ieta2,iphi2,p2.E());
794             fhTowerDecayPhotonAsymmetry[iSupMod2]->Fill(ieta2,iphi2,asym);
795             
796             if(!mask1)fhTowerDecayPhotonHitMaskFrame[iSupMod1]->Fill(ieta1,iphi1);
797             if(!mask2)fhTowerDecayPhotonHitMaskFrame[iSupMod2]->Fill(ieta2,iphi2);
798             
799           }// pair in mass of pi0
800         }       
801         else  {
802           //printf("Regroup N %d, eta1 %d, phi1 %d, eta2 %d, phi2 %d \n",fGroupNCells, ieta1, iphi1, ieta2, iphi2);
803           for (Int_t i = -fGroupNCells; i < fGroupNCells+1; i++) {
804             for (Int_t j = -fGroupNCells; j < fGroupNCells+1; j++) {
805               //printf("\t i %d, j %d\n",i,j);
806               if((ieta1+i >= 0) && (iphi1+j >= 0) && (ieta1+i < 48) && (iphi1+j < 24)){
807                 //printf("\t \t eta1+i %d, phi1+j %d\n", ieta1+i, iphi1+j);
808                 fHmpi0[iSupMod1][ieta1+i][iphi1+j]->Fill(invmass);
809               }
810               if((ieta2+i >= 0) && (iphi2+j >= 0) && (ieta2+i < 48) && (iphi2+j < 24)){
811                 //printf("\t \t eta2+i %d, phi2+j %d\n", ieta2+i, iphi2+j);
812                 fHmpi0[iSupMod2][ieta2+i][iphi2+j]->Fill(invmass);
813               }
814             }// j loop
815           }//i loop
816         }//group cells
817         
818         if(DebugLevel() > 1) printf("AliAnalysisTaskEMCALPi0CalibSelection 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",
819                                     iSupMod1,iphi1,ieta1,iSupMod2,iphi2,ieta2,p12.M(),e1i,c1->E(),e2i,c2->E());
820       }
821       
822     }
823     
824   } // end of loop over EMCAL clusters
825   
826   delete caloClustersArr;
827   
828   PostData(1,fOutputContainer);
829   
830 }
831
832 //_____________________________________________________
833 void AliAnalysisTaskEMCALPi0CalibSelection::PrintInfo(){
834   
835   //Print settings
836   printf("Cluster cuts: %2.2f < E < %2.2f GeV; number of cells > %d; Assymetry < %1.2f, pair time diff < %2.2f\n", 
837          fEmin,fEmax, fMinNCells, fAsyCut, fDTimeCut) ;
838   printf("Group %d cells\n", fGroupNCells) ;
839   printf("Cluster maximal cell away from border at least %d cells\n", fRecoUtils->GetNumberOfCellsFromEMCALBorder()) ;
840   printf("Histograms: bins %d; energy range: %2.2f < E < %2.2f GeV\n",fNbins,fMinBin,fMaxBin) ;
841   printf("Switchs:\n \t Remove Bad Channels? %d; Use filtered input? %d;  Correct Clusters? %d, \n \t Mass per channel same SM clusters? %d\n",
842          fRecoUtils->IsBadChannelsRemovalSwitchedOn(),fFilteredInput,fCorrectClusters, fSameSM) ;
843   printf("EMCAL Geometry name: < %s >, Load Matrices %d\n",fEMCALGeoName.Data(), fLoadMatrices) ;
844   if(fLoadMatrices) {for(Int_t ism = 0; ism < AliEMCALGeoParams::fgkEMCALModules; ism++) fMatrix[ism]->Print() ; }
845   
846   
847 }
848