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