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