]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - PWGHF/hfe/AliAnalysisTaskHFECal.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
on Proof dont delete the output objects in the destructor
[u/mrichter/AliRoot.git] / PWGHF / hfe / AliAnalysisTaskHFECal.cxx
... / ...
CommitLineData
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
83ClassImp(AliAnalysisTaskHFECal)
84//________________________________________________________________________
85AliAnalysisTaskHFECal::AliAnalysisTaskHFECal(const char *name)
86 : AliAnalysisTaskSE(name)
87 ,fESD(0)
88 ,fMC(0)
89 ,fGeom(0)
90 ,fOutputList(0)
91 ,fqahist(0)
92 ,fTrackCuts(0)
93 ,fCuts(0)
94 ,fIdentifiedAsOutInz(kFALSE)
95 ,fPassTheEventCut(kFALSE)
96 ,fRejectKinkMother(kFALSE)
97 ,fmcData(kFALSE)
98 ,fVz(0.0)
99 ,fCFM(0)
100 ,fPID(0)
101 ,fPIDqa(0)
102 ,fOpeningAngleCut(0.1)
103 ,fInvmassCut(0.01)
104 ,fNoEvents(0)
105 ,fEMCAccE(0)
106 ,fTrkpt(0)
107 ,fTrkEovPBef(0)
108 ,fTrkEovPAft(0)
109 ,fdEdxBef(0)
110 ,fdEdxAft(0)
111 ,fIncpT(0)
112 ,fIncpTM20(0)
113 ,fInvmassLS(0)
114 ,fInvmassULS(0)
115 ,fInvmassLSmc(0)
116 ,fInvmassULSmc(0)
117 ,fInvmassLSmc0(0)
118 ,fInvmassLSmc1(0)
119 ,fInvmassLSmc2(0)
120 ,fInvmassLSmc3(0)
121 ,fInvmassULSmc0(0)
122 ,fInvmassULSmc1(0)
123 ,fInvmassULSmc2(0)
124 ,fInvmassULSmc3(0)
125 ,fOpeningAngleLS(0)
126 ,fOpeningAngleULS(0)
127 ,fPhotoElecPt(0)
128 ,fPhoElecPt(0)
129 ,fPhoElecPtM20(0)
130 ,fSameElecPt(0)
131 ,fSameElecPtM20(0)
132 ,fTrackPtBefTrkCuts(0)
133 ,fTrackPtAftTrkCuts(0)
134 ,fTPCnsigma(0)
135 ,fCent(0)
136 ,fEleInfo(0)
137 /*
138 ,fClsEBftTrigCut(0)
139 ,fClsEAftTrigCut(0)
140 ,fClsEAftTrigCut1(0)
141 ,fClsEAftTrigCut2(0)
142 ,fClsEAftTrigCut3(0)
143 ,fClsEAftTrigCut4(0)
144 ,fClsETime(0)
145 ,fClsETime1(0)
146 ,fTrigTimes(0)
147 ,fCellCheck(0)
148 */
149 ,fInputHFEMC(0)
150 ,fInputAlle(0)
151 ,fIncpTMChfe(0)
152 ,fIncpTMChfeAll(0)
153 ,fIncpTMCM20hfe(0)
154 ,fIncpTMCM20hfeAll(0)
155 ,fIncpTMCM20hfeCheck(0)
156 ,fIncpTMCpho(0)
157 ,fIncpTMCM20pho(0)
158 ,fPhoElecPtMC(0)
159 ,fPhoElecPtMCM20(0)
160 ,fSameElecPtMC(0)
161 ,fSameElecPtMCM20(0)
162 ,fIncpTMCM20pho_pi0e(0)
163 ,fPhoElecPtMCM20_pi0e(0)
164 ,fSameElecPtMCM20_pi0e(0)
165 ,fIncpTMCM20pho_eta(0)
166 ,fPhoElecPtMCM20_eta(0)
167 ,fSameElecPtMCM20_eta(0)
168 ,CheckNclust(0)
169 ,CheckNits(0)
170 ,Hpi0pTcheck(0)
171 ,HETApTcheck(0)
172 ,HphopTcheck(0)
173 ,fpTCheck(0)
174 ,fMomDtoE(0)
175{
176 //Named constructor
177
178 fPID = new AliHFEpid("hfePid");
179 fTrackCuts = new AliESDtrackCuts();
180
181 // Define input and output slots here
182 // Input slot #0 works with a TChain
183 DefineInput(0, TChain::Class());
184 // Output slot #0 id reserved by the base class for AOD
185 // Output slot #1 writes into a TH1 container
186 // DefineOutput(1, TH1I::Class());
187 DefineOutput(1, TList::Class());
188 // DefineOutput(3, TTree::Class());
189}
190
191//________________________________________________________________________
192AliAnalysisTaskHFECal::AliAnalysisTaskHFECal()
193 : AliAnalysisTaskSE("DefaultAnalysis_AliAnalysisTaskHFECal")
194 ,fESD(0)
195 ,fMC(0)
196 ,fGeom(0)
197 ,fOutputList(0)
198 ,fqahist(0)
199 ,fTrackCuts(0)
200 ,fCuts(0)
201 ,fIdentifiedAsOutInz(kFALSE)
202 ,fPassTheEventCut(kFALSE)
203 ,fRejectKinkMother(kFALSE)
204 ,fmcData(kFALSE)
205 ,fVz(0.0)
206 ,fCFM(0)
207 ,fPID(0)
208 ,fPIDqa(0)
209 ,fOpeningAngleCut(0.1)
210 ,fInvmassCut(0.01)
211 ,fNoEvents(0)
212 ,fEMCAccE(0)
213 ,fTrkpt(0)
214 ,fTrkEovPBef(0)
215 ,fTrkEovPAft(0)
216 ,fdEdxBef(0)
217 ,fdEdxAft(0)
218 ,fIncpT(0)
219 ,fIncpTM20(0)
220 ,fInvmassLS(0)
221 ,fInvmassULS(0)
222 ,fInvmassLSmc(0)
223 ,fInvmassULSmc(0)
224 ,fInvmassLSmc0(0)
225 ,fInvmassLSmc1(0)
226 ,fInvmassLSmc2(0)
227 ,fInvmassLSmc3(0)
228 ,fInvmassULSmc0(0)
229 ,fInvmassULSmc1(0)
230 ,fInvmassULSmc2(0)
231 ,fInvmassULSmc3(0)
232 ,fOpeningAngleLS(0)
233 ,fOpeningAngleULS(0)
234 ,fPhotoElecPt(0)
235 ,fPhoElecPt(0)
236 ,fPhoElecPtM20(0)
237 ,fSameElecPt(0)
238 ,fSameElecPtM20(0)
239 ,fTrackPtBefTrkCuts(0)
240 ,fTrackPtAftTrkCuts(0)
241 ,fTPCnsigma(0)
242 ,fCent(0)
243 ,fEleInfo(0)
244 /*
245 ,fClsEBftTrigCut(0)
246 ,fClsEAftTrigCut(0)
247 ,fClsEAftTrigCut1(0)
248 ,fClsEAftTrigCut2(0)
249 ,fClsEAftTrigCut3(0)
250 ,fClsEAftTrigCut4(0)
251 ,fClsETime(0)
252 ,fClsETime1(0)
253 ,fTrigTimes(0)
254 ,fCellCheck(0)
255 */
256 ,fInputHFEMC(0)
257 ,fInputAlle(0)
258 ,fIncpTMChfe(0)
259 ,fIncpTMChfeAll(0)
260 ,fIncpTMCM20hfe(0)
261 ,fIncpTMCM20hfeAll(0)
262 ,fIncpTMCM20hfeCheck(0)
263 ,fIncpTMCpho(0)
264 ,fIncpTMCM20pho(0)
265 ,fPhoElecPtMC(0)
266 ,fPhoElecPtMCM20(0)
267 ,fSameElecPtMC(0)
268 ,fSameElecPtMCM20(0)
269 ,fIncpTMCM20pho_pi0e(0)
270 ,fPhoElecPtMCM20_pi0e(0)
271 ,fSameElecPtMCM20_pi0e(0)
272 ,fIncpTMCM20pho_eta(0)
273 ,fPhoElecPtMCM20_eta(0)
274 ,fSameElecPtMCM20_eta(0)
275 ,CheckNclust(0)
276 ,CheckNits(0)
277 ,Hpi0pTcheck(0)
278 ,HETApTcheck(0)
279 ,HphopTcheck(0)
280 ,fpTCheck(0)
281 ,fMomDtoE(0)
282{
283 //Default constructor
284 fPID = new AliHFEpid("hfePid");
285
286 fTrackCuts = new AliESDtrackCuts();
287
288 // Constructor
289 // Define input and output slots here
290 // Input slot #0 works with a TChain
291 DefineInput(0, TChain::Class());
292 // Output slot #0 id reserved by the base class for AOD
293 // Output slot #1 writes into a TH1 container
294 // DefineOutput(1, TH1I::Class());
295 DefineOutput(1, TList::Class());
296 //DefineOutput(3, TTree::Class());
297}
298//_________________________________________
299
300AliAnalysisTaskHFECal::~AliAnalysisTaskHFECal()
301{
302 //Destructor
303
304 delete fOutputList;
305 delete fGeom;
306 delete fPID;
307 delete fCuts;
308 delete fCFM;
309 delete fPIDqa;
310 delete fTrackCuts;
311}
312//_________________________________________
313
314void AliAnalysisTaskHFECal::UserExec(Option_t*)
315{
316 //Main loop
317 //Called for each event
318
319 // create pointer to event
320 fESD = dynamic_cast<AliESDEvent*>(InputEvent());
321 if (!fESD) {
322 printf("ERROR: fESD not available\n");
323 return;
324 }
325
326 if(!fCuts){
327 AliError("HFE cuts not available");
328 return;
329 }
330
331 if(!fPID->IsInitialized()){
332 // Initialize PID with the given run number
333 AliWarning("PID not initialised, get from Run no");
334 fPID->InitializePID(fESD->GetRunNumber());
335 }
336
337 if(fmcData)fMC = MCEvent();
338 AliStack* stack = NULL;
339 if(fmcData && fMC)stack = fMC->Stack();
340
341 Float_t cent = -1.;
342 AliCentrality *centrality = fESD->GetCentrality();
343 cent = centrality->GetCentralityPercentile("V0M");
344
345 //---- fill MC track info
346 if(fmcData && fMC)
347 {
348 Int_t nParticles = stack->GetNtrack();
349 for (Int_t iParticle = 0; iParticle < nParticles; iParticle++) {
350 TParticle* particle = stack->Particle(iParticle);
351 int fPDG = particle->GetPdgCode();
352 double mcZvertex = fMC->GetPrimaryVertex()->GetZ();
353 double pTMC = particle->Pt();
354 double proR = particle->R();
355 double etaMC = particle->Eta();
356 if(fabs(etaMC)>0.6)continue;
357 Bool_t mcInDtoE= kFALSE;
358 Bool_t mcInBtoE= kFALSE;
359
360 Bool_t MChijing = fMC->IsFromBGEvent(iParticle);
361 //if(!MChijing)printf("not MC hijing");
362 int iHijing = 1;
363 if(!MChijing)iHijing = 0;
364 if(fPDG==111)Hpi0pTcheck->Fill(pTMC,iHijing);
365 if(fPDG==221)HETApTcheck->Fill(pTMC,iHijing);
366
367 if(particle->GetFirstMother()>-1 && fPDG==22)
368 {
369 int parentPID = stack->Particle(particle->GetFirstMother())->GetPdgCode();
370 if(parentPID==111 || parentPID==221)HphopTcheck->Fill(pTMC,iHijing); // pi0->g & eta->g
371 }
372
373 if(particle->GetFirstMother()>-1 && fabs(fPDG)==11)
374 {
375 int parentPID = stack->Particle(particle->GetFirstMother())->GetPdgCode();
376 if((fabs(parentPID)==411 || fabs(parentPID)==413 || fabs(parentPID)==421 || fabs(parentPID)==423 || fabs(parentPID)==431)&& fabs(fPDG)==11)mcInDtoE = kTRUE;
377 if((fabs(parentPID)==511 || fabs(parentPID)==513 || fabs(parentPID)==521 || fabs(parentPID)==523 || fabs(parentPID)==531)&& fabs(fPDG)==11)mcInBtoE = kTRUE;
378 if((mcInBtoE || mcInDtoE) && fabs(mcZvertex)<10.0)fInputHFEMC->Fill(cent,pTMC);
379 }
380
381 if(proR<7.0 && fabs(fPDG)==11)fInputAlle->Fill(cent,pTMC);
382
383 }
384 }
385
386 fNoEvents->Fill(0);
387
388 Int_t fNOtrks = fESD->GetNumberOfTracks();
389 const AliESDVertex *pVtx = fESD->GetPrimaryVertex();
390
391 Double_t pVtxZ = -999;
392 pVtxZ = pVtx->GetZ();
393
394 if(TMath::Abs(pVtxZ)>10) return;
395 fNoEvents->Fill(1);
396
397 if(fNOtrks<2) return;
398 fNoEvents->Fill(2);
399
400 AliPIDResponse *pidResponse = fInputHandler->GetPIDResponse();
401 if(!pidResponse){
402 AliDebug(1, "Using default PID Response");
403 pidResponse = AliHFEtools::GetDefaultPID(kFALSE, fInputEvent->IsA() == AliAODEvent::Class());
404 }
405
406 fPID->SetPIDResponse(pidResponse);
407
408 fCFM->SetRecEventInfo(fESD);
409
410 //Float_t cent = -1.;
411 //AliCentrality *centrality = fESD->GetCentrality();
412 //cent = centrality->GetCentralityPercentile("V0M");
413 fCent->Fill(cent);
414
415 //if(cent>90.) return;
416
417 // Calorimeter info.
418
419 FindTriggerClusters();
420
421 // make EMCAL array
422 for(Int_t iCluster=0; iCluster<fESD->GetNumberOfCaloClusters(); iCluster++)
423 {
424 AliESDCaloCluster *clust = fESD->GetCaloCluster(iCluster);
425 if(clust->IsEMCAL())
426 {
427 double clustE = clust->E();
428 float emcx[3]; // cluster pos
429 clust->GetPosition(emcx);
430 TVector3 clustpos(emcx[0],emcx[1],emcx[2]);
431 double emcphi = clustpos.Phi();
432 double emceta = clustpos.Eta();
433 double calInfo[5];
434 calInfo[0] = emcphi; calInfo[1] = emceta; calInfo[2] = clustE; calInfo[3] = cent; calInfo[4] = clust->Chi2();
435 //if(clustE>1.5)fEMCAccE->Fill(calInfo);
436 //if(fqahist==1 && clustE>1.5)fEMCAccE->Fill(calInfo);
437 }
438 }
439
440 // Track loop
441 for (Int_t iTracks = 0; iTracks < fESD->GetNumberOfTracks(); iTracks++) {
442 AliESDtrack* track = fESD->GetTrack(iTracks);
443 if (!track) {
444 printf("ERROR: Could not receive track %d\n", iTracks);
445 continue;
446 }
447
448 Bool_t mcPho = kFALSE;
449 Bool_t mcDtoE= kFALSE;
450 Bool_t mcBtoE= kFALSE;
451 Bool_t mcOrgPi0 = kFALSE;
452 Bool_t mcOrgEta = kFALSE;
453 double mcele = -1.;
454 double mcpT = 0.0;
455 double mcMompT = 0.0;
456 //double mcGrandMompT = 0.0;
457 double mcWeight = -10.0;
458
459 int iHijing = 1;
460
461 if(fmcData && fMC && stack)
462 {
463 Int_t label = TMath::Abs(track->GetLabel());
464
465 Bool_t MChijing = fMC->IsFromBGEvent(label);
466 if(!MChijing)iHijing = 0;
467
468 TParticle* particle = stack->Particle(label);
469 int mcpid = particle->GetPdgCode();
470 mcpT = particle->Pt();
471 //printf("MCpid = %d",mcpid);
472 if(particle->GetFirstMother()>-1)
473 {
474 int parentlabel = particle->GetFirstMother();
475 int parentPID = stack->Particle(particle->GetFirstMother())->GetPdgCode();
476 mcMompT = stack->Particle(particle->GetFirstMother())->Pt();
477 if((parentPID==22 || parentPID==111 || parentPID==221)&& fabs(mcpid)==11)mcPho = kTRUE;
478 if((fabs(parentPID)==411 || fabs(parentPID)==413 || fabs(parentPID)==421 || fabs(parentPID)==423 || fabs(parentPID)==431)&& fabs(mcpid)==11)mcDtoE = kTRUE;
479 if((fabs(parentPID)==511 || fabs(parentPID)==513 || fabs(parentPID)==521 || fabs(parentPID)==523 || fabs(parentPID)==531)&& fabs(mcpid)==11)mcBtoE = kTRUE;
480
481 // make D->e pT correlation
482 if(mcDtoE)fMomDtoE->Fill(mcpT,mcMompT);
483
484 // pi->e (Dalitz)
485 if(parentPID==111 && fabs(mcpid)==11 && mcMompT>0.0)
486 {
487 mcOrgPi0 = kTRUE;
488 mcWeight = GetMCweight(mcMompT);
489 }
490 // eta->e (Dalitz)
491 if(parentPID==221 && fabs(mcpid)==11 && mcMompT>0.0)
492 {
493 mcOrgEta = kTRUE;
494 mcWeight = GetMCweightEta(mcMompT);
495 }
496
497
498 // access grand parent pi0->g->e
499 TParticle* particle_parent = stack->Particle(parentlabel); // get parent pointer
500 if(particle_parent->GetFirstMother()>-1 && parentPID==22 && fabs(mcpid)==11) // get grand parent g->e
501 {
502 int grand_parentPID = stack->Particle(particle_parent->GetFirstMother())->GetPdgCode();
503 double pTtmp = stack->Particle(particle_parent->GetFirstMother())->Pt();
504 if(grand_parentPID==111 && pTtmp>0.0)
505 {
506 mcOrgPi0 = kTRUE;
507 mcWeight = GetMCweight(pTtmp);
508 }
509 if(grand_parentPID==221 && pTtmp>0.0)
510 {
511 mcOrgEta = kTRUE;
512 mcWeight = GetMCweightEta(pTtmp);
513 }
514 }
515 }
516 if(fabs(mcpid)==11 && mcDtoE)mcele= 1.;
517 if(fabs(mcpid)==11 && mcBtoE)mcele= 2.;
518 if(fabs(mcpid)==11 && mcPho)mcele= 3.;
519 }
520
521 //printf("weight = %f\n",mcWeight);
522
523 if(TMath::Abs(track->Eta())>0.6) continue;
524 if(TMath::Abs(track->Pt()<2.0)) continue;
525
526 fTrackPtBefTrkCuts->Fill(track->Pt());
527 // RecKine: ITSTPC cuts
528 if(!ProcessCutStep(AliHFEcuts::kStepRecKineITSTPC, track)) continue;
529
530 //RecKink
531 if(fRejectKinkMother) { // Quick and dirty fix to reject both kink mothers and daughters
532 if(track->GetKinkIndex(0) != 0) continue;
533 }
534
535 // RecPrim
536 if(!ProcessCutStep(AliHFEcuts::kStepRecPrim, track)) continue;
537
538 // HFEcuts: ITS layers cuts
539 if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsITS, track)) continue;
540
541 // HFE cuts: TPC PID cleanup
542 if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTPC, track)) continue;
543
544 int nTPCcl = track->GetTPCNcls();
545 int nTPCclF = track->GetTPCNclsF();
546 int nITS = track->GetNcls(0);
547
548 fTrackPtAftTrkCuts->Fill(track->Pt());
549
550 Double_t mom = -999., eop=-999., pt = -999., dEdx=-999., fTPCnSigma=-10, phi=-999., eta=-999.;
551 pt = track->Pt();
552 if(pt<2.0)continue;
553
554 // Track extrapolation
555
556 Int_t charge = track->Charge();
557 fTrkpt->Fill(pt);
558 mom = track->P();
559 phi = track->Phi();
560 eta = track->Eta();
561 dEdx = track->GetTPCsignal();
562 fTPCnSigma = fPID->GetPIDResponse() ? fPID->GetPIDResponse()->NumberOfSigmasTPC(track, AliPID::kElectron) : 1000;
563
564 double ncells = -1.0;
565 double m20 = -1.0;
566 double m02 = -1.0;
567 double disp = -1.0;
568 double rmatch = -1.0;
569 double nmatch = -1.0;
570 double oppstatus = 0.0;
571
572 Bool_t fFlagPhotonicElec = kFALSE;
573 Bool_t fFlagConvinatElec = kFALSE;
574
575 Int_t clsId = track->GetEMCALcluster();
576 if (clsId>0){
577 AliESDCaloCluster *clust = fESD->GetCaloCluster(clsId);
578 if(clust && clust->IsEMCAL()){
579
580 double clustE = clust->E();
581 eop = clustE/fabs(mom);
582 //double clustT = clust->GetTOF();
583 ncells = clust->GetNCells();
584 m02 = clust->GetM02();
585 m20 = clust->GetM20();
586 disp = clust->GetDispersion();
587 double delphi = clust->GetTrackDx();
588 double deleta = clust->GetTrackDz();
589 rmatch = sqrt(pow(delphi,2)+pow(deleta,2));
590 nmatch = clust->GetNTracksMatched();
591
592 if(fTPCnSigma>-1.5 && fTPCnSigma<3.0)
593 {
594 SelectPhotonicElectron(iTracks,cent,track,stack,fFlagPhotonicElec,fFlagConvinatElec,fTPCnSigma,m20,eop,mcele,mcWeight,iHijing);
595 }
596 if(fFlagPhotonicElec)oppstatus = 1.0;
597 if(fFlagConvinatElec)oppstatus = 2.0;
598 if(fFlagPhotonicElec && fFlagConvinatElec)oppstatus = 3.0;
599
600 double valdedx[16];
601 valdedx[0] = pt; valdedx[1] = nITS; valdedx[2] = phi; valdedx[3] = eta; valdedx[4] = fTPCnSigma;
602 valdedx[5] = eop; valdedx[6] = rmatch; valdedx[7] = ncells, valdedx[8] = nTPCclF; valdedx[9] = m20; valdedx[10] = mcpT;
603 valdedx[11] = cent; valdedx[12] = charge; valdedx[13] = oppstatus; valdedx[14] = nTPCcl;
604 valdedx[15] = mcele;
605 if(fqahist==1)fEleInfo->Fill(valdedx);
606
607
608 }
609 }
610
611 if(nITS<2.5)continue;
612 if(nTPCcl<100)continue;
613
614 CheckNclust->Fill(nTPCcl);
615 CheckNits->Fill(nITS);
616
617 fdEdxBef->Fill(mom,fTPCnSigma);
618 fTPCnsigma->Fill(mom,fTPCnSigma);
619 if(fTPCnSigma >= -1.0 && fTPCnSigma <= 3)fTrkEovPBef->Fill(pt,eop);
620
621 Int_t pidpassed = 1;
622
623
624 //--- track accepted
625 AliHFEpidObject hfetrack;
626 hfetrack.SetAnalysisType(AliHFEpidObject::kESDanalysis);
627 hfetrack.SetRecTrack(track);
628 double binct = 10.5;
629 if((0.0< cent) && (cent<5.0)) binct = 0.5;
630 if((5.0< cent) && (cent<10.0)) binct = 1.5;
631 if((10.0< cent) && (cent<20.0)) binct = 2.5;
632 if((20.0< cent) && (cent<30.0)) binct = 3.5;
633 if((30.0< cent) && (cent<40.0)) binct = 4.5;
634 if((40.0< cent) && (cent<50.0)) binct = 5.5;
635 if((50.0< cent) && (cent<60.0)) binct = 6.5;
636 if((60.0< cent) && (cent<70.0)) binct = 7.5;
637 if((70.0< cent) && (cent<80.0)) binct = 8.5;
638 if((80.0< cent) && (cent<90.0)) binct = 9.5;
639 if((90.0< cent) && (cent<100.0)) binct = 10.5;
640
641 hfetrack.SetCentrality((int)binct); //added
642 hfetrack.SetPbPb();
643 if(!fPID->IsSelected(&hfetrack, NULL, "", fPIDqa)) pidpassed = 0;
644
645 if(pidpassed==0) continue;
646
647 fTrkEovPAft->Fill(pt,eop);
648 fdEdxAft->Fill(mom,fTPCnSigma);
649
650 // Fill real data
651 fIncpT->Fill(cent,pt);
652 if(fFlagPhotonicElec) fPhoElecPt->Fill(cent,pt);
653 if(fFlagConvinatElec) fSameElecPt->Fill(cent,pt);
654
655 if(m20>0.0 && m20<0.3)
656 {
657 fIncpTM20->Fill(cent,pt);
658 if(fFlagPhotonicElec) fPhoElecPtM20->Fill(cent,pt);
659 if(fFlagConvinatElec) fSameElecPtM20->Fill(cent,pt);
660 }
661
662 // MC
663 if(mcele>0) // select MC electrons
664 {
665
666 fIncpTMChfeAll->Fill(cent,pt);
667 if(m20>0.0 && m20<0.3)fIncpTMCM20hfeAll->Fill(cent,pt);
668
669 if(mcBtoE || mcDtoE) // select B->e & D->e
670 {
671 fIncpTMChfe->Fill(cent,pt);
672 if(m20>0.0 && m20<0.3)fIncpTMCM20hfe->Fill(cent,pt);
673 if(m20>0.0 && m20<0.3)fIncpTMCM20hfeCheck->Fill(cent,mcpT);
674 }
675
676 if(mcPho) // select photonic electrons
677 {
678 double phoval[5];
679 phoval[0] = cent;
680 phoval[1] = pt;
681 phoval[2] = fTPCnSigma;
682 phoval[3] = iHijing;
683 phoval[4] = mcMompT;
684
685 fIncpTMCpho->Fill(phoval);
686 if(fFlagPhotonicElec) fPhoElecPtMC->Fill(phoval);
687 if(fFlagConvinatElec) fSameElecPtMC->Fill(phoval);
688
689 if(m20>0.0 && m20<0.3)
690 {
691 fIncpTMCM20pho->Fill(phoval);
692 if(fFlagPhotonicElec) fPhoElecPtMCM20->Fill(phoval);
693 if(fFlagConvinatElec) fSameElecPtMCM20->Fill(phoval);
694 // pi0->g->e
695 if(mcWeight>-1)
696 {
697 if(iHijing==1)mcWeight = 1.0;
698 if(mcOrgPi0)
699 {
700 fIncpTMCM20pho_pi0e->Fill(phoval,mcWeight);
701 if(fFlagPhotonicElec) fPhoElecPtMCM20_pi0e->Fill(phoval,mcWeight);
702 if(fFlagConvinatElec) fSameElecPtMCM20_pi0e->Fill(phoval,mcWeight);
703 }
704 if(mcOrgEta)
705 {
706 fIncpTMCM20pho_eta->Fill(phoval,mcWeight);
707 if(fFlagPhotonicElec) fPhoElecPtMCM20_eta->Fill(phoval,mcWeight);
708 if(fFlagConvinatElec) fSameElecPtMCM20_eta->Fill(phoval,mcWeight);
709 }
710 // --- eta
711 }
712 }
713 }
714 }
715 }
716 PostData(1, fOutputList);
717}
718//_________________________________________
719void AliAnalysisTaskHFECal::UserCreateOutputObjects()
720{
721 //--- Check MC
722
723 //Bool_t mcData = kFALSE;
724 if(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler())
725 {
726 fmcData = kTRUE;
727 printf("+++++ MC Data available");
728 }
729 if(fmcData)
730 {
731 printf("++++++++= MC analysis \n");
732 }
733 else
734 {
735 printf("++++++++ real data analysis \n");
736 }
737
738 printf("+++++++ QA hist %d",fqahist);
739
740 //---- Geometry
741 fGeom = AliEMCALGeometry::GetInstance("EMCAL_COMPLETEV1");
742
743 //--------Initialize PID
744 //fPID->SetHasMCData(kFALSE);
745 fPID->SetHasMCData(fmcData);
746 if(!fPID->GetNumberOfPIDdetectors())
747 {
748 fPID->AddDetector("TPC", 0);
749 fPID->AddDetector("EMCAL", 1);
750 }
751
752 Double_t params[4];
753 const char *cutmodel;
754 cutmodel = "pol0";
755 params[0] = -1.0; //sigma min
756 double maxnSig = 3.0;
757 if(fmcData)
758 {
759 params[0] = -5.0; //sigma min
760 maxnSig = 5.0;
761 }
762
763 for(Int_t a=0;a<11;a++)fPID->ConfigureTPCcentralityCut(a,cutmodel,params,maxnSig);
764
765 fPID->SortDetectors();
766 fPIDqa = new AliHFEpidQAmanager();
767 fPIDqa->Initialize(fPID);
768
769 //------- fcut --------------
770 fCuts = new AliHFEcuts();
771 fCuts->CreateStandardCuts();
772 //fCuts->SetMinNClustersTPC(100);
773 fCuts->SetMinNClustersTPC(90);
774 fCuts->SetMinRatioTPCclusters(0.6);
775 fCuts->SetTPCmodes(AliHFEextraCuts::kFound, AliHFEextraCuts::kFoundOverFindable);
776 //fCuts->SetMinNClustersITS(3);
777 fCuts->SetMinNClustersITS(2);
778 fCuts->SetCutITSpixel(AliHFEextraCuts::kAny);
779 fCuts->SetCheckITSLayerStatus(kFALSE);
780 fCuts->SetVertexRange(10.);
781 fCuts->SetTOFPIDStep(kFALSE);
782 fCuts->SetPtRange(2, 50);
783 fCuts->SetMaxImpactParam(3.,3.);
784
785 //--------Initialize correction Framework and Cuts
786 fCFM = new AliCFManager;
787 const Int_t kNcutSteps = AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kNcutStepsRecTrack + AliHFEcuts::kNcutStepsDETrack;
788 fCFM->SetNStepParticle(kNcutSteps);
789 for(Int_t istep = 0; istep < kNcutSteps; istep++)
790 fCFM->SetParticleCutsList(istep, NULL);
791
792 if(!fCuts){
793 AliWarning("Cuts not available. Default cuts will be used");
794 fCuts = new AliHFEcuts;
795 fCuts->CreateStandardCuts();
796 }
797 fCuts->Initialize(fCFM);
798
799 //---------Output Tlist
800 fOutputList = new TList();
801 fOutputList->SetOwner();
802 fOutputList->Add(fPIDqa->MakeList("PIDQA"));
803
804 fNoEvents = new TH1F("fNoEvents","",4,-0.5,3.5) ;
805 fOutputList->Add(fNoEvents);
806
807 Int_t binsE[5] = {250, 100, 1000, 200, 10};
808 Double_t xminE[5] = {1.0, -1, 0.0, 0, -0.5};
809 Double_t xmaxE[5] = {3.5, 1, 100.0, 100, 9.5};
810 fEMCAccE = new THnSparseD("fEMCAccE","EMC acceptance & E;#eta;#phi;Energy;Centrality;trugCondition;",5,binsE,xminE,xmaxE);
811 if(fqahist==1)fOutputList->Add(fEMCAccE);
812
813 fTrkpt = new TH1F("fTrkpt","track pt",100,0,50);
814 fOutputList->Add(fTrkpt);
815
816 fTrackPtBefTrkCuts = new TH1F("fTrackPtBefTrkCuts","track pt before track cuts",100,0,50);
817 fOutputList->Add(fTrackPtBefTrkCuts);
818
819 fTrackPtAftTrkCuts = new TH1F("fTrackPtAftTrkCuts","track pt after track cuts",100,0,50);
820 fOutputList->Add(fTrackPtAftTrkCuts);
821
822 fTPCnsigma = new TH2F("fTPCnsigma", "TPC - n sigma",100,0,50,200,-10,10);
823 fOutputList->Add(fTPCnsigma);
824
825 fTrkEovPBef = new TH2F("fTrkEovPBef","track E/p before HFE pid",100,0,50,100,0,2);
826 fOutputList->Add(fTrkEovPBef);
827
828 fTrkEovPAft = new TH2F("fTrkEovPAft","track E/p after HFE pid",100,0,50,100,0,2);
829 fOutputList->Add(fTrkEovPAft);
830
831 fdEdxBef = new TH2F("fdEdxBef","track dEdx vs p before HFE pid",100,0,50,200,-10,10);
832 fOutputList->Add(fdEdxBef);
833
834 fdEdxAft = new TH2F("fdEdxAft","track dEdx vs p after HFE pid",100,0,50,200,-10,10);
835 fOutputList->Add(fdEdxAft);
836
837 fIncpT = new TH2F("fIncpT","HFE pid electro vs. centrality",200,0,100,100,0,50);
838 fOutputList->Add(fIncpT);
839
840 fIncpTM20 = new TH2F("fIncpTM20","HFE pid electro vs. centrality with M20",200,0,100,100,0,50);
841 fOutputList->Add(fIncpTM20);
842
843 Int_t nBinspho[9] = { 200, 100, 500, 12, 50, 4, 200, 8, 100};
844 Double_t minpho[9] = { 0., 0., 0., -2.5, 0, -0.5, 0,-1.5, 0};
845 Double_t maxpho[9] = {100., 50., 0.5, 3.5, 1, 3.5, 2, 6.5, 50};
846
847 fInvmassLS = new THnSparseD("fInvmassLS", "Inv mass of LS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2); nSigma; angle; m20cut; eop; Mcele;", 9, nBinspho,minpho, maxpho);
848 if(fqahist==1)fOutputList->Add(fInvmassLS);
849
850 fInvmassULS = new THnSparseD("fInvmassULS", "Inv mass of ULS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2); nSigma; angle; m20cut; eop; MCele", 9, nBinspho,minpho, maxpho);
851 if(fqahist==1)fOutputList->Add(fInvmassULS);
852
853 fInvmassLSmc = new THnSparseD("fInvmassLSmc", "Inv mass of LS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2); nSigma; angle; m20cut; eop; Mcele;", 9, nBinspho,minpho, maxpho);
854 if(fqahist==1)fOutputList->Add(fInvmassLSmc);
855
856 fInvmassULSmc = new THnSparseD("fInvmassULSmc", "Inv mass of ULS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2); nSigma; angle; m20cut; eop; MCele", 9, nBinspho,minpho, maxpho);
857 if(fqahist==1)fOutputList->Add(fInvmassULSmc);
858
859 fInvmassLSmc0 = new TH2D("fInvmassLSmc0", "Inv mass of LS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2)",20,0,20,500,0,0.5 );
860 fInvmassLSmc0->Sumw2();
861 fOutputList->Add(fInvmassLSmc0);
862
863 fInvmassLSmc1 = new TH2D("fInvmassLSmc1", "Inv mass of LS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2)",20,0,20,500,0,0.5 );
864 fInvmassLSmc1->Sumw2();
865 fOutputList->Add(fInvmassLSmc1);
866
867 fInvmassLSmc2 = new TH2D("fInvmassLSmc2", "Inv mass of LS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2)",20,0,20,500,0,0.5 );
868 fInvmassLSmc2->Sumw2();
869 fOutputList->Add(fInvmassLSmc2);
870
871 fInvmassLSmc3 = new TH2D("fInvmassLSmc3", "Inv mass of LS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2)",20,0,20,500,0,0.5 );
872 fInvmassLSmc3->Sumw2();
873 fOutputList->Add(fInvmassLSmc3);
874
875 fInvmassULSmc0 = new TH2D("fInvmassULSmc0", "Inv mass of ULS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2)",20,0,20,500,0,0.5 );
876 fInvmassULSmc0->Sumw2();
877 fOutputList->Add(fInvmassULSmc0);
878
879 fInvmassULSmc1 = new TH2D("fInvmassULSmc1", "Inv mass of ULS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2)",20,0,20,500,0,0.5 );
880 fInvmassULSmc1->Sumw2();
881 fOutputList->Add(fInvmassULSmc1);
882
883 fInvmassULSmc2 = new TH2D("fInvmassULSmc2", "Inv mass of ULS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2)",20,0,20,500,0,0.5 );
884 fInvmassULSmc2->Sumw2();
885 fOutputList->Add(fInvmassULSmc2);
886
887 fInvmassULSmc3 = new TH2D("fInvmassULSmc3", "Inv mass of ULS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2)",20,0,20,500,0,0.5 );
888 fInvmassULSmc3->Sumw2();
889 fOutputList->Add(fInvmassULSmc3);
890
891 fOpeningAngleLS = new TH1F("fOpeningAngleLS","Opening angle for LS pairs",100,0,1);
892 fOutputList->Add(fOpeningAngleLS);
893
894 fOpeningAngleULS = new TH1F("fOpeningAngleULS","Opening angle for ULS pairs",100,0,1);
895 fOutputList->Add(fOpeningAngleULS);
896
897 fPhotoElecPt = new TH1F("fPhotoElecPt", "photonic electron pt",100,0,50);
898 fOutputList->Add(fPhotoElecPt);
899
900 fPhoElecPt = new TH2F("fPhoElecPt", "Pho-inclusive electron pt",200,0,100,100,0,50);
901 fOutputList->Add(fPhoElecPt);
902
903 fPhoElecPtM20 = new TH2F("fPhoElecPtM20", "Pho-inclusive electron pt with M20",200,0,100,100,0,50);
904 fOutputList->Add(fPhoElecPtM20);
905
906 fSameElecPt = new TH2F("fSameElecPt", "Same-inclusive electron pt",200,0,100,100,0,50);
907 fOutputList->Add(fSameElecPt);
908
909 fSameElecPtM20 = new TH2F("fSameElecPtM20", "Same-inclusive electron pt with M20",200,0,100,100,0,50);
910 fOutputList->Add(fSameElecPtM20);
911
912 fCent = new TH1F("fCent","Centrality",200,0,100) ;
913 fOutputList->Add(fCent);
914
915 // Make common binning
916 const Double_t kMinP = 0.;
917 const Double_t kMaxP = 50.;
918 //const Double_t kTPCSigMim = 40.;
919 //const Double_t kTPCSigMax = 140.;
920
921 // 1st histogram: TPC dEdx with/without EMCAL (p, pT, TPC Signal, phi, eta, Sig, e/p, ,match, cell, M02, M20, Disp, Centrality, select)
922 Int_t nBins[16] = { 250, 10, 60, 20, 100, 300, 50, 40, 200, 200, 250, 200, 3, 5, 100, 8};
923 Double_t min[16] = {kMinP, -0.5, 1.0, -1.0, -6.0, 0, 0, 0, 0.0, 0.0, 0.0, 0, -1.5, -0.5, 80, -1.5};
924 Double_t max[16] = {kMaxP, 9.5, 4.0, 1.0, 4.0, 3.0, 0.1, 40, 200, 2.0, 50.0, 100, 1.5, 4.5, 180, 6.5};
925 fEleInfo = new THnSparseD("fEleInfo", "Electron Info; pT [GeV/c]; TPC signal;phi;eta;nSig; E/p;Rmatch;Ncell;clsF;M20;mcpT;Centrality;charge;opp;same;trigCond;MCele", 16, nBins, min, max);
926 if(fqahist==1)fOutputList->Add(fEleInfo);
927
928 //<--- Trigger info
929 /*
930 fClsEBftTrigCut = new TH1F("fClsEBftTrigCut","cluster E before trigger selection",1000,0,100);
931 fOutputList->Add(fClsEBftTrigCut);
932
933 fClsEAftTrigCut = new TH1F("fClsEAftTrigCut","cluster E if cls has 0 trigcut channel",1000,0,100);
934 fOutputList->Add(fClsEAftTrigCut);
935
936 fClsEAftTrigCut1 = new TH1F("fClsEAftTrigCut1","cluster E if cls with trig channel",1000,0,100);
937 fOutputList->Add(fClsEAftTrigCut1);
938
939 fClsEAftTrigCut2 = new TH1F("fClsEAftTrigCut2","cluster E if cls with trigcut channel",1000,0,100);
940 fOutputList->Add(fClsEAftTrigCut2);
941
942 fClsEAftTrigCut3 = new TH1F("fClsEAftTrigCut3","cluster E if cls with trigcut channel + nCell>Ecorrect",1000,0,100);
943 fOutputList->Add(fClsEAftTrigCut3);
944
945 fClsEAftTrigCut4 = new TH1F("fClsEAftTrigCut4","cluster E if cls with trigcut channel + nCell>Ecorrect + cls time cut",1000,0,100);
946 fOutputList->Add(fClsEAftTrigCut4);
947
948 fClsETime = new TH2F("fClsETime", "Cls time vs E; E; time;",1000,0,100,1000,-0.0000002,0.0000009);
949 fOutputList->Add(fClsETime);
950
951 fClsETime1 = new TH2F("fClsETime1", "Cls time vs E if cls contains trigger channel; E; time;",1000,0,100,1000,-0.0000002,0.0000009);
952 fOutputList->Add(fClsETime1);
953
954 fTrigTimes = new TH1F("fTrigTimes", "Trigger time; time; N;",25,0,25);
955 fOutputList->Add(fTrigTimes);
956
957 fCellCheck = new TH2F("fCellCheck", "Cell vs E; E GeV; Cell ID",10,6,26,12000,0,12000);
958 fOutputList->Add(fCellCheck);
959 */
960 //<---------- MC
961
962 fInputHFEMC = new TH2F("fInputHFEMC","Input MC HFE pid electro vs. centrality",200,0,100,100,0,50);
963 fOutputList->Add(fInputHFEMC);
964
965 fInputAlle = new TH2F("fInputAlle","Input MC electro vs. centrality",200,0,100,100,0,50);
966 fOutputList->Add(fInputAlle);
967
968 fIncpTMChfe = new TH2F("fIncpTMChfe","MC HFE pid electro vs. centrality",200,0,100,100,0,50);
969 fOutputList->Add(fIncpTMChfe);
970
971 fIncpTMChfeAll = new TH2F("fIncpTMChfeAll","MC Alle pid electro vs. centrality",200,0,100,100,0,50);
972 fOutputList->Add(fIncpTMChfeAll);
973
974 fIncpTMCM20hfe = new TH2F("fIncpTMCM20hfe","MC HFE pid electro vs. centrality with M20",200,0,100,100,0,50);
975 fOutputList->Add(fIncpTMCM20hfe);
976
977 fIncpTMCM20hfeAll = new TH2F("fIncpTMCM20hfeAll","MC Alle pid electro vs. centrality with M20",200,0,100,100,0,50);
978 fOutputList->Add(fIncpTMCM20hfeAll);
979
980 fIncpTMCM20hfeCheck = new TH2F("fIncpTMCM20hfeCheck","MC HFE pid electro vs. centrality with M20 Check",200,0,100,100,0,50);
981 fOutputList->Add(fIncpTMCM20hfeCheck);
982
983 Int_t nBinspho2[5] = { 200, 100, 7, 3, 100};
984 Double_t minpho2[5] = { 0., 0., -2.5, -0.5, 0.};
985 Double_t maxpho2[5] = {100., 50., 4.5, 2.5, 50.};
986
987 fIncpTMCpho = new THnSparseD("fIncpTMCpho","MC Pho pid electro vs. centrality",5,nBinspho2,minpho2,maxpho2);
988 fOutputList->Add(fIncpTMCpho);
989
990 fIncpTMCM20pho = new THnSparseD("fIncpTMCM20pho","MC Pho pid electro vs. centrality with M20",5,nBinspho2,minpho2,maxpho2);
991 fOutputList->Add(fIncpTMCM20pho);
992
993 fPhoElecPtMC = new THnSparseD("fPhoElecPtMC", "MC Pho-inclusive electron pt",5,nBinspho2,minpho2,maxpho2);
994 fOutputList->Add(fPhoElecPtMC);
995
996 fPhoElecPtMCM20 = new THnSparseD("fPhoElecPtMCM20", "MC Pho-inclusive electron pt with M20",5,nBinspho2,minpho2,maxpho2);
997 fOutputList->Add(fPhoElecPtMCM20);
998
999 fSameElecPtMC = new THnSparseD("fSameElecPtMC", "MC Same-inclusive electron pt",5,nBinspho2,minpho2,maxpho2);
1000 fOutputList->Add(fSameElecPtMC);
1001
1002 fSameElecPtMCM20 = new THnSparseD("fSameElecPtMCM20", "MC Same-inclusive electron pt with M20",5,nBinspho2,minpho2,maxpho2);
1003 fOutputList->Add(fSameElecPtMCM20);
1004
1005 fIncpTMCM20pho_pi0e = new THnSparseD("fIncpTMCM20pho_pi0e","MC Pho pi0->e pid electro vs. centrality with M20",5,nBinspho2,minpho2,maxpho2);
1006 fIncpTMCM20pho_pi0e->Sumw2();
1007 fOutputList->Add(fIncpTMCM20pho_pi0e);
1008
1009 fPhoElecPtMCM20_pi0e = new THnSparseD("fPhoElecPtMCM20_pi0e", "MC Pho-inclusive electron pt with M20 pi0->e",5,nBinspho2,minpho2,maxpho2);
1010 fPhoElecPtMCM20_pi0e->Sumw2();
1011 fOutputList->Add(fPhoElecPtMCM20_pi0e);
1012
1013 fSameElecPtMCM20_pi0e = new THnSparseD("fSameElecPtMCM20_pi0e", "MC Same-inclusive electron pt pi0->e",5,nBinspho2,minpho2,maxpho2);
1014 fSameElecPtMCM20_pi0e->Sumw2();
1015 fOutputList->Add(fSameElecPtMCM20_pi0e);
1016 //
1017 fIncpTMCM20pho_eta = new THnSparseD("fIncpTMCM20pho_eta","MC Pho pi0->e pid electro vs. centrality with M20",5,nBinspho2,minpho2,maxpho2);
1018 fIncpTMCM20pho_eta->Sumw2();
1019 fOutputList->Add(fIncpTMCM20pho_eta);
1020
1021 fPhoElecPtMCM20_eta = new THnSparseD("fPhoElecPtMCM20_eta", "MC Pho-inclusive electron pt with M20 pi0->e",5,nBinspho2,minpho2,maxpho2);
1022 fPhoElecPtMCM20_eta->Sumw2();
1023 fOutputList->Add(fPhoElecPtMCM20_eta);
1024
1025 fSameElecPtMCM20_eta = new THnSparseD("fSameElecPtMCM20_eta", "MC Same-inclusive electron pt pi0->e",5,nBinspho2,minpho2,maxpho2);
1026 fSameElecPtMCM20_eta->Sumw2();
1027 fOutputList->Add(fSameElecPtMCM20_eta);
1028
1029
1030 CheckNclust = new TH1D("CheckNclust","cluster check",200,0,200);
1031 fOutputList->Add(CheckNclust);
1032
1033 CheckNits = new TH1D("CheckNits","ITS cluster check",8,-0.5,7.5);
1034 fOutputList->Add(CheckNits);
1035
1036 Hpi0pTcheck = new TH2D("Hpi0pTcheck","Pi0 pT from Hijing",100,0,50,3,-0.5,2.5);
1037 fOutputList->Add(Hpi0pTcheck);
1038
1039 HETApTcheck = new TH2D("HETApTcheck","Eta pT from Hijing",100,0,50,3,-0.5,2.5);
1040 fOutputList->Add(HETApTcheck);
1041
1042 HphopTcheck = new TH2D("HphopTcheck","Pho pT from Hijing",100,0,50,3,-0.5,2.5);
1043 fOutputList->Add(HphopTcheck);
1044
1045 fpTCheck = new TH1D("fpTCheck","pT check",500,0,50);
1046 fOutputList->Add(fpTCheck);
1047
1048 fMomDtoE = new TH2D("fMomDtoE","D->E pT correlations;e p_{T} GeV/c;D p_{T} GeV/c",400,0,40,400,0,40);
1049 fOutputList->Add(fMomDtoE);
1050
1051 PostData(1,fOutputList);
1052}
1053
1054//________________________________________________________________________
1055void AliAnalysisTaskHFECal::Terminate(Option_t *)
1056{
1057 // Info("Terminate");
1058 AliAnalysisTaskSE::Terminate();
1059}
1060
1061//________________________________________________________________________
1062Bool_t AliAnalysisTaskHFECal::ProcessCutStep(Int_t cutStep, AliVParticle *track)
1063{
1064 // Check single track cuts for a given cut step
1065 const Int_t kMCOffset = AliHFEcuts::kNcutStepsMCTrack;
1066 if(!fCFM->CheckParticleCuts(cutStep + kMCOffset, track)) return kFALSE;
1067 return kTRUE;
1068}
1069//_________________________________________
1070//void AliAnalysisTaskHFECal::SelectPhotonicElectron(Int_t itrack, Double_t cent, AliESDtrack *track, Bool_t &fFlagPhotonicElec)
1071//void AliAnalysisTaskHFECal::SelectPhotonicElectron(Int_t itrack, Double_t cent, AliESDtrack *track, Bool_t &fFlagPhotonicElec, Bool_t &fFlagConvinatElec, Double_t nSig)
1072void AliAnalysisTaskHFECal::SelectPhotonicElectron(Int_t itrack, Double_t cent, AliESDtrack *track, AliStack *stack ,Bool_t &fFlagPhotonicElec, Bool_t &fFlagConvinatElec, Double_t nSig, Double_t shower, Double_t ep, Double_t mce, Double_t w, Int_t ibgevent)
1073{
1074 //Identify non-heavy flavour electrons using Invariant mass method
1075
1076 fTrackCuts->SetAcceptKinkDaughters(kFALSE);
1077 fTrackCuts->SetRequireTPCRefit(kTRUE);
1078 fTrackCuts->SetRequireITSRefit(kTRUE);
1079 fTrackCuts->SetEtaRange(-0.9,0.9);
1080 //fTrackCuts->SetRequireSigmaToVertex(kTRUE);
1081 fTrackCuts->SetMaxChi2PerClusterTPC(3.5);
1082 fTrackCuts->SetMinNClustersTPC(90);
1083
1084 const AliESDVertex *pVtx = fESD->GetPrimaryVertex();
1085
1086 double ptEle = track->Pt(); //add
1087 if(ibgevent==0 && w > 0.0)
1088 {
1089 fpTCheck->Fill(ptEle,w);
1090 }
1091
1092 Bool_t flagPhotonicElec = kFALSE;
1093 Bool_t flagConvinatElec = kFALSE;
1094
1095 int p1 = 0;
1096 if(mce==3)
1097 {
1098 Int_t label = TMath::Abs(track->GetLabel());
1099 TParticle* particle = stack->Particle(label);
1100 p1 = particle->GetFirstMother();
1101 }
1102
1103 //for(Int_t jTracks = itrack+1; jTracks<fESD->GetNumberOfTracks(); jTracks++){
1104 for(Int_t jTracks = 0; jTracks<fESD->GetNumberOfTracks(); jTracks++){
1105 AliESDtrack* trackAsso = fESD->GetTrack(jTracks);
1106 if (!trackAsso) {
1107 printf("ERROR: Could not receive track %d\n", jTracks);
1108 continue;
1109 }
1110 if(itrack==jTracks)continue;
1111 int jbgevent = 0;
1112
1113 int p2 = 0;
1114 if(mce==3)
1115 {
1116 Int_t label2 = TMath::Abs(trackAsso->GetLabel());
1117 TParticle* particle2 = stack->Particle(label2);
1118 Bool_t MChijing_ass = fMC->IsFromBGEvent(label2);
1119 if(MChijing_ass)jbgevent =1;
1120 if(particle2->GetFirstMother()>-1)
1121 p2 = particle2->GetFirstMother();
1122 }
1123
1124 Double_t dEdxAsso = -999., ptPrim=-999., ptAsso=-999., openingAngle = -999.;
1125 Double_t mass=999., width = -999;
1126 Bool_t fFlagLS=kFALSE, fFlagULS=kFALSE;
1127
1128 //ptPrim = track->Pt();
1129 ptPrim = ptEle;
1130
1131 dEdxAsso = trackAsso->GetTPCsignal();
1132 ptAsso = trackAsso->Pt();
1133 Int_t chargeAsso = trackAsso->Charge();
1134 Int_t charge = track->Charge();
1135
1136
1137 if(ptAsso <0.5) continue;
1138 if(!fTrackCuts->AcceptTrack(trackAsso)) continue;
1139 if(dEdxAsso <65 || dEdxAsso>100) continue; //11a pass1
1140
1141 Int_t fPDGe1 = 11; Int_t fPDGe2 = 11;
1142 if(charge>0) fPDGe1 = -11;
1143 if(chargeAsso>0) fPDGe2 = -11;
1144
1145 //printf("chargeAsso = %d\n",chargeAsso);
1146 //printf("charge = %d\n",charge);
1147 if(charge == chargeAsso) fFlagLS = kTRUE;
1148 if(charge != chargeAsso) fFlagULS = kTRUE;
1149
1150 //printf("fFlagLS = %d\n",fFlagLS);
1151 //printf("fFlagULS = %d\n",fFlagULS);
1152 //printf("\n");
1153
1154 AliKFParticle ge1(*track, fPDGe1);
1155 AliKFParticle ge2(*trackAsso, fPDGe2);
1156 AliKFParticle recg(ge1, ge2);
1157
1158 if(recg.GetNDF()<1) continue;
1159 Double_t chi2recg = recg.GetChi2()/recg.GetNDF();
1160 if(TMath::Sqrt(TMath::Abs(chi2recg))>3.) continue;
1161
1162 // for v5-03-70-AN comment
1163 AliKFVertex primV(*pVtx);
1164 primV += recg;
1165 recg.SetProductionVertex(primV);
1166
1167 recg.SetMassConstraint(0,0.0001);
1168
1169 openingAngle = ge1.GetAngle(ge2);
1170 if(fFlagLS) fOpeningAngleLS->Fill(openingAngle);
1171 if(fFlagULS) fOpeningAngleULS->Fill(openingAngle);
1172
1173
1174 recg.GetMass(mass,width);
1175
1176 double ishower = 0;
1177 if(shower>0.0 && shower<0.3)ishower = 1;
1178
1179 double phoinfo[9];
1180 phoinfo[0] = cent;
1181 phoinfo[1] = ptPrim;
1182 phoinfo[2] = mass;
1183 phoinfo[3] = nSig;
1184 phoinfo[4] = openingAngle;
1185 phoinfo[5] = ishower;
1186 phoinfo[6] = ep;
1187 phoinfo[7] = mce;
1188 phoinfo[8] = ptAsso;
1189
1190 if(fFlagLS) fInvmassLS->Fill(phoinfo);
1191 if(fFlagULS) fInvmassULS->Fill(phoinfo);
1192 if(fFlagLS && ibgevent==0 && jbgevent==0) fInvmassLSmc->Fill(phoinfo,w);
1193 if(fFlagULS && ibgevent==0 && jbgevent==0)
1194 {
1195 fInvmassULSmc->Fill(phoinfo,w);
1196 }
1197 //printf("fInvmassCut %f\n",fInvmassCut);
1198 //printf("openingAngle %f\n",fOpeningAngleCut);
1199
1200 if(openingAngle > fOpeningAngleCut) continue;
1201
1202 // for real data
1203 //printf("mce =%f\n",mce);
1204 if(mce<-0.5) // mce==-1. is real
1205 {
1206 //printf("Real data\n");
1207 if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec){
1208 //if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec && (p1==p2)){ <--- only MC train (55,56) v5-03-68-AN & 69 for check
1209 flagPhotonicElec = kTRUE;
1210 }
1211 if(mass<fInvmassCut && fFlagLS && !flagConvinatElec){
1212 flagConvinatElec = kTRUE;
1213 }
1214 }
1215 // for MC data
1216 else
1217 {
1218 //printf("MC data\n");
1219
1220 if(w>0.0)
1221 {
1222 if(fFlagLS && ibgevent==0 && jbgevent==0) fInvmassLSmc0->Fill(ptPrim,mass,w);
1223 if(fFlagULS && ibgevent==0 && jbgevent==0) fInvmassULSmc0->Fill(ptPrim,mass,w);
1224 if(fFlagLS && ibgevent==0 && jbgevent==0) fInvmassLSmc1->Fill(ptPrim,mass);
1225 if(fFlagULS && ibgevent==0 && jbgevent==0) fInvmassULSmc1->Fill(ptPrim,mass);
1226 if(fFlagLS && ibgevent==0 && jbgevent==0 && (p1==p2)) fInvmassLSmc2->Fill(ptPrim,mass,w);
1227 if(fFlagULS && ibgevent==0 && jbgevent==0 && (p1==p2)) fInvmassULSmc2->Fill(ptPrim,mass,w);
1228 if(fFlagLS && ibgevent==0 && jbgevent==0 && (p1==p2)) fInvmassLSmc3->Fill(ptPrim,mass);
1229 if(fFlagULS && ibgevent==0 && jbgevent==0 && (p1==p2)) fInvmassULSmc3->Fill(ptPrim,mass);
1230 }
1231 if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec && (ibgevent==jbgevent)){
1232 //if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec && (p1==p2)){ <--- only MC train (55,56) v5-03-68-AN & 69 for check
1233 flagPhotonicElec = kTRUE;
1234 }
1235 if(mass<fInvmassCut && fFlagLS && !flagConvinatElec && (ibgevent==jbgevent)){
1236 flagConvinatElec = kTRUE;
1237 }
1238 }
1239
1240 }
1241 fFlagPhotonicElec = flagPhotonicElec;
1242 fFlagConvinatElec = flagConvinatElec;
1243
1244}
1245//-------------------------------------------
1246//void AliAnalysisTaskHFECal::GetMCweight(double mcPi0pT)
1247double AliAnalysisTaskHFECal::GetMCweight(double mcPi0pT)
1248{
1249 double weight = 1.0;
1250
1251 if(mcPi0pT>0.0 && mcPi0pT<5.0)
1252 {
1253 weight = 0.323*mcPi0pT/(TMath::Exp(-1.6+0.767*mcPi0pT+0.0285*mcPi0pT*mcPi0pT));
1254 }
1255 else
1256 {
1257 weight = 115.0/(0.718*mcPi0pT*TMath::Power(mcPi0pT,3.65));
1258 }
1259 return weight;
1260}
1261
1262double AliAnalysisTaskHFECal::GetMCweightEta(double mcEtapT)
1263{
1264 double weight = 1.0;
1265
1266 weight = 223.3/TMath::Power((TMath::Exp(-0.17*mcEtapT-0.0322*mcEtapT*mcEtapT)+mcEtapT/1.69),5.65);
1267 return weight;
1268}
1269
1270
1271//_________________________________________
1272void AliAnalysisTaskHFECal::FindTriggerClusters()
1273{
1274 // constants
1275 const int nModuleCols = 2;
1276 const int nModuleRows = 5;
1277 const int nColsFeeModule = 48;
1278 const int nRowsFeeModule = 24;
1279 const int nColsFaltroModule = 24;
1280 const int nRowsFaltroModule = 12;
1281 //const int faltroWidthMax = 20;
1282
1283 // part 1, trigger extraction -------------------------------------
1284 Int_t globCol, globRow;
1285 //Int_t ntimes=0, nTrigChannel=0, nTrigChannelCut=0, trigInCut=0;
1286 Int_t ntimes=0, nTrigChannel=0, nTrigChannelCut=0;
1287
1288 //Int_t trigtimes[faltroWidthMax];
1289 Double_t cellTime[nColsFeeModule*nModuleCols][nRowsFeeModule*nModuleRows];
1290 Double_t cellEnergy[nColsFeeModule*nModuleCols][nRowsFeeModule*nModuleRows];
1291 //Double_t fTrigCutLow = 6;
1292 //Double_t fTrigCutHigh = 10;
1293 Double_t fTimeCutLow = 469e-09;
1294 Double_t fTimeCutHigh = 715e-09;
1295
1296 AliESDCaloTrigger * fCaloTrigger = fESD->GetCaloTrigger( "EMCAL" );
1297
1298 // erase trigger maps
1299 for(Int_t i = 0; i < nColsFaltroModule*nModuleCols; i++ )
1300 {
1301 for(Int_t j = 0; j < nRowsFaltroModule*nModuleRows; j++ )
1302 {
1303 ftriggersCut[i][j] = 0;
1304 ftriggers[i][j] = 0;
1305 ftriggersTime[i][j] = 0;
1306 }
1307 }
1308
1309 Int_t iglobCol=0, iglobRow=0;
1310 // go through triggers
1311 if( fCaloTrigger->GetEntries() > 0 )
1312 {
1313 // needs reset
1314 fCaloTrigger->Reset();
1315 while( fCaloTrigger->Next() )
1316 {
1317 fCaloTrigger->GetPosition( globCol, globRow );
1318 fCaloTrigger->GetNL0Times( ntimes );
1319 /*
1320 // no L0s
1321 if( ntimes < 1 ) continue;
1322 // get precise timings
1323 fCaloTrigger->GetL0Times( trigtimes );
1324 trigInCut = 0;
1325 for(Int_t i = 0; i < ntimes; i++ )
1326 {
1327 // save the first trigger time in channel
1328 if( i == 0 || triggersTime[globCol][globRow] > trigtimes[i] )
1329 triggersTime[globCol][globRow] = trigtimes[i];
1330 //printf("trigger times: %d\n",trigtimes[i]);
1331 // check if in cut
1332 if(trigtimes[i] > fTrigCutLow && trigtimes[i] < fTrigCutHigh )
1333 trigInCut = 1;
1334
1335 fTrigTimes->Fill(trigtimes[i]);
1336 }
1337 */
1338
1339 //L1 analysis from AliAnalysisTaskEMCALTriggerQA
1340 Int_t bit = 0;
1341 fCaloTrigger->GetTriggerBits(bit);
1342
1343 Int_t ts = 0;
1344 fCaloTrigger->GetL1TimeSum(ts);
1345 if (ts > 0)ftriggers[globCol][globRow] = 1;
1346 // number of triggered channels in event
1347 nTrigChannel++;
1348 // ... inside cut
1349 if(ts>0 && (bit >> 6 & 0x1))
1350 {
1351 iglobCol = globCol;
1352 iglobRow = globRow;
1353 nTrigChannelCut++;
1354 ftriggersCut[globCol][globRow] = 1;
1355 }
1356
1357 } // calo trigger entries
1358 } // has calo trigger entries
1359
1360 // part 2 go through the clusters here -----------------------------------
1361 Int_t nCluster=0, nCell=0, iCell=0, gCell=0;
1362 Short_t cellAddr, nSACell, mclabel;
1363 //Int_t nSACell, iSACell, mclabel;
1364 Int_t iSACell;
1365 Double_t cellAmp=0, cellTimeT=0, clusterTime=0, efrac=0;
1366 Int_t nSupMod, nModule, nIphi, nIeta, iphi, ieta, gphi, geta, feta, fphi;
1367
1368 TRefArray *fCaloClusters = new TRefArray();
1369 fESD->GetEMCALClusters( fCaloClusters );
1370 nCluster = fCaloClusters->GetEntries();
1371
1372
1373 // save all cells times if there are clusters
1374 if( nCluster > 0 ){
1375 // erase time map
1376 for(Int_t i = 0; i < nColsFeeModule*nModuleCols; i++ ){
1377 for(Int_t j = 0; j < nRowsFeeModule*nModuleRows; j++ ){
1378 cellTime[i][j] = 0.;
1379 cellEnergy[i][j] = 0.;
1380 }
1381 }
1382
1383 // get cells
1384 AliESDCaloCells *fCaloCells = fESD->GetEMCALCells();
1385 //AliVCaloCells fCaloCells = fESD->GetEMCALCells();
1386 nSACell = fCaloCells->GetNumberOfCells();
1387 for(iSACell = 0; iSACell < nSACell; iSACell++ ){
1388 // get the cell info *fCal
1389 fCaloCells->GetCell( iSACell, cellAddr, cellAmp, cellTimeT , mclabel, efrac);
1390
1391 // get cell position
1392 fGeom->GetCellIndex( cellAddr, nSupMod, nModule, nIphi, nIeta );
1393 fGeom->GetCellPhiEtaIndexInSModule( nSupMod,nModule, nIphi, nIeta, iphi, ieta);
1394
1395 // convert co global phi eta
1396 gphi = iphi + nRowsFeeModule*(nSupMod/2);
1397 geta = ieta + nColsFeeModule*(nSupMod%2);
1398
1399 // save cell time and energy
1400 cellTime[geta][gphi] = cellTimeT;
1401 cellEnergy[geta][gphi] = cellAmp;
1402
1403 }
1404 }
1405
1406 Int_t nClusterTrig, nClusterTrigCut;
1407 UShort_t *cellAddrs;
1408 Double_t clsE=-999, clsEta=-999, clsPhi=-999;
1409 Float_t clsPos[3] = {0.,0.,0.};
1410
1411 for(Int_t icl=0; icl<fESD->GetNumberOfCaloClusters(); icl++)
1412 {
1413 AliESDCaloCluster *cluster = fESD->GetCaloCluster(icl);
1414 if(!cluster || !cluster->IsEMCAL()) continue;
1415
1416 // get cluster cells
1417 nCell = cluster->GetNCells();
1418
1419 // get cluster energy
1420 clsE = cluster->E();
1421
1422 // get cluster position
1423 cluster->GetPosition(clsPos);
1424 TVector3 clsPosVec(clsPos[0],clsPos[1],clsPos[2]);
1425 clsEta = clsPosVec.Eta();
1426 clsPhi = clsPosVec.Phi();
1427
1428 // get the cell addresses
1429 cellAddrs = cluster->GetCellsAbsId();
1430
1431 // check if the cluster contains cell, that was marked as triggered
1432 nClusterTrig = 0;
1433 nClusterTrigCut = 0;
1434
1435 // loop the cells to check, if cluser in acceptance
1436 // any cluster with a cell outside acceptance is not considered
1437 for( iCell = 0; iCell < nCell; iCell++ )
1438 {
1439 // check hot cell
1440 //if(clsE>6.0)fCellCheck->Fill(clsE,cellAddrs[iCell]);
1441
1442 // get cell position
1443 fGeom->GetCellIndex( cellAddrs[iCell], nSupMod, nModule, nIphi, nIeta );
1444 fGeom->GetCellPhiEtaIndexInSModule( nSupMod,nModule, nIphi, nIeta, iphi, ieta);
1445
1446 // convert co global phi eta
1447 gphi = iphi + nRowsFeeModule*(nSupMod/2);
1448 geta = ieta + nColsFeeModule*(nSupMod%2);
1449
1450 if( cellTime[geta][gphi] > 0. ){
1451 clusterTime += cellTime[geta][gphi];
1452 gCell++;
1453 }
1454
1455 // get corresponding FALTRO
1456 fphi = gphi / 2;
1457 feta = geta / 2;
1458
1459 // try to match with a triggered
1460 if( ftriggers[feta][fphi]==1)
1461 { nClusterTrig++;
1462 }
1463 if( ftriggersCut[feta][fphi]==1)
1464 { nClusterTrigCut++;
1465 }
1466
1467 } // cells
1468
1469
1470 if( gCell > 0 )
1471 clusterTime = clusterTime / (Double_t)gCell;
1472 // fix the reconstriction code time 100ns jumps
1473 if( fESD->GetBunchCrossNumber() % 4 < 2 )
1474 clusterTime -= 0.0000001;
1475
1476 //fClsETime->Fill(clsE,clusterTime);
1477 //fClsEBftTrigCut->Fill(clsE);
1478
1479 if(nClusterTrig>0){
1480 //fClsETime1->Fill(clsE,clusterTime);
1481 }
1482
1483 if(nClusterTrig>0){
1484 cluster->SetChi2(1);
1485 //fClsEAftTrigCut1->Fill(clsE);
1486 }
1487
1488 if(nClusterTrigCut>0){
1489 cluster->SetChi2(2);
1490 //fClsEAftTrigCut2->Fill(clsE);
1491 }
1492
1493 if(nClusterTrigCut>0 && ( nCell > (1 + clsE / 3)))
1494 {
1495 cluster->SetChi2(3);
1496 //fClsEAftTrigCut3->Fill(clsE);
1497 }
1498
1499 if(nClusterTrigCut>0 && (nCell > (1 + clsE / 3) )&&( clusterTime > fTimeCutLow && clusterTime < fTimeCutHigh ))
1500 {
1501 // cluster->SetChi2(4);
1502 //fClsEAftTrigCut4->Fill(clsE);
1503 }
1504 if(nClusterTrigCut<1)
1505 {
1506 cluster->SetChi2(0);
1507
1508 //fClsEAftTrigCut->Fill(clsE);
1509 }
1510
1511 } // clusters
1512}
1513
1514
1515
1516
1517
1518
1519