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),//fCalibData(0x0),
49 fEmin(0.5), fEmax(15.), fDTimeCut(20.), fAsyCut(1.),fMinNCells(2), fGroupNCells(0),
50 fLogWeight(4.5), fSameSM(kFALSE), fFilteredInput(kFALSE),
51 fCorrectClusters(kFALSE), fEMCALGeoName("EMCAL_COMPLETEV1"),
52 fRecoUtils(new AliEMCALRecoUtils),
53 fNbins(300), fMinBin(0.), fMaxBin(300.),fOutputContainer(0x0),
54 fHmgg(0x0), fHmggDifferentSM(0x0),
55 fHmggMaskFrame(0x0), fHmggDifferentSMMaskFrame(0x0),
56 fHOpeningAngle(0x0), fHOpeningAngleDifferentSM(0x0),
57 fHIncidentAngle(0x0), fHIncidentAngleDifferentSM(0x0),
58 fHAsymmetry(0x0), fHAsymmetryDifferentSM(0x0),
59 fhNEvents(0x0),fCuts(0x0),fLoadMatrices(0),
60 fNMaskCellColumns(11), fMaskCellColumns(0x0),
61 fhClusterTime(0x0), fhClusterPairDiffTime(0x0)
63 //Named constructor which should be used.
65 for(Int_t iMod=0; iMod < AliEMCALGeoParams::fgkEMCALModules; iMod++) {
66 for(Int_t iX=0; iX<24; iX++) {
67 for(Int_t iZ=0; iZ<48; iZ++) {
68 fHmpi0[iMod][iZ][iX]=0;
73 fMaskCellColumns = new Int_t[fNMaskCellColumns];
74 fMaskCellColumns[0] = 6 ; fMaskCellColumns[1] = 7 ; fMaskCellColumns[2] = 8 ;
75 fMaskCellColumns[3] = 35; fMaskCellColumns[4] = 36; fMaskCellColumns[5] = 37;
76 fMaskCellColumns[6] = 12+AliEMCALGeoParams::fgkEMCALCols; fMaskCellColumns[7] = 13+AliEMCALGeoParams::fgkEMCALCols;
77 fMaskCellColumns[8] = 40+AliEMCALGeoParams::fgkEMCALCols; fMaskCellColumns[9] = 41+AliEMCALGeoParams::fgkEMCALCols;
78 fMaskCellColumns[10]= 42+AliEMCALGeoParams::fgkEMCALCols;
80 for(Int_t iSMPair = 0; iSMPair < AliEMCALGeoParams::fgkEMCALModules/2; iSMPair++) {
81 fHmggPairSameSectorSM[iSMPair] = 0;
82 fHmggPairSameSectorSMMaskFrame[iSMPair] = 0;
83 fhClusterPairDiffTimeSameSector[iSMPair]= 0;
85 for(Int_t iSMPair = 0; iSMPair < AliEMCALGeoParams::fgkEMCALModules-2; iSMPair++){
86 fHmggPairSameSideSM[iSMPair] = 0;
87 fHmggPairSameSideSMMaskFrame[iSMPair] = 0;
88 fhClusterPairDiffTimeSameSide[iSMPair] = 0;
91 for(Int_t iSM = 0; iSM < AliEMCALGeoParams::fgkEMCALModules; iSM++) {
93 fHmggSMMaskFrame[iSM] = 0;
94 fHOpeningAngleSM[iSM] = 0;
95 fHOpeningAnglePairSM[iSM] = 0;
96 fHAsymmetrySM[iSM] = 0;
97 fHAsymmetryPairSM[iSM] = 0;
98 fHIncidentAngleSM[iSM] = 0;
99 fHIncidentAnglePairSM[iSM] = 0;
100 fhTowerDecayPhotonHit[iSM] = 0;
101 fhTowerDecayPhotonEnergy[iSM] = 0;
102 fhTowerDecayPhotonAsymmetry[iSM] = 0;
103 fhTowerDecayPhotonHitMaskFrame[iSM]= 0;
105 fhClusterTimeSM[iSM] = 0;
106 fhClusterPairDiffTimeSameSM[iSM] = 0;
109 DefineOutput(1, TList::Class());
110 DefineOutput(2, TList::Class()); // will contain cuts or local params
114 //__________________________________________________
115 AliAnalysisTaskEMCALPi0CalibSelection::~AliAnalysisTaskEMCALPi0CalibSelection()
119 if(fOutputContainer){
120 fOutputContainer->Delete() ;
121 delete fOutputContainer ;
124 if(fEMCALGeo) delete fEMCALGeo ;
125 if(fRecoUtils) delete fRecoUtils ;
126 if(fNMaskCellColumns) delete [] fMaskCellColumns;
130 //_____________________________________________________
131 void AliAnalysisTaskEMCALPi0CalibSelection::LocalInit()
133 // Local Initialization
135 // Create cuts/param objects and publish to slot
136 const Int_t buffersize = 255;
137 char onePar[buffersize] ;
140 snprintf(onePar,buffersize, "Custer cuts: %2.2f < E < %2.2f GeV; min number of cells %d; Assymetry cut %1.2f", fEmin,fEmax, fMinNCells, fAsyCut) ;
141 fCuts->Add(new TObjString(onePar));
142 snprintf(onePar,buffersize, "Group %d cells;", fGroupNCells) ;
143 fCuts->Add(new TObjString(onePar));
144 snprintf(onePar,buffersize, "Cluster maximal cell away from border at least %d cells;", fRecoUtils->GetNumberOfCellsFromEMCALBorder()) ;
145 fCuts->Add(new TObjString(onePar));
146 snprintf(onePar,buffersize, "Histograms: bins %d; energy range: %2.2f < E < %2.2f GeV;",fNbins,fMinBin,fMaxBin) ;
147 fCuts->Add(new TObjString(onePar));
148 snprintf(onePar,buffersize, "Switchs: Remove Bad Channels? %d; Use filtered input? %d; Correct Clusters? %d, Mass per channel same SM clusters? %d ",
149 fRecoUtils->IsBadChannelsRemovalSwitchedOn(),fFilteredInput,fCorrectClusters, fSameSM) ;
150 fCuts->Add(new TObjString(onePar));
151 snprintf(onePar,buffersize, "EMCAL Geometry name: < %s >, Load Matrices? %d",fEMCALGeoName.Data(),fLoadMatrices) ;
152 fCuts->Add(new TObjString(onePar));
154 fCuts ->SetOwner(kTRUE);
161 //__________________________________________________
162 void AliAnalysisTaskEMCALPi0CalibSelection::UserCreateOutputObjects()
164 //Create output container, init geometry
166 fEMCALGeo = AliEMCALGeometry::GetInstance(fEMCALGeoName) ;
167 Int_t nSM = (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules();
169 fOutputContainer = new TList();
170 const Int_t buffersize = 255;
171 char hname[buffersize], htitl[buffersize];
173 for(Int_t iMod=0; iMod < nSM; iMod++) {
174 for(Int_t iRow=0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
175 for(Int_t iCol=0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++) {
176 snprintf(hname,buffersize, "%d_%d_%d",iMod,iCol,iRow);
177 snprintf(htitl,buffersize, "Two-gamma inv. mass for super mod %d, cell(col,row)=(%d,%d)",iMod,iCol,iRow);
178 fHmpi0[iMod][iCol][iRow] = new TH1F(hname,htitl,fNbins,fMinBin,fMaxBin);
179 fOutputContainer->Add(fHmpi0[iMod][iCol][iRow]);
184 fHmgg = new TH2F("hmgg","2-cluster invariant mass",fNbins,fMinBin,fMaxBin,100,0,10);
185 fHmgg->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
186 fHmgg->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
187 fOutputContainer->Add(fHmgg);
189 fHmggDifferentSM = new TH2F("hmggDifferentSM","2-cluster invariant mass, different SM",fNbins,fMinBin,fMaxBin,100,0,10);
190 fHmggDifferentSM->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
191 fHmggDifferentSM->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
192 fOutputContainer->Add(fHmggDifferentSM);
194 fHOpeningAngle = new TH2F("hopang","2-cluster opening angle",100,0.,50.,100,0,10);
195 fHOpeningAngle->SetXTitle("#alpha_{#gamma #gamma}");
196 fHOpeningAngle->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
197 fOutputContainer->Add(fHOpeningAngle);
199 fHOpeningAngleDifferentSM = new TH2F("hopangDifferentSM","2-cluster opening angle, different SM",100,0,50.,100,0,10);
200 fHOpeningAngleDifferentSM->SetXTitle("#alpha_{#gamma #gamma}");
201 fHOpeningAngleDifferentSM->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
202 fOutputContainer->Add(fHOpeningAngleDifferentSM);
204 fHIncidentAngle = new TH2F("hinang","#gamma incident angle in SM",100,0.,20.,100,0,10);
205 fHIncidentAngle->SetXTitle("#alpha_{#gamma SM center}");
206 fHIncidentAngle->SetYTitle("p_{T #gamma} (GeV/c)");
207 fOutputContainer->Add(fHIncidentAngle);
209 fHIncidentAngleDifferentSM = new TH2F("hinangDifferentSM","#gamma incident angle in SM, different SM pair",100,0,20.,100,0,10);
210 fHIncidentAngleDifferentSM->SetXTitle("#alpha_{#gamma - SM center}");
211 fHIncidentAngleDifferentSM->SetYTitle("p_{T #gamma} (GeV/c)");
212 fOutputContainer->Add(fHIncidentAngleDifferentSM);
214 fHAsymmetry = new TH2F("hasym","2-cluster opening angle",100,0.,1.,100,0,10);
215 fHAsymmetry->SetXTitle("a");
216 fHAsymmetry->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
217 fOutputContainer->Add(fHAsymmetry);
219 fHAsymmetryDifferentSM = new TH2F("hasymDifferentSM","2-cluster opening angle, different SM",100,0,1.,100,0,10);
220 fHAsymmetryDifferentSM->SetXTitle("a");
221 fHAsymmetryDifferentSM->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
222 fOutputContainer->Add(fHAsymmetryDifferentSM);
225 //TString pairname[] = {"A side (0-2)", "C side (1-3)","Row 0 (0-1)", "Row 1 (2-3)"};
227 fHmggMaskFrame = new TH2F("hmggMaskFrame","2-cluster invariant mass, frame masked",fNbins,fMinBin,fMaxBin,100,0,10);
228 fHmggMaskFrame->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
229 fHmggMaskFrame->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
230 fOutputContainer->Add(fHmggMaskFrame);
232 fHmggDifferentSMMaskFrame = new TH2F("hmggDifferentSMMaskFrame","2-cluster invariant mass, different SM, frame masked",
233 fNbins,fMinBin,fMaxBin,100,0,10);
234 fHmggDifferentSMMaskFrame->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
235 fHmggDifferentSMMaskFrame->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
236 fOutputContainer->Add(fHmggDifferentSMMaskFrame);
239 for(Int_t iSM = 0; iSM < nSM; iSM++) {
241 snprintf(hname, buffersize, "hmgg_SM%d",iSM);
242 snprintf(htitl, buffersize, "Two-gamma inv. mass for super mod %d",iSM);
243 fHmggSM[iSM] = new TH2F(hname,htitl,fNbins,fMinBin,fMaxBin,100,0,10);
244 fHmggSM[iSM]->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
245 fHmggSM[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
246 fOutputContainer->Add(fHmggSM[iSM]);
248 snprintf(hname, buffersize, "hmgg_SM%d_MaskFrame",iSM);
249 snprintf(htitl, buffersize, "Two-gamma inv. mass for super mod %d",iSM);
250 fHmggSMMaskFrame[iSM] = new TH2F(hname,htitl,fNbins,fMinBin,fMaxBin,100,0,10);
251 fHmggSMMaskFrame[iSM]->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
252 fHmggSMMaskFrame[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
253 fOutputContainer->Add(fHmggSMMaskFrame[iSM]);
257 snprintf(hname,buffersize, "hmgg_PairSameSectorSM%d",iSM);
258 snprintf(htitl,buffersize, "Two-gamma inv. mass for SM pair Sector: %d",iSM);
259 fHmggPairSameSectorSM[iSM] = new TH2F(hname,htitl,fNbins,fMinBin,fMaxBin,100,0,10);
260 fHmggPairSameSectorSM[iSM]->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
261 fHmggPairSameSectorSM[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
262 fOutputContainer->Add(fHmggPairSameSectorSM[iSM]);
264 snprintf(hname,buffersize, "hmgg_PairSameSectorSM%d_MaskFrame",iSM);
265 snprintf(htitl,buffersize, "Two-gamma inv. mass for SM pair Sector: %d",iSM);
266 fHmggPairSameSectorSMMaskFrame[iSM] = new TH2F(hname,htitl,fNbins,fMinBin,fMaxBin,100,0,10);
267 fHmggPairSameSectorSMMaskFrame[iSM]->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
268 fHmggPairSameSectorSMMaskFrame[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
269 fOutputContainer->Add(fHmggPairSameSectorSMMaskFrame[iSM]);
271 fhClusterPairDiffTimeSameSector[iSM] = new TH2F(Form("hClusterPairDiffTimeSameSector%d",iSM),
272 Form("cluster pair time difference vs E, Sector %d",iSM),
273 100,0,10, 200,-100,100);
274 fhClusterPairDiffTimeSameSector[iSM]->SetXTitle("E_{pair} (GeV)");
275 fhClusterPairDiffTimeSameSector[iSM]->SetYTitle("#Delta t (ns)");
276 fOutputContainer->Add(fhClusterPairDiffTimeSameSector[iSM]);
282 snprintf(hname,buffersize, "hmgg_PairSameSideSM%d",iSM);
283 snprintf(htitl,buffersize, "Two-gamma inv. mass for SM pair Sector: %d",iSM);
284 fHmggPairSameSideSM[iSM] = new TH2F(hname,htitl,fNbins,fMinBin,fMaxBin,100,0,10);
285 fHmggPairSameSideSM[iSM]->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
286 fHmggPairSameSideSM[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
287 fOutputContainer->Add(fHmggPairSameSideSM[iSM]);
289 snprintf(hname,buffersize, "hmgg_PairSameSideSM%d_MaskFrame",iSM);
290 snprintf(htitl,buffersize, "Two-gamma inv. mass for SM pair Sector: %d",iSM);
291 fHmggPairSameSideSMMaskFrame[iSM] = new TH2F(hname,htitl,fNbins,fMinBin,fMaxBin,100,0,10);
292 fHmggPairSameSideSMMaskFrame[iSM]->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
293 fHmggPairSameSideSMMaskFrame[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
294 fOutputContainer->Add(fHmggPairSameSideSMMaskFrame[iSM]);
296 fhClusterPairDiffTimeSameSide[iSM] = new TH2F(Form("hClusterPairDiffTimeSameSide%d",iSM),
297 Form("cluster pair time difference vs E, Side %d",iSM),
298 100,0,10, 200,-100,100);
299 fhClusterPairDiffTimeSameSide[iSM]->SetXTitle("E_{pair} (GeV)");
300 fhClusterPairDiffTimeSameSide[iSM]->SetYTitle("#Delta t (ns)");
301 fOutputContainer->Add(fhClusterPairDiffTimeSameSide[iSM]);
305 snprintf(hname, buffersize, "hopang_SM%d",iSM);
306 snprintf(htitl, buffersize, "Opening angle for super mod %d",iSM);
307 fHOpeningAngleSM[iSM] = new TH2F(hname,htitl,100,0.,50.,100,0,10);
308 fHOpeningAngleSM[iSM]->SetXTitle("#alpha_{#gamma #gamma} (deg)");
309 fHOpeningAngleSM[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
310 fOutputContainer->Add(fHOpeningAngleSM[iSM]);
312 snprintf(hname,buffersize, "hopang_PairSM%d",iSM);
313 snprintf(htitl,buffersize, "Opening angle for SM pair: %d",iSM);
314 fHOpeningAnglePairSM[iSM] = new TH2F(hname,htitl,100,0.,50.,100,0,10);
315 fHOpeningAnglePairSM[iSM]->SetXTitle("#alpha_{#gamma #gamma} (deg)");
316 fHOpeningAnglePairSM[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
317 fOutputContainer->Add(fHOpeningAnglePairSM[iSM]);
319 snprintf(hname, buffersize, "hinang_SM%d",iSM);
320 snprintf(htitl, buffersize, "Incident angle for super mod %d",iSM);
321 fHIncidentAngleSM[iSM] = new TH2F(hname,htitl,100,0.,20.,100,0,10);
322 fHIncidentAngleSM[iSM]->SetXTitle("#alpha_{#gamma - SM center} (deg)");
323 fHIncidentAngleSM[iSM]->SetYTitle("p_{T #gamma} (GeV/c)");
324 fOutputContainer->Add(fHIncidentAngleSM[iSM]);
326 snprintf(hname,buffersize, "hinang_PairSM%d",iSM);
327 snprintf(htitl,buffersize, "Incident angle for SM pair: %d",iSM);
328 fHIncidentAnglePairSM[iSM] = new TH2F(hname,htitl,100,0.,20.,100,0,10);
329 fHIncidentAnglePairSM[iSM]->SetXTitle("#alpha_{#gamma - SM center} (deg)");
330 fHIncidentAnglePairSM[iSM]->SetYTitle("p_{T #gamma} (GeV/c)");
331 fOutputContainer->Add(fHIncidentAnglePairSM[iSM]);
333 snprintf(hname, buffersize, "hasym_SM%d",iSM);
334 snprintf(htitl, buffersize, "Asymmetry for super mod %d",iSM);
335 fHAsymmetrySM[iSM] = new TH2F(hname,htitl,100,0.,1.,100,0,10);
336 fHAsymmetrySM[iSM]->SetXTitle("a");
337 fHAsymmetrySM[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
338 fOutputContainer->Add(fHAsymmetrySM[iSM]);
340 snprintf(hname,buffersize, "hasym_PairSM%d",iSM);
341 snprintf(htitl,buffersize, "Asymmetry for SM pair: %d",iSM);
342 fHAsymmetryPairSM[iSM] = new TH2F(hname,htitl,100,0.,1.,100,0,10);
343 fHAsymmetryPairSM[iSM]->SetXTitle("a");
344 fHAsymmetryPairSM[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
345 fOutputContainer->Add(fHAsymmetryPairSM[iSM]);
350 fhTowerDecayPhotonHit[iSM] = new TH2F (Form("hTowerDecPhotonHit_Mod%d",iSM),
351 Form("Entries in grid of cells in Module %d",iSM),
352 colmax+2,-1.5,colmax+0.5, rowmax+2,-1.5,rowmax+0.5);
353 fhTowerDecayPhotonHit[iSM]->SetYTitle("row (phi direction)");
354 fhTowerDecayPhotonHit[iSM]->SetXTitle("column (eta direction)");
355 fOutputContainer->Add(fhTowerDecayPhotonHit[iSM]);
357 fhTowerDecayPhotonEnergy[iSM] = new TH2F (Form("hTowerDecPhotonEnergy_Mod%d",iSM),
358 Form("Accumulated energy in grid of cells in Module %d",iSM),
359 colmax+2,-1.5,colmax+0.5, rowmax+2,-1.5,rowmax+0.5);
360 fhTowerDecayPhotonEnergy[iSM]->SetYTitle("row (phi direction)");
361 fhTowerDecayPhotonEnergy[iSM]->SetXTitle("column (eta direction)");
362 fOutputContainer->Add(fhTowerDecayPhotonEnergy[iSM]);
364 fhTowerDecayPhotonAsymmetry[iSM] = new TH2F (Form("hTowerDecPhotonAsymmetry_Mod%d",iSM),
365 Form("Accumulated asymmetry in grid of cells in Module %d",iSM),
366 colmax+2,-1.5,colmax+0.5, rowmax+2,-1.5,rowmax+0.5);
367 fhTowerDecayPhotonAsymmetry[iSM]->SetYTitle("row (phi direction)");
368 fhTowerDecayPhotonAsymmetry[iSM]->SetXTitle("column (eta direction)");
369 fOutputContainer->Add(fhTowerDecayPhotonAsymmetry[iSM]);
371 fhTowerDecayPhotonHitMaskFrame[iSM] = new TH2F (Form("hTowerDecPhotonHit_Mod%d_MaskFrame",iSM),Form("Entries in grid of cells in Module %d",iSM),
372 colmax+2,-1.5,colmax+0.5, rowmax+2,-1.5,rowmax+0.5);
373 fhTowerDecayPhotonHitMaskFrame[iSM]->SetYTitle("row (phi direction)");
374 fhTowerDecayPhotonHitMaskFrame[iSM]->SetXTitle("column (eta direction)");
375 fOutputContainer->Add(fhTowerDecayPhotonHitMaskFrame[iSM]);
377 fhClusterTimeSM[iSM] = new TH2F(Form("hClusterTime_SM%d",iSM),"cluster time vs E",100,0,10, 100,0,1000);
378 fhClusterTimeSM[iSM]->SetXTitle("E (GeV)");
379 fhClusterTimeSM[iSM]->SetYTitle("t (ns)");
380 fOutputContainer->Add(fhClusterTimeSM[iSM]);
382 fhClusterPairDiffTimeSameSM[iSM] = new TH2F(Form("hClusterPairDiffTimeSameSM%d",iSM),
383 Form("cluster pair time difference vs E, SM %d",iSM),
384 100,0,10, 200,-100,100);
385 fhClusterPairDiffTimeSameSM[iSM]->SetXTitle("E (GeV)");
386 fhClusterPairDiffTimeSameSM[iSM]->SetYTitle("#Delta t (ns)");
387 fOutputContainer->Add(fhClusterPairDiffTimeSameSM[iSM]);
391 fhClusterTime = new TH2F("hClusterTime","cluster time vs E",100,0,10, 100,0,1000);
392 fhClusterTime->SetXTitle("E (GeV)");
393 fhClusterTime->SetYTitle("t (ns)");
394 fOutputContainer->Add(fhClusterTime);
396 fhClusterPairDiffTime = new TH2F("hClusterPairDiffTime","cluster pair time difference vs E",100,0,10, 200,-100,100);
397 fhClusterPairDiffTime->SetXTitle("E_{pair} (GeV)");
398 fhClusterPairDiffTime->SetYTitle("#Delta t (ns)");
399 fOutputContainer->Add(fhClusterPairDiffTime);
402 fhNEvents = new TH1I("hNEvents", "Number of analyzed events" , 1 , 0 , 1 ) ;
403 fOutputContainer->Add(fhNEvents);
405 fOutputContainer->SetOwner(kTRUE);
407 // fCalibData = new AliEMCALCalibData();
409 PostData(1,fOutputContainer);
413 //__________________________________________________
414 Bool_t AliAnalysisTaskEMCALPi0CalibSelection::MaskFrameCluster(const Int_t iSM, const Int_t ieta) const {
415 //Check if cell is in one of the regions where we have significant amount of material in front of EMCAL
418 if(iSM%2) icol+=48; // Impair SM, shift index [0-47] to [48-96]
420 if (fNMaskCellColumns && fMaskCellColumns) {
421 for (Int_t imask = 0; imask < fNMaskCellColumns; imask++) {
422 if(icol==fMaskCellColumns[imask]) return kTRUE;
430 //__________________________________________________
431 void AliAnalysisTaskEMCALPi0CalibSelection::UserExec(Option_t* /* option */)
433 //Analysis per event.
435 if(fRecoUtils->GetParticleType()!=AliEMCALRecoUtils::kPhoton){
436 printf("Wrong particle type for cluster position recalculation! = %d\n", fRecoUtils->GetParticleType());
440 if(!(((AliESDEvent*)InputEvent())->GetFiredTriggerClasses()).Contains("EMC")) {
441 //printf("Reject Event %d, FiredClass %s\n",(Int_t)Entry(),(((AliESDEvent*)InputEvent())->GetFiredTriggerClasses()).Data());
445 fhNEvents->Fill(0); //Event analyzed
447 //Get the input event
448 AliVEvent* event = 0;
449 if(fFilteredInput) event = AODEvent();
450 else event = InputEvent();
453 printf("Input event not available!\n");
458 printf("AliAnalysisTaskEMCALPi0CalibSelection <<< %s: Event %d >>>\n",event->GetName(), (Int_t)Entry());
461 //Get the primary vertex
463 event->GetPrimaryVertex()->GetXYZ(v) ;
465 if(DebugLevel() > 1) printf("AliAnalysisTaskEMCALPi0CalibSelection Vertex: (%.3f,%.3f,%.3f)\n",v[0],v[1],v[2]);
467 //Int_t runNum = aod->GetRunNumber();
468 //if(DebugLevel() > 1) printf("Run number: %d\n",runNum);
470 Int_t nSM = (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules();
471 //Get the matrix with geometry information
472 if(fhNEvents->GetEntries()==1){
474 printf("AliAnalysisTaskEMCALPi0CalibSelection::UserExec() - Load user defined geometry matrices\n");
475 for(Int_t mod=0; mod < nSM ; mod++){
478 fMatrix[mod]->Print();
479 fEMCALGeo->SetMisalMatrix(fMatrix[mod],mod) ;
483 else if(!gGeoManager){
484 printf("AliAnalysisTaskEMCALPi0CalibSelection::UserExec() - Get geo matrices from data\n");
485 //Still not implemented in AOD, just a workaround to be able to work at least with ESDs
486 if(!strcmp(event->GetName(),"AliAODEvent")) {
488 printf("AliAnalysisTaskEMCALPi0CalibSelection Use ideal geometry, values geometry matrix not kept in AODs.\n");
491 if(DebugLevel() > 1) printf("AliAnalysisTaskEMCALPi0CalibSelection Load Misaligned matrices. \n");
492 AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event) ;
494 printf("AliAnalysisTaskEMCALPi0CalibSelection::UserExec() - This event does not contain ESDs?");
497 for(Int_t mod=0; mod < nSM; mod++){
499 esd->GetEMCALMatrix(mod)->Print();
500 if(esd->GetEMCALMatrix(mod)) fEMCALGeo->SetMisalMatrix(esd->GetEMCALMatrix(mod),mod) ;
503 }//Load matrices from Data
506 if(DebugLevel() > 1) printf("AliAnalysisTaskEMCALPi0CalibSelection Will use fLogWeight %.3f .\n",fLogWeight);
515 Bool_t shared = kFALSE;
521 //Get the list of clusters
522 TRefArray * caloClustersArr = new TRefArray();
523 event->GetEMCALClusters(caloClustersArr);
524 const Int_t kNumberOfEMCALClusters = caloClustersArr->GetEntries() ;
525 if(DebugLevel() > 1) printf("AliAnalysisTaskEMCALPi0CalibSelection - N CaloClusters: %d \n", kNumberOfEMCALClusters);
528 AliVCaloCells *emCells = event->GetEMCALCells();
530 // loop over EMCAL clusters
531 //----------------------------------------------------------
532 // First recalibrate and recalculate energy and position
533 Float_t pos[]={0,0,0};
534 if(fCorrectClusters){
535 for(Int_t iClu=0; iClu<kNumberOfEMCALClusters; iClu++) {
536 AliVCluster *c1 = (AliVCluster *) caloClustersArr->At(iClu);
538 Float_t e1i = c1->E(); // cluster energy before correction
539 if (e1i < fEmin) continue;
540 else if (e1i > fEmax) continue;
541 else if (c1->GetNCells() < fMinNCells) continue;
543 if(fRecoUtils->ClusterContainsBadChannel(fEMCALGeo, c1->GetCellsAbsId(), c1->GetNCells())) continue;
547 printf("Std : i %d, E %f, dispersion %f, m02 %f, m20 %f\n",c1->GetID(),c1->E(),c1->GetDispersion(),c1->GetM02(),c1->GetM20());
548 c1->GetPosition(pos);
549 printf("Std : i %d, x %f, y %f, z %f\n",c1->GetID(), pos[0], pos[1], pos[2]);
552 //Correct cluster energy and position if requested, and not corrected previously, by default Off
553 if(fRecoUtils->IsRecalibrationOn()) {
554 fRecoUtils->RecalibrateClusterEnergy(fEMCALGeo, c1, emCells);
555 fRecoUtils->RecalculateClusterShowerShapeParameters(fEMCALGeo, emCells,c1);
556 fRecoUtils->RecalculateClusterPID(c1);
559 printf("Energy: after recalibration %f; ",c1->E());
561 // Recalculate cluster position
562 fRecoUtils->RecalculateClusterPosition(fEMCALGeo, emCells,c1);
564 // Correct Non-Linearity
566 c1->SetE(fRecoUtils->CorrectClusterEnergyLinearity(c1));
568 printf("after linearity correction %f\n",c1->E());
572 printf("Cor : i %d, E %f, dispersion %f, m02 %f, m20 %f\n",c1->GetID(),c1->E(),c1->GetDispersion(),c1->GetM02(),c1->GetM20());
573 c1->GetPosition(pos);
574 printf("Cor : i %d, x %f, y %f, z %f\n",c1->GetID(), pos[0], pos[1], pos[2]);
579 //----------------------------------------------------------
580 //Now the invariant mass analysis with the corrected clusters
581 for(Int_t iClu=0; iClu<kNumberOfEMCALClusters-1; iClu++) {
583 AliVCluster *c1 = (AliVCluster *) caloClustersArr->At(iClu);
584 if(fRecoUtils->ClusterContainsBadChannel(fEMCALGeo, c1->GetCellsAbsId(), c1->GetNCells())) continue;
586 Float_t e1i = c1->E(); // cluster energy before correction
587 if (e1i < fEmin) continue;
588 else if (e1i > fEmax) continue;
589 else if (c1->GetNCells() < fMinNCells) continue;
593 printf("IMA : i %d, E %f, dispersion %f, m02 %f, m20 %f\n",c1->GetID(),e1i,c1->GetDispersion(),c1->GetM02(),c1->GetM20());
594 c1->GetPosition(pos);
595 printf("IMA : i %d, x %f, y %f, z %f\n",c1->GetID(), pos[0], pos[1], pos[2]);
598 fRecoUtils->GetMaxEnergyCell(fEMCALGeo, emCells,c1,absId1,iSupMod1,ieta1,iphi1,shared);
599 c1->GetMomentum(p1,v);
601 //Check if cluster is in fidutial region, not too close to borders
602 Bool_t in1 = fRecoUtils->CheckCellFiducialRegion(fEMCALGeo, c1, emCells);
603 // Clusters not facing frame structures
604 Bool_t mask1 = MaskFrameCluster(iSupMod1, ieta1);
605 //if(mask1) printf("Reject eta %d SM %d\n",ieta1, iSupMod1);
607 Double_t time1 = c1->GetTOF()*1.e9;
608 fhClusterTime ->Fill(c1->E(),time1);
609 fhClusterTimeSM[iSupMod1]->Fill(c1->E(),time1);
611 // Combine cluster with other clusters and get the invariant mass
612 for (Int_t jClu=iClu+1; jClu<kNumberOfEMCALClusters; jClu++) {
613 AliAODCaloCluster *c2 = (AliAODCaloCluster *) caloClustersArr->At(jClu);
615 if(fRecoUtils->ClusterContainsBadChannel(fEMCALGeo, c2->GetCellsAbsId(), c2->GetNCells())) continue;
617 Float_t e2i = c2->E();
618 if (e2i < fEmin) continue;
619 else if (e2i > fEmax) continue;
620 else if (c2->GetNCells() < fMinNCells) continue;
623 fRecoUtils->GetMaxEnergyCell(fEMCALGeo, emCells,c2,absId2,iSupMod2,ieta2,iphi2,shared);
624 c2->GetMomentum(p2,v);
627 Float_t invmass = p12.M()*1000;
628 //printf("*** mass %f\n",invmass);
631 Float_t asym = TMath::Abs(p1.E()-p2.E())/(p1.E()+p2.E());
632 //printf("asymmetry %f\n",asym);
633 if(asym > fAsyCut) continue;
636 Double_t time2 = c2->GetTOF()*1.e9;
637 fhClusterPairDiffTime->Fill(p12.E(),time1-time2);
638 if(TMath::Abs(time1-time2) > fDTimeCut) continue;
640 if(invmass < fMaxBin && invmass > fMinBin ){
642 //Check if cluster is in fidutial region, not too close to borders
643 Bool_t in2 = fRecoUtils->CheckCellFiducialRegion(fEMCALGeo, c2, emCells);
645 // Clusters not facing frame structures
646 Bool_t mask2 = MaskFrameCluster(iSupMod2, ieta2);
647 //if(mask2) printf("Reject eta %d SM %d\n",ieta2, iSupMod2);
651 fHmgg->Fill(invmass,p12.Pt());
653 if(iSupMod1==iSupMod2) {
654 fHmggSM[iSupMod1]->Fill(invmass,p12.Pt());
655 fhClusterPairDiffTimeSameSM[iSupMod1]->Fill(p12.E(),time1-time2);
658 fHmggDifferentSM ->Fill(invmass,p12.Pt());
662 for(Int_t i = 0; i < nSM/2; i++){
664 if((iSupMod1==j && iSupMod2==j+1) || (iSupMod1==j+1 && iSupMod2==j)) {
665 fHmggPairSameSectorSM[i]->Fill(invmass,p12.Pt());
666 fhClusterPairDiffTimeSameSector[i]->Fill(p12.E(),time1-time2);
671 for(Int_t i = 0; i < nSM-2; i++){
672 if((iSupMod1==i && iSupMod2==i+2) || (iSupMod1==i+2 && iSupMod2==i)) {
673 fHmggPairSameSideSM[i]->Fill(invmass,p12.Pt());
674 fhClusterPairDiffTimeSameSide[i]->Fill(p12.E(),time1-time2);
679 if(!mask1 && !mask2){
681 fHmggMaskFrame->Fill(invmass,p12.Pt());
683 if(iSupMod1==iSupMod2) fHmggSMMaskFrame[iSupMod1]->Fill(invmass,p12.Pt());
684 else fHmggDifferentSMMaskFrame ->Fill(invmass,p12.Pt());
688 for(Int_t i = 0; i < nSM/2; i++){
690 if((iSupMod1==j && iSupMod2==j+1) || (iSupMod1==j+1 && iSupMod2==j)) fHmggPairSameSectorSMMaskFrame[i]->Fill(invmass,p12.Pt());
694 for(Int_t i = 0; i < nSM-2; i++){
695 if((iSupMod1==i && iSupMod2==i+2) || (iSupMod1==i+2 && iSupMod2==i)) fHmggPairSameSideSMMaskFrame[i]->Fill(invmass,p12.Pt());
698 }// Pair not facing frame
701 if(invmass > 100. && invmass < 160.){//restrict to clusters really close to pi0 peak
703 //Opening angle of 2 photons
704 Float_t opangle = p1.Angle(p2.Vect())*TMath::RadToDeg();
705 //printf("*******>>>>>>>> In PEAK pt %f, angle %f \n",p12.Pt(),opangle);
707 //Incident angle of each photon
708 Float_t inangle1 =0., inangle2=0.;
710 Float_t posSM1cen[3]={0.,0.,0.};
711 Float_t depth = fRecoUtils->GetDepth(p1.Energy(),fRecoUtils->GetParticleType(),iSupMod1);
712 fEMCALGeo->RecalculateTowerPosition(11.5, 23.5, iSupMod1, depth, fRecoUtils->GetMisalTransShiftArray(),fRecoUtils->GetMisalRotShiftArray(),posSM1cen);
713 Float_t posSM2cen[3]={0.,0.,0.};
714 depth = fRecoUtils->GetDepth(p2.Energy(),fRecoUtils->GetParticleType(),iSupMod2);
715 fEMCALGeo->RecalculateTowerPosition(11.5, 23.5, iSupMod2, depth, fRecoUtils->GetMisalTransShiftArray(),fRecoUtils->GetMisalRotShiftArray(),posSM2cen);
716 //printf("SM1 %d pos (%2.3f,%2.3f,%2.3f) \n",iSupMod1,posSM1cen[0],posSM1cen[1],posSM1cen[2]);
717 //printf("SM2 %d pos (%2.3f,%2.3f,%2.3f) \n",iSupMod2,posSM2cen[0],posSM2cen[1],posSM2cen[2]);
719 TVector3 vecSM1cen(posSM1cen[0]-v[0],posSM1cen[1]-v[1],posSM1cen[2]-v[2]);
720 TVector3 vecSM2cen(posSM2cen[0]-v[0],posSM2cen[1]-v[1],posSM2cen[2]-v[2]);
721 inangle1 = p1.Angle(vecSM1cen)*TMath::RadToDeg();
722 inangle2 = p2.Angle(vecSM2cen)*TMath::RadToDeg();
723 //printf("Incident angle: cluster 1 %2.3f; cluster 2 %2.3f\n",inangle1,inangle2);
726 fHOpeningAngle ->Fill(opangle,p12.Pt());
727 fHIncidentAngle->Fill(inangle1,p1.Pt());
728 fHIncidentAngle->Fill(inangle2,p2.Pt());
729 fHAsymmetry ->Fill(asym,p12.Pt());
731 if(iSupMod1==iSupMod2) {
732 fHOpeningAngleSM[iSupMod1] ->Fill(opangle,p12.Pt());
733 fHIncidentAngleSM[iSupMod1]->Fill(inangle1,p1.Pt());
734 fHIncidentAngleSM[iSupMod1]->Fill(inangle2,p2.Pt());
735 fHAsymmetrySM[iSupMod1] ->Fill(asym,p12.Pt());
738 fHOpeningAngleDifferentSM ->Fill(opangle,p12.Pt());
739 fHIncidentAngleDifferentSM ->Fill(inangle1,p1.Pt());
740 fHIncidentAngleDifferentSM ->Fill(inangle2,p2.Pt());
741 fHAsymmetryDifferentSM ->Fill(asym,p12.Pt());
744 if((iSupMod1==0 && iSupMod2==2) || (iSupMod1==2 && iSupMod2==0)) {
745 fHOpeningAnglePairSM[0] ->Fill(opangle,p12.Pt());
746 fHIncidentAnglePairSM[0]->Fill(inangle1,p1.Pt());
747 fHIncidentAnglePairSM[0]->Fill(inangle2,p2.Pt());
748 fHAsymmetryPairSM[0] ->Fill(asym,p12.Pt());
751 if((iSupMod1==1 && iSupMod2==3) || (iSupMod1==3 && iSupMod2==1)) {
752 fHOpeningAnglePairSM[1] ->Fill(opangle,p12.Pt());
753 fHIncidentAnglePairSM[1]->Fill(inangle1,p1.Pt());
754 fHIncidentAnglePairSM[1]->Fill(inangle2,p2.Pt());
755 fHAsymmetryPairSM[1] ->Fill(asym,p12.Pt());
759 if((iSupMod1==0 && iSupMod2==1) || (iSupMod1==1 && iSupMod2==0)) {
760 fHOpeningAnglePairSM[2] ->Fill(opangle,p12.Pt());
761 fHIncidentAnglePairSM[2]->Fill(inangle1,p1.Pt());
762 fHIncidentAnglePairSM[2]->Fill(inangle2,p2.Pt());
763 fHAsymmetryPairSM[2] ->Fill(asym,p12.Pt());
767 if((iSupMod1==2 && iSupMod2==3) || (iSupMod1==3 && iSupMod2==2)) {
768 fHOpeningAnglePairSM[3] ->Fill(opangle,p12.Pt());
769 fHIncidentAnglePairSM[3]->Fill(inangle1,p1.Pt());
770 fHIncidentAnglePairSM[3]->Fill(inangle2,p2.Pt());
771 fHAsymmetryPairSM[3] ->Fill(asym,p12.Pt());
774 }// pair in 100 < mass < 160
776 }//in acceptance cuts
778 //In case of filling only channels with second cluster in same SM
779 if(fSameSM && iSupMod1!=iSupMod2) continue;
781 if (fGroupNCells == 0){
782 fHmpi0[iSupMod1][ieta1][iphi1]->Fill(invmass);
783 fHmpi0[iSupMod2][ieta2][iphi2]->Fill(invmass);
785 if(invmass > 100. && invmass < 160.){//restrict to clusters really close to pi0 peak
786 fhTowerDecayPhotonHit [iSupMod1]->Fill(ieta1,iphi1);
787 fhTowerDecayPhotonEnergy [iSupMod1]->Fill(ieta1,iphi1,p1.E());
788 fhTowerDecayPhotonAsymmetry[iSupMod1]->Fill(ieta1,iphi1,asym);
790 fhTowerDecayPhotonHit [iSupMod2]->Fill(ieta2,iphi2);
791 fhTowerDecayPhotonEnergy [iSupMod2]->Fill(ieta2,iphi2,p2.E());
792 fhTowerDecayPhotonAsymmetry[iSupMod2]->Fill(ieta2,iphi2,asym);
794 if(!mask1)fhTowerDecayPhotonHitMaskFrame[iSupMod1]->Fill(ieta1,iphi1);
795 if(!mask2)fhTowerDecayPhotonHitMaskFrame[iSupMod2]->Fill(ieta2,iphi2);
797 }// pair in mass of pi0
800 //printf("Regroup N %d, eta1 %d, phi1 %d, eta2 %d, phi2 %d \n",fGroupNCells, ieta1, iphi1, ieta2, iphi2);
801 for (Int_t i = -fGroupNCells; i < fGroupNCells+1; i++) {
802 for (Int_t j = -fGroupNCells; j < fGroupNCells+1; j++) {
803 //printf("\t i %d, j %d\n",i,j);
804 if((ieta1+i >= 0) && (iphi1+j >= 0) && (ieta1+i < 48) && (iphi1+j < 24)){
805 //printf("\t \t eta1+i %d, phi1+j %d\n", ieta1+i, iphi1+j);
806 fHmpi0[iSupMod1][ieta1+i][iphi1+j]->Fill(invmass);
808 if((ieta2+i >= 0) && (iphi2+j >= 0) && (ieta2+i < 48) && (iphi2+j < 24)){
809 //printf("\t \t eta2+i %d, phi2+j %d\n", ieta2+i, iphi2+j);
810 fHmpi0[iSupMod2][ieta2+i][iphi2+j]->Fill(invmass);
816 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",
817 iSupMod1,iphi1,ieta1,iSupMod2,iphi2,ieta2,p12.M(),e1i,c1->E(),e2i,c2->E());
822 } // end of loop over EMCAL clusters
824 delete caloClustersArr;
826 PostData(1,fOutputContainer);
830 //_____________________________________________________
831 void AliAnalysisTaskEMCALPi0CalibSelection::PrintInfo(){
834 printf("Cluster cuts: %2.2f < E < %2.2f GeV; number of cells > %d; Assymetry < %1.2f, pair time diff < %2.2f\n",
835 fEmin,fEmax, fMinNCells, fAsyCut, fDTimeCut) ;
836 printf("Group %d cells\n", fGroupNCells) ;
837 printf("Cluster maximal cell away from border at least %d cells\n", fRecoUtils->GetNumberOfCellsFromEMCALBorder()) ;
838 printf("Histograms: bins %d; energy range: %2.2f < E < %2.2f GeV\n",fNbins,fMinBin,fMaxBin) ;
839 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",
840 fRecoUtils->IsBadChannelsRemovalSwitchedOn(),fFilteredInput,fCorrectClusters, fSameSM) ;
841 printf("EMCAL Geometry name: < %s >, Load Matrices %d\n",fEMCALGeoName.Data(), fLoadMatrices) ;
842 if(fLoadMatrices) {for(Int_t ism = 0; ism < AliEMCALGeoParams::fgkEMCALModules; ism++) fMatrix[ism]->Print() ; }