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