]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGGA/EMCALTasks/AliAnalysisTaskEMCALPi0CalibSelection.cxx
submitted AddTasks for GG pPb
[u/mrichter/AliRoot.git] / PWGGA / EMCALTasks / AliAnalysisTaskEMCALPi0CalibSelection.cxx
CommitLineData
375cec9b 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
cd231d42 3
375cec9b 4 * Permission to use, copy, modify and distribute this software and its *
5 * documentation strictly for non-commercial purposes is hereby granted *
6 * without fee, provided that the above copyright notice appears in all *
7 * copies and that both the copyright notice and this permission notice *
8 * appear in the supporting documentation. The authors make no claims *
9 * about the suitability of this software for any purpose. It is *
10 * provided "as is" without express or implied warranty. *
11 **************************************************************************/
12
a8dc7d71 13// $Id$
477c5cd2 14
375cec9b 15//---------------------------------------------------------------------------//
16// //
17// Fill histograms (one per cell) with two-cluster invariant mass //
18// using calibration coefficients of the previous iteration. //
19// Histogram for a given cell is filled if the most energy of one cluster //
20// is deposited in this cell and the other cluster could be anywherein EMCAL.//
21// //
cd231d42 22// //
23// Author: Boris Polishchuk //
24// Adapted to AOD reading by Gustavo Conesa //
25// //
cd231d42 26// //
375cec9b 27//---------------------------------------------------------------------------//
28
375cec9b 29// Root
30#include "TLorentzVector.h"
375cec9b 31#include "TRefArray.h"
32#include "TList.h"
33#include "TH1F.h"
247abff4 34#include <TGeoManager.h>
375cec9b 35
36// AliRoot
37#include "AliAnalysisTaskEMCALPi0CalibSelection.h"
38#include "AliAODEvent.h"
39#include "AliESDEvent.h"
375cec9b 40#include "AliEMCALGeometry.h"
c8fe2783 41#include "AliVCluster.h"
42#include "AliVCaloCells.h"
9584c261 43#include "AliEMCALRecoUtils.h"
a8dc7d71 44#include "AliOADBContainer.h"
375cec9b 45
46ClassImp(AliAnalysisTaskEMCALPi0CalibSelection)
47
375cec9b 48
477c5cd2 49//______________________________________________________________________________________________
375cec9b 50AliAnalysisTaskEMCALPi0CalibSelection::AliAnalysisTaskEMCALPi0CalibSelection(const char* name) :
7b2d541a 51AliAnalysisTaskSE(name),
52fEMCALGeo(0x0), fLoadMatrices(0),
53fEMCALGeoName("EMCAL_COMPLETE12SMV1"),
54fTriggerName("EMC"),
55fRecoUtils(new AliEMCALRecoUtils),
56fOADBFilePath(""), fCorrectClusters(kFALSE),
57fCaloClustersArr(0x0), fEMCALCells(0x0),
58fCuts(0x0), fOutputContainer(0x0),
59fVertex(), fFilteredInput(kFALSE),
49b53920 60fEmin(0.5), fEmax(15.),
477c5cd2 61fL0min(0.01), fL0max(0.5),
62fDTimeCut(100.), fTimeMax(1000000), fTimeMin(-1000000),
a7e5a381 63fAsyCut(1.), fMinNCells(2), fGroupNCells(0),
7b2d541a 64fLogWeight(4.5), fSameSM(kFALSE),
49b53920 65fNMaskCellColumns(11), fMaskCellColumns(0x0),
a7e5a381 66fInvMassCutMin(110.), fInvMassCutMax(160.),
49b53920 67//Histograms
7b2d541a 68fNbins(300),
a7e5a381 69fMinBin(0.), fMaxBin(300.),
70fNTimeBins(1000), fMinTimeBin(0.), fMaxTimeBin(1000.),
49b53920 71fHmgg(0x0), fHmggDifferentSM(0x0),
72fHmggMaskFrame(0x0), fHmggDifferentSMMaskFrame(0x0),
73fHOpeningAngle(0x0), fHOpeningAngleDifferentSM(0x0),
49b53920 74fHAsymmetry(0x0), fHAsymmetryDifferentSM(0x0),
75fhNEvents(0x0),
76fhClusterTime(0x0), fhClusterPairDiffTime(0x0)
375cec9b 77{
78 //Named constructor which should be used.
79
44cf05d7 80 for(Int_t iMod=0; iMod < AliEMCALGeoParams::fgkEMCALModules; iMod++) {
bdd2a262 81 for(Int_t iX=0; iX<24; iX++) {
82 for(Int_t iZ=0; iZ<48; iZ++) {
a7e5a381 83 fHmpi0[iMod][iZ][iX] = 0 ;
375cec9b 84 }
85 }
86 }
6eb2a715 87
7b2d541a 88 fVertex[0]=fVertex[1]=fVertex[2]=-1000;
89
a7e5a381 90 fHTpi0[0]= 0 ;
91 fHTpi0[1]= 0 ;
92 fHTpi0[2]= 0 ;
93 fHTpi0[3]= 0 ;
94
42b19289 95 fMaskCellColumns = new Int_t[fNMaskCellColumns];
42b19289 96 fMaskCellColumns[0] = 6 ; fMaskCellColumns[1] = 7 ; fMaskCellColumns[2] = 8 ;
97 fMaskCellColumns[3] = 35; fMaskCellColumns[4] = 36; fMaskCellColumns[5] = 37;
98 fMaskCellColumns[6] = 12+AliEMCALGeoParams::fgkEMCALCols; fMaskCellColumns[7] = 13+AliEMCALGeoParams::fgkEMCALCols;
99 fMaskCellColumns[8] = 40+AliEMCALGeoParams::fgkEMCALCols; fMaskCellColumns[9] = 41+AliEMCALGeoParams::fgkEMCALCols;
100 fMaskCellColumns[10]= 42+AliEMCALGeoParams::fgkEMCALCols;
101
477c5cd2 102 for(Int_t iSMPair = 0; iSMPair < AliEMCALGeoParams::fgkEMCALModules/2; iSMPair++)
103 {
42b19289 104 fHmggPairSameSectorSM[iSMPair] = 0;
105 fHmggPairSameSectorSMMaskFrame[iSMPair] = 0;
af2d7c9b 106 fhClusterPairDiffTimeSameSector[iSMPair]= 0;
42b19289 107 }
477c5cd2 108
109 for(Int_t iSMPair = 0; iSMPair < AliEMCALGeoParams::fgkEMCALModules-2; iSMPair++)
110 {
42b19289 111 fHmggPairSameSideSM[iSMPair] = 0;
112 fHmggPairSameSideSMMaskFrame[iSMPair] = 0;
af2d7c9b 113 fhClusterPairDiffTimeSameSide[iSMPair] = 0;
42b19289 114 }
1dabc151 115
477c5cd2 116 for(Int_t iSM = 0; iSM < AliEMCALGeoParams::fgkEMCALModules; iSM++)
117 {
44cf05d7 118 fHmggSM[iSM] = 0;
42b19289 119 fHmggSMMaskFrame[iSM] = 0;
44cf05d7 120 fHOpeningAngleSM[iSM] = 0;
121 fHOpeningAnglePairSM[iSM] = 0;
122 fHAsymmetrySM[iSM] = 0;
123 fHAsymmetryPairSM[iSM] = 0;
44cf05d7 124 fhTowerDecayPhotonHit[iSM] = 0;
125 fhTowerDecayPhotonEnergy[iSM] = 0;
126 fhTowerDecayPhotonAsymmetry[iSM] = 0;
42b19289 127 fhTowerDecayPhotonHitMaskFrame[iSM]= 0;
44cf05d7 128 fMatrix[iSM] = 0x0;
af2d7c9b 129 fhClusterTimeSM[iSM] = 0;
130 fhClusterPairDiffTimeSameSM[iSM] = 0;
2dfb1428 131 }
132
cf028690 133 DefineOutput(1, TList::Class());
6eb2a715 134 DefineOutput(2, TList::Class()); // will contain cuts or local params
af2d7c9b 135
375cec9b 136}
137
477c5cd2 138//_____________________________________________________________________________
375cec9b 139AliAnalysisTaskEMCALPi0CalibSelection::~AliAnalysisTaskEMCALPi0CalibSelection()
140{
141 //Destructor.
142
477c5cd2 143 if(fOutputContainer)
144 {
375cec9b 145 fOutputContainer->Delete() ;
146 delete fOutputContainer ;
147 }
af2d7c9b 148
42b19289 149 if(fEMCALGeo) delete fEMCALGeo ;
150 if(fRecoUtils) delete fRecoUtils ;
151 if(fNMaskCellColumns) delete [] fMaskCellColumns;
152
375cec9b 153}
154
7b2d541a 155//____________________________________________________________
156void AliAnalysisTaskEMCALPi0CalibSelection::CorrectClusters()
a8dc7d71 157{
7b2d541a 158 // loop over EMCAL clusters
159 //----------------------------------------------------------
160 // First recalibrate and recalculate energy and position
a8dc7d71 161
a8dc7d71 162
7b2d541a 163 if(fCorrectClusters)
a8dc7d71 164 {
a8dc7d71 165
7b2d541a 166 if(fRecoUtils->GetParticleType()!=AliEMCALRecoUtils::kPhoton)
167 {
168 printf("Wrong particle type for cluster position recalculation! = %d\n", fRecoUtils->GetParticleType());
169 abort();
170 }
43dcae1f 171
7b2d541a 172 if(DebugLevel() > 1) printf("AliAnalysisTaskEMCALPi0CalibSelection Will use fLogWeight %.3f .\n",fLogWeight);
a8dc7d71 173
7b2d541a 174 Float_t pos[]={0,0,0};
a8dc7d71 175
7b2d541a 176 for(Int_t iClu=0; iClu < fCaloClustersArr->GetEntriesFast(); iClu++)
a8dc7d71 177 {
7b2d541a 178 AliVCluster *c1 = (AliVCluster *) fCaloClustersArr->At(iClu);
a8dc7d71 179
7b2d541a 180 Float_t e1i = c1->E(); // cluster energy before correction
181 if (e1i < fEmin) continue;
182 else if (e1i > fEmax) continue;
a8dc7d71 183
7b2d541a 184 else if (c1->GetNCells() < fMinNCells) continue;
185
186 else if (c1->GetM02() < fL0min || c1->GetM02() > fL0max) continue;
187
188 if(fRecoUtils->ClusterContainsBadChannel(fEMCALGeo, c1->GetCellsAbsId(), c1->GetNCells())) continue;
189
190 if(DebugLevel() > 2)
191 {
192 printf("Std : i %d, E %f, dispersion %f, m02 %f, m20 %f\n",c1->GetID(),c1->E(),c1->GetDispersion(),c1->GetM02(),c1->GetM20());
193 c1->GetPosition(pos);
194 printf("Std : i %d, x %f, y %f, z %f\n",c1->GetID(), pos[0], pos[1], pos[2]);
a8dc7d71 195 }
a8dc7d71 196
7b2d541a 197 //Correct cluster energy and position if requested, and not corrected previously, by default Off
198 if(fRecoUtils->IsRecalibrationOn())
a8dc7d71 199 {
7b2d541a 200 fRecoUtils->RecalibrateClusterEnergy(fEMCALGeo, c1, fEMCALCells);
201 fRecoUtils->RecalculateClusterShowerShapeParameters(fEMCALGeo, fEMCALCells,c1);
202 fRecoUtils->RecalculateClusterPID(c1);
203 }
204
205 if(DebugLevel() > 2)
206 printf("Energy: after recalibration %f; \n",c1->E());
207
208 // Recalculate cluster position
209 fRecoUtils->RecalculateClusterPosition(fEMCALGeo, fEMCALCells,c1);
210
211 // Correct Non-Linearity
212 c1->SetE(fRecoUtils->CorrectClusterEnergyLinearity(c1));
213
214 if(DebugLevel() > 2)
215 printf("\t after linearity correction %f\n",c1->E());
216
217 //In case of MC analysis, to match resolution/calibration in real data
218 c1->SetE(fRecoUtils->SmearClusterEnergy(c1));
219
220 if(DebugLevel() > 2)
221 printf("\t after smearing %f\n",c1->E());
222
223 if(DebugLevel() > 2)
224 {
225 printf("Cor : i %d, E %f, dispersion %f, m02 %f, m20 %f\n",c1->GetID(),c1->E(),c1->GetDispersion(),c1->GetM02(),c1->GetM20());
226 c1->GetPosition(pos);
227 printf("Cor : i %d, x %f, y %f, z %f\n",c1->GetID(), pos[0], pos[1], pos[2]);
228 }
229 }
230 }
a8dc7d71 231}
232
7b2d541a 233//__________________________________________________________
234void AliAnalysisTaskEMCALPi0CalibSelection::FillHistograms()
375cec9b 235{
7b2d541a 236 // Now fill the invariant mass analysis with the corrected clusters, and other general histograms
237
238 Int_t absId1 = -1;
239 Int_t iSupMod1 = -1;
240 Int_t iphi1 = -1;
241 Int_t ieta1 = -1;
242 Int_t absId2 = -1;
243 Int_t iSupMod2 = -1;
244 Int_t iphi2 = -1;
245 Int_t ieta2 = -1;
246 Bool_t shared = kFALSE;
42b19289 247
7b2d541a 248 TLorentzVector p1;
249 TLorentzVector p2;
250 TLorentzVector p12;
42b19289 251
7b2d541a 252 Float_t pos[]={0,0,0};
42b19289 253
7b2d541a 254 Int_t bc = InputEvent()->GetBunchCrossNumber();
255 Int_t nSM = (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules();
256
257 for(Int_t iClu=0; iClu<fCaloClustersArr->GetEntriesFast()-1; iClu++)
a8dc7d71 258 {
7b2d541a 259 AliVCluster *c1 = (AliVCluster *) fCaloClustersArr->At(iClu);
2dfb1428 260
7b2d541a 261 if(fRecoUtils->ClusterContainsBadChannel(fEMCALGeo, c1->GetCellsAbsId(), c1->GetNCells())) continue;
42b19289 262
7b2d541a 263 Float_t e1i = c1->E(); // cluster energy before correction
42b19289 264
7b2d541a 265 if (e1i < fEmin) continue;
266 else if (e1i > fEmax) continue;
267
268 else if (!fRecoUtils->IsGoodCluster(c1,fEMCALGeo,fEMCALCells,bc)) continue;
269
270 else if (c1->GetNCells() < fMinNCells) continue;
271
272 else if (c1->GetM02() < fL0min || c1->GetM02() > fL0max) continue;
273
274 if(DebugLevel() > 2)
275 {
276 printf("IMA : i %d, E %f, dispersion %f, m02 %f, m20 %f\n",c1->GetID(),e1i,c1->GetDispersion(),c1->GetM02(),c1->GetM20());
277 c1->GetPosition(pos);
278 printf("IMA : i %d, x %f, y %f, z %f\n",c1->GetID(), pos[0], pos[1], pos[2]);
279 }
280
281 fRecoUtils->GetMaxEnergyCell(fEMCALGeo, fEMCALCells,c1,absId1,iSupMod1,ieta1,iphi1,shared);
282 c1->GetMomentum(p1,fVertex);
283
284 //Check if cluster is in fidutial region, not too close to borders
285 Bool_t in1 = fRecoUtils->CheckCellFiducialRegion(fEMCALGeo, c1, fEMCALCells);
286
287 // Clusters not facing frame structures
288 Bool_t mask1 = MaskFrameCluster(iSupMod1, ieta1);
289 //if(mask1) printf("Reject eta %d SM %d\n",ieta1, iSupMod1);
290
291 Double_t time1 = c1->GetTOF()*1.e9;
292
293 if(time1 > fTimeMax || time1 < fTimeMin) continue;
294
295 fhClusterTime ->Fill(c1->E(),time1);
296 fhClusterTimeSM[iSupMod1]->Fill(c1->E(),time1);
297
298 // Combine cluster with other clusters and get the invariant mass
299 for (Int_t jClu=iClu+1; jClu < fCaloClustersArr->GetEntriesFast(); jClu++)
a8dc7d71 300 {
7b2d541a 301 AliAODCaloCluster *c2 = (AliAODCaloCluster *) fCaloClustersArr->At(jClu);
42b19289 302
7b2d541a 303 Float_t e2i = c2->E();
304 if (e2i < fEmin) continue;
305 else if (e2i > fEmax) continue;
49b53920 306
7b2d541a 307 else if (!fRecoUtils->IsGoodCluster(c2,fEMCALGeo,fEMCALCells,bc))continue;
49b53920 308
7b2d541a 309 else if (c2->GetNCells() < fMinNCells) continue;
49b53920 310
7b2d541a 311 else if (c2->GetM02() < fL0min || c2->GetM02() > fL0max) continue;
312
313
314 fRecoUtils->GetMaxEnergyCell(fEMCALGeo, fEMCALCells,c2,absId2,iSupMod2,ieta2,iphi2,shared);
315 c2->GetMomentum(p2,fVertex);
316
317 p12 = p1+p2;
318 Float_t invmass = p12.M()*1000;
319
320 //Asimetry cut
321 Float_t asym = TMath::Abs(p1.E()-p2.E())/(p1.E()+p2.E());
322
323 if(asym > fAsyCut) continue;
324
325 //Time cut
326 Double_t time2 = c2->GetTOF()*1.e9;
327
328 if(time2 > fTimeMax || time2 < fTimeMin) continue;
329
330 fhClusterPairDiffTime->Fill(p12.E(),time1-time2);
331 if(TMath::Abs(time1-time2) > fDTimeCut) continue;
332
333 if(invmass < fMaxBin && invmass > fMinBin )
334 {
335 //Check if cluster is in fidutial region, not too close to borders
336 Bool_t in2 = fRecoUtils->CheckCellFiducialRegion(fEMCALGeo, c2, fEMCALCells);
337
338 // Clusters not facing frame structures
339 Bool_t mask2 = MaskFrameCluster(iSupMod2, ieta2);
340 //if(mask2) printf("Reject eta %d SM %d\n",ieta2, iSupMod2);
341
342 if(in1 && in2)
343 {
344 fHmgg->Fill(invmass,p12.Pt());
345
346 if(iSupMod1==iSupMod2)
347 {
348 fHmggSM[iSupMod1]->Fill(invmass,p12.Pt());
349 fhClusterPairDiffTimeSameSM[iSupMod1]->Fill(p12.E(),time1-time2);
350 }
351 else
352 fHmggDifferentSM ->Fill(invmass,p12.Pt());
353
354 // Same sector
355 Int_t j=0;
356 for(Int_t i = 0; i < nSM/2; i++)
357 {
358 j=2*i;
359 if((iSupMod1==j && iSupMod2==j+1) || (iSupMod1==j+1 && iSupMod2==j))
360 {
361 fHmggPairSameSectorSM[i]->Fill(invmass,p12.Pt());
362 fhClusterPairDiffTimeSameSector[i]->Fill(p12.E(),time1-time2);
363 }
364 }
365
366 // Same side
367 for(Int_t i = 0; i < nSM-2; i++)
368 {
369 if((iSupMod1==i && iSupMod2==i+2) || (iSupMod1==i+2 && iSupMod2==i))
370 {
371 fHmggPairSameSideSM[i]->Fill(invmass,p12.Pt());
372 fhClusterPairDiffTimeSameSide[i]->Fill(p12.E(),time1-time2);
373 }
374 }
375
376
377 if(!mask1 && !mask2)
378 {
379 fHmggMaskFrame->Fill(invmass,p12.Pt());
380
381 if(iSupMod1==iSupMod2) fHmggSMMaskFrame[iSupMod1]->Fill(invmass,p12.Pt());
382 else fHmggDifferentSMMaskFrame ->Fill(invmass,p12.Pt());
383
384 // Same sector
385 j=0;
386 for(Int_t i = 0; i < nSM/2; i++)
387 {
388 j=2*i;
389 if((iSupMod1==j && iSupMod2==j+1) || (iSupMod1==j+1 && iSupMod2==j)) fHmggPairSameSectorSMMaskFrame[i]->Fill(invmass,p12.Pt());
390 }
391
392 // Same side
393 for(Int_t i = 0; i < nSM-2; i++)
394 {
395 if((iSupMod1==i && iSupMod2==i+2) || (iSupMod1==i+2 && iSupMod2==i)) fHmggPairSameSideSMMaskFrame[i]->Fill(invmass,p12.Pt());
396 }
397
398 }// Pair not facing frame
399
400
401 if(invmass > fInvMassCutMin && invmass < fInvMassCutMax) //restrict to clusters really close to pi0 peak
402 {
403
404 // Check time of cells in both clusters, and fill time histogram
405 for(Int_t icell = 0; icell < c1->GetNCells(); icell++)
406 {
407 Int_t absID = c1->GetCellAbsId(icell);
408 fHTpi0[bc%4]->Fill(absID, fEMCALCells->GetCellTime(absID)*1.e9);
409 }
410
411 for(Int_t icell = 0; icell < c2->GetNCells(); icell++)
412 {
413 Int_t absID = c2->GetCellAbsId(icell);
414 fHTpi0[bc%4]->Fill(absID, fEMCALCells->GetCellTime(absID)*1.e9);
415 }
416
417 //Opening angle of 2 photons
418 Float_t opangle = p1.Angle(p2.Vect())*TMath::RadToDeg();
419 //printf("*******>>>>>>>> In PEAK pt %f, angle %f \n",p12.Pt(),opangle);
420
421
422 fHOpeningAngle ->Fill(opangle,p12.Pt());
423 fHAsymmetry ->Fill(asym,p12.Pt());
424
425 if(iSupMod1==iSupMod2)
426 {
427 fHOpeningAngleSM[iSupMod1] ->Fill(opangle,p12.Pt());
428 fHAsymmetrySM[iSupMod1] ->Fill(asym,p12.Pt());
429 }
430 else
431 {
432 fHOpeningAngleDifferentSM ->Fill(opangle,p12.Pt());
433 fHAsymmetryDifferentSM ->Fill(asym,p12.Pt());
434 }
435
436 if((iSupMod1==0 && iSupMod2==2) || (iSupMod1==2 && iSupMod2==0))
437 {
438 fHOpeningAnglePairSM[0] ->Fill(opangle,p12.Pt());
439 fHAsymmetryPairSM[0] ->Fill(asym,p12.Pt());
440
441 }
442 if((iSupMod1==1 && iSupMod2==3) || (iSupMod1==3 && iSupMod2==1))
443 {
444 fHOpeningAnglePairSM[1] ->Fill(opangle,p12.Pt());
445 fHAsymmetryPairSM[1] ->Fill(asym,p12.Pt());
446 }
447
448 if((iSupMod1==0 && iSupMod2==1) || (iSupMod1==1 && iSupMod2==0))
449 {
450 fHOpeningAnglePairSM[2] ->Fill(opangle,p12.Pt());
451 fHAsymmetryPairSM[2] ->Fill(asym,p12.Pt());
452 }
453 if((iSupMod1==2 && iSupMod2==3) || (iSupMod1==3 && iSupMod2==2))
454 {
455 fHOpeningAnglePairSM[3] ->Fill(opangle,p12.Pt());
456 fHAsymmetryPairSM[3] ->Fill(asym,p12.Pt());
457 }
458
459 }// pair in 100 < mass < 160
460
461 }//in acceptance cuts
462
463 //In case of filling only channels with second cluster in same SM
464 if(fSameSM && iSupMod1!=iSupMod2) continue;
465
466 if (fGroupNCells == 0)
467 {
468 fHmpi0[iSupMod1][ieta1][iphi1]->Fill(invmass);
469 fHmpi0[iSupMod2][ieta2][iphi2]->Fill(invmass);
470
471 if(invmass > fInvMassCutMin && invmass < fInvMassCutMax)//restrict to clusters really close to pi0 peak
472 {
473 fhTowerDecayPhotonHit [iSupMod1]->Fill(ieta1,iphi1);
474 fhTowerDecayPhotonEnergy [iSupMod1]->Fill(ieta1,iphi1,p1.E());
475 fhTowerDecayPhotonAsymmetry[iSupMod1]->Fill(ieta1,iphi1,asym);
476
477 fhTowerDecayPhotonHit [iSupMod2]->Fill(ieta2,iphi2);
478 fhTowerDecayPhotonEnergy [iSupMod2]->Fill(ieta2,iphi2,p2.E());
479 fhTowerDecayPhotonAsymmetry[iSupMod2]->Fill(ieta2,iphi2,asym);
480
481 if(!mask1)fhTowerDecayPhotonHitMaskFrame[iSupMod1]->Fill(ieta1,iphi1);
482 if(!mask2)fhTowerDecayPhotonHitMaskFrame[iSupMod2]->Fill(ieta2,iphi2);
483
484 }// pair in mass of pi0
485 }
486 else {
487 //printf("Regroup N %d, eta1 %d, phi1 %d, eta2 %d, phi2 %d \n",fGroupNCells, ieta1, iphi1, ieta2, iphi2);
488 for (Int_t i = -fGroupNCells; i < fGroupNCells+1; i++)
489 {
490 for (Int_t j = -fGroupNCells; j < fGroupNCells+1; j++)
491 {
492 Int_t absId11 = fEMCALGeo->GetAbsCellIdFromCellIndexes(iSupMod1, iphi1+j, ieta1+i);
493 Int_t absId22 = fEMCALGeo->GetAbsCellIdFromCellIndexes(iSupMod2, iphi2+j, ieta2+i);
494 Bool_t ok1 = kFALSE;
495 Bool_t ok2 = kFALSE;
496 for(Int_t icell = 0; icell < c1->GetNCells(); icell++){
497 if(c1->GetCellsAbsId()[icell] == absId11) ok1=kTRUE;
498 }
499 for(Int_t icell = 0; icell < c2->GetNCells(); icell++){
500 if(c2->GetCellsAbsId()[icell] == absId22) ok2=kTRUE;
501 }
502
503 if(ok1 && (ieta1+i >= 0) && (iphi1+j >= 0) && (ieta1+i < 48) && (iphi1+j < 24))
504 {
505 fHmpi0[iSupMod1][ieta1+i][iphi1+j]->Fill(invmass);
506 }
507 if(ok2 && (ieta2+i >= 0) && (iphi2+j >= 0) && (ieta2+i < 48) && (iphi2+j < 24))
508 {
509 fHmpi0[iSupMod2][ieta2+i][iphi2+j]->Fill(invmass);
510 }
511 }// j loop
512 }//i loop
513 }//group cells
514
515 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",
516 iSupMod1,iphi1,ieta1,iSupMod2,iphi2,ieta2,p12.M(),e1i,c1->E(),e2i,c2->E());
517 }
518
519 }
520
521 } // end of loop over EMCAL clusters
522}
523
524//________________________________________________________________
525void AliAnalysisTaskEMCALPi0CalibSelection::InitGeometryMatrices()
526{
527 // Init geometry and set the geometry matrix, for the first event, skip the rest
528 // Also set once the run dependent calibrations
529
530
531 Int_t runnumber = InputEvent()->GetRunNumber() ;
532
533 if(fLoadMatrices)
534 {
535 printf("AliAnalysisTaskEMCALPi0CalibSelection::InitGeometryMatrices() - Load user defined EMCAL geometry matrices\n");
536
537 // OADB if available
538 AliOADBContainer emcGeoMat("AliEMCALgeo");
539
540 if(fOADBFilePath=="") fOADBFilePath = "$ALICE_ROOT/OADB/EMCAL" ;
541
542 emcGeoMat.InitFromFile(Form("%s/EMCALlocal2master.root",fOADBFilePath.Data()),"AliEMCALgeo");
543
544 TObjArray *matEMCAL=(TObjArray*)emcGeoMat.GetObject(runnumber,"EmcalMatrices");
545
546 for(Int_t mod=0; mod < (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules(); mod++)
547 {
548
549 if (!fMatrix[mod]) // Get it from OADB
550 {
551 if(fDebug > 1 )
552 printf("AliAnalysisTaskEMCALPi0CalibSelection::InitGeometryMatrices() - EMCAL matrices SM %d, %p\n",
553 mod,((TGeoHMatrix*) matEMCAL->At(mod)));
554 //((TGeoHMatrix*) matEMCAL->At(mod))->Print();
555
556 fMatrix[mod] = (TGeoHMatrix*) matEMCAL->At(mod) ;
557 }
558
559 if(fMatrix[mod])
560 {
561 if(DebugLevel() > 1)
562 fMatrix[mod]->Print();
563
564 fEMCALGeo->SetMisalMatrix(fMatrix[mod],mod) ;
565 }
566
567 }//SM loop
568 }//Load matrices
569 else if(!gGeoManager)
570 {
571 printf("AliAnalysisTaskEMCALPi0CalibSelection::InitGeometryMatrices() - Get geo matrices from data");
572 //Still not implemented in AOD, just a workaround to be able to work at least with ESDs
573 if(!strcmp(InputEvent()->GetName(),"AliAODEvent"))
574 {
575 if(DebugLevel() > 1)
576 Warning("UserExec","Use ideal geometry, values geometry matrix not kept in AODs.");
577 }//AOD
578 else
579 {
580 if(DebugLevel() > 1)
581 printf("AliAnalysisTaskEMCALPi0CalibSelection::InitGeometryMatrices() - AliAnalysisTaskEMCALClusterize Load Misaligned matrices.");
582
583 for(Int_t mod=0; mod < (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules(); mod++)
584 {
585 if(DebugLevel() > 1)
586 InputEvent()->GetEMCALMatrix(mod)->Print();
587
588 if(InputEvent()->GetEMCALMatrix(mod)) fEMCALGeo->SetMisalMatrix(InputEvent()->GetEMCALMatrix(mod),mod) ;
589
590 }
591
592 }//ESD
593 }//Load matrices from Data
594
595}
596
597//______________________________________________________________________
598void AliAnalysisTaskEMCALPi0CalibSelection::InitTemperatureCorrections()
599{
600 // Apply run dependent calibration correction
601
602 if(!fRecoUtils->IsRunDepRecalibrationOn()) return;
603
604 AliOADBContainer *contRFTD=new AliOADBContainer("");
605
606 contRFTD->InitFromFile(Form("%s/EMCALTemperatureCorrCalib.root",fOADBFilePath.Data()),"AliEMCALRunDepTempCalibCorrections");
607
608 Int_t runnumber = InputEvent()->GetRunNumber() ;
609
610 TH1S *htd=(TH1S*)contRFTD->GetObject(runnumber);
611
92e834ad 612 //If it did not exist for this run, get closes one
613 if (!htd)
614 {
615 AliWarning(Form("No TemperatureCorrCalib Objects for run: %d",runnumber));
616 // let's get the closest runnumber instead then..
617 Int_t lower = 0;
618 Int_t ic = 0;
619 Int_t maxEntry = contRFTD->GetNumberOfEntries();
620
621 while ( (ic < maxEntry) && (contRFTD->UpperLimit(ic) < runnumber) ) {
622 lower = ic;
623 ic++;
624 }
625
626 Int_t closest = lower;
627 if ( (ic<maxEntry) &&
628 (contRFTD->LowerLimit(ic)-runnumber) < (runnumber - contRFTD->UpperLimit(lower)) ) {
629 closest = ic;
630 }
631
632 AliWarning(Form("TemperatureCorrCalib Objects found closest id %d from run: %d", closest, contRFTD->LowerLimit(closest)));
633 htd = (TH1S*) contRFTD->GetObjectByIndex(closest);
634 }
635
636 // Fill parameters
7b2d541a 637 if(htd)
638 {
639 printf("AliAnalysisTaskEMCALPi0CalibSelection::SetOADBParameters() - Recalibrate (Temperature) EMCAL \n");
640
641 Int_t nSM = fEMCALGeo->GetNumberOfSuperModules();
642
643 for (Int_t ism = 0; ism < nSM; ++ism)
644 {
645 for (Int_t icol = 0; icol < 48; ++icol)
646 {
647 for (Int_t irow = 0; irow < 24; ++irow)
648 {
649 Float_t factor = fRecoUtils->GetEMCALChannelRecalibrationFactor(ism,icol,irow);
650
651 Int_t absID = fEMCALGeo->GetAbsCellIdFromCellIndexes(ism, irow, icol); // original calibration factor
652
653 if(DebugLevel() > 3)
654 printf(" ism %d, icol %d, irow %d,absID %d - Calib factor %1.5f - ",ism, icol, irow, absID, factor);
655
656 factor *= htd->GetBinContent(absID) / 10000. ; // correction dependent on T
657
658 fRecoUtils->SetEMCALChannelRecalibrationFactor(ism,icol,irow,factor);
659
660 if(DebugLevel() > 3)
661 printf(" T factor %1.5f - final factor %1.5f \n",htd->GetBinContent(absID) / 10000.,
662 fRecoUtils->GetEMCALChannelRecalibrationFactor(ism,icol,irow));
663
664 } // columns
665 } // rows
666 } // SM loop
667 }else printf("AliAnalysisTaskEMCALPi0CalibSelection::SetOADBParameters() - Do NOT recalibrate EMCAL with T variations, no params TH1 \n");
668
669}
670
671//___________________________________________________________________
672void AliAnalysisTaskEMCALPi0CalibSelection::UserCreateOutputObjects()
673{
674 //Create output container, init geometry
675
676 fEMCALGeo = AliEMCALGeometry::GetInstance(fEMCALGeoName) ;
677 Int_t nSM = (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules();
678
679 fOutputContainer = new TList();
680 const Int_t buffersize = 255;
681 char hname[buffersize], htitl[buffersize];
682
683 fhNEvents = new TH1I("hNEvents", "Number of analyzed events" , 1 , 0 , 1 ) ;
684 fOutputContainer->Add(fhNEvents);
685
686 fHmgg = new TH2F("hmgg","2-cluster invariant mass",fNbins,fMinBin,fMaxBin,100,0,10);
687 fHmgg->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
688 fHmgg->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
689 fOutputContainer->Add(fHmgg);
690
691 fHmggDifferentSM = new TH2F("hmggDifferentSM","2-cluster invariant mass, different SM",fNbins,fMinBin,fMaxBin,100,0,10);
692 fHmggDifferentSM->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
693 fHmggDifferentSM->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
694 fOutputContainer->Add(fHmggDifferentSM);
695
696 fHOpeningAngle = new TH2F("hopang","2-cluster opening angle",100,0.,50.,100,0,10);
697 fHOpeningAngle->SetXTitle("#alpha_{#gamma #gamma}");
698 fHOpeningAngle->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
699 fOutputContainer->Add(fHOpeningAngle);
700
701 fHOpeningAngleDifferentSM = new TH2F("hopangDifferentSM","2-cluster opening angle, different SM",100,0,50.,100,0,10);
702 fHOpeningAngleDifferentSM->SetXTitle("#alpha_{#gamma #gamma}");
703 fHOpeningAngleDifferentSM->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
704 fOutputContainer->Add(fHOpeningAngleDifferentSM);
705
706 fHAsymmetry = new TH2F("hasym","2-cluster opening angle",100,0.,1.,100,0,10);
707 fHAsymmetry->SetXTitle("a");
708 fHAsymmetry->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
709 fOutputContainer->Add(fHAsymmetry);
710
711 fHAsymmetryDifferentSM = new TH2F("hasymDifferentSM","2-cluster opening angle, different SM",100,0,1.,100,0,10);
712 fHAsymmetryDifferentSM->SetXTitle("a");
713 fHAsymmetryDifferentSM->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
714 fOutputContainer->Add(fHAsymmetryDifferentSM);
715
716
717 //TString pairname[] = {"A side (0-2)", "C side (1-3)","Row 0 (0-1)", "Row 1 (2-3)"};
718
719 fHmggMaskFrame = new TH2F("hmggMaskFrame","2-cluster invariant mass, frame masked",fNbins,fMinBin,fMaxBin,100,0,10);
720 fHmggMaskFrame->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
721 fHmggMaskFrame->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
722 fOutputContainer->Add(fHmggMaskFrame);
723
724 fHmggDifferentSMMaskFrame = new TH2F("hmggDifferentSMMaskFrame","2-cluster invariant mass, different SM, frame masked",
725 fNbins,fMinBin,fMaxBin,100,0,10);
726 fHmggDifferentSMMaskFrame->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
727 fHmggDifferentSMMaskFrame->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
728 fOutputContainer->Add(fHmggDifferentSMMaskFrame);
729
730
731 for(Int_t iSM = 0; iSM < nSM; iSM++)
732 {
733 snprintf(hname, buffersize, "hmgg_SM%d",iSM);
734 snprintf(htitl, buffersize, "Two-gamma inv. mass for super mod %d",iSM);
735 fHmggSM[iSM] = new TH2F(hname,htitl,fNbins,fMinBin,fMaxBin,100,0,10);
736 fHmggSM[iSM]->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
737 fHmggSM[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
738 fOutputContainer->Add(fHmggSM[iSM]);
739
740 snprintf(hname, buffersize, "hmgg_SM%d_MaskFrame",iSM);
741 snprintf(htitl, buffersize, "Two-gamma inv. mass for super mod %d",iSM);
742 fHmggSMMaskFrame[iSM] = new TH2F(hname,htitl,fNbins,fMinBin,fMaxBin,100,0,10);
743 fHmggSMMaskFrame[iSM]->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
744 fHmggSMMaskFrame[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
745 fOutputContainer->Add(fHmggSMMaskFrame[iSM]);
746
747
748 if(iSM < nSM/2)
749 {
750 snprintf(hname,buffersize, "hmgg_PairSameSectorSM%d",iSM);
751 snprintf(htitl,buffersize, "Two-gamma inv. mass for SM pair Sector: %d",iSM);
752 fHmggPairSameSectorSM[iSM] = new TH2F(hname,htitl,fNbins,fMinBin,fMaxBin,100,0,10);
753 fHmggPairSameSectorSM[iSM]->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
754 fHmggPairSameSectorSM[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
755 fOutputContainer->Add(fHmggPairSameSectorSM[iSM]);
756
757 snprintf(hname,buffersize, "hmgg_PairSameSectorSM%d_MaskFrame",iSM);
758 snprintf(htitl,buffersize, "Two-gamma inv. mass for SM pair Sector: %d",iSM);
759 fHmggPairSameSectorSMMaskFrame[iSM] = new TH2F(hname,htitl,fNbins,fMinBin,fMaxBin,100,0,10);
760 fHmggPairSameSectorSMMaskFrame[iSM]->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
761 fHmggPairSameSectorSMMaskFrame[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
762 fOutputContainer->Add(fHmggPairSameSectorSMMaskFrame[iSM]);
763
764 fhClusterPairDiffTimeSameSector[iSM] = new TH2F(Form("hClusterPairDiffTimeSameSector%d",iSM),
765 Form("cluster pair time difference vs E, Sector %d",iSM),
766 100,0,10, 200,-100,100);
767 fhClusterPairDiffTimeSameSector[iSM]->SetXTitle("E_{pair} (GeV)");
768 fhClusterPairDiffTimeSameSector[iSM]->SetYTitle("#Delta t (ns)");
769 fOutputContainer->Add(fhClusterPairDiffTimeSameSector[iSM]);
770
771
772 }
773
774 if(iSM < nSM-2)
775 {
776 snprintf(hname,buffersize, "hmgg_PairSameSideSM%d",iSM);
1dabc151 777 snprintf(htitl,buffersize, "Two-gamma inv. mass for SM pair Sector: %d",iSM);
778 fHmggPairSameSideSM[iSM] = new TH2F(hname,htitl,fNbins,fMinBin,fMaxBin,100,0,10);
779 fHmggPairSameSideSM[iSM]->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
780 fHmggPairSameSideSM[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
781 fOutputContainer->Add(fHmggPairSameSideSM[iSM]);
42b19289 782
783 snprintf(hname,buffersize, "hmgg_PairSameSideSM%d_MaskFrame",iSM);
784 snprintf(htitl,buffersize, "Two-gamma inv. mass for SM pair Sector: %d",iSM);
785 fHmggPairSameSideSMMaskFrame[iSM] = new TH2F(hname,htitl,fNbins,fMinBin,fMaxBin,100,0,10);
786 fHmggPairSameSideSMMaskFrame[iSM]->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
787 fHmggPairSameSideSMMaskFrame[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
af2d7c9b 788 fOutputContainer->Add(fHmggPairSameSideSMMaskFrame[iSM]);
49b53920 789
af2d7c9b 790 fhClusterPairDiffTimeSameSide[iSM] = new TH2F(Form("hClusterPairDiffTimeSameSide%d",iSM),
49b53920 791 Form("cluster pair time difference vs E, Side %d",iSM),
792 100,0,10, 200,-100,100);
af2d7c9b 793 fhClusterPairDiffTimeSameSide[iSM]->SetXTitle("E_{pair} (GeV)");
794 fhClusterPairDiffTimeSameSide[iSM]->SetYTitle("#Delta t (ns)");
795 fOutputContainer->Add(fhClusterPairDiffTimeSameSide[iSM]);
49b53920 796
1dabc151 797 }
9584c261 798
799 snprintf(hname, buffersize, "hopang_SM%d",iSM);
800 snprintf(htitl, buffersize, "Opening angle for super mod %d",iSM);
801 fHOpeningAngleSM[iSM] = new TH2F(hname,htitl,100,0.,50.,100,0,10);
802 fHOpeningAngleSM[iSM]->SetXTitle("#alpha_{#gamma #gamma} (deg)");
803 fHOpeningAngleSM[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
804 fOutputContainer->Add(fHOpeningAngleSM[iSM]);
805
806 snprintf(hname,buffersize, "hopang_PairSM%d",iSM);
1dabc151 807 snprintf(htitl,buffersize, "Opening angle for SM pair: %d",iSM);
9584c261 808 fHOpeningAnglePairSM[iSM] = new TH2F(hname,htitl,100,0.,50.,100,0,10);
809 fHOpeningAnglePairSM[iSM]->SetXTitle("#alpha_{#gamma #gamma} (deg)");
810 fHOpeningAnglePairSM[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
811 fOutputContainer->Add(fHOpeningAnglePairSM[iSM]);
812
9584c261 813 snprintf(hname, buffersize, "hasym_SM%d",iSM);
1dabc151 814 snprintf(htitl, buffersize, "Asymmetry for super mod %d",iSM);
9584c261 815 fHAsymmetrySM[iSM] = new TH2F(hname,htitl,100,0.,1.,100,0,10);
816 fHAsymmetrySM[iSM]->SetXTitle("a");
817 fHAsymmetrySM[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
818 fOutputContainer->Add(fHAsymmetrySM[iSM]);
819
820 snprintf(hname,buffersize, "hasym_PairSM%d",iSM);
1dabc151 821 snprintf(htitl,buffersize, "Asymmetry for SM pair: %d",iSM);
9584c261 822 fHAsymmetryPairSM[iSM] = new TH2F(hname,htitl,100,0.,1.,100,0,10);
823 fHAsymmetryPairSM[iSM]->SetXTitle("a");
824 fHAsymmetryPairSM[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
825 fOutputContainer->Add(fHAsymmetryPairSM[iSM]);
826
9584c261 827 Int_t colmax = 48;
828 Int_t rowmax = 24;
829
af2d7c9b 830 fhTowerDecayPhotonHit[iSM] = new TH2F (Form("hTowerDecPhotonHit_Mod%d",iSM),
49b53920 831 Form("Entries in grid of cells in Module %d",iSM),
832 colmax+2,-1.5,colmax+0.5, rowmax+2,-1.5,rowmax+0.5);
9584c261 833 fhTowerDecayPhotonHit[iSM]->SetYTitle("row (phi direction)");
834 fhTowerDecayPhotonHit[iSM]->SetXTitle("column (eta direction)");
835 fOutputContainer->Add(fhTowerDecayPhotonHit[iSM]);
836
af2d7c9b 837 fhTowerDecayPhotonEnergy[iSM] = new TH2F (Form("hTowerDecPhotonEnergy_Mod%d",iSM),
49b53920 838 Form("Accumulated energy in grid of cells in Module %d",iSM),
839 colmax+2,-1.5,colmax+0.5, rowmax+2,-1.5,rowmax+0.5);
9584c261 840 fhTowerDecayPhotonEnergy[iSM]->SetYTitle("row (phi direction)");
841 fhTowerDecayPhotonEnergy[iSM]->SetXTitle("column (eta direction)");
842 fOutputContainer->Add(fhTowerDecayPhotonEnergy[iSM]);
843
af2d7c9b 844 fhTowerDecayPhotonAsymmetry[iSM] = new TH2F (Form("hTowerDecPhotonAsymmetry_Mod%d",iSM),
49b53920 845 Form("Accumulated asymmetry in grid of cells in Module %d",iSM),
846 colmax+2,-1.5,colmax+0.5, rowmax+2,-1.5,rowmax+0.5);
9584c261 847 fhTowerDecayPhotonAsymmetry[iSM]->SetYTitle("row (phi direction)");
848 fhTowerDecayPhotonAsymmetry[iSM]->SetXTitle("column (eta direction)");
849 fOutputContainer->Add(fhTowerDecayPhotonAsymmetry[iSM]);
850
42b19289 851 fhTowerDecayPhotonHitMaskFrame[iSM] = new TH2F (Form("hTowerDecPhotonHit_Mod%d_MaskFrame",iSM),Form("Entries in grid of cells in Module %d",iSM),
49b53920 852 colmax+2,-1.5,colmax+0.5, rowmax+2,-1.5,rowmax+0.5);
42b19289 853 fhTowerDecayPhotonHitMaskFrame[iSM]->SetYTitle("row (phi direction)");
854 fhTowerDecayPhotonHitMaskFrame[iSM]->SetXTitle("column (eta direction)");
855 fOutputContainer->Add(fhTowerDecayPhotonHitMaskFrame[iSM]);
49b53920 856
af2d7c9b 857 fhClusterTimeSM[iSM] = new TH2F(Form("hClusterTime_SM%d",iSM),"cluster time vs E",100,0,10, 100,0,1000);
858 fhClusterTimeSM[iSM]->SetXTitle("E (GeV)");
859 fhClusterTimeSM[iSM]->SetYTitle("t (ns)");
860 fOutputContainer->Add(fhClusterTimeSM[iSM]);
42b19289 861
af2d7c9b 862 fhClusterPairDiffTimeSameSM[iSM] = new TH2F(Form("hClusterPairDiffTimeSameSM%d",iSM),
49b53920 863 Form("cluster pair time difference vs E, SM %d",iSM),
864 100,0,10, 200,-100,100);
af2d7c9b 865 fhClusterPairDiffTimeSameSM[iSM]->SetXTitle("E (GeV)");
866 fhClusterPairDiffTimeSameSM[iSM]->SetYTitle("#Delta t (ns)");
867 fOutputContainer->Add(fhClusterPairDiffTimeSameSM[iSM]);
49b53920 868
2dfb1428 869 }
6eb2a715 870
a8dc7d71 871 Int_t nchannels = nSM*AliEMCALGeoParams::fgkEMCALRows*AliEMCALGeoParams::fgkEMCALCols;
872 for(Int_t ibc = 0; ibc < 4; ibc++)
873 {
874 fHTpi0[ibc] = new TH2F(Form("hTime_BC%d",ibc),Form("Time of cell clusters under pi0 peak, bunch crossing %d",ibc),
875 nchannels,0,nchannels, fNTimeBins,fMinTimeBin,fMaxTimeBin);
876 fOutputContainer->Add(fHTpi0[ibc]);
877 fHTpi0[ibc]->SetYTitle("time (ns)");
878 fHTpi0[ibc]->SetXTitle("abs. Id. ");
879 }
880
881
af2d7c9b 882 fhClusterTime = new TH2F("hClusterTime","cluster time vs E",100,0,10, 100,0,1000);
883 fhClusterTime->SetXTitle("E (GeV)");
884 fhClusterTime->SetYTitle("t (ns)");
885 fOutputContainer->Add(fhClusterTime);
49b53920 886
a7e5a381 887 fhClusterPairDiffTime = new TH2F("hClusterPairDiffTime","cluster pair time difference vs E",100,0,10, 800,-400,400);
af2d7c9b 888 fhClusterPairDiffTime->SetXTitle("E_{pair} (GeV)");
889 fhClusterPairDiffTime->SetYTitle("#Delta t (ns)");
890 fOutputContainer->Add(fhClusterPairDiffTime);
49b53920 891
a8dc7d71 892 for(Int_t iMod=0; iMod < nSM; iMod++)
893 {
894 for(Int_t iRow=0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++)
895 {
896 for(Int_t iCol=0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++)
897 {
898 snprintf(hname,buffersize, "%d_%d_%d",iMod,iCol,iRow);
899 snprintf(htitl,buffersize, "Two-gamma inv. mass for super mod %d, cell(col,row)=(%d,%d)",iMod,iCol,iRow);
900 fHmpi0[iMod][iCol][iRow] = new TH1F(hname,htitl,fNbins,fMinBin,fMaxBin);
901 fHmpi0[iMod][iCol][iRow]->SetXTitle("mass (MeV/c^{2})");
902 fOutputContainer->Add(fHmpi0[iMod][iCol][iRow]);
903 }
904 }
905 }
247abff4 906
907 fOutputContainer->SetOwner(kTRUE);
a8dc7d71 908
cf028690 909 PostData(1,fOutputContainer);
49b53920 910
477c5cd2 911 // cuts container, set in terminate but init and post here
912
913 fCuts = new TList();
914
915 fCuts ->SetOwner(kTRUE);
916
917 PostData(2, fCuts);
918
375cec9b 919}
920
477c5cd2 921//______________________________________________________________________________________________________
922Bool_t AliAnalysisTaskEMCALPi0CalibSelection::MaskFrameCluster(const Int_t iSM, const Int_t ieta) const
923{
af2d7c9b 924 //Check if cell is in one of the regions where we have significant amount of material in front of EMCAL
42b19289 925
926 Int_t icol = ieta;
927 if(iSM%2) icol+=48; // Impair SM, shift index [0-47] to [48-96]
928
477c5cd2 929 if (fNMaskCellColumns && fMaskCellColumns)
930 {
931 for (Int_t imask = 0; imask < fNMaskCellColumns; imask++)
932 {
42b19289 933 if(icol==fMaskCellColumns[imask]) return kTRUE;
934 }
935 }
af2d7c9b 936
42b19289 937 return kFALSE;
938
939}
940
477c5cd2 941//__________________________________________________________________________
375cec9b 942void AliAnalysisTaskEMCALPi0CalibSelection::UserExec(Option_t* /* option */)
943{
7b2d541a 944 // Do analysis, first select the events, then correct the clusters if needed
945 // and finally fill the histograms per channel after recalibration
247abff4 946
7b2d541a 947 //Event selection
a8dc7d71 948 if(fTriggerName!="")
477c5cd2 949 {
7b2d541a 950 AliESDEvent* esdevent = dynamic_cast<AliESDEvent*> (InputEvent());
951 AliAODEvent* aodevent = dynamic_cast<AliAODEvent*> (InputEvent());
477c5cd2 952
7b2d541a 953 TString triggerClass = "";
954 if (esdevent) triggerClass = esdevent->GetFiredTriggerClasses();
955 else if(aodevent) triggerClass = aodevent->GetFiredTriggerClasses();
42b19289 956
ccd9df97 957 if(DebugLevel() > 1)
958 printf("AliAnalysisTaskEMCALPi0CalibSelection::UserExec() - Event %d, FiredClass %s",
959 (Int_t)Entry(),(((AliESDEvent*)InputEvent())->GetFiredTriggerClasses()).Data());
960
961 if(!triggerClass.Contains(fTriggerName))
477c5cd2 962 {
ccd9df97 963 if(DebugLevel() > 1) printf("Reject event! \n");
964 return;
965 }
966 else
967 if(DebugLevel() > 1) printf("Accept Event! \n");
7b2d541a 968 }
ccd9df97 969
7b2d541a 970 //Get the input event
971 AliVEvent* event = 0;
972 if(fFilteredInput) event = AODEvent();
973 else event = InputEvent();
974
975 if(!event)
976 {
977 printf("Input event not available!\n");
978 return;
979 }
980
981 if(DebugLevel() > 1)
982 printf("AliAnalysisTaskEMCALPi0CalibSelection <<< %s: Event %d >>>\n",event->GetName(), (Int_t)Entry());
983
984 //Get the primary vertex
985 event->GetPrimaryVertex()->GetXYZ(fVertex) ;
986
987 if(DebugLevel() > 1) printf("AliAnalysisTaskEMCALPi0CalibSelection Vertex: (%.3f,%.3f,%.3f)\n",fVertex[0],fVertex[1],fVertex[2]);
988
989 //Int_t runNum = aod->GetRunNumber();
990 //if(DebugLevel() > 1) printf("Run number: %d\n",runNum);
991
992 fhNEvents->Fill(0); //Count the events to be analyzed
993
994 // Acccess once the geometry matrix and temperature corrections
995 if(fhNEvents->GetEntries()==1)
996 {
997 InitGeometryMatrices();
998
999 InitTemperatureCorrections();
1000 }
1001
1002 //Get the list of clusters and cells
1003 fEMCALCells = event->GetEMCALCells();
1004
1005 fCaloClustersArr = new TRefArray();
1006 event->GetEMCALClusters(fCaloClustersArr);
375cec9b 1007
7b2d541a 1008 if(DebugLevel() > 1) printf("AliAnalysisTaskEMCALPi0CalibSelection - N CaloClusters: %d - N CaloCells %d \n",
1009 fCaloClustersArr->GetEntriesFast(), fEMCALCells->GetNumberOfCells());
1010
1011 CorrectClusters(); // Non linearity, new calibration, T calibration
1012
1013 FillHistograms();
1014
1015 delete fCaloClustersArr;
6eb2a715 1016
375cec9b 1017 PostData(1,fOutputContainer);
6eb2a715 1018
375cec9b 1019}
cfce8d44 1020
5ef94e1b 1021//_____________________________________________________
477c5cd2 1022void AliAnalysisTaskEMCALPi0CalibSelection::PrintInfo()
1023{
5ef94e1b 1024 //Print settings
477c5cd2 1025
1026 printf("Cluster cuts: %2.2f < E < %2.2f GeV; number of cells > %d; Assymetry < %1.2f, pair time diff < %2.2f, %2.2f < t < %2.2f ns\n",
1027 fEmin,fEmax, fMinNCells, fAsyCut, fDTimeCut,fTimeMin,fTimeMax) ;
1028
af2d7c9b 1029 printf("Group %d cells\n", fGroupNCells) ;
477c5cd2 1030
5ef94e1b 1031 printf("Cluster maximal cell away from border at least %d cells\n", fRecoUtils->GetNumberOfCellsFromEMCALBorder()) ;
477c5cd2 1032
af2d7c9b 1033 printf("Histograms: bins %d; energy range: %2.2f < E < %2.2f GeV\n",fNbins,fMinBin,fMaxBin) ;
477c5cd2 1034
af2d7c9b 1035 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 1036 fRecoUtils->IsBadChannelsRemovalSwitchedOn(),fFilteredInput,fCorrectClusters, fSameSM) ;
477c5cd2 1037
af2d7c9b 1038 printf("EMCAL Geometry name: < %s >, Load Matrices %d\n",fEMCALGeoName.Data(), fLoadMatrices) ;
a8dc7d71 1039 if(fLoadMatrices) {for(Int_t ism = 0; ism < AliEMCALGeoParams::fgkEMCALModules; ism++) if(fMatrix[ism]) fMatrix[ism]->Print() ; }
af2d7c9b 1040
5ef94e1b 1041}
1042
477c5cd2 1043//____________________________________________________________________
1044void AliAnalysisTaskEMCALPi0CalibSelection::Terminate(Option_t*)
1045{
1046 // Create cuts/param objects and publish to slot
1047 const Int_t buffersize = 255;
1048 char onePar[buffersize] ;
1049
1050 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; %2.2f < T < %2.2f ns; %3.1f < Mass < %3.1f",
1051 fEmin,fEmax, fL0min, fL0max, fMinNCells, fAsyCut, fDTimeCut, fTimeMin, fTimeMax, fInvMassCutMin, fInvMassCutMax) ;
1052 fCuts->Add(new TObjString(onePar));
1053 snprintf(onePar,buffersize, "Group %d cells;", fGroupNCells) ;
1054 fCuts->Add(new TObjString(onePar));
1055 snprintf(onePar,buffersize, "Cluster maximal cell away from border at least %d cells;", fRecoUtils->GetNumberOfCellsFromEMCALBorder()) ;
1056 fCuts->Add(new TObjString(onePar));
1057 snprintf(onePar,buffersize, "Histograms, Mass bins %d; energy range: %2.2f < E < %2.2f GeV;",fNbins,fMinBin,fMaxBin) ;
1058 fCuts->Add(new TObjString(onePar));
1059 snprintf(onePar,buffersize, "Histograms, Time bins %d; energy range: %2.2f < E < %2.2f GeV;",fNTimeBins,fMinTimeBin,fMaxTimeBin) ;
1060 fCuts->Add(new TObjString(onePar));
1061 snprintf(onePar,buffersize, "Switchs: Remove Bad Channels? %d; Use filtered input? %d; Correct Clusters? %d, Mass per channel same SM clusters? %d ",
1062 fRecoUtils->IsBadChannelsRemovalSwitchedOn(),fFilteredInput,fCorrectClusters, fSameSM) ;
1063 fCuts->Add(new TObjString(onePar));
1064 snprintf(onePar,buffersize, "EMCAL Geometry name: < %s >, Load Matrices? %d",fEMCALGeoName.Data(),fLoadMatrices) ;
1065 fCuts->Add(new TObjString(onePar));
1066
1067 // Post Data
1068 PostData(2, fCuts);
1069
1070}
1071