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