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