bfc7c23b |
1 | /************************************************************************** |
2 | * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * |
3 | * * |
4 | * Author: The ALICE Off-line Project. * |
5 | * Contributors are mentioned in the code where appropriate. * |
6 | * * |
7 | * Permission to use, copy, modify and distribute this software and its * |
8 | * documentation strictly for non-commercial purposes is hereby granted * |
9 | * without fee, provided that the above copyright notice appears in all * |
10 | * copies and that both the copyright notice and this permission notice * |
11 | * appear in the supporting documentation. The authors make no claims * |
12 | * about the suitability of this software for any purpose. It is * |
13 | * provided "as is" without express or implied warranty. * |
14 | **************************************************************************/ |
15 | |
16 | // Class for heavy-flavour electron with EMCal triggered events |
17 | // Author: Shingo Sakai |
18 | |
19 | |
20 | #include "TChain.h" |
21 | #include "TTree.h" |
22 | #include "TH2F.h" |
23 | #include "TMath.h" |
24 | #include "TCanvas.h" |
25 | #include "THnSparse.h" |
26 | #include "TLorentzVector.h" |
27 | #include "TString.h" |
28 | #include "TFile.h" |
29 | |
30 | #include "AliAnalysisTask.h" |
31 | #include "AliAnalysisManager.h" |
32 | |
33 | #include "AliESDEvent.h" |
34 | #include "AliESDHandler.h" |
35 | #include "AliAODEvent.h" |
36 | #include "AliAODHandler.h" |
37 | |
e7d87aef |
38 | #include "AliMCEventHandler.h" |
39 | #include "AliMCEvent.h" |
40 | #include "AliMCParticle.h" |
41 | |
bfc7c23b |
42 | #include "AliAnalysisTaskHFECal.h" |
43 | #include "TGeoGlobalMagField.h" |
44 | #include "AliLog.h" |
45 | #include "AliAnalysisTaskSE.h" |
46 | #include "TRefArray.h" |
47 | #include "TVector.h" |
48 | #include "AliESDInputHandler.h" |
49 | #include "AliESDpid.h" |
50 | #include "AliESDtrackCuts.h" |
51 | #include "AliPhysicsSelection.h" |
52 | #include "AliESDCaloCluster.h" |
53 | #include "AliAODCaloCluster.h" |
54 | #include "AliEMCALRecoUtils.h" |
55 | #include "AliEMCALGeometry.h" |
56 | #include "AliGeomManager.h" |
57 | #include "stdio.h" |
58 | #include "TGeoManager.h" |
59 | #include "iostream" |
60 | #include "fstream" |
61 | |
62 | #include "AliEMCALTrack.h" |
63 | #include "AliMagF.h" |
64 | |
65 | #include "AliKFParticle.h" |
66 | #include "AliKFVertex.h" |
67 | |
68 | #include "AliPID.h" |
69 | #include "AliPIDResponse.h" |
70 | #include "AliHFEcontainer.h" |
71 | #include "AliHFEcuts.h" |
72 | #include "AliHFEpid.h" |
73 | #include "AliHFEpidBase.h" |
74 | #include "AliHFEpidQAmanager.h" |
75 | #include "AliHFEtools.h" |
76 | #include "AliCFContainer.h" |
77 | #include "AliCFManager.h" |
78 | |
e7d87aef |
79 | #include "AliStack.h" |
80 | |
bfc7c23b |
81 | #include "AliCentrality.h" |
82 | |
83 | ClassImp(AliAnalysisTaskHFECal) |
84 | //________________________________________________________________________ |
85 | AliAnalysisTaskHFECal::AliAnalysisTaskHFECal(const char *name) |
86 | : AliAnalysisTaskSE(name) |
87 | ,fESD(0) |
e7d87aef |
88 | ,fMC(0) |
fb87d707 |
89 | ,fGeom(0) |
bfc7c23b |
90 | ,fOutputList(0) |
91 | ,fTrackCuts(0) |
92 | ,fCuts(0) |
93 | ,fIdentifiedAsOutInz(kFALSE) |
94 | ,fPassTheEventCut(kFALSE) |
95 | ,fRejectKinkMother(kFALSE) |
e7d87aef |
96 | ,fmcData(kFALSE) |
bfc7c23b |
97 | ,fVz(0.0) |
98 | ,fCFM(0) |
99 | ,fPID(0) |
100 | ,fPIDqa(0) |
101 | ,fOpeningAngleCut(0.1) |
102 | ,fInvmassCut(0.01) |
103 | ,fNoEvents(0) |
104 | ,fEMCAccE(0) |
105 | ,fTrkpt(0) |
106 | ,fTrkEovPBef(0) |
107 | ,fTrkEovPAft(0) |
108 | ,fdEdxBef(0) |
109 | ,fdEdxAft(0) |
110 | ,fIncpT(0) |
fb87d707 |
111 | ,fIncpTM20(0) |
bfc7c23b |
112 | ,fInvmassLS(0) |
113 | ,fInvmassULS(0) |
114 | ,fOpeningAngleLS(0) |
115 | ,fOpeningAngleULS(0) |
116 | ,fPhotoElecPt(0) |
117 | ,fPhoElecPt(0) |
fb87d707 |
118 | ,fPhoElecPtM20(0) |
f09766dd |
119 | ,fSameElecPt(0) |
fb87d707 |
120 | ,fSameElecPtM20(0) |
bfc7c23b |
121 | ,fTrackPtBefTrkCuts(0) |
122 | ,fTrackPtAftTrkCuts(0) |
123 | ,fTPCnsigma(0) |
124 | ,fCent(0) |
125 | ,fEleInfo(0) |
fb87d707 |
126 | ,fClsEBftTrigCut(0) |
127 | ,fClsEAftTrigCut(0) |
128 | ,fClsEAftTrigCut1(0) |
129 | ,fClsEAftTrigCut2(0) |
130 | ,fClsEAftTrigCut3(0) |
131 | ,fClsEAftTrigCut4(0) |
132 | ,fClsETime(0) |
133 | ,fClsETime1(0) |
134 | ,fTrigTimes(0) |
42c75692 |
135 | ,fCellCheck(0) |
e7d87aef |
136 | ,fInputHFEMC(0) |
137 | ,fIncpTMChfe(0) |
138 | ,fIncpTMCM20hfe(0) |
139 | ,fIncpTMCpho(0) |
140 | ,fIncpTMCM20pho(0) |
141 | ,fPhoElecPtMC(0) |
142 | ,fPhoElecPtMCM20(0) |
143 | ,fSameElecPtMC(0) |
144 | ,fSameElecPtMCM20(0) |
bfc7c23b |
145 | { |
146 | //Named constructor |
147 | |
148 | fPID = new AliHFEpid("hfePid"); |
149 | fTrackCuts = new AliESDtrackCuts(); |
150 | |
151 | // Define input and output slots here |
152 | // Input slot #0 works with a TChain |
153 | DefineInput(0, TChain::Class()); |
154 | // Output slot #0 id reserved by the base class for AOD |
155 | // Output slot #1 writes into a TH1 container |
156 | // DefineOutput(1, TH1I::Class()); |
157 | DefineOutput(1, TList::Class()); |
158 | // DefineOutput(3, TTree::Class()); |
159 | } |
160 | |
161 | //________________________________________________________________________ |
162 | AliAnalysisTaskHFECal::AliAnalysisTaskHFECal() |
163 | : AliAnalysisTaskSE("DefaultAnalysis_AliAnalysisTaskHFECal") |
164 | ,fESD(0) |
e7d87aef |
165 | ,fMC(0) |
fb87d707 |
166 | ,fGeom(0) |
bfc7c23b |
167 | ,fOutputList(0) |
168 | ,fTrackCuts(0) |
169 | ,fCuts(0) |
170 | ,fIdentifiedAsOutInz(kFALSE) |
171 | ,fPassTheEventCut(kFALSE) |
172 | ,fRejectKinkMother(kFALSE) |
e7d87aef |
173 | ,fmcData(kFALSE) |
bfc7c23b |
174 | ,fVz(0.0) |
175 | ,fCFM(0) |
176 | ,fPID(0) |
177 | ,fPIDqa(0) |
178 | ,fOpeningAngleCut(0.1) |
179 | ,fInvmassCut(0.01) |
180 | ,fNoEvents(0) |
181 | ,fEMCAccE(0) |
182 | ,fTrkpt(0) |
183 | ,fTrkEovPBef(0) |
184 | ,fTrkEovPAft(0) |
185 | ,fdEdxBef(0) |
186 | ,fdEdxAft(0) |
187 | ,fIncpT(0) |
fb87d707 |
188 | ,fIncpTM20(0) |
bfc7c23b |
189 | ,fInvmassLS(0) |
190 | ,fInvmassULS(0) |
191 | ,fOpeningAngleLS(0) |
192 | ,fOpeningAngleULS(0) |
193 | ,fPhotoElecPt(0) |
194 | ,fPhoElecPt(0) |
fb87d707 |
195 | ,fPhoElecPtM20(0) |
f09766dd |
196 | ,fSameElecPt(0) |
fb87d707 |
197 | ,fSameElecPtM20(0) |
bfc7c23b |
198 | ,fTrackPtBefTrkCuts(0) |
199 | ,fTrackPtAftTrkCuts(0) |
200 | ,fTPCnsigma(0) |
201 | ,fCent(0) |
202 | ,fEleInfo(0) |
fb87d707 |
203 | ,fClsEBftTrigCut(0) |
204 | ,fClsEAftTrigCut(0) |
205 | ,fClsEAftTrigCut1(0) |
206 | ,fClsEAftTrigCut2(0) |
207 | ,fClsEAftTrigCut3(0) |
208 | ,fClsEAftTrigCut4(0) |
209 | ,fClsETime(0) |
210 | ,fClsETime1(0) |
211 | ,fTrigTimes(0) |
42c75692 |
212 | ,fCellCheck(0) |
e7d87aef |
213 | ,fInputHFEMC(0) |
214 | ,fIncpTMChfe(0) |
215 | ,fIncpTMCM20hfe(0) |
216 | ,fIncpTMCpho(0) |
217 | ,fIncpTMCM20pho(0) |
218 | ,fPhoElecPtMC(0) |
219 | ,fPhoElecPtMCM20(0) |
220 | ,fSameElecPtMC(0) |
221 | ,fSameElecPtMCM20(0) |
bfc7c23b |
222 | { |
223 | //Default constructor |
224 | fPID = new AliHFEpid("hfePid"); |
225 | |
226 | fTrackCuts = new AliESDtrackCuts(); |
227 | |
228 | // Constructor |
229 | // Define input and output slots here |
230 | // Input slot #0 works with a TChain |
231 | DefineInput(0, TChain::Class()); |
232 | // Output slot #0 id reserved by the base class for AOD |
233 | // Output slot #1 writes into a TH1 container |
234 | // DefineOutput(1, TH1I::Class()); |
235 | DefineOutput(1, TList::Class()); |
236 | //DefineOutput(3, TTree::Class()); |
237 | } |
238 | //_________________________________________ |
239 | |
240 | AliAnalysisTaskHFECal::~AliAnalysisTaskHFECal() |
241 | { |
242 | //Destructor |
243 | |
244 | delete fOutputList; |
fb87d707 |
245 | delete fGeom; |
bfc7c23b |
246 | delete fPID; |
247 | delete fCFM; |
248 | delete fPIDqa; |
249 | delete fTrackCuts; |
250 | } |
251 | //_________________________________________ |
252 | |
253 | void AliAnalysisTaskHFECal::UserExec(Option_t*) |
254 | { |
255 | //Main loop |
256 | //Called for each event |
257 | |
258 | // create pointer to event |
259 | fESD = dynamic_cast<AliESDEvent*>(InputEvent()); |
260 | if (!fESD) { |
261 | printf("ERROR: fESD not available\n"); |
262 | return; |
263 | } |
264 | |
265 | if(!fCuts){ |
266 | AliError("HFE cuts not available"); |
267 | return; |
268 | } |
269 | |
270 | if(!fPID->IsInitialized()){ |
271 | // Initialize PID with the given run number |
272 | AliWarning("PID not initialised, get from Run no"); |
273 | fPID->InitializePID(fESD->GetRunNumber()); |
274 | } |
275 | |
e7d87aef |
276 | if(fmcData)fMC = MCEvent(); |
277 | AliStack* stack = NULL; |
278 | if(fmcData && fMC)stack = fMC->Stack(); |
279 | |
280 | Float_t cent = -1.; |
281 | AliCentrality *centrality = fESD->GetCentrality(); |
282 | cent = centrality->GetCentralityPercentile("V0M"); |
283 | |
284 | //---- fill MC track info |
285 | if(fmcData && fMC) |
286 | { |
287 | Int_t nParticles = stack->GetNtrack(); |
288 | for (Int_t iParticle = 0; iParticle < nParticles; iParticle++) { |
289 | TParticle* particle = stack->Particle(iParticle); |
290 | int fPDG = particle->GetPdgCode(); |
291 | double mcZvertex = fMC->GetPrimaryVertex()->GetZ(); |
292 | double pTMC = particle->Pt(); |
293 | Bool_t mcInDtoE= kFALSE; |
294 | Bool_t mcInBtoE= kFALSE; |
295 | |
296 | if(particle->GetFirstMother()>-1 && fabs(fPDG)==11) |
297 | { |
298 | int parentPID = stack->Particle(particle->GetFirstMother())->GetPdgCode(); |
299 | if((fabs(parentPID)==411 || fabs(parentPID)==413 || fabs(parentPID)==421 || fabs(parentPID)==423 || fabs(parentPID)==431)&& fabs(fPDG)==11)mcInDtoE = kTRUE; |
300 | if((fabs(parentPID)==511 || fabs(parentPID)==513 || fabs(parentPID)==521 || fabs(parentPID)==523 || fabs(parentPID)==531)&& fabs(fPDG)==11)mcInBtoE = kTRUE; |
301 | if((mcInBtoE || mcInDtoE) && fabs(mcZvertex)<10.0)fInputHFEMC->Fill(cent,pTMC); |
302 | } |
303 | |
304 | } |
305 | } |
306 | |
bfc7c23b |
307 | fNoEvents->Fill(0); |
308 | |
309 | Int_t fNOtrks = fESD->GetNumberOfTracks(); |
310 | const AliESDVertex *pVtx = fESD->GetPrimaryVertex(); |
311 | |
312 | Double_t pVtxZ = -999; |
313 | pVtxZ = pVtx->GetZ(); |
314 | |
315 | if(TMath::Abs(pVtxZ)>10) return; |
316 | fNoEvents->Fill(1); |
317 | |
318 | if(fNOtrks<2) return; |
319 | fNoEvents->Fill(2); |
320 | |
321 | AliPIDResponse *pidResponse = fInputHandler->GetPIDResponse(); |
322 | if(!pidResponse){ |
323 | AliDebug(1, "Using default PID Response"); |
324 | pidResponse = AliHFEtools::GetDefaultPID(kFALSE, fInputEvent->IsA() == AliAODEvent::Class()); |
325 | } |
326 | |
327 | fPID->SetPIDResponse(pidResponse); |
328 | |
329 | fCFM->SetRecEventInfo(fESD); |
330 | |
e7d87aef |
331 | //Float_t cent = -1.; |
332 | //AliCentrality *centrality = fESD->GetCentrality(); |
333 | //cent = centrality->GetCentralityPercentile("V0M"); |
bfc7c23b |
334 | fCent->Fill(cent); |
335 | |
336 | if(cent>90.) return; |
337 | |
338 | // Calorimeter info. |
339 | |
42c75692 |
340 | FindTriggerClusters(); |
341 | |
bfc7c23b |
342 | // make EMCAL array |
343 | for(Int_t iCluster=0; iCluster<fESD->GetNumberOfCaloClusters(); iCluster++) |
344 | { |
345 | AliESDCaloCluster *clust = fESD->GetCaloCluster(iCluster); |
346 | if(clust->IsEMCAL()) |
347 | { |
348 | double clustE = clust->E(); |
349 | float emcx[3]; // cluster pos |
350 | clust->GetPosition(emcx); |
351 | TVector3 clustpos(emcx[0],emcx[1],emcx[2]); |
352 | double emcphi = clustpos.Phi(); |
353 | double emceta = clustpos.Eta(); |
fb87d707 |
354 | double calInfo[5]; |
355 | calInfo[0] = emcphi; calInfo[1] = emceta; calInfo[2] = clustE; calInfo[3] = cent; calInfo[4] = clust->Chi2(); |
f09766dd |
356 | if(clustE>1.5)fEMCAccE->Fill(calInfo); |
bfc7c23b |
357 | } |
358 | } |
359 | |
360 | // Track loop |
361 | for (Int_t iTracks = 0; iTracks < fESD->GetNumberOfTracks(); iTracks++) { |
362 | AliESDtrack* track = fESD->GetTrack(iTracks); |
363 | if (!track) { |
364 | printf("ERROR: Could not receive track %d\n", iTracks); |
365 | continue; |
366 | } |
e7d87aef |
367 | |
368 | Bool_t mcPho = kFALSE; |
369 | Bool_t mcDtoE= kFALSE; |
370 | Bool_t mcBtoE= kFALSE; |
371 | double mcele = -1.; |
372 | if(fmcData && fMC && stack) |
373 | { |
374 | Int_t label = track->GetLabel(); |
375 | TParticle* particle = stack->Particle(label); |
376 | int mcpid = particle->GetPdgCode(); |
377 | printf("MCpid = %d",mcpid); |
378 | if(particle->GetFirstMother()>-1) |
379 | { |
380 | int parentPID = stack->Particle(particle->GetFirstMother())->GetPdgCode(); |
381 | if((parentPID==22 || parentPID==111 || parentPID==221)&& fabs(mcpid)==11)mcPho = kTRUE; |
382 | if((fabs(parentPID)==411 || fabs(parentPID)==413 || fabs(parentPID)==421 || fabs(parentPID)==423 || fabs(parentPID)==431)&& fabs(mcpid)==11)mcDtoE = kTRUE; |
383 | if((fabs(parentPID)==511 || fabs(parentPID)==513 || fabs(parentPID)==521 || fabs(parentPID)==523 || fabs(parentPID)==531)&& fabs(mcpid)==11)mcBtoE = kTRUE; |
384 | } |
385 | if(fabs(mcpid)==11)mcele= 1.; |
386 | } |
387 | |
bfc7c23b |
388 | if(TMath::Abs(track->Eta())>0.7) continue; |
f09766dd |
389 | if(TMath::Abs(track->Pt()<2.0)) continue; |
bfc7c23b |
390 | |
391 | fTrackPtBefTrkCuts->Fill(track->Pt()); |
392 | // RecKine: ITSTPC cuts |
393 | if(!ProcessCutStep(AliHFEcuts::kStepRecKineITSTPC, track)) continue; |
394 | |
395 | //RecKink |
396 | if(fRejectKinkMother) { // Quick and dirty fix to reject both kink mothers and daughters |
397 | if(track->GetKinkIndex(0) != 0) continue; |
398 | } |
399 | |
400 | // RecPrim |
401 | if(!ProcessCutStep(AliHFEcuts::kStepRecPrim, track)) continue; |
402 | |
403 | // HFEcuts: ITS layers cuts |
404 | if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsITS, track)) continue; |
405 | |
f09766dd |
406 | // HFE cuts: TPC PID cleanup |
407 | if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTPC, track)) continue; |
408 | |
bfc7c23b |
409 | |
410 | fTrackPtAftTrkCuts->Fill(track->Pt()); |
411 | |
412 | Double_t mom = -999., eop=-999., pt = -999., dEdx=-999., fTPCnSigma=-10, phi=-999., eta=-999.; |
fb87d707 |
413 | //pt = track->Pt(); |
414 | //if(pt<2.0)continue; |
bfc7c23b |
415 | |
416 | // Track extrapolation |
417 | |
bfc7c23b |
418 | Int_t charge = track->Charge(); |
419 | fTrkpt->Fill(pt); |
420 | mom = track->P(); |
421 | phi = track->Phi(); |
422 | eta = track->Eta(); |
423 | dEdx = track->GetTPCsignal(); |
424 | fTPCnSigma = fPID->GetPIDResponse() ? fPID->GetPIDResponse()->NumberOfSigmasTPC(track, AliPID::kElectron) : 1000; |
425 | |
426 | double ncells = -1.0; |
427 | double m20 = -1.0; |
428 | double m02 = -1.0; |
429 | double disp = -1.0; |
430 | double rmatch = -1.0; |
431 | double nmatch = -1.0; |
f09766dd |
432 | double oppstatus = 0.0; |
433 | double samestatus = 0.0; |
bfc7c23b |
434 | |
f09766dd |
435 | Bool_t fFlagPhotonicElec = kFALSE; |
436 | Bool_t fFlagConvinatElec = kFALSE; |
437 | SelectPhotonicElectron(iTracks,cent,track,fFlagPhotonicElec,fFlagConvinatElec); |
438 | if(fFlagPhotonicElec)oppstatus = 1.0; |
439 | if(fFlagConvinatElec)samestatus = 1.0; |
bfc7c23b |
440 | |
441 | Int_t clsId = track->GetEMCALcluster(); |
442 | if (clsId>0){ |
443 | AliESDCaloCluster *clust = fESD->GetCaloCluster(clsId); |
444 | if(clust && clust->IsEMCAL()){ |
445 | |
446 | double clustE = clust->E(); |
447 | eop = clustE/fabs(mom); |
448 | //double clustT = clust->GetTOF(); |
449 | ncells = clust->GetNCells(); |
450 | m02 = clust->GetM02(); |
451 | m20 = clust->GetM20(); |
452 | disp = clust->GetDispersion(); |
453 | double delphi = clust->GetTrackDx(); |
454 | double deleta = clust->GetTrackDz(); |
455 | rmatch = sqrt(pow(delphi,2)+pow(deleta,2)); |
456 | nmatch = clust->GetNTracksMatched(); |
457 | |
e7d87aef |
458 | double valdedx[17]; |
fb87d707 |
459 | valdedx[0] = pt; valdedx[1] = dEdx; valdedx[2] = phi; valdedx[3] = eta; valdedx[4] = fTPCnSigma; |
460 | valdedx[5] = eop; valdedx[6] = rmatch; valdedx[7] = ncells, valdedx[8] = m02; valdedx[9] = m20; valdedx[10] = disp; |
461 | valdedx[11] = cent; valdedx[12] = charge; valdedx[13] = oppstatus; valdedx[14] = samestatus; valdedx[15] = clust->Chi2(); |
e7d87aef |
462 | valdedx[16] = mcele; |
bfc7c23b |
463 | fEleInfo->Fill(valdedx); |
464 | |
465 | |
466 | } |
467 | } |
468 | |
42c75692 |
469 | fdEdxBef->Fill(mom,fTPCnSigma); |
bfc7c23b |
470 | fTPCnsigma->Fill(mom,fTPCnSigma); |
f09766dd |
471 | if(fTPCnSigma >= -1.0 && fTPCnSigma <= 3)fTrkEovPBef->Fill(pt,eop); |
bfc7c23b |
472 | |
bfc7c23b |
473 | Int_t pidpassed = 1; |
474 | |
475 | |
476 | //--- track accepted |
477 | AliHFEpidObject hfetrack; |
478 | hfetrack.SetAnalysisType(AliHFEpidObject::kESDanalysis); |
479 | hfetrack.SetRecTrack(track); |
fb87d707 |
480 | int centf = (int)cent; |
481 | hfetrack.SetCentrality(centf); //added |
bfc7c23b |
482 | hfetrack.SetPbPb(); |
483 | if(!fPID->IsSelected(&hfetrack, NULL, "", fPIDqa)) pidpassed = 0; |
484 | |
485 | if(pidpassed==0) continue; |
486 | |
487 | fTrkEovPAft->Fill(pt,eop); |
42c75692 |
488 | fdEdxAft->Fill(mom,fTPCnSigma); |
bfc7c23b |
489 | |
e7d87aef |
490 | // Fill real data |
fb87d707 |
491 | fIncpT->Fill(cent,pt); |
bfc7c23b |
492 | if(fFlagPhotonicElec) fPhoElecPt->Fill(cent,pt); |
f09766dd |
493 | if(fFlagConvinatElec) fSameElecPt->Fill(cent,pt); |
fb87d707 |
494 | |
495 | if(m20>0.0 && m20<0.3) |
496 | { |
497 | fIncpTM20->Fill(cent,pt); |
498 | if(fFlagPhotonicElec) fPhoElecPtM20->Fill(cent,pt); |
499 | if(fFlagConvinatElec) fSameElecPtM20->Fill(cent,pt); |
500 | } |
e7d87aef |
501 | |
502 | // MC |
503 | if(mcele==1) |
504 | { |
505 | if(mcBtoE || mcDtoE) |
506 | { |
507 | fIncpTMChfe->Fill(cent,pt); |
508 | if(m20>0.0 && m20<0.3)fIncpTMCM20hfe->Fill(cent,pt); |
509 | } |
510 | if(mcPho) |
511 | { |
512 | fIncpTMCpho->Fill(cent,pt); |
513 | if(fFlagPhotonicElec) fPhoElecPtMC->Fill(cent,pt); |
514 | if(fFlagConvinatElec) fSameElecPtMC->Fill(cent,pt); |
515 | |
516 | if(m20>0.0 && m20<0.3) |
517 | { |
518 | fIncpTMCM20pho->Fill(cent,pt); |
519 | if(fFlagPhotonicElec) fPhoElecPtMCM20->Fill(cent,pt); |
520 | if(fFlagConvinatElec) fSameElecPtMCM20->Fill(cent,pt); |
521 | } |
522 | } |
523 | } |
bfc7c23b |
524 | } |
525 | PostData(1, fOutputList); |
526 | } |
527 | //_________________________________________ |
528 | void AliAnalysisTaskHFECal::UserCreateOutputObjects() |
529 | { |
42c75692 |
530 | //--- Check MC |
531 | |
e7d87aef |
532 | //Bool_t mcData = kFALSE; |
42c75692 |
533 | if(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler()) |
534 | { |
e7d87aef |
535 | fmcData = kTRUE; |
42c75692 |
536 | printf("+++++ MC Data available"); |
537 | } |
e7d87aef |
538 | if(fmcData) |
42c75692 |
539 | { |
540 | printf("++++++++= MC analysis \n"); |
541 | } |
542 | else |
543 | { |
544 | printf("++++++++ real data analysis \n"); |
545 | } |
546 | |
fb87d707 |
547 | //---- Geometry |
548 | fGeom = AliEMCALGeometry::GetInstance("EMCAL_COMPLETEV1"); |
549 | |
bfc7c23b |
550 | //--------Initialize PID |
42c75692 |
551 | //fPID->SetHasMCData(kFALSE); |
e7d87aef |
552 | fPID->SetHasMCData(fmcData); |
bfc7c23b |
553 | if(!fPID->GetNumberOfPIDdetectors()) |
554 | { |
555 | fPID->AddDetector("TPC", 0); |
556 | fPID->AddDetector("EMCAL", 1); |
557 | } |
558 | |
559 | fPID->SortDetectors(); |
560 | fPIDqa = new AliHFEpidQAmanager(); |
561 | fPIDqa->Initialize(fPID); |
562 | |
563 | //--------Initialize correction Framework and Cuts |
564 | fCFM = new AliCFManager; |
565 | const Int_t kNcutSteps = AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kNcutStepsRecTrack + AliHFEcuts::kNcutStepsDETrack; |
566 | fCFM->SetNStepParticle(kNcutSteps); |
567 | for(Int_t istep = 0; istep < kNcutSteps; istep++) |
568 | fCFM->SetParticleCutsList(istep, NULL); |
569 | |
570 | if(!fCuts){ |
571 | AliWarning("Cuts not available. Default cuts will be used"); |
572 | fCuts = new AliHFEcuts; |
573 | fCuts->CreateStandardCuts(); |
574 | } |
575 | fCuts->Initialize(fCFM); |
576 | |
577 | //---------Output Tlist |
578 | fOutputList = new TList(); |
579 | fOutputList->SetOwner(); |
580 | fOutputList->Add(fPIDqa->MakeList("PIDQA")); |
581 | |
582 | fNoEvents = new TH1F("fNoEvents","",4,-0.5,3.5) ; |
583 | fOutputList->Add(fNoEvents); |
584 | |
fb87d707 |
585 | Int_t binsE[5] = {250, 100, 1000, 100, 10}; |
586 | Double_t xminE[5] = {1.0, -1, 0.0, 0, -0.5}; |
587 | Double_t xmaxE[5] = {3.5, 1, 100.0, 100, 9.5}; |
588 | fEMCAccE = new THnSparseD("fEMCAccE","EMC acceptance & E;#eta;#phi;Energy;Centrality;trugCondition;",5,binsE,xminE,xmaxE); |
bfc7c23b |
589 | fOutputList->Add(fEMCAccE); |
590 | |
591 | fTrkpt = new TH1F("fTrkpt","track pt",100,0,50); |
592 | fOutputList->Add(fTrkpt); |
593 | |
594 | fTrackPtBefTrkCuts = new TH1F("fTrackPtBefTrkCuts","track pt before track cuts",100,0,50); |
595 | fOutputList->Add(fTrackPtBefTrkCuts); |
596 | |
597 | fTrackPtAftTrkCuts = new TH1F("fTrackPtAftTrkCuts","track pt after track cuts",100,0,50); |
598 | fOutputList->Add(fTrackPtAftTrkCuts); |
599 | |
600 | fTPCnsigma = new TH2F("fTPCnsigma", "TPC - n sigma",100,0,50,200,-10,10); |
601 | fOutputList->Add(fTPCnsigma); |
602 | |
603 | fTrkEovPBef = new TH2F("fTrkEovPBef","track E/p before HFE pid",100,0,50,100,0,2); |
604 | fOutputList->Add(fTrkEovPBef); |
605 | |
606 | fTrkEovPAft = new TH2F("fTrkEovPAft","track E/p after HFE pid",100,0,50,100,0,2); |
607 | fOutputList->Add(fTrkEovPAft); |
608 | |
42c75692 |
609 | fdEdxBef = new TH2F("fdEdxBef","track dEdx vs p before HFE pid",100,0,50,200,-10,10); |
bfc7c23b |
610 | fOutputList->Add(fdEdxBef); |
611 | |
42c75692 |
612 | fdEdxAft = new TH2F("fdEdxAft","track dEdx vs p after HFE pid",100,0,50,200,-10,10); |
bfc7c23b |
613 | fOutputList->Add(fdEdxAft); |
614 | |
615 | fIncpT = new TH2F("fIncpT","HFE pid electro vs. centrality",100,0,100,100,0,50); |
616 | fOutputList->Add(fIncpT); |
617 | |
fb87d707 |
618 | fIncpTM20 = new TH2F("fIncpTM20","HFE pid electro vs. centrality with M20",100,0,100,100,0,50); |
619 | fOutputList->Add(fIncpTM20); |
f09766dd |
620 | |
621 | Int_t nBinspho[3] = { 100, 100, 500}; |
622 | Double_t minpho[3] = { 0., 0., 0.}; |
623 | Double_t maxpho[3] = {100., 50., 0.5}; |
624 | |
625 | fInvmassLS = new THnSparseD("fInvmassLS", "Inv mass of LS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2);", 3, nBinspho,minpho, maxpho); |
bfc7c23b |
626 | fOutputList->Add(fInvmassLS); |
627 | |
f09766dd |
628 | fInvmassULS = new THnSparseD("fInvmassULS", "Inv mass of ULS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2);", 3, nBinspho,minpho, maxpho); |
bfc7c23b |
629 | fOutputList->Add(fInvmassULS); |
630 | |
631 | fOpeningAngleLS = new TH1F("fOpeningAngleLS","Opening angle for LS pairs",100,0,1); |
632 | fOutputList->Add(fOpeningAngleLS); |
633 | |
634 | fOpeningAngleULS = new TH1F("fOpeningAngleULS","Opening angle for ULS pairs",100,0,1); |
635 | fOutputList->Add(fOpeningAngleULS); |
636 | |
637 | fPhotoElecPt = new TH1F("fPhotoElecPt", "photonic electron pt",100,0,50); |
638 | fOutputList->Add(fPhotoElecPt); |
639 | |
f09766dd |
640 | fPhoElecPt = new TH2F("fPhoElecPt", "Pho-inclusive electron pt",100,0,100,100,0,50); |
bfc7c23b |
641 | fOutputList->Add(fPhoElecPt); |
642 | |
fb87d707 |
643 | fPhoElecPtM20 = new TH2F("fPhoElecPtM20", "Pho-inclusive electron pt with M20",100,0,100,100,0,50); |
644 | fOutputList->Add(fPhoElecPtM20); |
645 | |
f09766dd |
646 | fSameElecPt = new TH2F("fSameElecPt", "Same-inclusive electron pt",100,0,100,100,0,50); |
647 | fOutputList->Add(fSameElecPt); |
648 | |
fb87d707 |
649 | fSameElecPtM20 = new TH2F("fSameElecPtM20", "Same-inclusive electron pt with M20",100,0,100,100,0,50); |
650 | fOutputList->Add(fSameElecPtM20); |
651 | |
bfc7c23b |
652 | fCent = new TH1F("fCent","Centrality",100,0,100) ; |
653 | fOutputList->Add(fCent); |
654 | |
655 | // Make common binning |
f09766dd |
656 | const Double_t kMinP = 2.; |
bfc7c23b |
657 | const Double_t kMaxP = 50.; |
658 | const Double_t kTPCSigMim = 40.; |
659 | const Double_t kTPCSigMax = 140.; |
660 | |
661 | // 1st histogram: TPC dEdx with/without EMCAL (p, pT, TPC Signal, phi, eta, Sig, e/p, ,match, cell, M02, M20, Disp, Centrality, select) |
e7d87aef |
662 | Int_t nBins[17] = { 480, 200, 60, 20, 600, 300, 100, 40, 200, 200, 200, 100, 3, 3, 3, 10, 3}; |
663 | Double_t min[17] = {kMinP, kTPCSigMim, 1.0, -1.0, -8.0, 0, 0, 0, 0.0, 0.0, 0.0, 0, -1.5, -0.5, -0.5, -0.5, -1.5}; |
664 | Double_t max[17] = {kMaxP, kTPCSigMax, 4.0, 1.0, 4.0, 3.0, 0.1, 40, 2.0, 2.0, 2.0, 100, 1.5, 2.5, 2.5, 9.5, 1.5}; |
665 | fEleInfo = new THnSparseD("fEleInfo", "Electron Info; pT [GeV/c]; TPC signal;phi;eta;nSig; E/p;Rmatch;Ncell;M02;M20;Disp;Centrality;charge;opp;same;trigCond;", 17, nBins, min, max); |
bfc7c23b |
666 | fOutputList->Add(fEleInfo); |
667 | |
fb87d707 |
668 | //<--- Trigger info |
669 | |
670 | fClsEBftTrigCut = new TH1F("fClsEBftTrigCut","cluster E before trigger selection",1000,0,100); |
671 | fOutputList->Add(fClsEBftTrigCut); |
672 | |
673 | fClsEAftTrigCut = new TH1F("fClsEAftTrigCut","cluster E if cls has 0 trigcut channel",1000,0,100); |
674 | fOutputList->Add(fClsEAftTrigCut); |
675 | |
676 | fClsEAftTrigCut1 = new TH1F("fClsEAftTrigCut1","cluster E if cls with trig channel",1000,0,100); |
677 | fOutputList->Add(fClsEAftTrigCut1); |
678 | |
679 | fClsEAftTrigCut2 = new TH1F("fClsEAftTrigCut2","cluster E if cls with trigcut channel",1000,0,100); |
680 | fOutputList->Add(fClsEAftTrigCut2); |
681 | |
682 | fClsEAftTrigCut3 = new TH1F("fClsEAftTrigCut3","cluster E if cls with trigcut channel + nCell>Ecorrect",1000,0,100); |
683 | fOutputList->Add(fClsEAftTrigCut3); |
684 | |
685 | fClsEAftTrigCut4 = new TH1F("fClsEAftTrigCut4","cluster E if cls with trigcut channel + nCell>Ecorrect + cls time cut",1000,0,100); |
686 | fOutputList->Add(fClsEAftTrigCut4); |
687 | |
688 | fClsETime = new TH2F("fClsETime", "Cls time vs E; E; time;",1000,0,100,1000,-0.0000002,0.0000009); |
689 | fOutputList->Add(fClsETime); |
690 | |
691 | fClsETime1 = new TH2F("fClsETime1", "Cls time vs E if cls contains trigger channel; E; time;",1000,0,100,1000,-0.0000002,0.0000009); |
692 | fOutputList->Add(fClsETime1); |
693 | |
694 | fTrigTimes = new TH1F("fTrigTimes", "Trigger time; time; N;",25,0,25); |
695 | fOutputList->Add(fTrigTimes); |
696 | |
42c75692 |
697 | fCellCheck = new TH2F("fCellCheck", "Cell vs E; E GeV; Cell ID",10,6,26,12000,0,12000); |
698 | fOutputList->Add(fCellCheck); |
699 | |
e7d87aef |
700 | //<---------- MC |
701 | |
702 | fInputHFEMC = new TH2F("fInputHFEMC","Input MC HFE pid electro vs. centrality",100,0,100,100,0,50); |
703 | fOutputList->Add(fInputHFEMC); |
704 | |
705 | fIncpTMChfe = new TH2F("fIncpTMChfe","MC HFE pid electro vs. centrality",100,0,100,100,0,50); |
706 | fOutputList->Add(fIncpTMChfe); |
707 | |
708 | fIncpTMCM20hfe = new TH2F("fIncpTMCM20hfe","MC HFE pid electro vs. centrality with M20",100,0,100,100,0,50); |
709 | fOutputList->Add(fIncpTMCM20hfe); |
710 | |
711 | fIncpTMCpho = new TH2F("fIncpTMCpho","MC Pho pid electro vs. centrality",100,0,100,100,0,50); |
712 | fOutputList->Add(fIncpTMCpho); |
713 | |
714 | fIncpTMCM20pho = new TH2F("fIncpTMCM20pho","MC Pho pid electro vs. centrality with M20",100,0,100,100,0,50); |
715 | fOutputList->Add(fIncpTMCM20pho); |
716 | |
717 | fPhoElecPtMC = new TH2F("fPhoElecPtMC", "MC Pho-inclusive electron pt",100,0,100,100,0,50); |
718 | fOutputList->Add(fPhoElecPtMC); |
719 | |
720 | fPhoElecPtMCM20 = new TH2F("fPhoElecPtMCM20", "MC Pho-inclusive electron pt with M20",100,0,100,100,0,50); |
721 | fOutputList->Add(fPhoElecPtM20); |
722 | |
723 | fSameElecPtMC = new TH2F("fSameElecPtMC", "MC Same-inclusive electron pt",100,0,100,100,0,50); |
724 | fOutputList->Add(fSameElecPtMC); |
725 | |
726 | fSameElecPtMCM20 = new TH2F("fSameElecPtMCM20", "MC Same-inclusive electron pt with M20",100,0,100,100,0,50); |
727 | fOutputList->Add(fSameElecPtMCM20); |
728 | |
729 | |
bfc7c23b |
730 | PostData(1,fOutputList); |
731 | } |
732 | |
733 | //________________________________________________________________________ |
734 | void AliAnalysisTaskHFECal::Terminate(Option_t *) |
735 | { |
736 | // Info("Terminate"); |
737 | AliAnalysisTaskSE::Terminate(); |
738 | } |
739 | |
740 | //________________________________________________________________________ |
741 | Bool_t AliAnalysisTaskHFECal::ProcessCutStep(Int_t cutStep, AliVParticle *track) |
742 | { |
743 | // Check single track cuts for a given cut step |
744 | const Int_t kMCOffset = AliHFEcuts::kNcutStepsMCTrack; |
745 | if(!fCFM->CheckParticleCuts(cutStep + kMCOffset, track)) return kFALSE; |
746 | return kTRUE; |
747 | } |
748 | //_________________________________________ |
f09766dd |
749 | //void AliAnalysisTaskHFECal::SelectPhotonicElectron(Int_t itrack, Double_t cent, AliESDtrack *track, Bool_t &fFlagPhotonicElec) |
750 | void AliAnalysisTaskHFECal::SelectPhotonicElectron(Int_t itrack, Double_t cent, AliESDtrack *track, Bool_t &fFlagPhotonicElec, Bool_t &fFlagConvinatElec) |
bfc7c23b |
751 | { |
752 | //Identify non-heavy flavour electrons using Invariant mass method |
753 | |
fb87d707 |
754 | fTrackCuts->SetAcceptKinkDaughters(kFALSE); |
f09766dd |
755 | //fTrackCuts->SetRequireTPCRefit(kTRUE); |
756 | //fTrackCuts->SetEtaRange(-0.7,0.7); |
757 | //fTrackCuts->SetRequireSigmaToVertex(kTRUE); |
bfc7c23b |
758 | fTrackCuts->SetMaxChi2PerClusterTPC(3.5); |
f09766dd |
759 | fTrackCuts->SetMinNClustersTPC(70); |
bfc7c23b |
760 | |
761 | const AliESDVertex *pVtx = fESD->GetPrimaryVertex(); |
762 | |
763 | Bool_t flagPhotonicElec = kFALSE; |
f09766dd |
764 | Bool_t flagConvinatElec = kFALSE; |
bfc7c23b |
765 | |
5e65985c |
766 | //for(Int_t jTracks = itrack+1; jTracks<fESD->GetNumberOfTracks(); jTracks++){ |
767 | for(Int_t jTracks = 0; jTracks<fESD->GetNumberOfTracks(); jTracks++){ |
bfc7c23b |
768 | AliESDtrack* trackAsso = fESD->GetTrack(jTracks); |
769 | if (!trackAsso) { |
770 | printf("ERROR: Could not receive track %d\n", jTracks); |
771 | continue; |
772 | } |
5e65985c |
773 | if(itrack==jTracks)continue; |
774 | |
f09766dd |
775 | Double_t dEdxAsso = -999., ptPrim=-999., ptAsso=-999., openingAngle = -999.; |
bfc7c23b |
776 | Double_t mass=999., width = -999; |
777 | Bool_t fFlagLS=kFALSE, fFlagULS=kFALSE; |
778 | |
f09766dd |
779 | ptPrim = track->Pt(); |
780 | |
bfc7c23b |
781 | dEdxAsso = trackAsso->GetTPCsignal(); |
782 | ptAsso = trackAsso->Pt(); |
783 | Int_t chargeAsso = trackAsso->Charge(); |
784 | Int_t charge = track->Charge(); |
785 | |
5e65985c |
786 | //if(ptAsso <0.3) continue; |
787 | if(ptAsso <0.5) continue; |
bfc7c23b |
788 | if(!fTrackCuts->AcceptTrack(trackAsso)) continue; |
f09766dd |
789 | if(dEdxAsso <65 || dEdxAsso>100) continue; //11a pass1 |
bfc7c23b |
790 | |
791 | Int_t fPDGe1 = 11; Int_t fPDGe2 = 11; |
792 | if(charge>0) fPDGe1 = -11; |
793 | if(chargeAsso>0) fPDGe2 = -11; |
794 | |
795 | if(charge == chargeAsso) fFlagLS = kTRUE; |
796 | if(charge != chargeAsso) fFlagULS = kTRUE; |
797 | |
798 | AliKFParticle ge1(*track, fPDGe1); |
799 | AliKFParticle ge2(*trackAsso, fPDGe2); |
800 | AliKFParticle recg(ge1, ge2); |
801 | |
802 | if(recg.GetNDF()<1) continue; |
803 | Double_t chi2recg = recg.GetChi2()/recg.GetNDF(); |
804 | if(TMath::Sqrt(TMath::Abs(chi2recg))>3.) continue; |
805 | |
806 | AliKFVertex primV(*pVtx); |
807 | primV += recg; |
808 | recg.SetProductionVertex(primV); |
809 | |
810 | recg.SetMassConstraint(0,0.0001); |
811 | |
812 | openingAngle = ge1.GetAngle(ge2); |
813 | if(fFlagLS) fOpeningAngleLS->Fill(openingAngle); |
814 | if(fFlagULS) fOpeningAngleULS->Fill(openingAngle); |
815 | |
816 | if(openingAngle > fOpeningAngleCut) continue; |
817 | |
818 | recg.GetMass(mass,width); |
819 | |
f09766dd |
820 | double phoinfo[3]; |
821 | phoinfo[0] = cent; |
822 | phoinfo[1] = ptPrim; |
823 | phoinfo[2] = mass; |
824 | |
825 | if(fFlagLS) fInvmassLS->Fill(phoinfo); |
826 | if(fFlagULS) fInvmassULS->Fill(phoinfo); |
bfc7c23b |
827 | |
828 | if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec){ |
829 | flagPhotonicElec = kTRUE; |
830 | } |
f09766dd |
831 | if(mass<fInvmassCut && fFlagLS && !flagConvinatElec){ |
832 | flagConvinatElec = kTRUE; |
833 | } |
bfc7c23b |
834 | |
835 | } |
836 | fFlagPhotonicElec = flagPhotonicElec; |
f09766dd |
837 | fFlagConvinatElec = flagConvinatElec; |
bfc7c23b |
838 | |
839 | } |
840 | |
fb87d707 |
841 | |
842 | //_________________________________________ |
843 | void AliAnalysisTaskHFECal::FindTriggerClusters() |
844 | { |
845 | // constants |
846 | const int nModuleCols = 2; |
847 | const int nModuleRows = 5; |
848 | const int nColsFeeModule = 48; |
849 | const int nRowsFeeModule = 24; |
850 | const int nColsFaltroModule = 24; |
851 | const int nRowsFaltroModule = 12; |
852 | //const int faltroWidthMax = 20; |
853 | |
854 | // part 1, trigger extraction ------------------------------------- |
855 | Int_t globCol, globRow; |
856 | //Int_t ntimes=0, nTrigChannel=0, nTrigChannelCut=0, trigInCut=0; |
857 | Int_t ntimes=0, nTrigChannel=0, nTrigChannelCut=0; |
858 | |
859 | //Int_t trigtimes[faltroWidthMax]; |
860 | Double_t cellTime[nColsFeeModule*nModuleCols][nRowsFeeModule*nModuleRows]; |
861 | Double_t cellEnergy[nColsFeeModule*nModuleCols][nRowsFeeModule*nModuleRows]; |
862 | //Double_t fTrigCutLow = 6; |
863 | //Double_t fTrigCutHigh = 10; |
864 | Double_t fTimeCutLow = 469e-09; |
865 | Double_t fTimeCutHigh = 715e-09; |
866 | |
867 | AliESDCaloTrigger * fCaloTrigger = fESD->GetCaloTrigger( "EMCAL" ); |
868 | |
869 | // erase trigger maps |
870 | for(Int_t i = 0; i < nColsFaltroModule*nModuleCols; i++ ) |
871 | { |
872 | for(Int_t j = 0; j < nRowsFaltroModule*nModuleRows; j++ ) |
873 | { |
874 | ftriggersCut[i][j] = 0; |
875 | ftriggers[i][j] = 0; |
876 | ftriggersTime[i][j] = 0; |
877 | } |
878 | } |
879 | |
880 | Int_t iglobCol=0, iglobRow=0; |
881 | // go through triggers |
882 | if( fCaloTrigger->GetEntries() > 0 ) |
883 | { |
884 | // needs reset |
885 | fCaloTrigger->Reset(); |
886 | while( fCaloTrigger->Next() ) |
887 | { |
888 | fCaloTrigger->GetPosition( globCol, globRow ); |
889 | fCaloTrigger->GetNL0Times( ntimes ); |
890 | /* |
891 | // no L0s |
892 | if( ntimes < 1 ) continue; |
893 | // get precise timings |
894 | fCaloTrigger->GetL0Times( trigtimes ); |
895 | trigInCut = 0; |
896 | for(Int_t i = 0; i < ntimes; i++ ) |
897 | { |
898 | // save the first trigger time in channel |
899 | if( i == 0 || triggersTime[globCol][globRow] > trigtimes[i] ) |
900 | triggersTime[globCol][globRow] = trigtimes[i]; |
901 | //printf("trigger times: %d\n",trigtimes[i]); |
902 | // check if in cut |
903 | if(trigtimes[i] > fTrigCutLow && trigtimes[i] < fTrigCutHigh ) |
904 | trigInCut = 1; |
905 | |
906 | fTrigTimes->Fill(trigtimes[i]); |
907 | } |
908 | */ |
909 | |
910 | //L1 analysis from AliAnalysisTaskEMCALTriggerQA |
911 | Int_t bit = 0; |
912 | fCaloTrigger->GetTriggerBits(bit); |
913 | |
914 | Int_t ts = 0; |
915 | fCaloTrigger->GetL1TimeSum(ts); |
916 | if (ts > 0)ftriggers[globCol][globRow] = 1; |
917 | // number of triggered channels in event |
918 | nTrigChannel++; |
919 | // ... inside cut |
920 | if(ts>0 && (bit >> 6 & 0x1)) |
921 | { |
922 | iglobCol = globCol; |
923 | iglobRow = globRow; |
924 | nTrigChannelCut++; |
925 | ftriggersCut[globCol][globRow] = 1; |
926 | } |
927 | |
928 | } // calo trigger entries |
929 | } // has calo trigger entries |
930 | |
931 | // part 2 go through the clusters here ----------------------------------- |
932 | Int_t nCluster=0, nCell=0, iCell=0, gCell=0; |
933 | Short_t cellAddr, nSACell, mclabel; |
934 | //Int_t nSACell, iSACell, mclabel; |
935 | Int_t iSACell; |
42c75692 |
936 | Double_t cellAmp=0, cellTimeT=0, clusterTime=0, efrac=0; |
fb87d707 |
937 | Int_t nSupMod, nModule, nIphi, nIeta, iphi, ieta, gphi, geta, feta, fphi; |
938 | |
939 | TRefArray *fCaloClusters = new TRefArray(); |
940 | fESD->GetEMCALClusters( fCaloClusters ); |
941 | nCluster = fCaloClusters->GetEntries(); |
942 | |
943 | |
944 | // save all cells times if there are clusters |
945 | if( nCluster > 0 ){ |
946 | // erase time map |
947 | for(Int_t i = 0; i < nColsFeeModule*nModuleCols; i++ ){ |
948 | for(Int_t j = 0; j < nRowsFeeModule*nModuleRows; j++ ){ |
949 | cellTime[i][j] = 0.; |
950 | cellEnergy[i][j] = 0.; |
951 | } |
952 | } |
953 | |
954 | // get cells |
955 | AliESDCaloCells *fCaloCells = fESD->GetEMCALCells(); |
956 | //AliVCaloCells fCaloCells = fESD->GetEMCALCells(); |
957 | nSACell = fCaloCells->GetNumberOfCells(); |
958 | for(iSACell = 0; iSACell < nSACell; iSACell++ ){ |
959 | // get the cell info *fCal |
960 | fCaloCells->GetCell( iSACell, cellAddr, cellAmp, cellTimeT , mclabel, efrac); |
961 | |
962 | // get cell position |
963 | fGeom->GetCellIndex( cellAddr, nSupMod, nModule, nIphi, nIeta ); |
964 | fGeom->GetCellPhiEtaIndexInSModule( nSupMod,nModule, nIphi, nIeta, iphi, ieta); |
965 | |
966 | // convert co global phi eta |
967 | gphi = iphi + nRowsFeeModule*(nSupMod/2); |
968 | geta = ieta + nColsFeeModule*(nSupMod%2); |
969 | |
970 | // save cell time and energy |
971 | cellTime[geta][gphi] = cellTimeT; |
972 | cellEnergy[geta][gphi] = cellAmp; |
973 | |
974 | } |
975 | } |
976 | |
977 | Int_t nClusterTrig, nClusterTrigCut; |
978 | UShort_t *cellAddrs; |
979 | Double_t clsE=-999, clsEta=-999, clsPhi=-999; |
980 | Float_t clsPos[3] = {0.,0.,0.}; |
981 | |
982 | for(Int_t icl=0; icl<fESD->GetNumberOfCaloClusters(); icl++) |
983 | { |
984 | AliESDCaloCluster *cluster = fESD->GetCaloCluster(icl); |
985 | if(!cluster || !cluster->IsEMCAL()) continue; |
986 | |
987 | // get cluster cells |
988 | nCell = cluster->GetNCells(); |
989 | |
990 | // get cluster energy |
991 | clsE = cluster->E(); |
992 | |
993 | // get cluster position |
994 | cluster->GetPosition(clsPos); |
995 | TVector3 clsPosVec(clsPos[0],clsPos[1],clsPos[2]); |
996 | clsEta = clsPosVec.Eta(); |
997 | clsPhi = clsPosVec.Phi(); |
998 | |
999 | // get the cell addresses |
1000 | cellAddrs = cluster->GetCellsAbsId(); |
1001 | |
1002 | // check if the cluster contains cell, that was marked as triggered |
1003 | nClusterTrig = 0; |
1004 | nClusterTrigCut = 0; |
1005 | |
1006 | // loop the cells to check, if cluser in acceptance |
1007 | // any cluster with a cell outside acceptance is not considered |
1008 | for( iCell = 0; iCell < nCell; iCell++ ) |
1009 | { |
42c75692 |
1010 | // check hot cell |
1011 | if(clsE>6.0)fCellCheck->Fill(clsE,cellAddrs[iCell]); |
1012 | |
fb87d707 |
1013 | // get cell position |
1014 | fGeom->GetCellIndex( cellAddrs[iCell], nSupMod, nModule, nIphi, nIeta ); |
1015 | fGeom->GetCellPhiEtaIndexInSModule( nSupMod,nModule, nIphi, nIeta, iphi, ieta); |
1016 | |
1017 | // convert co global phi eta |
1018 | gphi = iphi + nRowsFeeModule*(nSupMod/2); |
1019 | geta = ieta + nColsFeeModule*(nSupMod%2); |
1020 | |
1021 | if( cellTime[geta][gphi] > 0. ){ |
1022 | clusterTime += cellTime[geta][gphi]; |
1023 | gCell++; |
1024 | } |
1025 | |
1026 | // get corresponding FALTRO |
1027 | fphi = gphi / 2; |
1028 | feta = geta / 2; |
1029 | |
1030 | // try to match with a triggered |
1031 | if( ftriggers[feta][fphi]==1) |
1032 | { nClusterTrig++; |
1033 | } |
1034 | if( ftriggersCut[feta][fphi]==1) |
1035 | { nClusterTrigCut++; |
1036 | } |
1037 | |
1038 | } // cells |
1039 | |
1040 | |
1041 | if( gCell > 0 ) |
1042 | clusterTime = clusterTime / (Double_t)gCell; |
1043 | // fix the reconstriction code time 100ns jumps |
1044 | if( fESD->GetBunchCrossNumber() % 4 < 2 ) |
1045 | clusterTime -= 0.0000001; |
1046 | |
1047 | fClsETime->Fill(clsE,clusterTime); |
1048 | fClsEBftTrigCut->Fill(clsE); |
1049 | |
1050 | if(nClusterTrig>0){ |
1051 | fClsETime1->Fill(clsE,clusterTime); |
1052 | } |
1053 | |
1054 | if(nClusterTrig>0){ |
1055 | cluster->SetChi2(1); |
1056 | fClsEAftTrigCut1->Fill(clsE); |
1057 | } |
1058 | |
1059 | if(nClusterTrigCut>0){ |
1060 | cluster->SetChi2(2); |
1061 | fClsEAftTrigCut2->Fill(clsE); |
1062 | } |
1063 | |
1064 | if(nClusterTrigCut>0 && ( nCell > (1 + clsE / 3))) |
1065 | { |
1066 | cluster->SetChi2(3); |
1067 | fClsEAftTrigCut3->Fill(clsE); |
1068 | } |
1069 | |
1070 | if(nClusterTrigCut>0 && (nCell > (1 + clsE / 3) )&&( clusterTime > fTimeCutLow && clusterTime < fTimeCutHigh )) |
1071 | { |
1072 | // cluster->SetChi2(4); |
1073 | fClsEAftTrigCut4->Fill(clsE); |
1074 | } |
1075 | if(nClusterTrigCut<1) |
1076 | { |
1077 | cluster->SetChi2(0); |
1078 | |
1079 | fClsEAftTrigCut->Fill(clsE); |
1080 | } |
1081 | |
1082 | } // clusters |
1083 | } |
1084 | |
1085 | |
1086 | |
1087 | |
1088 | |
1089 | |
1090 | |