Fix Coverity
[u/mrichter/AliRoot.git] / PWG4 / CaloCalib / AliAnalysisTaskEMCALPi0CalibSelection.cxx
CommitLineData
375cec9b 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: Boris Polishchuk *
44cf05d7 5 * Adapted to AOD reading by Gustavo Conesa *
375cec9b 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"
375cec9b 29#include "TRefArray.h"
30#include "TList.h"
31#include "TH1F.h"
247abff4 32#include <TGeoManager.h>
375cec9b 33
34// AliRoot
35#include "AliAnalysisTaskEMCALPi0CalibSelection.h"
36#include "AliAODEvent.h"
37#include "AliESDEvent.h"
375cec9b 38#include "AliEMCALGeometry.h"
c8fe2783 39#include "AliVCluster.h"
40#include "AliVCaloCells.h"
9584c261 41#include "AliEMCALRecoUtils.h"
375cec9b 42
43ClassImp(AliAnalysisTaskEMCALPi0CalibSelection)
44
375cec9b 45
46//__________________________________________________
47AliAnalysisTaskEMCALPi0CalibSelection::AliAnalysisTaskEMCALPi0CalibSelection(const char* name) :
49b53920 48AliAnalysisTaskSE(name),fEMCALGeo(0x0),
49fEmin(0.5), fEmax(15.),
a7e5a381 50fL0min(0.01), fL0max(0.5), fDTimeCut(100.),
51fAsyCut(1.), fMinNCells(2), fGroupNCells(0),
52fLogWeight(4.5), fSameSM(kFALSE), fFilteredInput(kFALSE),
49b53920 53fCorrectClusters(kFALSE), fEMCALGeoName("EMCAL_COMPLETEV1"),
54fTriggerName("EMC"), fRecoUtils(new AliEMCALRecoUtils),
55fCuts(0x0), fLoadMatrices(0),
56fNMaskCellColumns(11), fMaskCellColumns(0x0),
a7e5a381 57fInvMassCutMin(110.), fInvMassCutMax(160.),
49b53920 58//Histograms
a7e5a381 59fOutputContainer(0x0), fNbins(300),
60fMinBin(0.), fMaxBin(300.),
61fNTimeBins(1000), fMinTimeBin(0.), fMaxTimeBin(1000.),
49b53920 62fHmgg(0x0), fHmggDifferentSM(0x0),
63fHmggMaskFrame(0x0), fHmggDifferentSMMaskFrame(0x0),
64fHOpeningAngle(0x0), fHOpeningAngleDifferentSM(0x0),
65fHIncidentAngle(0x0), fHIncidentAngleDifferentSM(0x0),
66fHAsymmetry(0x0), fHAsymmetryDifferentSM(0x0),
67fhNEvents(0x0),
68fhClusterTime(0x0), fhClusterPairDiffTime(0x0)
375cec9b 69{
70 //Named constructor which should be used.
71
44cf05d7 72 for(Int_t iMod=0; iMod < AliEMCALGeoParams::fgkEMCALModules; iMod++) {
bdd2a262 73 for(Int_t iX=0; iX<24; iX++) {
74 for(Int_t iZ=0; iZ<48; iZ++) {
a7e5a381 75 fHmpi0[iMod][iZ][iX] = 0 ;
375cec9b 76 }
77 }
78 }
6eb2a715 79
a7e5a381 80 fHTpi0[0]= 0 ;
81 fHTpi0[1]= 0 ;
82 fHTpi0[2]= 0 ;
83 fHTpi0[3]= 0 ;
84
85
42b19289 86 fMaskCellColumns = new Int_t[fNMaskCellColumns];
42b19289 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;
af2d7c9b 96 fhClusterPairDiffTimeSameSector[iSMPair]= 0;
42b19289 97 }
98 for(Int_t iSMPair = 0; iSMPair < AliEMCALGeoParams::fgkEMCALModules-2; iSMPair++){
99 fHmggPairSameSideSM[iSMPair] = 0;
100 fHmggPairSameSideSMMaskFrame[iSMPair] = 0;
af2d7c9b 101 fhClusterPairDiffTimeSameSide[iSMPair] = 0;
42b19289 102 }
1dabc151 103
44cf05d7 104 for(Int_t iSM = 0; iSM < AliEMCALGeoParams::fgkEMCALModules; iSM++) {
105 fHmggSM[iSM] = 0;
42b19289 106 fHmggSMMaskFrame[iSM] = 0;
44cf05d7 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;
42b19289 116 fhTowerDecayPhotonHitMaskFrame[iSM]= 0;
44cf05d7 117 fMatrix[iSM] = 0x0;
af2d7c9b 118 fhClusterTimeSM[iSM] = 0;
119 fhClusterPairDiffTimeSameSM[iSM] = 0;
2dfb1428 120 }
121
cf028690 122 DefineOutput(1, TList::Class());
6eb2a715 123 DefineOutput(2, TList::Class()); // will contain cuts or local params
af2d7c9b 124
375cec9b 125}
126
127//__________________________________________________
128AliAnalysisTaskEMCALPi0CalibSelection::~AliAnalysisTaskEMCALPi0CalibSelection()
129{
130 //Destructor.
131
132 if(fOutputContainer){
133 fOutputContainer->Delete() ;
134 delete fOutputContainer ;
135 }
af2d7c9b 136
42b19289 137 if(fEMCALGeo) delete fEMCALGeo ;
138 if(fRecoUtils) delete fRecoUtils ;
139 if(fNMaskCellColumns) delete [] fMaskCellColumns;
140
375cec9b 141}
142
6eb2a715 143//_____________________________________________________
144void AliAnalysisTaskEMCALPi0CalibSelection::LocalInit()
145{
af2d7c9b 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
a7e5a381 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) ;
af2d7c9b 155 fCuts->Add(new TObjString(onePar));
156 snprintf(onePar,buffersize, "Group %d cells;", fGroupNCells) ;
157 fCuts->Add(new TObjString(onePar));
247abff4 158 snprintf(onePar,buffersize, "Cluster maximal cell away from border at least %d cells;", fRecoUtils->GetNumberOfCellsFromEMCALBorder()) ;
af2d7c9b 159 fCuts->Add(new TObjString(onePar));
a7e5a381 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) ;
af2d7c9b 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 ",
49b53920 165 fRecoUtils->IsBadChannelsRemovalSwitchedOn(),fFilteredInput,fCorrectClusters, fSameSM) ;
af2d7c9b 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
6eb2a715 175}
375cec9b 176
375cec9b 177//__________________________________________________
178void AliAnalysisTaskEMCALPi0CalibSelection::UserCreateOutputObjects()
179{
247abff4 180 //Create output container, init geometry
af2d7c9b 181
cf028690 182 fEMCALGeo = AliEMCALGeometry::GetInstance(fEMCALGeoName) ;
44cf05d7 183 Int_t nSM = (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules();
9584c261 184
375cec9b 185 fOutputContainer = new TList();
2dfb1428 186 const Int_t buffersize = 255;
187 char hname[buffersize], htitl[buffersize];
375cec9b 188
44cf05d7 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++) {
2dfb1428 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);
70ae4900 194 fHmpi0[iMod][iCol][iRow] = new TH1F(hname,htitl,fNbins,fMinBin,fMaxBin);
a7e5a381 195 fHmpi0[iMod][iCol][iRow]->SetXTitle("mass (MeV/c^{2})");
70ae4900 196 fOutputContainer->Add(fHmpi0[iMod][iCol][iRow]);
375cec9b 197 }
198 }
199 }
af2d7c9b 200
a7e5a381 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
70ae4900 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)");
375cec9b 213 fOutputContainer->Add(fHmgg);
49b53920 214
2dfb1428 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);
49b53920 219
9584c261 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
2dfb1428 250
1dabc151 251 //TString pairname[] = {"A side (0-2)", "C side (1-3)","Row 0 (0-1)", "Row 1 (2-3)"};
2dfb1428 252
42b19289 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
44cf05d7 265 for(Int_t iSM = 0; iSM < nSM; iSM++) {
2dfb1428 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)");
2dfb1428 272 fOutputContainer->Add(fHmggSM[iSM]);
273
42b19289 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
1dabc151 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]);
42b19289 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]);
49b53920 296
af2d7c9b 297 fhClusterPairDiffTimeSameSector[iSM] = new TH2F(Form("hClusterPairDiffTimeSameSector%d",iSM),
49b53920 298 Form("cluster pair time difference vs E, Sector %d",iSM),
299 100,0,10, 200,-100,100);
af2d7c9b 300 fhClusterPairDiffTimeSameSector[iSM]->SetXTitle("E_{pair} (GeV)");
301 fhClusterPairDiffTimeSameSector[iSM]->SetYTitle("#Delta t (ns)");
302 fOutputContainer->Add(fhClusterPairDiffTimeSameSector[iSM]);
49b53920 303
304
1dabc151 305 }
9584c261 306
1dabc151 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]);
42b19289 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)");
af2d7c9b 320 fOutputContainer->Add(fHmggPairSameSideSMMaskFrame[iSM]);
49b53920 321
af2d7c9b 322 fhClusterPairDiffTimeSameSide[iSM] = new TH2F(Form("hClusterPairDiffTimeSameSide%d",iSM),
49b53920 323 Form("cluster pair time difference vs E, Side %d",iSM),
324 100,0,10, 200,-100,100);
af2d7c9b 325 fhClusterPairDiffTimeSameSide[iSM]->SetXTitle("E_{pair} (GeV)");
326 fhClusterPairDiffTimeSameSide[iSM]->SetYTitle("#Delta t (ns)");
327 fOutputContainer->Add(fhClusterPairDiffTimeSameSide[iSM]);
49b53920 328
1dabc151 329 }
9584c261 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);
1dabc151 339 snprintf(htitl,buffersize, "Opening angle for SM pair: %d",iSM);
9584c261 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);
1dabc151 353 snprintf(htitl,buffersize, "Incident angle for SM pair: %d",iSM);
9584c261 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);
1dabc151 360 snprintf(htitl, buffersize, "Asymmetry for super mod %d",iSM);
9584c261 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);
1dabc151 367 snprintf(htitl,buffersize, "Asymmetry for SM pair: %d",iSM);
9584c261 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
9584c261 373 Int_t colmax = 48;
374 Int_t rowmax = 24;
375
af2d7c9b 376 fhTowerDecayPhotonHit[iSM] = new TH2F (Form("hTowerDecPhotonHit_Mod%d",iSM),
49b53920 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);
9584c261 379 fhTowerDecayPhotonHit[iSM]->SetYTitle("row (phi direction)");
380 fhTowerDecayPhotonHit[iSM]->SetXTitle("column (eta direction)");
381 fOutputContainer->Add(fhTowerDecayPhotonHit[iSM]);
382
af2d7c9b 383 fhTowerDecayPhotonEnergy[iSM] = new TH2F (Form("hTowerDecPhotonEnergy_Mod%d",iSM),
49b53920 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);
9584c261 386 fhTowerDecayPhotonEnergy[iSM]->SetYTitle("row (phi direction)");
387 fhTowerDecayPhotonEnergy[iSM]->SetXTitle("column (eta direction)");
388 fOutputContainer->Add(fhTowerDecayPhotonEnergy[iSM]);
389
af2d7c9b 390 fhTowerDecayPhotonAsymmetry[iSM] = new TH2F (Form("hTowerDecPhotonAsymmetry_Mod%d",iSM),
49b53920 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);
9584c261 393 fhTowerDecayPhotonAsymmetry[iSM]->SetYTitle("row (phi direction)");
394 fhTowerDecayPhotonAsymmetry[iSM]->SetXTitle("column (eta direction)");
395 fOutputContainer->Add(fhTowerDecayPhotonAsymmetry[iSM]);
396
42b19289 397 fhTowerDecayPhotonHitMaskFrame[iSM] = new TH2F (Form("hTowerDecPhotonHit_Mod%d_MaskFrame",iSM),Form("Entries in grid of cells in Module %d",iSM),
49b53920 398 colmax+2,-1.5,colmax+0.5, rowmax+2,-1.5,rowmax+0.5);
42b19289 399 fhTowerDecayPhotonHitMaskFrame[iSM]->SetYTitle("row (phi direction)");
400 fhTowerDecayPhotonHitMaskFrame[iSM]->SetXTitle("column (eta direction)");
401 fOutputContainer->Add(fhTowerDecayPhotonHitMaskFrame[iSM]);
49b53920 402
af2d7c9b 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]);
42b19289 407
af2d7c9b 408 fhClusterPairDiffTimeSameSM[iSM] = new TH2F(Form("hClusterPairDiffTimeSameSM%d",iSM),
49b53920 409 Form("cluster pair time difference vs E, SM %d",iSM),
410 100,0,10, 200,-100,100);
af2d7c9b 411 fhClusterPairDiffTimeSameSM[iSM]->SetXTitle("E (GeV)");
412 fhClusterPairDiffTimeSameSM[iSM]->SetYTitle("#Delta t (ns)");
413 fOutputContainer->Add(fhClusterPairDiffTimeSameSM[iSM]);
49b53920 414
2dfb1428 415 }
6eb2a715 416
af2d7c9b 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);
49b53920 421
a7e5a381 422 fhClusterPairDiffTime = new TH2F("hClusterPairDiffTime","cluster pair time difference vs E",100,0,10, 800,-400,400);
af2d7c9b 423 fhClusterPairDiffTime->SetXTitle("E_{pair} (GeV)");
424 fhClusterPairDiffTime->SetYTitle("#Delta t (ns)");
425 fOutputContainer->Add(fhClusterPairDiffTime);
49b53920 426
af2d7c9b 427
6eb2a715 428 fhNEvents = new TH1I("hNEvents", "Number of analyzed events" , 1 , 0 , 1 ) ;
429 fOutputContainer->Add(fhNEvents);
247abff4 430
431 fOutputContainer->SetOwner(kTRUE);
432
af2d7c9b 433 // fCalibData = new AliEMCALCalibData();
434
cf028690 435 PostData(1,fOutputContainer);
49b53920 436
375cec9b 437}
438
42b19289 439//__________________________________________________
440Bool_t AliAnalysisTaskEMCALPi0CalibSelection::MaskFrameCluster(const Int_t iSM, const Int_t ieta) const {
af2d7c9b 441 //Check if cell is in one of the regions where we have significant amount of material in front of EMCAL
42b19289 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 }
af2d7c9b 451
42b19289 452 return kFALSE;
453
454}
455
375cec9b 456//__________________________________________________
457void AliAnalysisTaskEMCALPi0CalibSelection::UserExec(Option_t* /* option */)
458{
459 //Analysis per event.
375cec9b 460
19db8f8c 461 if(fRecoUtils->GetParticleType()!=AliEMCALRecoUtils::kPhoton){
462 printf("Wrong particle type for cluster position recalculation! = %d\n", fRecoUtils->GetParticleType());
463 abort();
464 }
247abff4 465
afaaef51 466 if(!(((AliESDEvent*)InputEvent())->GetFiredTriggerClasses()).Contains(fTriggerName)) {
af2d7c9b 467 //printf("Reject Event %d, FiredClass %s\n",(Int_t)Entry(),(((AliESDEvent*)InputEvent())->GetFiredTriggerClasses()).Data());
468 return;
469 }
49b53920 470
6eb2a715 471 fhNEvents->Fill(0); //Event analyzed
472
247abff4 473 //Get the input event
474 AliVEvent* event = 0;
475 if(fFilteredInput) event = AODEvent();
476 else event = InputEvent();
70ae4900 477
247abff4 478 if(!event) {
479 printf("Input event not available!\n");
480 return;
375cec9b 481 }
af2d7c9b 482
247abff4 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) ;
375cec9b 490
375cec9b 491 if(DebugLevel() > 1) printf("AliAnalysisTaskEMCALPi0CalibSelection Vertex: (%.3f,%.3f,%.3f)\n",v[0],v[1],v[2]);
492
247abff4 493 //Int_t runNum = aod->GetRunNumber();
494 //if(DebugLevel() > 1) printf("Run number: %d\n",runNum);
375cec9b 495
1dabc151 496 Int_t nSM = (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules();
375cec9b 497 //Get the matrix with geometry information
3b13c34c 498 if(fhNEvents->GetEntries()==1){
499 if(fLoadMatrices){
500 printf("AliAnalysisTaskEMCALPi0CalibSelection::UserExec() - Load user defined geometry matrices\n");
1dabc151 501 for(Int_t mod=0; mod < nSM ; mod++){
3b13c34c 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")) {
f2ccb5b8 513 if(DebugLevel() > 1)
514 printf("AliAnalysisTaskEMCALPi0CalibSelection Use ideal geometry, values geometry matrix not kept in AODs.\n");
3b13c34c 515 }//AOD
516 else {
517 if(DebugLevel() > 1) printf("AliAnalysisTaskEMCALPi0CalibSelection Load Misaligned matrices. \n");
f2ccb5b8 518 AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event) ;
519 if(!esd) {
520 printf("AliAnalysisTaskEMCALPi0CalibSelection::UserExec() - This event does not contain ESDs?");
3b13c34c 521 return;
f2ccb5b8 522 }
1dabc151 523 for(Int_t mod=0; mod < nSM; mod++){
af2d7c9b 524 if(DebugLevel() > 1)
3b13c34c 525 esd->GetEMCALMatrix(mod)->Print();
f2ccb5b8 526 if(esd->GetEMCALMatrix(mod)) fEMCALGeo->SetMisalMatrix(esd->GetEMCALMatrix(mod),mod) ;
527 }
3b13c34c 528 }//ESD
529 }//Load matrices from Data
f2ccb5b8 530 }//first event
375cec9b 531
532 if(DebugLevel() > 1) printf("AliAnalysisTaskEMCALPi0CalibSelection Will use fLogWeight %.3f .\n",fLogWeight);
19db8f8c 533 Int_t absId1 = -1;
6eb2a715 534 Int_t iSupMod1 = -1;
535 Int_t iphi1 = -1;
536 Int_t ieta1 = -1;
19db8f8c 537 Int_t absId2 = -1;
6eb2a715 538 Int_t iSupMod2 = -1;
539 Int_t iphi2 = -1;
540 Int_t ieta2 = -1;
af2d7c9b 541 Bool_t shared = kFALSE;
3b13c34c 542
375cec9b 543 TLorentzVector p1;
544 TLorentzVector p2;
545 TLorentzVector p12;
546
247abff4 547 //Get the list of clusters
375cec9b 548 TRefArray * caloClustersArr = new TRefArray();
42b19289 549 event->GetEMCALClusters(caloClustersArr);
375cec9b 550 const Int_t kNumberOfEMCALClusters = caloClustersArr->GetEntries() ;
cf028690 551 if(DebugLevel() > 1) printf("AliAnalysisTaskEMCALPi0CalibSelection - N CaloClusters: %d \n", kNumberOfEMCALClusters);
375cec9b 552
247abff4 553 // Get EMCAL cells
554 AliVCaloCells *emCells = event->GetEMCALCells();
555
375cec9b 556 // loop over EMCAL clusters
247abff4 557 //----------------------------------------------------------
558 // First recalibrate and recalculate energy and position
559 Float_t pos[]={0,0,0};
560 if(fCorrectClusters){
9a2c875d 561 for(Int_t iClu=0; iClu<kNumberOfEMCALClusters; iClu++) {
247abff4 562 AliVCluster *c1 = (AliVCluster *) caloClustersArr->At(iClu);
563
af2d7c9b 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;
49b53920 568 else if (c1->GetM02() < fL0min || c1->GetM02() > fL0max) continue;
569
247abff4 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
5ef94e1b 580 if(fRecoUtils->IsRecalibrationOn()) {
581 fRecoUtils->RecalibrateClusterEnergy(fEMCALGeo, c1, emCells);
582 fRecoUtils->RecalculateClusterShowerShapeParameters(fEMCALGeo, emCells,c1);
583 fRecoUtils->RecalculateClusterPID(c1);
584 }
247abff4 585 if(DebugLevel() > 2)
49b53920 586 printf("Energy: after recalibration %f; \n",c1->E());
247abff4 587
44907916 588 // Recalculate cluster position
589 fRecoUtils->RecalculateClusterPosition(fEMCALGeo, emCells,c1);
590
247abff4 591 // Correct Non-Linearity
592 c1->SetE(fRecoUtils->CorrectClusterEnergyLinearity(c1));
593 if(DebugLevel() > 2)
49b53920 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());
44907916 600
247abff4 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 }
5ef94e1b 609
247abff4 610 //----------------------------------------------------------
a7e5a381 611 //Now the invariant mass analysis with the corrected clusters
612 Int_t bc = event->GetBunchCrossNumber();
247abff4 613 for(Int_t iClu=0; iClu<kNumberOfEMCALClusters-1; iClu++) {
375cec9b 614
247abff4 615 AliVCluster *c1 = (AliVCluster *) caloClustersArr->At(iClu);
616 if(fRecoUtils->ClusterContainsBadChannel(fEMCALGeo, c1->GetCellsAbsId(), c1->GetNCells())) continue;
6eb2a715 617
618 Float_t e1i = c1->E(); // cluster energy before correction
1dabc151 619 if (e1i < fEmin) continue;
620 else if (e1i > fEmax) continue;
a7e5a381 621 else if (!fRecoUtils->IsGoodCluster(c1,fEMCALGeo,emCells,bc));
6eb2a715 622 else if (c1->GetNCells() < fMinNCells) continue;
49b53920 623 else if (c1->GetM02() < fL0min || c1->GetM02() > fL0max) continue;
624
6eb2a715 625 if(DebugLevel() > 2)
70ae4900 626 {
247abff4 627 printf("IMA : i %d, E %f, dispersion %f, m02 %f, m20 %f\n",c1->GetID(),e1i,c1->GetDispersion(),c1->GetM02(),c1->GetM20());
70ae4900 628 c1->GetPosition(pos);
247abff4 629 printf("IMA : i %d, x %f, y %f, z %f\n",c1->GetID(), pos[0], pos[1], pos[2]);
70ae4900 630 }
6eb2a715 631
3b13c34c 632 fRecoUtils->GetMaxEnergyCell(fEMCALGeo, emCells,c1,absId1,iSupMod1,ieta1,iphi1,shared);
9584c261 633 c1->GetMomentum(p1,v);
247abff4 634
42b19289 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);
49b53920 640
af2d7c9b 641 Double_t time1 = c1->GetTOF()*1.e9;
642 fhClusterTime ->Fill(c1->E(),time1);
643 fhClusterTimeSM[iSupMod1]->Fill(c1->E(),time1);
42b19289 644
247abff4 645 // Combine cluster with other clusters and get the invariant mass
646 for (Int_t jClu=iClu+1; jClu<kNumberOfEMCALClusters; jClu++) {
375cec9b 647 AliAODCaloCluster *c2 = (AliAODCaloCluster *) caloClustersArr->At(jClu);
af2d7c9b 648
a7e5a381 649
375cec9b 650 Float_t e2i = c2->E();
1dabc151 651 if (e2i < fEmin) continue;
6eb2a715 652 else if (e2i > fEmax) continue;
a7e5a381 653 else if (!fRecoUtils->IsGoodCluster(c2,fEMCALGeo,emCells,bc))continue;
6eb2a715 654 else if (c2->GetNCells() < fMinNCells) continue;
49b53920 655 else if (c2->GetM02() < fL0min || c2->GetM02() > fL0max) continue;
656
19db8f8c 657
3b13c34c 658 fRecoUtils->GetMaxEnergyCell(fEMCALGeo, emCells,c2,absId2,iSupMod2,ieta2,iphi2,shared);
6eb2a715 659 c2->GetMomentum(p2,v);
49b53920 660
375cec9b 661 p12 = p1+p2;
662 Float_t invmass = p12.M()*1000;
9584c261 663 //printf("*** mass %f\n",invmass);
49b53920 664
af2d7c9b 665 //Asimetry cut
9584c261 666 Float_t asym = TMath::Abs(p1.E()-p2.E())/(p1.E()+p2.E());
667 //printf("asymmetry %f\n",asym);
9fdaff9a 668 if(asym > fAsyCut) continue;
49b53920 669
af2d7c9b 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;
49b53920 674
af2d7c9b 675 if(invmass < fMaxBin && invmass > fMinBin ){
70ae4900 676
cfce8d44 677 //Check if cluster is in fidutial region, not too close to borders
247abff4 678 Bool_t in2 = fRecoUtils->CheckCellFiducialRegion(fEMCALGeo, c2, emCells);
679
42b19289 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);
49b53920 683
cfce8d44 684 if(in1 && in2){
685
686 fHmgg->Fill(invmass,p12.Pt());
247abff4 687
af2d7c9b 688 if(iSupMod1==iSupMod2) {
49b53920 689 fHmggSM[iSupMod1]->Fill(invmass,p12.Pt());
690 fhClusterPairDiffTimeSameSM[iSupMod1]->Fill(p12.E(),time1-time2);
691 }
af2d7c9b 692 else
49b53920 693 fHmggDifferentSM ->Fill(invmass,p12.Pt());
247abff4 694
1dabc151 695 // Same sector
696 Int_t j=0;
697 for(Int_t i = 0; i < nSM/2; i++){
698 j=2*i;
af2d7c9b 699 if((iSupMod1==j && iSupMod2==j+1) || (iSupMod1==j+1 && iSupMod2==j)) {
49b53920 700 fHmggPairSameSectorSM[i]->Fill(invmass,p12.Pt());
701 fhClusterPairDiffTimeSameSector[i]->Fill(p12.E(),time1-time2);
702 }
1dabc151 703 }
9584c261 704
1dabc151 705 // Same side
706 for(Int_t i = 0; i < nSM-2; i++){
af2d7c9b 707 if((iSupMod1==i && iSupMod2==i+2) || (iSupMod1==i+2 && iSupMod2==i)) {
49b53920 708 fHmggPairSameSideSM[i]->Fill(invmass,p12.Pt());
709 fhClusterPairDiffTimeSameSide[i]->Fill(p12.E(),time1-time2);
710 }
1dabc151 711 }
42b19289 712
42b19289 713
49b53920 714 if(!mask1 && !mask2){
715
42b19289 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
a7e5a381 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 }
9584c261 749
a7e5a381 750 //Opening angle of 2 photons
9584c261 751 Float_t opangle = p1.Angle(p2.Vect())*TMath::RadToDeg();
752 //printf("*******>>>>>>>> In PEAK pt %f, angle %f \n",p12.Pt(),opangle);
247abff4 753
f2ccb5b8 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 }
42b19289 772
9584c261 773 fHOpeningAngle ->Fill(opangle,p12.Pt());
774 fHIncidentAngle->Fill(inangle1,p1.Pt());
775 fHIncidentAngle->Fill(inangle2,p2.Pt());
776 fHAsymmetry ->Fill(asym,p12.Pt());
247abff4 777
9584c261 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());
247abff4 796
9584c261 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());
247abff4 803
9584c261 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());
247abff4 811
812
9584c261 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 }
247abff4 820
9584c261 821 }// pair in 100 < mass < 160
247abff4 822
9584c261 823 }//in acceptance cuts
2dfb1428 824
825 //In case of filling only channels with second cluster in same SM
826 if(fSameSM && iSupMod1!=iSupMod2) continue;
827
70ae4900 828 if (fGroupNCells == 0){
247abff4 829 fHmpi0[iSupMod1][ieta1][iphi1]->Fill(invmass);
830 fHmpi0[iSupMod2][ieta2][iphi2]->Fill(invmass);
9584c261 831
247abff4 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
42b19289 841 if(!mask1)fhTowerDecayPhotonHitMaskFrame[iSupMod1]->Fill(ieta1,iphi1);
842 if(!mask2)fhTowerDecayPhotonHitMaskFrame[iSupMod2]->Fill(ieta2,iphi2);
843
247abff4 844 }// pair in mass of pi0
70ae4900 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);
49b53920 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)){
70ae4900 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 }
49b53920 867 if(ok2 && (ieta2+i >= 0) && (iphi2+j >= 0) && (ieta2+i < 48) && (iphi2+j < 24)){
70ae4900 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",
247abff4 876 iSupMod1,iphi1,ieta1,iSupMod2,iphi2,ieta2,p12.M(),e1i,c1->E(),e2i,c2->E());
6eb2a715 877 }
878
375cec9b 879 }
880
881 } // end of loop over EMCAL clusters
882
883 delete caloClustersArr;
6eb2a715 884
375cec9b 885 PostData(1,fOutputContainer);
6eb2a715 886
375cec9b 887}
cfce8d44 888
5ef94e1b 889//_____________________________________________________
890void AliAnalysisTaskEMCALPi0CalibSelection::PrintInfo(){
891
892 //Print settings
af2d7c9b 893 printf("Cluster cuts: %2.2f < E < %2.2f GeV; number of cells > %d; Assymetry < %1.2f, pair time diff < %2.2f\n",
49b53920 894 fEmin,fEmax, fMinNCells, fAsyCut, fDTimeCut) ;
af2d7c9b 895 printf("Group %d cells\n", fGroupNCells) ;
5ef94e1b 896 printf("Cluster maximal cell away from border at least %d cells\n", fRecoUtils->GetNumberOfCellsFromEMCALBorder()) ;
af2d7c9b 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",
49b53920 899 fRecoUtils->IsBadChannelsRemovalSwitchedOn(),fFilteredInput,fCorrectClusters, fSameSM) ;
af2d7c9b 900 printf("EMCAL Geometry name: < %s >, Load Matrices %d\n",fEMCALGeoName.Data(), fLoadMatrices) ;
42b19289 901 if(fLoadMatrices) {for(Int_t ism = 0; ism < AliEMCALGeoParams::fgkEMCALModules; ism++) fMatrix[ism]->Print() ; }
af2d7c9b 902
5ef94e1b 903
904}
905