1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: Boris Polishchuk *
5 * Adapted to AOD reading by Gustavo Conesa *
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 **************************************************************************/
16 //---------------------------------------------------------------------------//
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.//
23 //---------------------------------------------------------------------------//
26 //#include <Riostream.h>
28 #include "TLorentzVector.h"
29 #include "TRefArray.h"
32 #include <TGeoManager.h>
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"
43 ClassImp(AliAnalysisTaskEMCALPi0CalibSelection)
46 //__________________________________________________
47 AliAnalysisTaskEMCALPi0CalibSelection::AliAnalysisTaskEMCALPi0CalibSelection(const char* name) :
48 AliAnalysisTaskSE(name),fEMCALGeo(0x0),
49 fEmin(0.5), fEmax(15.), fDTimeCut(20.),
50 fAsyCut(1.), fMinNCells(2), fGroupNCells(0),
51 fLogWeight(4.5), fSameSM(kFALSE), fFilteredInput(kFALSE),
52 fCorrectClusters(kFALSE), fEMCALGeoName("EMCAL_COMPLETEV1"),
53 fTriggerName("EMC"), fRecoUtils(new AliEMCALRecoUtils),
54 fCuts(0x0), fLoadMatrices(0),
55 fNMaskCellColumns(11), fMaskCellColumns(0x0),
57 fNbins(300), fMinBin(0.), fMaxBin(300.), fOutputContainer(0x0),
58 fHmgg(0x0), fHmggDifferentSM(0x0),
59 fHmggMaskFrame(0x0), fHmggDifferentSMMaskFrame(0x0),
60 fHOpeningAngle(0x0), fHOpeningAngleDifferentSM(0x0),
61 fHIncidentAngle(0x0), fHIncidentAngleDifferentSM(0x0),
62 fHAsymmetry(0x0), fHAsymmetryDifferentSM(0x0),
64 fhClusterTime(0x0), fhClusterPairDiffTime(0x0)
66 //Named constructor which should be used.
68 for(Int_t iMod=0; iMod < AliEMCALGeoParams::fgkEMCALModules; iMod++) {
69 for(Int_t iX=0; iX<24; iX++) {
70 for(Int_t iZ=0; iZ<48; iZ++) {
71 fHmpi0[iMod][iZ][iX]=0;
76 fMaskCellColumns = new Int_t[fNMaskCellColumns];
77 fMaskCellColumns[0] = 6 ; fMaskCellColumns[1] = 7 ; fMaskCellColumns[2] = 8 ;
78 fMaskCellColumns[3] = 35; fMaskCellColumns[4] = 36; fMaskCellColumns[5] = 37;
79 fMaskCellColumns[6] = 12+AliEMCALGeoParams::fgkEMCALCols; fMaskCellColumns[7] = 13+AliEMCALGeoParams::fgkEMCALCols;
80 fMaskCellColumns[8] = 40+AliEMCALGeoParams::fgkEMCALCols; fMaskCellColumns[9] = 41+AliEMCALGeoParams::fgkEMCALCols;
81 fMaskCellColumns[10]= 42+AliEMCALGeoParams::fgkEMCALCols;
83 for(Int_t iSMPair = 0; iSMPair < AliEMCALGeoParams::fgkEMCALModules/2; iSMPair++) {
84 fHmggPairSameSectorSM[iSMPair] = 0;
85 fHmggPairSameSectorSMMaskFrame[iSMPair] = 0;
86 fhClusterPairDiffTimeSameSector[iSMPair]= 0;
88 for(Int_t iSMPair = 0; iSMPair < AliEMCALGeoParams::fgkEMCALModules-2; iSMPair++){
89 fHmggPairSameSideSM[iSMPair] = 0;
90 fHmggPairSameSideSMMaskFrame[iSMPair] = 0;
91 fhClusterPairDiffTimeSameSide[iSMPair] = 0;
94 for(Int_t iSM = 0; iSM < AliEMCALGeoParams::fgkEMCALModules; iSM++) {
96 fHmggSMMaskFrame[iSM] = 0;
97 fHOpeningAngleSM[iSM] = 0;
98 fHOpeningAnglePairSM[iSM] = 0;
99 fHAsymmetrySM[iSM] = 0;
100 fHAsymmetryPairSM[iSM] = 0;
101 fHIncidentAngleSM[iSM] = 0;
102 fHIncidentAnglePairSM[iSM] = 0;
103 fhTowerDecayPhotonHit[iSM] = 0;
104 fhTowerDecayPhotonEnergy[iSM] = 0;
105 fhTowerDecayPhotonAsymmetry[iSM] = 0;
106 fhTowerDecayPhotonHitMaskFrame[iSM]= 0;
108 fhClusterTimeSM[iSM] = 0;
109 fhClusterPairDiffTimeSameSM[iSM] = 0;
112 DefineOutput(1, TList::Class());
113 DefineOutput(2, TList::Class()); // will contain cuts or local params
117 //__________________________________________________
118 AliAnalysisTaskEMCALPi0CalibSelection::~AliAnalysisTaskEMCALPi0CalibSelection()
122 if(fOutputContainer){
123 fOutputContainer->Delete() ;
124 delete fOutputContainer ;
127 if(fEMCALGeo) delete fEMCALGeo ;
128 if(fRecoUtils) delete fRecoUtils ;
129 if(fNMaskCellColumns) delete [] fMaskCellColumns;
133 //_____________________________________________________
134 void AliAnalysisTaskEMCALPi0CalibSelection::LocalInit()
136 // Local Initialization
138 // Create cuts/param objects and publish to slot
139 const Int_t buffersize = 255;
140 char onePar[buffersize] ;
143 snprintf(onePar,buffersize, "Custer cuts: %2.2f < E < %2.2f GeV; min number of cells %d; Assymetry cut %1.2f", fEmin,fEmax, fMinNCells, fAsyCut) ;
144 fCuts->Add(new TObjString(onePar));
145 snprintf(onePar,buffersize, "Group %d cells;", fGroupNCells) ;
146 fCuts->Add(new TObjString(onePar));
147 snprintf(onePar,buffersize, "Cluster maximal cell away from border at least %d cells;", fRecoUtils->GetNumberOfCellsFromEMCALBorder()) ;
148 fCuts->Add(new TObjString(onePar));
149 snprintf(onePar,buffersize, "Histograms: bins %d; energy range: %2.2f < E < %2.2f GeV;",fNbins,fMinBin,fMaxBin) ;
150 fCuts->Add(new TObjString(onePar));
151 snprintf(onePar,buffersize, "Switchs: Remove Bad Channels? %d; Use filtered input? %d; Correct Clusters? %d, Mass per channel same SM clusters? %d ",
152 fRecoUtils->IsBadChannelsRemovalSwitchedOn(),fFilteredInput,fCorrectClusters, fSameSM) ;
153 fCuts->Add(new TObjString(onePar));
154 snprintf(onePar,buffersize, "EMCAL Geometry name: < %s >, Load Matrices? %d",fEMCALGeoName.Data(),fLoadMatrices) ;
155 fCuts->Add(new TObjString(onePar));
157 fCuts ->SetOwner(kTRUE);
164 //__________________________________________________
165 void AliAnalysisTaskEMCALPi0CalibSelection::UserCreateOutputObjects()
167 //Create output container, init geometry
169 fEMCALGeo = AliEMCALGeometry::GetInstance(fEMCALGeoName) ;
170 Int_t nSM = (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules();
172 fOutputContainer = new TList();
173 const Int_t buffersize = 255;
174 char hname[buffersize], htitl[buffersize];
176 for(Int_t iMod=0; iMod < nSM; iMod++) {
177 for(Int_t iRow=0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
178 for(Int_t iCol=0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++) {
179 snprintf(hname,buffersize, "%d_%d_%d",iMod,iCol,iRow);
180 snprintf(htitl,buffersize, "Two-gamma inv. mass for super mod %d, cell(col,row)=(%d,%d)",iMod,iCol,iRow);
181 fHmpi0[iMod][iCol][iRow] = new TH1F(hname,htitl,fNbins,fMinBin,fMaxBin);
182 fOutputContainer->Add(fHmpi0[iMod][iCol][iRow]);
187 fHmgg = new TH2F("hmgg","2-cluster invariant mass",fNbins,fMinBin,fMaxBin,100,0,10);
188 fHmgg->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
189 fHmgg->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
190 fOutputContainer->Add(fHmgg);
192 fHmggDifferentSM = new TH2F("hmggDifferentSM","2-cluster invariant mass, different SM",fNbins,fMinBin,fMaxBin,100,0,10);
193 fHmggDifferentSM->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
194 fHmggDifferentSM->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
195 fOutputContainer->Add(fHmggDifferentSM);
197 fHOpeningAngle = new TH2F("hopang","2-cluster opening angle",100,0.,50.,100,0,10);
198 fHOpeningAngle->SetXTitle("#alpha_{#gamma #gamma}");
199 fHOpeningAngle->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
200 fOutputContainer->Add(fHOpeningAngle);
202 fHOpeningAngleDifferentSM = new TH2F("hopangDifferentSM","2-cluster opening angle, different SM",100,0,50.,100,0,10);
203 fHOpeningAngleDifferentSM->SetXTitle("#alpha_{#gamma #gamma}");
204 fHOpeningAngleDifferentSM->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
205 fOutputContainer->Add(fHOpeningAngleDifferentSM);
207 fHIncidentAngle = new TH2F("hinang","#gamma incident angle in SM",100,0.,20.,100,0,10);
208 fHIncidentAngle->SetXTitle("#alpha_{#gamma SM center}");
209 fHIncidentAngle->SetYTitle("p_{T #gamma} (GeV/c)");
210 fOutputContainer->Add(fHIncidentAngle);
212 fHIncidentAngleDifferentSM = new TH2F("hinangDifferentSM","#gamma incident angle in SM, different SM pair",100,0,20.,100,0,10);
213 fHIncidentAngleDifferentSM->SetXTitle("#alpha_{#gamma - SM center}");
214 fHIncidentAngleDifferentSM->SetYTitle("p_{T #gamma} (GeV/c)");
215 fOutputContainer->Add(fHIncidentAngleDifferentSM);
217 fHAsymmetry = new TH2F("hasym","2-cluster opening angle",100,0.,1.,100,0,10);
218 fHAsymmetry->SetXTitle("a");
219 fHAsymmetry->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
220 fOutputContainer->Add(fHAsymmetry);
222 fHAsymmetryDifferentSM = new TH2F("hasymDifferentSM","2-cluster opening angle, different SM",100,0,1.,100,0,10);
223 fHAsymmetryDifferentSM->SetXTitle("a");
224 fHAsymmetryDifferentSM->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
225 fOutputContainer->Add(fHAsymmetryDifferentSM);
228 //TString pairname[] = {"A side (0-2)", "C side (1-3)","Row 0 (0-1)", "Row 1 (2-3)"};
230 fHmggMaskFrame = new TH2F("hmggMaskFrame","2-cluster invariant mass, frame masked",fNbins,fMinBin,fMaxBin,100,0,10);
231 fHmggMaskFrame->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
232 fHmggMaskFrame->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
233 fOutputContainer->Add(fHmggMaskFrame);
235 fHmggDifferentSMMaskFrame = new TH2F("hmggDifferentSMMaskFrame","2-cluster invariant mass, different SM, frame masked",
236 fNbins,fMinBin,fMaxBin,100,0,10);
237 fHmggDifferentSMMaskFrame->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
238 fHmggDifferentSMMaskFrame->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
239 fOutputContainer->Add(fHmggDifferentSMMaskFrame);
242 for(Int_t iSM = 0; iSM < nSM; iSM++) {
244 snprintf(hname, buffersize, "hmgg_SM%d",iSM);
245 snprintf(htitl, buffersize, "Two-gamma inv. mass for super mod %d",iSM);
246 fHmggSM[iSM] = new TH2F(hname,htitl,fNbins,fMinBin,fMaxBin,100,0,10);
247 fHmggSM[iSM]->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
248 fHmggSM[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
249 fOutputContainer->Add(fHmggSM[iSM]);
251 snprintf(hname, buffersize, "hmgg_SM%d_MaskFrame",iSM);
252 snprintf(htitl, buffersize, "Two-gamma inv. mass for super mod %d",iSM);
253 fHmggSMMaskFrame[iSM] = new TH2F(hname,htitl,fNbins,fMinBin,fMaxBin,100,0,10);
254 fHmggSMMaskFrame[iSM]->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
255 fHmggSMMaskFrame[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
256 fOutputContainer->Add(fHmggSMMaskFrame[iSM]);
260 snprintf(hname,buffersize, "hmgg_PairSameSectorSM%d",iSM);
261 snprintf(htitl,buffersize, "Two-gamma inv. mass for SM pair Sector: %d",iSM);
262 fHmggPairSameSectorSM[iSM] = new TH2F(hname,htitl,fNbins,fMinBin,fMaxBin,100,0,10);
263 fHmggPairSameSectorSM[iSM]->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
264 fHmggPairSameSectorSM[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
265 fOutputContainer->Add(fHmggPairSameSectorSM[iSM]);
267 snprintf(hname,buffersize, "hmgg_PairSameSectorSM%d_MaskFrame",iSM);
268 snprintf(htitl,buffersize, "Two-gamma inv. mass for SM pair Sector: %d",iSM);
269 fHmggPairSameSectorSMMaskFrame[iSM] = new TH2F(hname,htitl,fNbins,fMinBin,fMaxBin,100,0,10);
270 fHmggPairSameSectorSMMaskFrame[iSM]->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
271 fHmggPairSameSectorSMMaskFrame[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
272 fOutputContainer->Add(fHmggPairSameSectorSMMaskFrame[iSM]);
274 fhClusterPairDiffTimeSameSector[iSM] = new TH2F(Form("hClusterPairDiffTimeSameSector%d",iSM),
275 Form("cluster pair time difference vs E, Sector %d",iSM),
276 100,0,10, 200,-100,100);
277 fhClusterPairDiffTimeSameSector[iSM]->SetXTitle("E_{pair} (GeV)");
278 fhClusterPairDiffTimeSameSector[iSM]->SetYTitle("#Delta t (ns)");
279 fOutputContainer->Add(fhClusterPairDiffTimeSameSector[iSM]);
285 snprintf(hname,buffersize, "hmgg_PairSameSideSM%d",iSM);
286 snprintf(htitl,buffersize, "Two-gamma inv. mass for SM pair Sector: %d",iSM);
287 fHmggPairSameSideSM[iSM] = new TH2F(hname,htitl,fNbins,fMinBin,fMaxBin,100,0,10);
288 fHmggPairSameSideSM[iSM]->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
289 fHmggPairSameSideSM[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
290 fOutputContainer->Add(fHmggPairSameSideSM[iSM]);
292 snprintf(hname,buffersize, "hmgg_PairSameSideSM%d_MaskFrame",iSM);
293 snprintf(htitl,buffersize, "Two-gamma inv. mass for SM pair Sector: %d",iSM);
294 fHmggPairSameSideSMMaskFrame[iSM] = new TH2F(hname,htitl,fNbins,fMinBin,fMaxBin,100,0,10);
295 fHmggPairSameSideSMMaskFrame[iSM]->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
296 fHmggPairSameSideSMMaskFrame[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
297 fOutputContainer->Add(fHmggPairSameSideSMMaskFrame[iSM]);
299 fhClusterPairDiffTimeSameSide[iSM] = new TH2F(Form("hClusterPairDiffTimeSameSide%d",iSM),
300 Form("cluster pair time difference vs E, Side %d",iSM),
301 100,0,10, 200,-100,100);
302 fhClusterPairDiffTimeSameSide[iSM]->SetXTitle("E_{pair} (GeV)");
303 fhClusterPairDiffTimeSameSide[iSM]->SetYTitle("#Delta t (ns)");
304 fOutputContainer->Add(fhClusterPairDiffTimeSameSide[iSM]);
308 snprintf(hname, buffersize, "hopang_SM%d",iSM);
309 snprintf(htitl, buffersize, "Opening angle for super mod %d",iSM);
310 fHOpeningAngleSM[iSM] = new TH2F(hname,htitl,100,0.,50.,100,0,10);
311 fHOpeningAngleSM[iSM]->SetXTitle("#alpha_{#gamma #gamma} (deg)");
312 fHOpeningAngleSM[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
313 fOutputContainer->Add(fHOpeningAngleSM[iSM]);
315 snprintf(hname,buffersize, "hopang_PairSM%d",iSM);
316 snprintf(htitl,buffersize, "Opening angle for SM pair: %d",iSM);
317 fHOpeningAnglePairSM[iSM] = new TH2F(hname,htitl,100,0.,50.,100,0,10);
318 fHOpeningAnglePairSM[iSM]->SetXTitle("#alpha_{#gamma #gamma} (deg)");
319 fHOpeningAnglePairSM[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
320 fOutputContainer->Add(fHOpeningAnglePairSM[iSM]);
322 snprintf(hname, buffersize, "hinang_SM%d",iSM);
323 snprintf(htitl, buffersize, "Incident angle for super mod %d",iSM);
324 fHIncidentAngleSM[iSM] = new TH2F(hname,htitl,100,0.,20.,100,0,10);
325 fHIncidentAngleSM[iSM]->SetXTitle("#alpha_{#gamma - SM center} (deg)");
326 fHIncidentAngleSM[iSM]->SetYTitle("p_{T #gamma} (GeV/c)");
327 fOutputContainer->Add(fHIncidentAngleSM[iSM]);
329 snprintf(hname,buffersize, "hinang_PairSM%d",iSM);
330 snprintf(htitl,buffersize, "Incident angle for SM pair: %d",iSM);
331 fHIncidentAnglePairSM[iSM] = new TH2F(hname,htitl,100,0.,20.,100,0,10);
332 fHIncidentAnglePairSM[iSM]->SetXTitle("#alpha_{#gamma - SM center} (deg)");
333 fHIncidentAnglePairSM[iSM]->SetYTitle("p_{T #gamma} (GeV/c)");
334 fOutputContainer->Add(fHIncidentAnglePairSM[iSM]);
336 snprintf(hname, buffersize, "hasym_SM%d",iSM);
337 snprintf(htitl, buffersize, "Asymmetry for super mod %d",iSM);
338 fHAsymmetrySM[iSM] = new TH2F(hname,htitl,100,0.,1.,100,0,10);
339 fHAsymmetrySM[iSM]->SetXTitle("a");
340 fHAsymmetrySM[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
341 fOutputContainer->Add(fHAsymmetrySM[iSM]);
343 snprintf(hname,buffersize, "hasym_PairSM%d",iSM);
344 snprintf(htitl,buffersize, "Asymmetry for SM pair: %d",iSM);
345 fHAsymmetryPairSM[iSM] = new TH2F(hname,htitl,100,0.,1.,100,0,10);
346 fHAsymmetryPairSM[iSM]->SetXTitle("a");
347 fHAsymmetryPairSM[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
348 fOutputContainer->Add(fHAsymmetryPairSM[iSM]);
353 fhTowerDecayPhotonHit[iSM] = new TH2F (Form("hTowerDecPhotonHit_Mod%d",iSM),
354 Form("Entries in grid of cells in Module %d",iSM),
355 colmax+2,-1.5,colmax+0.5, rowmax+2,-1.5,rowmax+0.5);
356 fhTowerDecayPhotonHit[iSM]->SetYTitle("row (phi direction)");
357 fhTowerDecayPhotonHit[iSM]->SetXTitle("column (eta direction)");
358 fOutputContainer->Add(fhTowerDecayPhotonHit[iSM]);
360 fhTowerDecayPhotonEnergy[iSM] = new TH2F (Form("hTowerDecPhotonEnergy_Mod%d",iSM),
361 Form("Accumulated energy in grid of cells in Module %d",iSM),
362 colmax+2,-1.5,colmax+0.5, rowmax+2,-1.5,rowmax+0.5);
363 fhTowerDecayPhotonEnergy[iSM]->SetYTitle("row (phi direction)");
364 fhTowerDecayPhotonEnergy[iSM]->SetXTitle("column (eta direction)");
365 fOutputContainer->Add(fhTowerDecayPhotonEnergy[iSM]);
367 fhTowerDecayPhotonAsymmetry[iSM] = new TH2F (Form("hTowerDecPhotonAsymmetry_Mod%d",iSM),
368 Form("Accumulated asymmetry in grid of cells in Module %d",iSM),
369 colmax+2,-1.5,colmax+0.5, rowmax+2,-1.5,rowmax+0.5);
370 fhTowerDecayPhotonAsymmetry[iSM]->SetYTitle("row (phi direction)");
371 fhTowerDecayPhotonAsymmetry[iSM]->SetXTitle("column (eta direction)");
372 fOutputContainer->Add(fhTowerDecayPhotonAsymmetry[iSM]);
374 fhTowerDecayPhotonHitMaskFrame[iSM] = new TH2F (Form("hTowerDecPhotonHit_Mod%d_MaskFrame",iSM),Form("Entries in grid of cells in Module %d",iSM),
375 colmax+2,-1.5,colmax+0.5, rowmax+2,-1.5,rowmax+0.5);
376 fhTowerDecayPhotonHitMaskFrame[iSM]->SetYTitle("row (phi direction)");
377 fhTowerDecayPhotonHitMaskFrame[iSM]->SetXTitle("column (eta direction)");
378 fOutputContainer->Add(fhTowerDecayPhotonHitMaskFrame[iSM]);
380 fhClusterTimeSM[iSM] = new TH2F(Form("hClusterTime_SM%d",iSM),"cluster time vs E",100,0,10, 100,0,1000);
381 fhClusterTimeSM[iSM]->SetXTitle("E (GeV)");
382 fhClusterTimeSM[iSM]->SetYTitle("t (ns)");
383 fOutputContainer->Add(fhClusterTimeSM[iSM]);
385 fhClusterPairDiffTimeSameSM[iSM] = new TH2F(Form("hClusterPairDiffTimeSameSM%d",iSM),
386 Form("cluster pair time difference vs E, SM %d",iSM),
387 100,0,10, 200,-100,100);
388 fhClusterPairDiffTimeSameSM[iSM]->SetXTitle("E (GeV)");
389 fhClusterPairDiffTimeSameSM[iSM]->SetYTitle("#Delta t (ns)");
390 fOutputContainer->Add(fhClusterPairDiffTimeSameSM[iSM]);
394 fhClusterTime = new TH2F("hClusterTime","cluster time vs E",100,0,10, 100,0,1000);
395 fhClusterTime->SetXTitle("E (GeV)");
396 fhClusterTime->SetYTitle("t (ns)");
397 fOutputContainer->Add(fhClusterTime);
399 fhClusterPairDiffTime = new TH2F("hClusterPairDiffTime","cluster pair time difference vs E",100,0,10, 200,-100,100);
400 fhClusterPairDiffTime->SetXTitle("E_{pair} (GeV)");
401 fhClusterPairDiffTime->SetYTitle("#Delta t (ns)");
402 fOutputContainer->Add(fhClusterPairDiffTime);
405 fhNEvents = new TH1I("hNEvents", "Number of analyzed events" , 1 , 0 , 1 ) ;
406 fOutputContainer->Add(fhNEvents);
408 fOutputContainer->SetOwner(kTRUE);
410 // fCalibData = new AliEMCALCalibData();
412 PostData(1,fOutputContainer);
416 //__________________________________________________
417 Bool_t AliAnalysisTaskEMCALPi0CalibSelection::MaskFrameCluster(const Int_t iSM, const Int_t ieta) const {
418 //Check if cell is in one of the regions where we have significant amount of material in front of EMCAL
421 if(iSM%2) icol+=48; // Impair SM, shift index [0-47] to [48-96]
423 if (fNMaskCellColumns && fMaskCellColumns) {
424 for (Int_t imask = 0; imask < fNMaskCellColumns; imask++) {
425 if(icol==fMaskCellColumns[imask]) return kTRUE;
433 //__________________________________________________
434 void AliAnalysisTaskEMCALPi0CalibSelection::UserExec(Option_t* /* option */)
436 //Analysis per event.
438 if(fRecoUtils->GetParticleType()!=AliEMCALRecoUtils::kPhoton){
439 printf("Wrong particle type for cluster position recalculation! = %d\n", fRecoUtils->GetParticleType());
443 if(!(((AliESDEvent*)InputEvent())->GetFiredTriggerClasses()).Contains(fTriggerName)) {
444 //printf("Reject Event %d, FiredClass %s\n",(Int_t)Entry(),(((AliESDEvent*)InputEvent())->GetFiredTriggerClasses()).Data());
448 fhNEvents->Fill(0); //Event analyzed
450 //Get the input event
451 AliVEvent* event = 0;
452 if(fFilteredInput) event = AODEvent();
453 else event = InputEvent();
456 printf("Input event not available!\n");
461 printf("AliAnalysisTaskEMCALPi0CalibSelection <<< %s: Event %d >>>\n",event->GetName(), (Int_t)Entry());
464 //Get the primary vertex
466 event->GetPrimaryVertex()->GetXYZ(v) ;
468 if(DebugLevel() > 1) printf("AliAnalysisTaskEMCALPi0CalibSelection Vertex: (%.3f,%.3f,%.3f)\n",v[0],v[1],v[2]);
470 //Int_t runNum = aod->GetRunNumber();
471 //if(DebugLevel() > 1) printf("Run number: %d\n",runNum);
473 Int_t nSM = (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules();
474 //Get the matrix with geometry information
475 if(fhNEvents->GetEntries()==1){
477 printf("AliAnalysisTaskEMCALPi0CalibSelection::UserExec() - Load user defined geometry matrices\n");
478 for(Int_t mod=0; mod < nSM ; mod++){
481 fMatrix[mod]->Print();
482 fEMCALGeo->SetMisalMatrix(fMatrix[mod],mod) ;
486 else if(!gGeoManager){
487 printf("AliAnalysisTaskEMCALPi0CalibSelection::UserExec() - Get geo matrices from data\n");
488 //Still not implemented in AOD, just a workaround to be able to work at least with ESDs
489 if(!strcmp(event->GetName(),"AliAODEvent")) {
491 printf("AliAnalysisTaskEMCALPi0CalibSelection Use ideal geometry, values geometry matrix not kept in AODs.\n");
494 if(DebugLevel() > 1) printf("AliAnalysisTaskEMCALPi0CalibSelection Load Misaligned matrices. \n");
495 AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event) ;
497 printf("AliAnalysisTaskEMCALPi0CalibSelection::UserExec() - This event does not contain ESDs?");
500 for(Int_t mod=0; mod < nSM; mod++){
502 esd->GetEMCALMatrix(mod)->Print();
503 if(esd->GetEMCALMatrix(mod)) fEMCALGeo->SetMisalMatrix(esd->GetEMCALMatrix(mod),mod) ;
506 }//Load matrices from Data
509 if(DebugLevel() > 1) printf("AliAnalysisTaskEMCALPi0CalibSelection Will use fLogWeight %.3f .\n",fLogWeight);
518 Bool_t shared = kFALSE;
524 //Get the list of clusters
525 TRefArray * caloClustersArr = new TRefArray();
526 event->GetEMCALClusters(caloClustersArr);
527 const Int_t kNumberOfEMCALClusters = caloClustersArr->GetEntries() ;
528 if(DebugLevel() > 1) printf("AliAnalysisTaskEMCALPi0CalibSelection - N CaloClusters: %d \n", kNumberOfEMCALClusters);
531 AliVCaloCells *emCells = event->GetEMCALCells();
533 // loop over EMCAL clusters
534 //----------------------------------------------------------
535 // First recalibrate and recalculate energy and position
536 Float_t pos[]={0,0,0};
537 if(fCorrectClusters){
538 for(Int_t iClu=0; iClu<kNumberOfEMCALClusters; iClu++) {
539 AliVCluster *c1 = (AliVCluster *) caloClustersArr->At(iClu);
541 Float_t e1i = c1->E(); // cluster energy before correction
542 if (e1i < fEmin) continue;
543 else if (e1i > fEmax) continue;
544 else if (c1->GetNCells() < fMinNCells) continue;
546 if(fRecoUtils->ClusterContainsBadChannel(fEMCALGeo, c1->GetCellsAbsId(), c1->GetNCells())) continue;
550 printf("Std : i %d, E %f, dispersion %f, m02 %f, m20 %f\n",c1->GetID(),c1->E(),c1->GetDispersion(),c1->GetM02(),c1->GetM20());
551 c1->GetPosition(pos);
552 printf("Std : i %d, x %f, y %f, z %f\n",c1->GetID(), pos[0], pos[1], pos[2]);
555 //Correct cluster energy and position if requested, and not corrected previously, by default Off
556 if(fRecoUtils->IsRecalibrationOn()) {
557 fRecoUtils->RecalibrateClusterEnergy(fEMCALGeo, c1, emCells);
558 fRecoUtils->RecalculateClusterShowerShapeParameters(fEMCALGeo, emCells,c1);
559 fRecoUtils->RecalculateClusterPID(c1);
562 printf("Energy: after recalibration %f; ",c1->E());
564 // Recalculate cluster position
565 fRecoUtils->RecalculateClusterPosition(fEMCALGeo, emCells,c1);
567 // Correct Non-Linearity
569 c1->SetE(fRecoUtils->CorrectClusterEnergyLinearity(c1));
571 printf("after linearity correction %f\n",c1->E());
575 printf("Cor : i %d, E %f, dispersion %f, m02 %f, m20 %f\n",c1->GetID(),c1->E(),c1->GetDispersion(),c1->GetM02(),c1->GetM20());
576 c1->GetPosition(pos);
577 printf("Cor : i %d, x %f, y %f, z %f\n",c1->GetID(), pos[0], pos[1], pos[2]);
582 //----------------------------------------------------------
583 //Now the invariant mass analysis with the corrected clusters
584 for(Int_t iClu=0; iClu<kNumberOfEMCALClusters-1; iClu++) {
586 AliVCluster *c1 = (AliVCluster *) caloClustersArr->At(iClu);
587 if(fRecoUtils->ClusterContainsBadChannel(fEMCALGeo, c1->GetCellsAbsId(), c1->GetNCells())) continue;
589 Float_t e1i = c1->E(); // cluster energy before correction
590 if (e1i < fEmin) continue;
591 else if (e1i > fEmax) continue;
592 else if (c1->GetNCells() < fMinNCells) continue;
596 printf("IMA : i %d, E %f, dispersion %f, m02 %f, m20 %f\n",c1->GetID(),e1i,c1->GetDispersion(),c1->GetM02(),c1->GetM20());
597 c1->GetPosition(pos);
598 printf("IMA : i %d, x %f, y %f, z %f\n",c1->GetID(), pos[0], pos[1], pos[2]);
601 fRecoUtils->GetMaxEnergyCell(fEMCALGeo, emCells,c1,absId1,iSupMod1,ieta1,iphi1,shared);
602 c1->GetMomentum(p1,v);
604 //Check if cluster is in fidutial region, not too close to borders
605 Bool_t in1 = fRecoUtils->CheckCellFiducialRegion(fEMCALGeo, c1, emCells);
606 // Clusters not facing frame structures
607 Bool_t mask1 = MaskFrameCluster(iSupMod1, ieta1);
608 //if(mask1) printf("Reject eta %d SM %d\n",ieta1, iSupMod1);
610 Double_t time1 = c1->GetTOF()*1.e9;
611 fhClusterTime ->Fill(c1->E(),time1);
612 fhClusterTimeSM[iSupMod1]->Fill(c1->E(),time1);
614 // Combine cluster with other clusters and get the invariant mass
615 for (Int_t jClu=iClu+1; jClu<kNumberOfEMCALClusters; jClu++) {
616 AliAODCaloCluster *c2 = (AliAODCaloCluster *) caloClustersArr->At(jClu);
618 if(fRecoUtils->ClusterContainsBadChannel(fEMCALGeo, c2->GetCellsAbsId(), c2->GetNCells())) continue;
620 Float_t e2i = c2->E();
621 if (e2i < fEmin) continue;
622 else if (e2i > fEmax) continue;
623 else if (c2->GetNCells() < fMinNCells) continue;
626 fRecoUtils->GetMaxEnergyCell(fEMCALGeo, emCells,c2,absId2,iSupMod2,ieta2,iphi2,shared);
627 c2->GetMomentum(p2,v);
630 Float_t invmass = p12.M()*1000;
631 //printf("*** mass %f\n",invmass);
634 Float_t asym = TMath::Abs(p1.E()-p2.E())/(p1.E()+p2.E());
635 //printf("asymmetry %f\n",asym);
636 if(asym > fAsyCut) continue;
639 Double_t time2 = c2->GetTOF()*1.e9;
640 fhClusterPairDiffTime->Fill(p12.E(),time1-time2);
641 if(TMath::Abs(time1-time2) > fDTimeCut) continue;
643 if(invmass < fMaxBin && invmass > fMinBin ){
645 //Check if cluster is in fidutial region, not too close to borders
646 Bool_t in2 = fRecoUtils->CheckCellFiducialRegion(fEMCALGeo, c2, emCells);
648 // Clusters not facing frame structures
649 Bool_t mask2 = MaskFrameCluster(iSupMod2, ieta2);
650 //if(mask2) printf("Reject eta %d SM %d\n",ieta2, iSupMod2);
654 fHmgg->Fill(invmass,p12.Pt());
656 if(iSupMod1==iSupMod2) {
657 fHmggSM[iSupMod1]->Fill(invmass,p12.Pt());
658 fhClusterPairDiffTimeSameSM[iSupMod1]->Fill(p12.E(),time1-time2);
661 fHmggDifferentSM ->Fill(invmass,p12.Pt());
665 for(Int_t i = 0; i < nSM/2; i++){
667 if((iSupMod1==j && iSupMod2==j+1) || (iSupMod1==j+1 && iSupMod2==j)) {
668 fHmggPairSameSectorSM[i]->Fill(invmass,p12.Pt());
669 fhClusterPairDiffTimeSameSector[i]->Fill(p12.E(),time1-time2);
674 for(Int_t i = 0; i < nSM-2; i++){
675 if((iSupMod1==i && iSupMod2==i+2) || (iSupMod1==i+2 && iSupMod2==i)) {
676 fHmggPairSameSideSM[i]->Fill(invmass,p12.Pt());
677 fhClusterPairDiffTimeSameSide[i]->Fill(p12.E(),time1-time2);
682 if(!mask1 && !mask2){
684 fHmggMaskFrame->Fill(invmass,p12.Pt());
686 if(iSupMod1==iSupMod2) fHmggSMMaskFrame[iSupMod1]->Fill(invmass,p12.Pt());
687 else fHmggDifferentSMMaskFrame ->Fill(invmass,p12.Pt());
691 for(Int_t i = 0; i < nSM/2; i++){
693 if((iSupMod1==j && iSupMod2==j+1) || (iSupMod1==j+1 && iSupMod2==j)) fHmggPairSameSectorSMMaskFrame[i]->Fill(invmass,p12.Pt());
697 for(Int_t i = 0; i < nSM-2; i++){
698 if((iSupMod1==i && iSupMod2==i+2) || (iSupMod1==i+2 && iSupMod2==i)) fHmggPairSameSideSMMaskFrame[i]->Fill(invmass,p12.Pt());
701 }// Pair not facing frame
704 if(invmass > 100. && invmass < 160.){//restrict to clusters really close to pi0 peak
706 //Opening angle of 2 photons
707 Float_t opangle = p1.Angle(p2.Vect())*TMath::RadToDeg();
708 //printf("*******>>>>>>>> In PEAK pt %f, angle %f \n",p12.Pt(),opangle);
710 //Incident angle of each photon
711 Float_t inangle1 =0., inangle2=0.;
713 Float_t posSM1cen[3]={0.,0.,0.};
714 Float_t depth = fRecoUtils->GetDepth(p1.Energy(),fRecoUtils->GetParticleType(),iSupMod1);
715 fEMCALGeo->RecalculateTowerPosition(11.5, 23.5, iSupMod1, depth, fRecoUtils->GetMisalTransShiftArray(),fRecoUtils->GetMisalRotShiftArray(),posSM1cen);
716 Float_t posSM2cen[3]={0.,0.,0.};
717 depth = fRecoUtils->GetDepth(p2.Energy(),fRecoUtils->GetParticleType(),iSupMod2);
718 fEMCALGeo->RecalculateTowerPosition(11.5, 23.5, iSupMod2, depth, fRecoUtils->GetMisalTransShiftArray(),fRecoUtils->GetMisalRotShiftArray(),posSM2cen);
719 //printf("SM1 %d pos (%2.3f,%2.3f,%2.3f) \n",iSupMod1,posSM1cen[0],posSM1cen[1],posSM1cen[2]);
720 //printf("SM2 %d pos (%2.3f,%2.3f,%2.3f) \n",iSupMod2,posSM2cen[0],posSM2cen[1],posSM2cen[2]);
722 TVector3 vecSM1cen(posSM1cen[0]-v[0],posSM1cen[1]-v[1],posSM1cen[2]-v[2]);
723 TVector3 vecSM2cen(posSM2cen[0]-v[0],posSM2cen[1]-v[1],posSM2cen[2]-v[2]);
724 inangle1 = p1.Angle(vecSM1cen)*TMath::RadToDeg();
725 inangle2 = p2.Angle(vecSM2cen)*TMath::RadToDeg();
726 //printf("Incident angle: cluster 1 %2.3f; cluster 2 %2.3f\n",inangle1,inangle2);
729 fHOpeningAngle ->Fill(opangle,p12.Pt());
730 fHIncidentAngle->Fill(inangle1,p1.Pt());
731 fHIncidentAngle->Fill(inangle2,p2.Pt());
732 fHAsymmetry ->Fill(asym,p12.Pt());
734 if(iSupMod1==iSupMod2) {
735 fHOpeningAngleSM[iSupMod1] ->Fill(opangle,p12.Pt());
736 fHIncidentAngleSM[iSupMod1]->Fill(inangle1,p1.Pt());
737 fHIncidentAngleSM[iSupMod1]->Fill(inangle2,p2.Pt());
738 fHAsymmetrySM[iSupMod1] ->Fill(asym,p12.Pt());
741 fHOpeningAngleDifferentSM ->Fill(opangle,p12.Pt());
742 fHIncidentAngleDifferentSM ->Fill(inangle1,p1.Pt());
743 fHIncidentAngleDifferentSM ->Fill(inangle2,p2.Pt());
744 fHAsymmetryDifferentSM ->Fill(asym,p12.Pt());
747 if((iSupMod1==0 && iSupMod2==2) || (iSupMod1==2 && iSupMod2==0)) {
748 fHOpeningAnglePairSM[0] ->Fill(opangle,p12.Pt());
749 fHIncidentAnglePairSM[0]->Fill(inangle1,p1.Pt());
750 fHIncidentAnglePairSM[0]->Fill(inangle2,p2.Pt());
751 fHAsymmetryPairSM[0] ->Fill(asym,p12.Pt());
754 if((iSupMod1==1 && iSupMod2==3) || (iSupMod1==3 && iSupMod2==1)) {
755 fHOpeningAnglePairSM[1] ->Fill(opangle,p12.Pt());
756 fHIncidentAnglePairSM[1]->Fill(inangle1,p1.Pt());
757 fHIncidentAnglePairSM[1]->Fill(inangle2,p2.Pt());
758 fHAsymmetryPairSM[1] ->Fill(asym,p12.Pt());
762 if((iSupMod1==0 && iSupMod2==1) || (iSupMod1==1 && iSupMod2==0)) {
763 fHOpeningAnglePairSM[2] ->Fill(opangle,p12.Pt());
764 fHIncidentAnglePairSM[2]->Fill(inangle1,p1.Pt());
765 fHIncidentAnglePairSM[2]->Fill(inangle2,p2.Pt());
766 fHAsymmetryPairSM[2] ->Fill(asym,p12.Pt());
770 if((iSupMod1==2 && iSupMod2==3) || (iSupMod1==3 && iSupMod2==2)) {
771 fHOpeningAnglePairSM[3] ->Fill(opangle,p12.Pt());
772 fHIncidentAnglePairSM[3]->Fill(inangle1,p1.Pt());
773 fHIncidentAnglePairSM[3]->Fill(inangle2,p2.Pt());
774 fHAsymmetryPairSM[3] ->Fill(asym,p12.Pt());
777 }// pair in 100 < mass < 160
779 }//in acceptance cuts
781 //In case of filling only channels with second cluster in same SM
782 if(fSameSM && iSupMod1!=iSupMod2) continue;
784 if (fGroupNCells == 0){
785 fHmpi0[iSupMod1][ieta1][iphi1]->Fill(invmass);
786 fHmpi0[iSupMod2][ieta2][iphi2]->Fill(invmass);
788 if(invmass > 100. && invmass < 160.){//restrict to clusters really close to pi0 peak
789 fhTowerDecayPhotonHit [iSupMod1]->Fill(ieta1,iphi1);
790 fhTowerDecayPhotonEnergy [iSupMod1]->Fill(ieta1,iphi1,p1.E());
791 fhTowerDecayPhotonAsymmetry[iSupMod1]->Fill(ieta1,iphi1,asym);
793 fhTowerDecayPhotonHit [iSupMod2]->Fill(ieta2,iphi2);
794 fhTowerDecayPhotonEnergy [iSupMod2]->Fill(ieta2,iphi2,p2.E());
795 fhTowerDecayPhotonAsymmetry[iSupMod2]->Fill(ieta2,iphi2,asym);
797 if(!mask1)fhTowerDecayPhotonHitMaskFrame[iSupMod1]->Fill(ieta1,iphi1);
798 if(!mask2)fhTowerDecayPhotonHitMaskFrame[iSupMod2]->Fill(ieta2,iphi2);
800 }// pair in mass of pi0
803 //printf("Regroup N %d, eta1 %d, phi1 %d, eta2 %d, phi2 %d \n",fGroupNCells, ieta1, iphi1, ieta2, iphi2);
804 for (Int_t i = -fGroupNCells; i < fGroupNCells+1; i++) {
805 for (Int_t j = -fGroupNCells; j < fGroupNCells+1; j++) {
806 //printf("\t i %d, j %d\n",i,j);
807 if((ieta1+i >= 0) && (iphi1+j >= 0) && (ieta1+i < 48) && (iphi1+j < 24)){
808 //printf("\t \t eta1+i %d, phi1+j %d\n", ieta1+i, iphi1+j);
809 fHmpi0[iSupMod1][ieta1+i][iphi1+j]->Fill(invmass);
811 if((ieta2+i >= 0) && (iphi2+j >= 0) && (ieta2+i < 48) && (iphi2+j < 24)){
812 //printf("\t \t eta2+i %d, phi2+j %d\n", ieta2+i, iphi2+j);
813 fHmpi0[iSupMod2][ieta2+i][iphi2+j]->Fill(invmass);
819 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",
820 iSupMod1,iphi1,ieta1,iSupMod2,iphi2,ieta2,p12.M(),e1i,c1->E(),e2i,c2->E());
825 } // end of loop over EMCAL clusters
827 delete caloClustersArr;
829 PostData(1,fOutputContainer);
833 //_____________________________________________________
834 void AliAnalysisTaskEMCALPi0CalibSelection::PrintInfo(){
837 printf("Cluster cuts: %2.2f < E < %2.2f GeV; number of cells > %d; Assymetry < %1.2f, pair time diff < %2.2f\n",
838 fEmin,fEmax, fMinNCells, fAsyCut, fDTimeCut) ;
839 printf("Group %d cells\n", fGroupNCells) ;
840 printf("Cluster maximal cell away from border at least %d cells\n", fRecoUtils->GetNumberOfCellsFromEMCALBorder()) ;
841 printf("Histograms: bins %d; energy range: %2.2f < E < %2.2f GeV\n",fNbins,fMinBin,fMaxBin) ;
842 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",
843 fRecoUtils->IsBadChannelsRemovalSwitchedOn(),fFilteredInput,fCorrectClusters, fSameSM) ;
844 printf("EMCAL Geometry name: < %s >, Load Matrices %d\n",fEMCALGeoName.Data(), fLoadMatrices) ;
845 if(fLoadMatrices) {for(Int_t ism = 0; ism < AliEMCALGeoParams::fgkEMCALModules; ism++) fMatrix[ism]->Print() ; }