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