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