]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGHF/hfe/AliAnalysisTaskHFECal.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Update HFE code
[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"
8efacc22 23#include "TF1.h"
bfc7c23b 24#include "TMath.h"
25#include "TCanvas.h"
26#include "THnSparse.h"
27#include "TLorentzVector.h"
28#include "TString.h"
29#include "TFile.h"
09022877 30#include "TGraphErrors.h"
bfc7c23b 31
6451f693 32#include "TDatabasePDG.h"
33
bfc7c23b 34#include "AliAnalysisTask.h"
35#include "AliAnalysisManager.h"
36
37#include "AliESDEvent.h"
38#include "AliESDHandler.h"
39#include "AliAODEvent.h"
40#include "AliAODHandler.h"
41
e7d87aef 42#include "AliMCEventHandler.h"
43#include "AliMCEvent.h"
44#include "AliMCParticle.h"
45
bfc7c23b 46#include "AliAnalysisTaskHFECal.h"
47#include "TGeoGlobalMagField.h"
48#include "AliLog.h"
49#include "AliAnalysisTaskSE.h"
50#include "TRefArray.h"
51#include "TVector.h"
52#include "AliESDInputHandler.h"
53#include "AliESDpid.h"
54#include "AliESDtrackCuts.h"
55#include "AliPhysicsSelection.h"
56#include "AliESDCaloCluster.h"
57#include "AliAODCaloCluster.h"
58#include "AliEMCALRecoUtils.h"
59#include "AliEMCALGeometry.h"
60#include "AliGeomManager.h"
61#include "stdio.h"
62#include "TGeoManager.h"
63#include "iostream"
64#include "fstream"
65
66#include "AliEMCALTrack.h"
67#include "AliMagF.h"
68
69#include "AliKFParticle.h"
70#include "AliKFVertex.h"
71
72#include "AliPID.h"
73#include "AliPIDResponse.h"
74#include "AliHFEcontainer.h"
75#include "AliHFEcuts.h"
76#include "AliHFEpid.h"
77#include "AliHFEpidBase.h"
78#include "AliHFEpidQAmanager.h"
79#include "AliHFEtools.h"
80#include "AliCFContainer.h"
81#include "AliCFManager.h"
82
e7d87aef 83#include "AliStack.h"
84
bfc7c23b 85#include "AliCentrality.h"
86
c2aad3ae 87using namespace std;
88
bfc7c23b 89ClassImp(AliAnalysisTaskHFECal)
90//________________________________________________________________________
91AliAnalysisTaskHFECal::AliAnalysisTaskHFECal(const char *name)
92 : AliAnalysisTaskSE(name)
93 ,fESD(0)
e7d87aef 94 ,fMC(0)
46305ed6 95 ,stack(0)
fb87d707 96 ,fGeom(0)
bfc7c23b 97 ,fOutputList(0)
c6bca3f2 98 ,fqahist(0)
bfc7c23b 99 ,fTrackCuts(0)
100 ,fCuts(0)
101 ,fIdentifiedAsOutInz(kFALSE)
102 ,fPassTheEventCut(kFALSE)
103 ,fRejectKinkMother(kFALSE)
e7d87aef 104 ,fmcData(kFALSE)
bfc7c23b 105 ,fVz(0.0)
106 ,fCFM(0)
107 ,fPID(0)
108 ,fPIDqa(0)
109 ,fOpeningAngleCut(0.1)
6451f693 110 ,fMimpTassCut(0.5)
6b2cee73 111 ,fMimNsigassCut(-3)
a009053a 112 ,fInvmassCut(0) // no mass
113 ,fSetMassConstraint(kTRUE)
3c475a9d 114 ,fSetMassWidthCut(kFALSE)
2f8ef315 115 ,fSetMassNonlinear(kFALSE)
6451f693 116 ,fSetKFpart(kTRUE)
bfc7c23b 117 ,fNoEvents(0)
118 ,fEMCAccE(0)
f4e0d2d5 119 ,hEMCAccE(0)
bfc7c23b 120 ,fTrkpt(0)
121 ,fTrkEovPBef(0)
122 ,fTrkEovPAft(0)
123 ,fdEdxBef(0)
124 ,fdEdxAft(0)
125 ,fIncpT(0)
fb87d707 126 ,fIncpTM20(0)
bfc7c23b 127 ,fInvmassLS(0)
128 ,fInvmassULS(0)
8806ce6c 129 ,fInvmassLSmc(0)
130 ,fInvmassULSmc(0)
6ec70f8c 131 ,fInvmassLSreco(0)
132 ,fInvmassULSreco(0)
93c189c5 133 ,fInvmassLSmc0(0)
134 ,fInvmassLSmc1(0)
135 ,fInvmassLSmc2(0)
136 ,fInvmassLSmc3(0)
137 ,fInvmassULSmc0(0)
138 ,fInvmassULSmc1(0)
139 ,fInvmassULSmc2(0)
140 ,fInvmassULSmc3(0)
bfc7c23b 141 ,fOpeningAngleLS(0)
142 ,fOpeningAngleULS(0)
143 ,fPhotoElecPt(0)
144 ,fPhoElecPt(0)
fb87d707 145 ,fPhoElecPtM20(0)
0f566acb 146 ,fPhoElecPtM20Mass(0)
f09766dd 147 ,fSameElecPt(0)
fb87d707 148 ,fSameElecPtM20(0)
0f566acb 149 ,fSameElecPtM20Mass(0)
ebef1da9 150 ,fSemiElecPtM20(0)
bfc7c23b 151 ,fTrackPtBefTrkCuts(0)
152 ,fTrackPtAftTrkCuts(0)
153 ,fTPCnsigma(0)
154 ,fCent(0)
155 ,fEleInfo(0)
e1f0fb74 156 ,fElenSigma(0)
feffe705 157 /*
fb87d707 158 ,fClsEBftTrigCut(0)
159 ,fClsEAftTrigCut(0)
160 ,fClsEAftTrigCut1(0)
161 ,fClsEAftTrigCut2(0)
162 ,fClsEAftTrigCut3(0)
163 ,fClsEAftTrigCut4(0)
164 ,fClsETime(0)
165 ,fClsETime1(0)
166 ,fTrigTimes(0)
42c75692 167 ,fCellCheck(0)
feffe705 168 */
e7d87aef 169 ,fInputHFEMC(0)
feffe705 170 ,fInputAlle(0)
e7d87aef 171 ,fIncpTMChfe(0)
3db00c72 172 ,fIncpTMChfeAll(0)
e7d87aef 173 ,fIncpTMCM20hfe(0)
3db00c72 174 ,fIncpTMCM20hfeAll(0)
60544aea 175 ,fIncpTMCM20hfeCheck(0)
bd6ee6dd 176 ,fInputHFEMC_weight(0)
177 ,fIncpTMCM20hfeCheck_weight(0)
e7d87aef 178 ,fIncpTMCpho(0)
179 ,fIncpTMCM20pho(0)
180 ,fPhoElecPtMC(0)
181 ,fPhoElecPtMCM20(0)
0f566acb 182 ,fPhoElecPtMCM20Mass(0)
e7d87aef 183 ,fSameElecPtMC(0)
184 ,fSameElecPtMCM20(0)
0f566acb 185 ,fSameElecPtMCM20Mass(0)
e4b0faf2 186 ,fIncpTMCM20pho_pi0e(0)
187 ,fPhoElecPtMCM20_pi0e(0)
188 ,fSameElecPtMCM20_pi0e(0)
93c189c5 189 ,fIncpTMCM20pho_eta(0)
190 ,fPhoElecPtMCM20_eta(0)
191 ,fSameElecPtMCM20_eta(0)
697ecf6b 192 ,fIncpTMCpho_pi0e_TPC(0)
193 ,fPhoElecPtMC_pi0e_TPC(0)
194 ,fSameElecPtMC_pi0e_TPC(0)
2f8ef315 195 ,fIncpTMCpho_eta_TPC(0)
196 ,fPhoElecPtMC_eta_TPC(0)
197 ,fSameElecPtMC_eta_TPC(0)
a6df418c 198 ,CheckNclust(0)
199 ,CheckNits(0)
ae88597c 200 ,CheckDCA(0)
d113d7cd 201 ,Hpi0pTcheck(0)
19325033 202 ,HETApTcheck(0)
bd6ee6dd 203 ,HphopTcheck(0)
5e0e45b3 204 ,HDpTcheck(0)
205 ,HBpTcheck(0)
93c189c5 206 ,fpTCheck(0)
50919258 207 ,fMomDtoE(0)
70448e3b 208 ,fLabelCheck(0)
209 ,fgeoFake(0)
d19af75d 210 ,fFakeTrk0(0)
211 ,fFakeTrk1(0)
5e38d973 212 ,ftimingEle(0)
6f4a9952 213 ,fIncMaxE(0)
5e38d973 214 ,fIncReco(0)
215 ,fPhoReco(0)
216 ,fSamReco(0)
c82f1b9d 217 ,fIncRecoMaxE(0)
218 ,fPhoRecoMaxE(0)
219 ,fSamRecoMaxE(0)
466d4f7a 220 ,fPhoVertexReco_TPC(0)
221 ,fPhoVertexReco_TPC_Invmass(0)
05ac9660 222 ,fPhoVertexReco_EMCal(0)
223 ,fPhoVertexReco_Invmass(0)
224 ,fPhoVertexReco_step0(0)
225 ,fPhoVertexReco_step1(0)
f2fa5a00 226 ,fMatchV0_0(0)
227 ,fMatchV0_1(0)
228 ,fMatchMC_0(0)
229 ,fMatchMC_1(0)
36e9774f 230 ,fpair(0)
c5a4ddbf 231 ,fFakeRejection0(0)
232 ,fFakeRejection1(0)
233 ,fFakeRejection2(0)
09022877 234 //,fnSigEtaCorr(NULL)
bfc7c23b 235{
236 //Named constructor
237
238 fPID = new AliHFEpid("hfePid");
239 fTrackCuts = new AliESDtrackCuts();
240
09022877 241 for(int i=0; i<7; i++)fnSigEtaCorr[i] = 0;
242
bfc7c23b 243 // Define input and output slots here
244 // Input slot #0 works with a TChain
245 DefineInput(0, TChain::Class());
246 // Output slot #0 id reserved by the base class for AOD
247 // Output slot #1 writes into a TH1 container
248 // DefineOutput(1, TH1I::Class());
249 DefineOutput(1, TList::Class());
250 // DefineOutput(3, TTree::Class());
251}
252
253//________________________________________________________________________
254AliAnalysisTaskHFECal::AliAnalysisTaskHFECal()
255 : AliAnalysisTaskSE("DefaultAnalysis_AliAnalysisTaskHFECal")
256 ,fESD(0)
e7d87aef 257 ,fMC(0)
46305ed6 258 ,stack(0)
fb87d707 259 ,fGeom(0)
bfc7c23b 260 ,fOutputList(0)
c6bca3f2 261 ,fqahist(0)
bfc7c23b 262 ,fTrackCuts(0)
263 ,fCuts(0)
264 ,fIdentifiedAsOutInz(kFALSE)
265 ,fPassTheEventCut(kFALSE)
266 ,fRejectKinkMother(kFALSE)
e7d87aef 267 ,fmcData(kFALSE)
bfc7c23b 268 ,fVz(0.0)
269 ,fCFM(0)
270 ,fPID(0)
271 ,fPIDqa(0)
272 ,fOpeningAngleCut(0.1)
6451f693 273 ,fMimpTassCut(0.5)
6b2cee73 274 ,fMimNsigassCut(-3)
a009053a 275 ,fInvmassCut(0) // no mass
276 ,fSetMassConstraint(kTRUE)
3c475a9d 277 ,fSetMassWidthCut(kFALSE)
2f8ef315 278 ,fSetMassNonlinear(kFALSE)
6451f693 279 ,fSetKFpart(kTRUE)
bfc7c23b 280 ,fNoEvents(0)
281 ,fEMCAccE(0)
f4e0d2d5 282 ,hEMCAccE(0)
bfc7c23b 283 ,fTrkpt(0)
284 ,fTrkEovPBef(0)
285 ,fTrkEovPAft(0)
286 ,fdEdxBef(0)
287 ,fdEdxAft(0)
288 ,fIncpT(0)
fb87d707 289 ,fIncpTM20(0)
bfc7c23b 290 ,fInvmassLS(0)
291 ,fInvmassULS(0)
8806ce6c 292 ,fInvmassLSmc(0)
293 ,fInvmassULSmc(0)
6ec70f8c 294 ,fInvmassLSreco(0)
295 ,fInvmassULSreco(0)
93c189c5 296 ,fInvmassLSmc0(0)
297 ,fInvmassLSmc1(0)
298 ,fInvmassLSmc2(0)
299 ,fInvmassLSmc3(0)
300 ,fInvmassULSmc0(0)
301 ,fInvmassULSmc1(0)
302 ,fInvmassULSmc2(0)
303 ,fInvmassULSmc3(0)
bfc7c23b 304 ,fOpeningAngleLS(0)
305 ,fOpeningAngleULS(0)
306 ,fPhotoElecPt(0)
307 ,fPhoElecPt(0)
fb87d707 308 ,fPhoElecPtM20(0)
0f566acb 309 ,fPhoElecPtM20Mass(0)
f09766dd 310 ,fSameElecPt(0)
fb87d707 311 ,fSameElecPtM20(0)
0f566acb 312 ,fSameElecPtM20Mass(0)
ebef1da9 313 ,fSemiElecPtM20(0)
bfc7c23b 314 ,fTrackPtBefTrkCuts(0)
315 ,fTrackPtAftTrkCuts(0)
316 ,fTPCnsigma(0)
317 ,fCent(0)
318 ,fEleInfo(0)
e1f0fb74 319 ,fElenSigma(0)
feffe705 320 /*
fb87d707 321 ,fClsEBftTrigCut(0)
322 ,fClsEAftTrigCut(0)
323 ,fClsEAftTrigCut1(0)
324 ,fClsEAftTrigCut2(0)
325 ,fClsEAftTrigCut3(0)
326 ,fClsEAftTrigCut4(0)
327 ,fClsETime(0)
328 ,fClsETime1(0)
329 ,fTrigTimes(0)
42c75692 330 ,fCellCheck(0)
feffe705 331 */
e7d87aef 332 ,fInputHFEMC(0)
feffe705 333 ,fInputAlle(0)
e7d87aef 334 ,fIncpTMChfe(0)
3db00c72 335 ,fIncpTMChfeAll(0)
e7d87aef 336 ,fIncpTMCM20hfe(0)
3db00c72 337 ,fIncpTMCM20hfeAll(0)
60544aea 338 ,fIncpTMCM20hfeCheck(0)
bd6ee6dd 339 ,fInputHFEMC_weight(0)
340 ,fIncpTMCM20hfeCheck_weight(0)
e7d87aef 341 ,fIncpTMCpho(0)
342 ,fIncpTMCM20pho(0)
343 ,fPhoElecPtMC(0)
344 ,fPhoElecPtMCM20(0)
0f566acb 345 ,fPhoElecPtMCM20Mass(0)
e7d87aef 346 ,fSameElecPtMC(0)
347 ,fSameElecPtMCM20(0)
0f566acb 348 ,fSameElecPtMCM20Mass(0)
e4b0faf2 349 ,fIncpTMCM20pho_pi0e(0)
350 ,fPhoElecPtMCM20_pi0e(0)
351 ,fSameElecPtMCM20_pi0e(0)
93c189c5 352 ,fIncpTMCM20pho_eta(0)
353 ,fPhoElecPtMCM20_eta(0)
354 ,fSameElecPtMCM20_eta(0)
697ecf6b 355 ,fIncpTMCpho_pi0e_TPC(0)
356 ,fPhoElecPtMC_pi0e_TPC(0)
357 ,fSameElecPtMC_pi0e_TPC(0)
2f8ef315 358 ,fIncpTMCpho_eta_TPC(0)
359 ,fPhoElecPtMC_eta_TPC(0)
360 ,fSameElecPtMC_eta_TPC(0)
a6df418c 361 ,CheckNclust(0)
362 ,CheckNits(0)
ae88597c 363 ,CheckDCA(0)
d113d7cd 364 ,Hpi0pTcheck(0)
19325033 365 ,HETApTcheck(0)
e4b0faf2 366 ,HphopTcheck(0)
5e0e45b3 367 ,HDpTcheck(0)
368 ,HBpTcheck(0)
93c189c5 369 ,fpTCheck(0)
50919258 370 ,fMomDtoE(0)
70448e3b 371 ,fLabelCheck(0)
372 ,fgeoFake(0)
d19af75d 373 ,fFakeTrk0(0)
374 ,fFakeTrk1(0)
59d998de 375 ,ftimingEle(0)
6f4a9952 376 ,fIncMaxE(0)
5e38d973 377 ,fIncReco(0)
378 ,fPhoReco(0)
379 ,fSamReco(0)
c82f1b9d 380 ,fIncRecoMaxE(0)
381 ,fPhoRecoMaxE(0)
382 ,fSamRecoMaxE(0)
466d4f7a 383 ,fPhoVertexReco_TPC(0)
384 ,fPhoVertexReco_TPC_Invmass(0)
05ac9660 385 ,fPhoVertexReco_EMCal(0)
386 ,fPhoVertexReco_Invmass(0)
387 ,fPhoVertexReco_step0(0)
388 ,fPhoVertexReco_step1(0)
f2fa5a00 389 ,fMatchV0_0(0)
390 ,fMatchV0_1(0)
391 ,fMatchMC_0(0)
392 ,fMatchMC_1(0)
36e9774f 393 ,fpair(0)
c5a4ddbf 394 ,fFakeRejection0(0)
395 ,fFakeRejection1(0)
396 ,fFakeRejection2(0)
09022877 397 //,fnSigEtaCorr(NULL)
bfc7c23b 398{
399 //Default constructor
400 fPID = new AliHFEpid("hfePid");
401
402 fTrackCuts = new AliESDtrackCuts();
403
09022877 404 for(int i=0; i<7; i++)fnSigEtaCorr[i] = 0;
405
bfc7c23b 406 // Constructor
407 // Define input and output slots here
408 // Input slot #0 works with a TChain
409 DefineInput(0, TChain::Class());
410 // Output slot #0 id reserved by the base class for AOD
411 // Output slot #1 writes into a TH1 container
412 // DefineOutput(1, TH1I::Class());
413 DefineOutput(1, TList::Class());
414 //DefineOutput(3, TTree::Class());
415}
416//_________________________________________
417
418AliAnalysisTaskHFECal::~AliAnalysisTaskHFECal()
419{
420 //Destructor
421
422 delete fOutputList;
fb87d707 423 delete fGeom;
bfc7c23b 424 delete fPID;
a6df418c 425 delete fCuts;
bfc7c23b 426 delete fCFM;
427 delete fPIDqa;
428 delete fTrackCuts;
429}
430//_________________________________________
431
432void AliAnalysisTaskHFECal::UserExec(Option_t*)
433{
434 //Main loop
435 //Called for each event
436
437 // create pointer to event
438 fESD = dynamic_cast<AliESDEvent*>(InputEvent());
439 if (!fESD) {
440 printf("ERROR: fESD not available\n");
441 return;
442 }
443
444 if(!fCuts){
445 AliError("HFE cuts not available");
446 return;
447 }
448
449 if(!fPID->IsInitialized()){
450 // Initialize PID with the given run number
451 AliWarning("PID not initialised, get from Run no");
452 fPID->InitializePID(fESD->GetRunNumber());
453 }
454
e7d87aef 455 if(fmcData)fMC = MCEvent();
46305ed6 456 //AliStack* stack = NULL;
e7d87aef 457 if(fmcData && fMC)stack = fMC->Stack();
458
459 Float_t cent = -1.;
460 AliCentrality *centrality = fESD->GetCentrality();
461 cent = centrality->GetCentralityPercentile("V0M");
462
463 //---- fill MC track info
464 if(fmcData && fMC)
465 {
466 Int_t nParticles = stack->GetNtrack();
467 for (Int_t iParticle = 0; iParticle < nParticles; iParticle++) {
468 TParticle* particle = stack->Particle(iParticle);
469 int fPDG = particle->GetPdgCode();
470 double mcZvertex = fMC->GetPrimaryVertex()->GetZ();
471 double pTMC = particle->Pt();
feffe705 472 double proR = particle->R();
9b3495ff 473 double etaMC = particle->Eta();
474 if(fabs(etaMC)>0.6)continue;
e7d87aef 475 Bool_t mcInDtoE= kFALSE;
476 Bool_t mcInBtoE= kFALSE;
477
d113d7cd 478 Bool_t MChijing = fMC->IsFromBGEvent(iParticle);
fd426741 479 //if(!MChijing)printf("not MC hijing");
d113d7cd 480 int iHijing = 1;
481 if(!MChijing)iHijing = 0;
68d718e7 482 double mcphoinfo[5];
483 mcphoinfo[0] = cent;
484 mcphoinfo[1] = pTMC;
485 mcphoinfo[2] = iHijing;
68d718e7 486 if(fPDG==111)Hpi0pTcheck->Fill(mcphoinfo);
487 if(fPDG==221)HETApTcheck->Fill(mcphoinfo);
bd6ee6dd 488 if(fabs(fPDG)==411 || fabs(fPDG)==413 || fabs(fPDG)==421 || fabs(fPDG)==423 || fabs(fPDG)==431)HDpTcheck->Fill(pTMC,iHijing);
489 if(fabs(fPDG)==511 || fabs(fPDG)==513 || fabs(fPDG)==521 || fabs(fPDG)==523 || fabs(fPDG)==531)HBpTcheck->Fill(pTMC,iHijing);
d113d7cd 490
e4b0faf2 491 if(particle->GetFirstMother()>-1 && fPDG==22)
492 {
493 int parentPID = stack->Particle(particle->GetFirstMother())->GetPdgCode();
494 if(parentPID==111 || parentPID==221)HphopTcheck->Fill(pTMC,iHijing); // pi0->g & eta->g
495 }
496
e7d87aef 497 if(particle->GetFirstMother()>-1 && fabs(fPDG)==11)
498 {
499 int parentPID = stack->Particle(particle->GetFirstMother())->GetPdgCode();
bd6ee6dd 500 double pTMCparent = stack->Particle(particle->GetFirstMother())->Pt();
e7d87aef 501 if((fabs(parentPID)==411 || fabs(parentPID)==413 || fabs(parentPID)==421 || fabs(parentPID)==423 || fabs(parentPID)==431)&& fabs(fPDG)==11)mcInDtoE = kTRUE;
502 if((fabs(parentPID)==511 || fabs(parentPID)==513 || fabs(parentPID)==521 || fabs(parentPID)==523 || fabs(parentPID)==531)&& fabs(fPDG)==11)mcInBtoE = kTRUE;
bd6ee6dd 503 if((mcInBtoE || mcInDtoE) && fabs(mcZvertex)<10.0)
504 {
505 fInputHFEMC->Fill(cent,pTMC);
506 double mcinfo[5];
507 mcinfo[0] = cent;
508 mcinfo[1] = pTMC;
509 mcinfo[2] = 0.0;
510 mcinfo[3] = iHijing;
511 mcinfo[4] = pTMCparent;
512 fInputHFEMC_weight->Fill(mcinfo);
513 }
e7d87aef 514 }
515
5e0e45b3 516
feffe705 517 if(proR<7.0 && fabs(fPDG)==11)fInputAlle->Fill(cent,pTMC);
518
e7d87aef 519 }
520 }
521
bfc7c23b 522 fNoEvents->Fill(0);
523
524 Int_t fNOtrks = fESD->GetNumberOfTracks();
525 const AliESDVertex *pVtx = fESD->GetPrimaryVertex();
526
527 Double_t pVtxZ = -999;
528 pVtxZ = pVtx->GetZ();
529
530 if(TMath::Abs(pVtxZ)>10) return;
531 fNoEvents->Fill(1);
532
f4e0d2d5 533 if(fNOtrks<1) return;
bfc7c23b 534 fNoEvents->Fill(2);
535
536 AliPIDResponse *pidResponse = fInputHandler->GetPIDResponse();
537 if(!pidResponse){
538 AliDebug(1, "Using default PID Response");
539 pidResponse = AliHFEtools::GetDefaultPID(kFALSE, fInputEvent->IsA() == AliAODEvent::Class());
540 }
541
542 fPID->SetPIDResponse(pidResponse);
543
544 fCFM->SetRecEventInfo(fESD);
545
bfc7c23b 546 fCent->Fill(cent);
547
b12bc641 548 //if(cent>90.) return;
bfc7c23b 549
550 // Calorimeter info.
551
42c75692 552 FindTriggerClusters();
553
bfc7c23b 554 // make EMCAL array
6f4a9952 555 double maxE = 0.0;
bfc7c23b 556 for(Int_t iCluster=0; iCluster<fESD->GetNumberOfCaloClusters(); iCluster++)
557 {
558 AliESDCaloCluster *clust = fESD->GetCaloCluster(iCluster);
68d718e7 559 if(clust && clust->IsEMCAL())
bfc7c23b 560 {
561 double clustE = clust->E();
562 float emcx[3]; // cluster pos
563 clust->GetPosition(emcx);
564 TVector3 clustpos(emcx[0],emcx[1],emcx[2]);
565 double emcphi = clustpos.Phi();
566 double emceta = clustpos.Eta();
fb87d707 567 double calInfo[5];
568 calInfo[0] = emcphi; calInfo[1] = emceta; calInfo[2] = clustE; calInfo[3] = cent; calInfo[4] = clust->Chi2();
f4e0d2d5 569 //fEMCAccE->Fill(calInfo);
f4e0d2d5 570 hEMCAccE->Fill(cent,clustE);
6f4a9952 571 if(clustE>maxE)maxE = clustE;
bfc7c23b 572 }
573 }
574
575 // Track loop
576 for (Int_t iTracks = 0; iTracks < fESD->GetNumberOfTracks(); iTracks++) {
577 AliESDtrack* track = fESD->GetTrack(iTracks);
578 if (!track) {
579 printf("ERROR: Could not receive track %d\n", iTracks);
580 continue;
581 }
e7d87aef 582
89f41a30 583 //--- Get MC informtion
584
46305ed6 585 int parentlabel = 99999;
586 int parentPID = 99999;
587 int grand_parentlabel = 99999;
588 int grand_parentPID = 99999;
e7d87aef 589 Bool_t mcPho = kFALSE;
590 Bool_t mcDtoE= kFALSE;
591 Bool_t mcBtoE= kFALSE;
93c189c5 592 Bool_t mcOrgPi0 = kFALSE;
593 Bool_t mcOrgEta = kFALSE;
e7d87aef 594 double mcele = -1.;
44be9e1d 595 double mcpT = 0.0;
b12bc641 596 double mcMompT = 0.0;
e4b0faf2 597 //double mcGrandMompT = 0.0;
598 double mcWeight = -10.0;
bd50bf24 599 double conv_proR = -1.0;
d113d7cd 600
601 int iHijing = 1;
70448e3b 602 int mcLabel = -1;
d113d7cd 603
e7d87aef 604 if(fmcData && fMC && stack)
605 {
feffe705 606 Int_t label = TMath::Abs(track->GetLabel());
7d495dc5 607 //mcLabel = track->GetLabel();
1dd85087 608 mcLabel = abs(track->GetLabel()); // check for conv. issue
bd6ee6dd 609
8efacc22 610 if(mcLabel>-1)
bd6ee6dd 611 {
612
613 Bool_t MChijing = fMC->IsFromBGEvent(label);
614 if(!MChijing)iHijing = 0;
615
616 TParticle* particle = stack->Particle(label);
617 int mcpid = particle->GetPdgCode();
618 mcpT = particle->Pt();
bd50bf24 619 conv_proR = particle->R();
bd6ee6dd 620 //printf("MCpid = %d",mcpid);
621 if(particle->GetFirstMother()>-1)
622 {
623 //int parentlabel = particle->GetFirstMother();
624 //int parentPID = stack->Particle(particle->GetFirstMother())->GetPdgCode();
625 //mcMompT = stack->Particle(particle->GetFirstMother())->Pt();
626 FindMother(particle, parentlabel, parentPID);
627 mcMompT = stack->Particle(parentlabel)->Pt();
628 if((parentPID==22 || parentPID==111 || parentPID==221)&& fabs(mcpid)==11)mcPho = kTRUE;
629 if((fabs(parentPID)==411 || fabs(parentPID)==413 || fabs(parentPID)==421 || fabs(parentPID)==423 || fabs(parentPID)==431)&& fabs(mcpid)==11)mcDtoE = kTRUE;
630 if((fabs(parentPID)==511 || fabs(parentPID)==513 || fabs(parentPID)==521 || fabs(parentPID)==523 || fabs(parentPID)==531)&& fabs(mcpid)==11)mcBtoE = kTRUE;
631
632 // make D->e pT correlation
633 if(mcDtoE)fMomDtoE->Fill(mcpT,mcMompT);
634
635 //cout << "check PID = " << parentPID << endl;
636 //cout << "check pho = " << mcPho << endl;
637 //cout << "check D or B = " << mcDtoE << endl;
638 // pi->e (Dalitz)
639 if(parentPID==111 && fabs(mcpid)==11 && mcMompT>0.0)
640 {
641 //cout << "find pi0->e " << endl;
642 mcOrgPi0 = kTRUE;
643 mcWeight = GetMCweight(mcMompT);
644 }
645 // eta->e (Dalitz)
646 if(parentPID==221 && fabs(mcpid)==11 && mcMompT>0.0)
647 {
648 //cout << "find Eta->e " << endl;
649 mcOrgEta = kTRUE;
650 mcWeight = GetMCweightEta(mcMompT);
651 }
652
653 // access grand parent
654 TParticle* particle_parent = stack->Particle(parentlabel); // get parent pointer
655 //if(particle_parent->GetFirstMother()>-1 && parentPID==22 && fabs(mcpid)==11) // get grand parent g->e
656 if(particle_parent->GetFirstMother()>-1 && (parentPID==22 || parentPID==111) && fabs(mcpid)==11) // get grand parent g->e
657 {
658 //int grand_parentPID = stack->Particle(particle_parent->GetFirstMother())->GetPdgCode();
659 //double pTtmp = stack->Particle(particle_parent->GetFirstMother())->Pt();
660 FindMother(particle_parent, grand_parentlabel, grand_parentPID);
661 double mcGrandpT = stack->Particle(grand_parentlabel)->Pt();
662 if(grand_parentPID==111 && mcGrandpT>0.0)
663 {
664 // check eta->pi0 decay !
665 int grand2_parentlabel = 99999; int grand2_parentPID = 99999;
666 TParticle* particle_grand = stack->Particle(grand_parentlabel); // get parent pointer
667 FindMother(particle_grand, grand2_parentlabel, grand2_parentPID);
668 if(grand2_parentPID==221)
669 {
670 //cout << "find Eta->e " << endl;
671 double mcGrandpT2 = stack->Particle(grand2_parentlabel)->Pt();
672 mcOrgEta = kTRUE;
673 mcWeight = GetMCweight(mcGrandpT2);
674 mcMompT = mcGrandpT2;
675 }
676 else
677 {
678 //cout << "find pi0->e " << endl;
679 mcOrgPi0 = kTRUE;
680 mcWeight = GetMCweight(mcGrandpT);
681 mcMompT = mcGrandpT;
682 }
683 }
684
685 if(grand_parentPID==221 && mcGrandpT>0.0)
686 {
687 //cout << "find Eta->e " << endl;
688 mcOrgEta = kTRUE;
689 mcOrgPi0 = kFALSE;
690 mcWeight = GetMCweightEta(mcGrandpT);
691 mcMompT = mcGrandpT;
692 }
693 }
694 }
695
cc1605fe 696 //cout << "===================="<<endl;
697 //cout << "mcDtoE : " << mcDtoE << endl;
698 //cout << "mcBtoE : " << mcBtoE << endl;
699 //cout << "mcPho : " << mcPho << endl;
700
701 if(fabs(mcpid)==11)mcele= 0.;
702 //cout << "check e: " << mcele << endl;
bd6ee6dd 703 if(fabs(mcpid)==11 && mcDtoE)mcele= 1.;
cc1605fe 704 //cout << "check D->e: " << mcele << endl;
bd6ee6dd 705 if(fabs(mcpid)==11 && mcBtoE)mcele= 2.;
cc1605fe 706 //cout << "check B->e: " << mcele << endl;
bd6ee6dd 707 if(fabs(mcpid)==11 && mcPho)mcele= 3.;
cc1605fe 708 //cout << "check Pho->e: " << mcele << endl;
bd6ee6dd 709
09022877 710 //cout << "check PID " << endl;
fd2126aa 711 if(fabs(mcpid)!=11)
712 {
09022877 713 //cout << "!= 11" << endl;
714 //cout << mcpid << endl;
fd2126aa 715 }
716 if(mcele==-1)
717 {
09022877 718 //cout << "mcele==-1" << endl;
719 //cout << mcele << endl;
720 //cout << mcpid << endl;
fd2126aa 721 }
722
bd6ee6dd 723 } // end of mcLabel>-1
724
89f41a30 725 } // <------ end of MC info.
e7d87aef 726
2b4460a6 727 //cout << "Pi0 = " << mcOrgPi0 << " ; Eta = " << mcOrgEta << endl;
50919258 728 //printf("weight = %f\n",mcWeight);
729
9b3495ff 730 if(TMath::Abs(track->Eta())>0.6) continue;
fbe402d7 731 //if(TMath::Abs(track->Pt()<2.5)) continue;
732 if(TMath::Abs(track->Pt()<0.1)) continue;
bfc7c23b 733
05ac9660 734 int nITS = track->GetNcls(0);
735
bfc7c23b 736 fTrackPtBefTrkCuts->Fill(track->Pt());
f2111233 737
05ac9660 738 /*
f2fa5a00 739 UChar_t itsPixel = track->GetITSClusterMap();
05ac9660 740 cout << "nITS = " << nITS << endl;
741 if(itsPixel & BIT(0))cout << "1st layer hit" << endl;
742 if(itsPixel & BIT(1))cout << "2nd layer hit" << endl;
743 if(itsPixel & BIT(2))cout << "3rd layer hit" << endl;
744 if(itsPixel & BIT(3))cout << "4th layer hit" << endl;
745 if(itsPixel & BIT(4))cout << "5th layer hit" << endl;
746 */
747
bfc7c23b 748 // RecKine: ITSTPC cuts
749 if(!ProcessCutStep(AliHFEcuts::kStepRecKineITSTPC, track)) continue;
750
751 //RecKink
752 if(fRejectKinkMother) { // Quick and dirty fix to reject both kink mothers and daughters
753 if(track->GetKinkIndex(0) != 0) continue;
754 }
755
756 // RecPrim
757 if(!ProcessCutStep(AliHFEcuts::kStepRecPrim, track)) continue;
758
759 // HFEcuts: ITS layers cuts
6e1d7a8a 760 if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsITS, track)) continue;
bfc7c23b 761
f09766dd 762 // HFE cuts: TPC PID cleanup
763 if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTPC, track)) continue;
764
36e9774f 765 if(mcPho && iHijing==0)fPhoVertexReco_step0->Fill(track->Pt(),conv_proR); // check MC vertex
766 if(mcPho && iHijing==1)fPhoVertexReco_step1->Fill(track->Pt(),conv_proR); // check MC vertex
767
a6df418c 768 int nTPCcl = track->GetTPCNcls();
96167a04 769 //int nTPCclF = track->GetTPCNclsF(); // warnings
05ac9660 770 //int nITS = track->GetNcls(0);
771
05ac9660 772
bfc7c23b 773 fTrackPtAftTrkCuts->Fill(track->Pt());
774
775 Double_t mom = -999., eop=-999., pt = -999., dEdx=-999., fTPCnSigma=-10, phi=-999., eta=-999.;
feffe705 776 pt = track->Pt();
d26226a1 777 if(pt<2.5)continue;
fbe402d7 778 if(pt<0.1)continue;
bfc7c23b 779
6f4a9952 780 //Int_t charge = track->Charge();
bfc7c23b 781 fTrkpt->Fill(pt);
782 mom = track->P();
783 phi = track->Phi();
784 eta = track->Eta();
ae88597c 785 float dca_xy;
786 float dca_z;
787 track->GetImpactParameters(dca_xy,dca_z);
788
bfc7c23b 789 dEdx = track->GetTPCsignal();
790 fTPCnSigma = fPID->GetPIDResponse() ? fPID->GetPIDResponse()->NumberOfSigmasTPC(track, AliPID::kElectron) : 1000;
8af58b0d 791
09022877 792 //cout << "nSigma correctoon-----" << endl;
793 //cout << "org = " << fTPCnSigma << endl;
9eae9ea5 794 /*
09022877 795 if(!fmcData) // nsigma eta correction
796 {
797 double nSigexpCorr = NsigmaCorrection(eta,cent);
798 fTPCnSigma -= nSigexpCorr;
799 }
9eae9ea5 800 */
09022877 801 //cout << "correction = " << fTPCnSigma << endl;
802
89f41a30 803 //--- track cluster match
804
bfc7c23b 805 double ncells = -1.0;
806 double m20 = -1.0;
807 double m02 = -1.0;
808 double disp = -1.0;
809 double rmatch = -1.0;
810 double nmatch = -1.0;
89f41a30 811 //double oppstatus = 0.0;
59d998de 812 double emctof = 0.0;
6f4a9952 813 Bool_t MaxEmatch = kFALSE;
bfc7c23b 814
bfc7c23b 815 Int_t clsId = track->GetEMCALcluster();
816 if (clsId>0){
817 AliESDCaloCluster *clust = fESD->GetCaloCluster(clsId);
818 if(clust && clust->IsEMCAL()){
819
820 double clustE = clust->E();
6f4a9952 821 if(clustE==maxE)MaxEmatch = kTRUE;
bfc7c23b 822 eop = clustE/fabs(mom);
8efacc22 823
bfc7c23b 824 //double clustT = clust->GetTOF();
825 ncells = clust->GetNCells();
826 m02 = clust->GetM02();
827 m20 = clust->GetM20();
828 disp = clust->GetDispersion();
829 double delphi = clust->GetTrackDx();
830 double deleta = clust->GetTrackDz();
831 rmatch = sqrt(pow(delphi,2)+pow(deleta,2));
832 nmatch = clust->GetNTracksMatched();
59d998de 833 emctof = clust->GetTOF();
834 //cout << "emctof = " << emctof << endl;
89f41a30 835 cout << "eop org = "<< eop << endl;
836 double eoporg = eop;
837 if(fmcData)
9eae9ea5 838 {
839 double mceopcorr = MCEopMeanCorrection(pt,cent);
840 eop += mceopcorr;
841 }
89f41a30 842 cout << "eop corr = " << eop << endl;
9eae9ea5 843
feffe705 844 double valdedx[16];
a6df418c 845 valdedx[0] = pt; valdedx[1] = nITS; valdedx[2] = phi; valdedx[3] = eta; valdedx[4] = fTPCnSigma;
05ac9660 846 valdedx[5] = eop; valdedx[6] = rmatch; valdedx[7] = dca_xy, valdedx[8] = dca_z; valdedx[9] = m20; valdedx[10] = mcpT;
89f41a30 847 valdedx[11] = cent; valdedx[12] = dEdx; valdedx[13] = eoporg; valdedx[14] = nTPCcl;
feffe705 848 valdedx[15] = mcele;
89f41a30 849 fEleInfo->Fill(valdedx);
bfc7c23b 850
851 }
852 }
e1f0fb74 853
854 //Get Cal info PID response
9eae9ea5 855 /*
e1f0fb74 856 double eop2;
857 double ss[4];
858 Double_t nSigmaEop = fPID->GetPIDResponse()->NumberOfSigmasEMCAL(track,AliPID::kElectron,eop2,ss);
859 if(fTPCnSigma>-1.5 && fTPCnSigma<3.0 && nITS>2.5 && nTPCcl>100)
860 {
861 double valEop[3];
862 valEop[0] = cent;
863 valEop[1] = pt;
864 valEop[2] = nSigmaEop;
865 fElenSigma->Fill(valEop);
866 }
9eae9ea5 867 */
89f41a30 868 // --- tarck cut & e ID
e1f0fb74 869
a6df418c 870 if(nITS<2.5)continue;
871 if(nTPCcl<100)continue;
7d495dc5 872
7d495dc5 873
a6df418c 874 CheckNclust->Fill(nTPCcl);
875 CheckNits->Fill(nITS);
ae88597c 876 CheckDCA->Fill(dca_xy,dca_z);
f2111233 877 // check production vertex of photons
f2111233 878
a6df418c 879
42c75692 880 fdEdxBef->Fill(mom,fTPCnSigma);
bfc7c23b 881 fTPCnsigma->Fill(mom,fTPCnSigma);
f09766dd 882 if(fTPCnSigma >= -1.0 && fTPCnSigma <= 3)fTrkEovPBef->Fill(pt,eop);
bfc7c23b 883
89f41a30 884 //--- track accepted by HFE
bfc7c23b 885
89f41a30 886 Int_t pidpassed = 1;
bfc7c23b 887 AliHFEpidObject hfetrack;
888 hfetrack.SetAnalysisType(AliHFEpidObject::kESDanalysis);
889 hfetrack.SetRecTrack(track);
3db00c72 890 double binct = 10.5;
891 if((0.0< cent) && (cent<5.0)) binct = 0.5;
892 if((5.0< cent) && (cent<10.0)) binct = 1.5;
893 if((10.0< cent) && (cent<20.0)) binct = 2.5;
894 if((20.0< cent) && (cent<30.0)) binct = 3.5;
895 if((30.0< cent) && (cent<40.0)) binct = 4.5;
896 if((40.0< cent) && (cent<50.0)) binct = 5.5;
897 if((50.0< cent) && (cent<60.0)) binct = 6.5;
898 if((60.0< cent) && (cent<70.0)) binct = 7.5;
899 if((70.0< cent) && (cent<80.0)) binct = 8.5;
900 if((80.0< cent) && (cent<90.0)) binct = 9.5;
901 if((90.0< cent) && (cent<100.0)) binct = 10.5;
902
903 hfetrack.SetCentrality((int)binct); //added
bfc7c23b 904 hfetrack.SetPbPb();
905 if(!fPID->IsSelected(&hfetrack, NULL, "", fPIDqa)) pidpassed = 0;
906
f2fa5a00 907 if(pidpassed==0) continue; // nSigma rejection
89f41a30 908
909 //--- photonic ID
910
0f566acb 911 Bool_t fFlagPhotonicTPC = kFALSE;
912 Bool_t fFlagConvinatTPC = kFALSE;
913 SelectPhotonicElectron(iTracks,cent,track,fFlagPhotonicTPC,fFlagConvinatTPC,fTPCnSigma,m20,eop,mcele,mcWeight,iHijing,mcOrgPi0,mcOrgEta,0);
89f41a30 914
915 //--- check reco eff. with TPC
916 double phoval[5];
917 phoval[0] = cent;
918 phoval[1] = pt;
919 phoval[2] = fTPCnSigma;
920 phoval[3] = iHijing;
921 phoval[4] = mcMompT;
922
2f8ef315 923 if((fTPCnSigma >= -5.0 && fTPCnSigma <= 5) && (mcOrgPi0 || mcOrgEta))
89f41a30 924 {
925 if(iHijing==1)mcWeight = 1.0;
2f8ef315 926 if(mcOrgPi0)
927 {
928 fIncpTMCpho_pi0e_TPC->Fill(phoval,mcWeight);
0f566acb 929 if(fFlagPhotonicTPC) fPhoElecPtMC_pi0e_TPC->Fill(phoval,mcWeight);
930 if(fFlagConvinatTPC) fSameElecPtMC_pi0e_TPC->Fill(phoval,mcWeight);
466d4f7a 931 fPhoVertexReco_TPC->Fill(track->Pt(),conv_proR,mcWeight); // check MC vertex
932 if(fFlagPhotonicTPC)fPhoVertexReco_TPC_Invmass->Fill(track->Pt(),conv_proR,mcWeight); // check MC vertex
2f8ef315 933 }
934 if(mcOrgEta)
935 {
936 fIncpTMCpho_eta_TPC->Fill(phoval,mcWeight);
0f566acb 937 if(fFlagPhotonicTPC) fPhoElecPtMC_eta_TPC->Fill(phoval,mcWeight);
938 if(fFlagConvinatTPC) fSameElecPtMC_eta_TPC->Fill(phoval,mcWeight);
2f8ef315 939 }
89f41a30 940 }
941
0f566acb 942 //--- matching check
f2fa5a00 943 double emcphimim = 1.396;
944 double emcphimax = 3.14;
945 if(phi>emcphimim && phi<emcphimax)
946 {
0f566acb 947 if(fFlagPhotonicTPC)fMatchV0_0->Fill(pt); // data
d26226a1 948 if(mcele>2.1)fMatchMC_0->Fill(pt,mcWeight); // MC
949
950 if(eop>=0.0) // have a match
f2fa5a00 951 {
0f566acb 952 if(fFlagPhotonicTPC)fMatchV0_1->Fill(pt);
d26226a1 953 if(mcele>2.1)fMatchMC_1->Fill(pt,mcWeight);
f2fa5a00 954 }
955 }
956
c5a4ddbf 957 // check fake rejection
958 if(mcOrgPi0 || mcOrgEta)
959 {
1c9a06c8 960 double phiacc0 = 80.0/180.0*acos(-1);
961 double phiacc1 = 180.0/180.0*acos(-1);
c5a4ddbf 962 int TrStat = 0;
1c9a06c8 963 if(phi>phiacc0 && phi<phiacc1)
964 {
965 if(track->GetLabel()>0)TrStat = 1;
966 fFakeRejection0->Fill(TrStat,pt,mcWeight);
967 if(eop>-1.0)fFakeRejection1->Fill(TrStat,pt,mcWeight); // have match
968 if(eop>0.9 && eop<1.3)fFakeRejection2->Fill(TrStat,pt,mcWeight); // have PID
969 }
c5a4ddbf 970 }
971
89f41a30 972 //+++++++ E/p cut ++++++++++++++++
973
974 if(eop<0.9 || eop>1.3)continue;
0f566acb 975
976 Bool_t fFlagPhotonicElec = kFALSE;
977 Bool_t fFlagConvinatElec = kFALSE;
978 SelectPhotonicElectron(iTracks,cent,track,fFlagPhotonicElec,fFlagConvinatElec,fTPCnSigma,m20,eop,mcele,mcWeight,iHijing,mcOrgPi0,mcOrgEta,1);
979
bfc7c23b 980 fTrkEovPAft->Fill(pt,eop);
42c75692 981 fdEdxAft->Fill(mom,fTPCnSigma);
bfc7c23b 982
e7d87aef 983 // Fill real data
fb87d707 984 fIncpT->Fill(cent,pt);
bfc7c23b 985 if(fFlagPhotonicElec) fPhoElecPt->Fill(cent,pt);
f09766dd 986 if(fFlagConvinatElec) fSameElecPt->Fill(cent,pt);
fb87d707 987
988 if(m20>0.0 && m20<0.3)
989 {
59d998de 990 fIncpTM20->Fill(cent,pt);
991 ftimingEle->Fill(pt,emctof);
fb87d707 992 if(fFlagPhotonicElec) fPhoElecPtM20->Fill(cent,pt);
993 if(fFlagConvinatElec) fSameElecPtM20->Fill(cent,pt);
ebef1da9 994 if(!fFlagPhotonicElec) fSemiElecPtM20->Fill(cent,pt);
fb87d707 995 }
feffe705 996
5e38d973 997
e7d87aef 998 // MC
70448e3b 999 // check label for electron candidiates
70448e3b 1000 int idlabel = 1;
1001 if(mcLabel==0)idlabel = 0;
1002 fLabelCheck->Fill(pt,idlabel);
1003 if(mcLabel==0)fgeoFake->Fill(phi,eta);
1004
d19af75d 1005 if(mcLabel<0 && m20>0.0 && m20<0.3 && fTPCnSigma>-1 && fTPCnSigma<3)
1006 {
1007 fFakeTrk0->Fill(cent,pt);
1008 }
1009
1010 if(mcele>-1) // select MC electrons
e7d87aef 1011 {
3db00c72 1012
1013 fIncpTMChfeAll->Fill(cent,pt);
1014 if(m20>0.0 && m20<0.3)fIncpTMCM20hfeAll->Fill(cent,pt);
d19af75d 1015 if(m20>0.0 && m20<0.3 && fTPCnSigma>-1 && fTPCnSigma<3)fFakeTrk1->Fill(cent,pt);
3db00c72 1016
e4b0faf2 1017 if(mcBtoE || mcDtoE) // select B->e & D->e
e7d87aef 1018 {
1019 fIncpTMChfe->Fill(cent,pt);
bd6ee6dd 1020 if(m20>0.0 && m20<0.3)
1021 {
09022877 1022 //cout << "MC label = " << mcLabel << endl;
bd6ee6dd 1023 fIncpTMCM20hfe->Fill(cent,pt);
1024 fIncpTMCM20hfeCheck->Fill(cent,mcpT);
1025 fIncpTMCM20hfeCheck_weight->Fill(phoval);
1026 }
e7d87aef 1027 }
e4b0faf2 1028
1029 if(mcPho) // select photonic electrons
e7d87aef 1030 {
44be9e1d 1031
1032 fIncpTMCpho->Fill(phoval);
1033 if(fFlagPhotonicElec) fPhoElecPtMC->Fill(phoval);
1034 if(fFlagConvinatElec) fSameElecPtMC->Fill(phoval);
e7d87aef 1035
1036 if(m20>0.0 && m20<0.3)
1037 {
44be9e1d 1038 fIncpTMCM20pho->Fill(phoval);
1039 if(fFlagPhotonicElec) fPhoElecPtMCM20->Fill(phoval);
1040 if(fFlagConvinatElec) fSameElecPtMCM20->Fill(phoval);
e4b0faf2 1041 // pi0->g->e
1042 if(mcWeight>-1)
1043 {
987053ce 1044 if(iHijing==1)mcWeight = 1.0;
93c189c5 1045 if(mcOrgPi0)
1046 {
697ecf6b 1047 fIncpTMCM20pho_pi0e->Fill(phoval,mcWeight);
1048 if(fFlagPhotonicElec) fPhoElecPtMCM20_pi0e->Fill(phoval,mcWeight);
1049 if(fFlagConvinatElec) fSameElecPtMCM20_pi0e->Fill(phoval,mcWeight);
05ac9660 1050
1051 // check production vertex
7d495dc5 1052 //fPhoVertexReco_EMCal->Fill(track->Pt(),conv_proR);
1053 //if(fFlagPhotonicElec)fPhoVertexReco_Invmass->Fill(track->Pt(),conv_proR);
1054 fPhoVertexReco_EMCal->Fill(track->Pt(),conv_proR,mcWeight);
1055 if(fFlagPhotonicElec)fPhoVertexReco_Invmass->Fill(track->Pt(),conv_proR,mcWeight);
05ac9660 1056
1057 }
bd50bf24 1058 // --- eta
93c189c5 1059 if(mcOrgEta)
1060 {
697ecf6b 1061 fIncpTMCM20pho_eta->Fill(phoval,mcWeight);
1062 if(fFlagPhotonicElec) fPhoElecPtMCM20_eta->Fill(phoval,mcWeight);
1063 if(fFlagConvinatElec) fSameElecPtMCM20_eta->Fill(phoval,mcWeight);
93c189c5 1064 }
bd50bf24 1065 // check production vertex
e4b0faf2 1066 }
e7d87aef 1067 }
1068 }
1069 }
bfc7c23b 1070 }
1071 PostData(1, fOutputList);
1072}
1073//_________________________________________
1074void AliAnalysisTaskHFECal::UserCreateOutputObjects()
1075{
42c75692 1076 //--- Check MC
1077
e7d87aef 1078 //Bool_t mcData = kFALSE;
42c75692 1079 if(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler())
1080 {
e7d87aef 1081 fmcData = kTRUE;
42c75692 1082 printf("+++++ MC Data available");
1083 }
e7d87aef 1084 if(fmcData)
42c75692 1085 {
1086 printf("++++++++= MC analysis \n");
1087 }
1088 else
1089 {
1090 printf("++++++++ real data analysis \n");
1091 }
1092
85985bb0 1093 printf("+++++++ QA hist %d",fqahist);
1094
fb87d707 1095 //---- Geometry
1096 fGeom = AliEMCALGeometry::GetInstance("EMCAL_COMPLETEV1");
1097
bfc7c23b 1098 //--------Initialize PID
42c75692 1099 //fPID->SetHasMCData(kFALSE);
e7d87aef 1100 fPID->SetHasMCData(fmcData);
bfc7c23b 1101 if(!fPID->GetNumberOfPIDdetectors())
1102 {
1103 fPID->AddDetector("TPC", 0);
89f41a30 1104 //fPID->AddDetector("EMCAL", 1); <--- apply PID selection in task
bfc7c23b 1105 }
1106
3db00c72 1107 Double_t params[4];
78ff1095 1108 const char *cutmodel;
3db00c72 1109 cutmodel = "pol0";
1110 params[0] = -1.0; //sigma min
9b3495ff 1111 double maxnSig = 3.0;
1112 if(fmcData)
1113 {
1114 params[0] = -5.0; //sigma min
1115 maxnSig = 5.0;
1116 }
1117
1118 for(Int_t a=0;a<11;a++)fPID->ConfigureTPCcentralityCut(a,cutmodel,params,maxnSig);
3db00c72 1119
bfc7c23b 1120 fPID->SortDetectors();
1121 fPIDqa = new AliHFEpidQAmanager();
1122 fPIDqa->Initialize(fPID);
a6df418c 1123
1124 //------- fcut --------------
1125 fCuts = new AliHFEcuts();
1126 fCuts->CreateStandardCuts();
1127 //fCuts->SetMinNClustersTPC(100);
1128 fCuts->SetMinNClustersTPC(90);
1129 fCuts->SetMinRatioTPCclusters(0.6);
1130 fCuts->SetTPCmodes(AliHFEextraCuts::kFound, AliHFEextraCuts::kFoundOverFindable);
1131 //fCuts->SetMinNClustersITS(3);
1132 fCuts->SetMinNClustersITS(2);
05ac9660 1133 fCuts->SetProductionVertex(0,50,0,50);
a6df418c 1134 fCuts->SetCutITSpixel(AliHFEextraCuts::kAny);
1135 fCuts->SetCheckITSLayerStatus(kFALSE);
1136 fCuts->SetVertexRange(10.);
1137 fCuts->SetTOFPIDStep(kFALSE);
fbe402d7 1138 //fCuts->SetPtRange(2, 50);
1139 fCuts->SetPtRange(0.1, 50);
0f566acb 1140 //fCuts->SetMaxImpactParam(3.,3.);
1141 fCuts->SetMaxImpactParam(2.4,3.2); // standard in 2011
a6df418c 1142
bfc7c23b 1143 //--------Initialize correction Framework and Cuts
1144 fCFM = new AliCFManager;
1145 const Int_t kNcutSteps = AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kNcutStepsRecTrack + AliHFEcuts::kNcutStepsDETrack;
1146 fCFM->SetNStepParticle(kNcutSteps);
1147 for(Int_t istep = 0; istep < kNcutSteps; istep++)
1148 fCFM->SetParticleCutsList(istep, NULL);
1149
1150 if(!fCuts){
1151 AliWarning("Cuts not available. Default cuts will be used");
1152 fCuts = new AliHFEcuts;
1153 fCuts->CreateStandardCuts();
1154 }
1155 fCuts->Initialize(fCFM);
1156
1157 //---------Output Tlist
1158 fOutputList = new TList();
1159 fOutputList->SetOwner();
1160 fOutputList->Add(fPIDqa->MakeList("PIDQA"));
1161
1162 fNoEvents = new TH1F("fNoEvents","",4,-0.5,3.5) ;
1163 fOutputList->Add(fNoEvents);
1164
02cd82a0 1165 Int_t binsE[5] = {250, 100, 1000, 200, 10};
fb87d707 1166 Double_t xminE[5] = {1.0, -1, 0.0, 0, -0.5};
1167 Double_t xmaxE[5] = {3.5, 1, 100.0, 100, 9.5};
1168 fEMCAccE = new THnSparseD("fEMCAccE","EMC acceptance & E;#eta;#phi;Energy;Centrality;trugCondition;",5,binsE,xminE,xmaxE);
c82f1b9d 1169 //if(fqahist==1)fOutputList->Add(fEMCAccE);
bfc7c23b 1170
f4e0d2d5 1171 hEMCAccE = new TH2F("hEMCAccE","Cluster Energy",200,0,100,100,0,20);
1172 fOutputList->Add(hEMCAccE);
1173
bfc7c23b 1174 fTrkpt = new TH1F("fTrkpt","track pt",100,0,50);
1175 fOutputList->Add(fTrkpt);
1176
1177 fTrackPtBefTrkCuts = new TH1F("fTrackPtBefTrkCuts","track pt before track cuts",100,0,50);
1178 fOutputList->Add(fTrackPtBefTrkCuts);
1179
1180 fTrackPtAftTrkCuts = new TH1F("fTrackPtAftTrkCuts","track pt after track cuts",100,0,50);
1181 fOutputList->Add(fTrackPtAftTrkCuts);
1182
1183 fTPCnsigma = new TH2F("fTPCnsigma", "TPC - n sigma",100,0,50,200,-10,10);
1184 fOutputList->Add(fTPCnsigma);
1185
1186 fTrkEovPBef = new TH2F("fTrkEovPBef","track E/p before HFE pid",100,0,50,100,0,2);
1187 fOutputList->Add(fTrkEovPBef);
1188
1189 fTrkEovPAft = new TH2F("fTrkEovPAft","track E/p after HFE pid",100,0,50,100,0,2);
1190 fOutputList->Add(fTrkEovPAft);
1191
42c75692 1192 fdEdxBef = new TH2F("fdEdxBef","track dEdx vs p before HFE pid",100,0,50,200,-10,10);
bfc7c23b 1193 fOutputList->Add(fdEdxBef);
1194
42c75692 1195 fdEdxAft = new TH2F("fdEdxAft","track dEdx vs p after HFE pid",100,0,50,200,-10,10);
bfc7c23b 1196 fOutputList->Add(fdEdxAft);
1197
02cd82a0 1198 fIncpT = new TH2F("fIncpT","HFE pid electro vs. centrality",200,0,100,100,0,50);
bfc7c23b 1199 fOutputList->Add(fIncpT);
1200
02cd82a0 1201 fIncpTM20 = new TH2F("fIncpTM20","HFE pid electro vs. centrality with M20",200,0,100,100,0,50);
fb87d707 1202 fOutputList->Add(fIncpTM20);
4e01b68c 1203
5e38d973 1204 Int_t nBinspho[9] = { 10, 30, 600, 60, 50, 4, 40, 8, 30};
1205 Double_t minpho[9] = { 0., 0., -0.1, 40, 0, -0.5, 0,-1.5, 0};
1206 Double_t maxpho[9] = {100., 30., 0.5, 100, 1, 3.5, 2, 6.5, 30};
f09766dd 1207
b12bc641 1208 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);
5cc5103f 1209 if(fqahist==1)fOutputList->Add(fInvmassLS);
bfc7c23b 1210
b12bc641 1211 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);
5cc5103f 1212 if(fqahist==1)fOutputList->Add(fInvmassULS);
bfc7c23b 1213
8806ce6c 1214 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);
5cc5103f 1215 if(fqahist==1)fOutputList->Add(fInvmassLSmc);
8806ce6c 1216
1217 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);
5cc5103f 1218 if(fqahist==1)fOutputList->Add(fInvmassULSmc);
8806ce6c 1219
d89866fd 1220 fInvmassLSreco = new TH2D("fInvmassLSreco", "Inv mass of LS (e,e) reco; cent; p_{T} (GeV/c); mass(GeV/c^2)",20,0,20,600,-0.1,0.5 );
6ec70f8c 1221 fInvmassLSreco->Sumw2();
1222 fOutputList->Add(fInvmassLSreco);
1223
d89866fd 1224 fInvmassULSreco = new TH2D("fInvmassULSreco", "Inv mass of ULS (e,e) reco; cent; p_{T} (GeV/c); mass(GeV/c^2)",20,0,20,600,-0.1,0.5 );
6ec70f8c 1225 fInvmassULSreco->Sumw2();
1226 fOutputList->Add(fInvmassULSreco);
1227
d89866fd 1228 fInvmassLSmc0 = new TH2D("fInvmassLSmc0", "Inv mass of LS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2)",20,0,20,600,-0.1,0.5 );
93c189c5 1229 fInvmassLSmc0->Sumw2();
1230 fOutputList->Add(fInvmassLSmc0);
1231
d89866fd 1232 fInvmassLSmc1 = new TH2D("fInvmassLSmc1", "Inv mass of LS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2)",20,0,20,600,-0.1,0.5 );
93c189c5 1233 fInvmassLSmc1->Sumw2();
1234 fOutputList->Add(fInvmassLSmc1);
1235
d89866fd 1236 fInvmassLSmc2 = new TH2D("fInvmassLSmc2", "Inv mass of LS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2)",20,0,20,600,-0.1,0.5 );
93c189c5 1237 fInvmassLSmc2->Sumw2();
1238 fOutputList->Add(fInvmassLSmc2);
1239
d89866fd 1240 fInvmassLSmc3 = new TH2D("fInvmassLSmc3", "Inv mass of LS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2)",20,0,20,600,-0.1,0.5 );
93c189c5 1241 fInvmassLSmc3->Sumw2();
1242 fOutputList->Add(fInvmassLSmc3);
1243
d89866fd 1244 fInvmassULSmc0 = new TH2D("fInvmassULSmc0", "Inv mass of ULS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2)",20,0,20,600,-0.1,0.5 );
93c189c5 1245 fInvmassULSmc0->Sumw2();
1246 fOutputList->Add(fInvmassULSmc0);
1247
d89866fd 1248 fInvmassULSmc1 = new TH2D("fInvmassULSmc1", "Inv mass of ULS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2)",20,0,20,600,-0.1,0.5 );
93c189c5 1249 fInvmassULSmc1->Sumw2();
1250 fOutputList->Add(fInvmassULSmc1);
1251
d89866fd 1252 fInvmassULSmc2 = new TH2D("fInvmassULSmc2", "Inv mass of ULS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2)",20,0,20,600,-0.1,0.5 );
93c189c5 1253 fInvmassULSmc2->Sumw2();
1254 fOutputList->Add(fInvmassULSmc2);
1255
d89866fd 1256 fInvmassULSmc3 = new TH2D("fInvmassULSmc3", "Inv mass of ULS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2)",20,0,20,600,-0.1,0.5 );
93c189c5 1257 fInvmassULSmc3->Sumw2();
93c189c5 1258 fOutputList->Add(fInvmassULSmc3);
1259
bfc7c23b 1260 fOpeningAngleLS = new TH1F("fOpeningAngleLS","Opening angle for LS pairs",100,0,1);
1261 fOutputList->Add(fOpeningAngleLS);
1262
1263 fOpeningAngleULS = new TH1F("fOpeningAngleULS","Opening angle for ULS pairs",100,0,1);
1264 fOutputList->Add(fOpeningAngleULS);
1265
1266 fPhotoElecPt = new TH1F("fPhotoElecPt", "photonic electron pt",100,0,50);
1267 fOutputList->Add(fPhotoElecPt);
1268
02cd82a0 1269 fPhoElecPt = new TH2F("fPhoElecPt", "Pho-inclusive electron pt",200,0,100,100,0,50);
bfc7c23b 1270 fOutputList->Add(fPhoElecPt);
1271
02cd82a0 1272 fPhoElecPtM20 = new TH2F("fPhoElecPtM20", "Pho-inclusive electron pt with M20",200,0,100,100,0,50);
fb87d707 1273 fOutputList->Add(fPhoElecPtM20);
1274
0f566acb 1275 fPhoElecPtM20Mass = new TH2F("fPhoElecPtM20Mass", "Pho-inclusive electron pt with M20 Mass",200,0,100,100,0,50);
1276 fPhoElecPtM20Mass->Sumw2();
1277 fOutputList->Add(fPhoElecPtM20Mass);
1278
02cd82a0 1279 fSameElecPt = new TH2F("fSameElecPt", "Same-inclusive electron pt",200,0,100,100,0,50);
f09766dd 1280 fOutputList->Add(fSameElecPt);
1281
02cd82a0 1282 fSameElecPtM20 = new TH2F("fSameElecPtM20", "Same-inclusive electron pt with M20",200,0,100,100,0,50);
fb87d707 1283 fOutputList->Add(fSameElecPtM20);
1284
0f566acb 1285 fSameElecPtM20Mass = new TH2F("fSameElecPtM20Mass", "Same-inclusive electron pt with M20 Mass",200,0,100,100,0,50);
1286 fSameElecPtM20Mass->Sumw2();
1287 fOutputList->Add(fSameElecPtM20Mass);
1288
ebef1da9 1289 fSemiElecPtM20 = new TH2F("fSemiElecPtM20", "Semi-inclusive electron pt with M20",200,0,100,100,0,50);
1290 fOutputList->Add(fSemiElecPtM20);
1291
02cd82a0 1292 fCent = new TH1F("fCent","Centrality",200,0,100) ;
bfc7c23b 1293 fOutputList->Add(fCent);
1294
1295 // Make common binning
44be9e1d 1296 const Double_t kMinP = 0.;
d19af75d 1297 const Double_t kMaxP = 20.;
bfc7c23b 1298
05ac9660 1299 //+++ 1st histogram: TPC dEdx with/without EMCAL (p, pT, TPC Signal, phi, eta, Sig, e/p, ,match, cell, M02, M20, Disp, Centrality, select)
1300 // 1st histogram: TPC dEdx with/without EMCAL (p, pT, TPC Signal, phi, eta, Sig, e/p, ,match, dca_xy, dca_z, M20, Disp, Centrality, select)
1301 Int_t nBins[16] = { 100, 7, 60, 20, 90, 100, 25, 60, 60, 100, 40, 10, 250, 100, 100, 8};
1302 Double_t min[16] = {kMinP, -0.5, 1.0, -1.0, -5.0, 0, 0, -3.0, -3.0, 0.0, 0.0, 0, 0, 0, 80, -1.5};
1303 Double_t max[16] = {kMaxP, 6.5, 4.0, 1.0, 4.0, 2.0, 0.05, 3.0, 3.0, 1.0, 20.0, 100, 100, 2.0, 180, 6.5};
1304 //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);
1305 fEleInfo = new THnSparseD("fEleInfo", "Electron Info; pT [GeV/c]; TPC signal;phi;eta;nSig; E/p;Rmatch;DCA_xy;DCA_z;M20;mcpT;Centrality;charge;opp;same;trigCond;MCele", 16, nBins, min, max);
f2fa5a00 1306 //if(fqahist==1)fOutputList->Add(fEleInfo);
bfc7c23b 1307
e1f0fb74 1308 // Make common binning
1309 Int_t nBinsEop[3] = { 10, 50, 100};
1310 Double_t minEop[3] = { 0, 0, -5};
1311 Double_t maxEop[3] = {100, 50, 5};
1312 fElenSigma= new THnSparseD("fElenSigma", "Electron nSigma; cent; pT [GeV/c]; nSigma", 3, nBinsEop, minEop, maxEop);
1313 fOutputList->Add(fElenSigma);
1314
1315
fb87d707 1316 //<--- Trigger info
feffe705 1317 /*
fb87d707 1318 fClsEBftTrigCut = new TH1F("fClsEBftTrigCut","cluster E before trigger selection",1000,0,100);
1319 fOutputList->Add(fClsEBftTrigCut);
1320
1321 fClsEAftTrigCut = new TH1F("fClsEAftTrigCut","cluster E if cls has 0 trigcut channel",1000,0,100);
1322 fOutputList->Add(fClsEAftTrigCut);
1323
1324 fClsEAftTrigCut1 = new TH1F("fClsEAftTrigCut1","cluster E if cls with trig channel",1000,0,100);
1325 fOutputList->Add(fClsEAftTrigCut1);
1326
1327 fClsEAftTrigCut2 = new TH1F("fClsEAftTrigCut2","cluster E if cls with trigcut channel",1000,0,100);
1328 fOutputList->Add(fClsEAftTrigCut2);
1329
1330 fClsEAftTrigCut3 = new TH1F("fClsEAftTrigCut3","cluster E if cls with trigcut channel + nCell>Ecorrect",1000,0,100);
1331 fOutputList->Add(fClsEAftTrigCut3);
1332
1333 fClsEAftTrigCut4 = new TH1F("fClsEAftTrigCut4","cluster E if cls with trigcut channel + nCell>Ecorrect + cls time cut",1000,0,100);
1334 fOutputList->Add(fClsEAftTrigCut4);
1335
1336 fClsETime = new TH2F("fClsETime", "Cls time vs E; E; time;",1000,0,100,1000,-0.0000002,0.0000009);
1337 fOutputList->Add(fClsETime);
1338
1339 fClsETime1 = new TH2F("fClsETime1", "Cls time vs E if cls contains trigger channel; E; time;",1000,0,100,1000,-0.0000002,0.0000009);
1340 fOutputList->Add(fClsETime1);
1341
1342 fTrigTimes = new TH1F("fTrigTimes", "Trigger time; time; N;",25,0,25);
1343 fOutputList->Add(fTrigTimes);
1344
42c75692 1345 fCellCheck = new TH2F("fCellCheck", "Cell vs E; E GeV; Cell ID",10,6,26,12000,0,12000);
1346 fOutputList->Add(fCellCheck);
feffe705 1347 */
e7d87aef 1348 //<---------- MC
1349
02cd82a0 1350 fInputHFEMC = new TH2F("fInputHFEMC","Input MC HFE pid electro vs. centrality",200,0,100,100,0,50);
e7d87aef 1351 fOutputList->Add(fInputHFEMC);
1352
02cd82a0 1353 fInputAlle = new TH2F("fInputAlle","Input MC electro vs. centrality",200,0,100,100,0,50);
feffe705 1354 fOutputList->Add(fInputAlle);
1355
02cd82a0 1356 fIncpTMChfe = new TH2F("fIncpTMChfe","MC HFE pid electro vs. centrality",200,0,100,100,0,50);
e7d87aef 1357 fOutputList->Add(fIncpTMChfe);
1358
02cd82a0 1359 fIncpTMChfeAll = new TH2F("fIncpTMChfeAll","MC Alle pid electro vs. centrality",200,0,100,100,0,50);
3db00c72 1360 fOutputList->Add(fIncpTMChfeAll);
1361
02cd82a0 1362 fIncpTMCM20hfe = new TH2F("fIncpTMCM20hfe","MC HFE pid electro vs. centrality with M20",200,0,100,100,0,50);
e7d87aef 1363 fOutputList->Add(fIncpTMCM20hfe);
1364
02cd82a0 1365 fIncpTMCM20hfeAll = new TH2F("fIncpTMCM20hfeAll","MC Alle pid electro vs. centrality with M20",200,0,100,100,0,50);
3db00c72 1366 fOutputList->Add(fIncpTMCM20hfeAll);
1367
60544aea 1368 fIncpTMCM20hfeCheck = new TH2F("fIncpTMCM20hfeCheck","MC HFE pid electro vs. centrality with M20 Check",200,0,100,100,0,50);
1369 fOutputList->Add(fIncpTMCM20hfeCheck);
44be9e1d 1370
2b4460a6 1371 Int_t nBinspho2[5] = { 200, 100, 7, 3, 700};
987053ce 1372 Double_t minpho2[5] = { 0., 0., -2.5, -0.5, 0.};
2b4460a6 1373 Double_t maxpho2[5] = {100., 50., 4.5, 2.5, 70.};
bd6ee6dd 1374
1375 fInputHFEMC_weight = new THnSparseD("fInputHFEMC_weight", "MC HFE electron pt",5,nBinspho2,minpho2,maxpho2);
1376 fOutputList->Add(fInputHFEMC_weight);
1377
1378 fIncpTMCM20hfeCheck_weight = new THnSparseD("fIncpTMCM20hfeCheck_weight", "HFE electron pt with M20",5,nBinspho2,minpho2,maxpho2);
1379 fOutputList->Add(fIncpTMCM20hfeCheck_weight);
1380
987053ce 1381 fIncpTMCpho = new THnSparseD("fIncpTMCpho","MC Pho pid electro vs. centrality",5,nBinspho2,minpho2,maxpho2);
e7d87aef 1382 fOutputList->Add(fIncpTMCpho);
1383
987053ce 1384 fIncpTMCM20pho = new THnSparseD("fIncpTMCM20pho","MC Pho pid electro vs. centrality with M20",5,nBinspho2,minpho2,maxpho2);
e7d87aef 1385 fOutputList->Add(fIncpTMCM20pho);
1386
987053ce 1387 fPhoElecPtMC = new THnSparseD("fPhoElecPtMC", "MC Pho-inclusive electron pt",5,nBinspho2,minpho2,maxpho2);
e7d87aef 1388 fOutputList->Add(fPhoElecPtMC);
1389
987053ce 1390 fPhoElecPtMCM20 = new THnSparseD("fPhoElecPtMCM20", "MC Pho-inclusive electron pt with M20",5,nBinspho2,minpho2,maxpho2);
feffe705 1391 fOutputList->Add(fPhoElecPtMCM20);
e7d87aef 1392
0f566acb 1393 fPhoElecPtMCM20Mass = new TH2D("fPhoElecPtMCM20Mass", "MC Pho-inclusive electron pt with M20 Mass",200,0,100,100,0,50);
1394 fOutputList->Add(fPhoElecPtMCM20Mass);
1395
987053ce 1396 fSameElecPtMC = new THnSparseD("fSameElecPtMC", "MC Same-inclusive electron pt",5,nBinspho2,minpho2,maxpho2);
e7d87aef 1397 fOutputList->Add(fSameElecPtMC);
1398
987053ce 1399 fSameElecPtMCM20 = new THnSparseD("fSameElecPtMCM20", "MC Same-inclusive electron pt with M20",5,nBinspho2,minpho2,maxpho2);
e7d87aef 1400 fOutputList->Add(fSameElecPtMCM20);
1401
0f566acb 1402 fSameElecPtMCM20Mass = new TH2D("fSameElecPtMCM20Mass", "MC Same-inclusive electron pt with M20 Mass",200,0,100,100,0,50);
1403 fOutputList->Add(fSameElecPtMCM20Mass);
1404
987053ce 1405 fIncpTMCM20pho_pi0e = new THnSparseD("fIncpTMCM20pho_pi0e","MC Pho pi0->e pid electro vs. centrality with M20",5,nBinspho2,minpho2,maxpho2);
1406 fIncpTMCM20pho_pi0e->Sumw2();
e4b0faf2 1407 fOutputList->Add(fIncpTMCM20pho_pi0e);
1408
987053ce 1409 fPhoElecPtMCM20_pi0e = new THnSparseD("fPhoElecPtMCM20_pi0e", "MC Pho-inclusive electron pt with M20 pi0->e",5,nBinspho2,minpho2,maxpho2);
1410 fPhoElecPtMCM20_pi0e->Sumw2();
e4b0faf2 1411 fOutputList->Add(fPhoElecPtMCM20_pi0e);
1412
987053ce 1413 fSameElecPtMCM20_pi0e = new THnSparseD("fSameElecPtMCM20_pi0e", "MC Same-inclusive electron pt pi0->e",5,nBinspho2,minpho2,maxpho2);
1414 fSameElecPtMCM20_pi0e->Sumw2();
e4b0faf2 1415 fOutputList->Add(fSameElecPtMCM20_pi0e);
697ecf6b 1416 //
93c189c5 1417 fIncpTMCM20pho_eta = new THnSparseD("fIncpTMCM20pho_eta","MC Pho pi0->e pid electro vs. centrality with M20",5,nBinspho2,minpho2,maxpho2);
1418 fIncpTMCM20pho_eta->Sumw2();
1419 fOutputList->Add(fIncpTMCM20pho_eta);
1420
1421 fPhoElecPtMCM20_eta = new THnSparseD("fPhoElecPtMCM20_eta", "MC Pho-inclusive electron pt with M20 pi0->e",5,nBinspho2,minpho2,maxpho2);
1422 fPhoElecPtMCM20_eta->Sumw2();
1423 fOutputList->Add(fPhoElecPtMCM20_eta);
1424
1425 fSameElecPtMCM20_eta = new THnSparseD("fSameElecPtMCM20_eta", "MC Same-inclusive electron pt pi0->e",5,nBinspho2,minpho2,maxpho2);
1426 fSameElecPtMCM20_eta->Sumw2();
1427 fOutputList->Add(fSameElecPtMCM20_eta);
697ecf6b 1428 // ------------
1429 fIncpTMCpho_pi0e_TPC = new THnSparseD("fIncpTMCpho_pi0e_TPC","MC Pho pi0->e pid electro vs. centrality with M20",5,nBinspho2,minpho2,maxpho2);
1430 fIncpTMCpho_pi0e_TPC->Sumw2();
1431 fOutputList->Add(fIncpTMCpho_pi0e_TPC);
1432
1433 fPhoElecPtMC_pi0e_TPC = new THnSparseD("fPhoElecPtMC_pi0e_TPC", "MC Pho-inclusive electron pt with pi0->e",5,nBinspho2,minpho2,maxpho2);
1434 fPhoElecPtMC_pi0e_TPC->Sumw2();
1435 fOutputList->Add(fPhoElecPtMC_pi0e_TPC);
1436
1437 fSameElecPtMC_pi0e_TPC = new THnSparseD("fSameElecPtMC_pi0e_TPC", "MC Same-inclusive electron pt pi0->e",5,nBinspho2,minpho2,maxpho2);
1438 fSameElecPtMC_pi0e_TPC->Sumw2();
1439 fOutputList->Add(fSameElecPtMC_pi0e_TPC);
2f8ef315 1440
1441 fIncpTMCpho_eta_TPC = new THnSparseD("fIncpTMCpho_eta_TPC","MC Pho pi0->e pid electro vs. centrality with M20",5,nBinspho2,minpho2,maxpho2);
1442 fIncpTMCpho_eta_TPC->Sumw2();
1443 fOutputList->Add(fIncpTMCpho_eta_TPC);
1444
1445 fPhoElecPtMC_eta_TPC = new THnSparseD("fPhoElecPtMC_eta_TPC", "MC Pho-inclusive electron pt with pi0->e",5,nBinspho2,minpho2,maxpho2);
1446 fPhoElecPtMC_eta_TPC->Sumw2();
1447 fOutputList->Add(fPhoElecPtMC_eta_TPC);
1448
1449 fSameElecPtMC_eta_TPC = new THnSparseD("fSameElecPtMC_eta_TPC", "MC Same-inclusive electron pt pi0->e",5,nBinspho2,minpho2,maxpho2);
1450 fSameElecPtMC_eta_TPC->Sumw2();
1451 fOutputList->Add(fSameElecPtMC_eta_TPC);
1452
1453
697ecf6b 1454 //-------------
e4b0faf2 1455
a6df418c 1456 CheckNclust = new TH1D("CheckNclust","cluster check",200,0,200);
1457 fOutputList->Add(CheckNclust);
1458
1459 CheckNits = new TH1D("CheckNits","ITS cluster check",8,-0.5,7.5);
1460 fOutputList->Add(CheckNits);
ae88597c 1461
1462 CheckDCA = new TH2D("CheckDCA","DCA check",200,-5,5,200,-5,5);
1463 fOutputList->Add(CheckDCA);
68d718e7 1464 /*
d113d7cd 1465 Hpi0pTcheck = new TH2D("Hpi0pTcheck","Pi0 pT from Hijing",100,0,50,3,-0.5,2.5);
1466 fOutputList->Add(Hpi0pTcheck);
1467
19325033 1468 HETApTcheck = new TH2D("HETApTcheck","Eta pT from Hijing",100,0,50,3,-0.5,2.5);
1469 fOutputList->Add(HETApTcheck);
68d718e7 1470 */
1471
1472 Int_t nBinspho3[3] = { 200, 100, 3};
1473 Double_t minpho3[3] = { 0., 0., -0.5};
1474 Double_t maxpho3[3] = {100., 50., 2.5};
19325033 1475
68d718e7 1476 Hpi0pTcheck = new THnSparseD("Hpi0pTcheck","Pi0 pT from Hijing",3,nBinspho3,minpho3,maxpho3);
1477 fOutputList->Add(Hpi0pTcheck);
1478
1479 HETApTcheck = new THnSparseD("HETApTcheck","Eta pT from Hijing",3,nBinspho3,minpho3,maxpho3);
1480 fOutputList->Add(HETApTcheck);
1481 //--
e4b0faf2 1482 HphopTcheck = new TH2D("HphopTcheck","Pho pT from Hijing",100,0,50,3,-0.5,2.5);
1483 fOutputList->Add(HphopTcheck);
5e0e45b3 1484 //
1485 HDpTcheck = new TH2D("HDpTcheck","D pT from Hijing",100,0,50,3,-0.5,2.5);
1486 fOutputList->Add(HDpTcheck);
1487
1488 HBpTcheck = new TH2D("HBpTcheck","B pT from Hijing",100,0,50,3,-0.5,2.5);
1489 fOutputList->Add(HBpTcheck);
1490 //
e4b0faf2 1491
93c189c5 1492 fpTCheck = new TH1D("fpTCheck","pT check",500,0,50);
1493 fOutputList->Add(fpTCheck);
1494
50919258 1495 fMomDtoE = new TH2D("fMomDtoE","D->E pT correlations;e p_{T} GeV/c;D p_{T} GeV/c",400,0,40,400,0,40);
1496 fOutputList->Add(fMomDtoE);
d113d7cd 1497
70448e3b 1498 fLabelCheck = new TH2D("fLabelCheck","MC label",50,0,50,5,-1.5,3.5);
1499 fOutputList->Add(fLabelCheck);
1500
1501 fgeoFake = new TH2D("fgeoFake","Label==0 eta and phi",628,0,6.28,200,-1,1);
1502 fOutputList->Add(fgeoFake);
1503
d19af75d 1504 fFakeTrk0 = new TH2D("fFakeTrk0","fake trakcs",10,0,100,20,0,20);
1505 fOutputList->Add(fFakeTrk0);
1506
1507 fFakeTrk1 = new TH2D("fFakeTrk1","true all e a.f. eID",10,0,100,20,0,20);
1508 fOutputList->Add(fFakeTrk1);
1509
59d998de 1510 ftimingEle = new TH2D("ftimingEle","electron TOF",100,0,20,100,1e-7,1e-6);
1511 fOutputList->Add(ftimingEle);
1512
09022877 1513 // eta correction
1514 // note: parameters 01/31new.TPCnSigmaEtaDep
1515 // 70-90 delta_eta = 0.2
1516
1517 double etaval[12] = {-0.55,-0.45,-0.35,-0.25,-0.15,-0.05,0.05,0.15,0.25,0.35,0.45,0.55};
1518 double corr0[12]= {-0.569177,-0.528844,-0.391979,-0.165494,0.0283495,0.156171,0.266353,0.13103,-0.0250842,-0.274089,-0.45481,-0.536291}; // 0-10 (done)
1519 double corr1[12]= {-0.404742,-0.278953,-0.218069,0.00139927,0.191412,0.354403,0.524594,0.341778,0.244199,-0.112146,-0.160692,-0.352832}; // 10-20 (done)
1520 double corr2[12] = {-0.306007,-0.16821,-0.0248635,0.202233,0.447051,0.497197,0.712251,0.433482,0.337907,0.168426,-0.0693229,-0.0728351}; // 20-30 (done)
1521 double corr3[12] = {-0.13884,-0.0503553,0.104403,0.389773,0.50697,0.539048,0.751642,0.655636,0.518563,0.308156,0.0361159,-0.0491439}; // 30-40 (done)
1522 double corr4[12] = {-0.0319431,0.0808711,0.208774,0.443217,0.557762,0.61453,0.889519,0.808282,0.620394,0.267092,0.15241,-0.0458664}; // 40-50 (done)
1523 double corr5[12] = {-0.130625,0.0189124,0.190344,0.467431,0.546353,0.672251,0.731541,0.802101,0.437108,0.294081,0.193682,0.159074}; // 50-70(done)
1524 double corr6[12] = {0.0600197,0.0600197,0.358366,0.358366,0.973734,0.973734,0.759812,0.759812,0.667861,0.667861,0.415635,0.415635}; // 70-90(done)
1525
1526 fnSigEtaCorr[0] = new TGraphErrors(12,etaval,corr0); // 0-10
1527 fnSigEtaCorr[1] = new TGraphErrors(12,etaval,corr1); // 10-20
1528 fnSigEtaCorr[2] = new TGraphErrors(12,etaval,corr2); // 20-30
1529 fnSigEtaCorr[3] = new TGraphErrors(12,etaval,corr3); // 30-40
1530 fnSigEtaCorr[4] = new TGraphErrors(12,etaval,corr4); // 40-50
1531 fnSigEtaCorr[5] = new TGraphErrors(12,etaval,corr5); // 50-70
1532 fnSigEtaCorr[6] = new TGraphErrors(12,etaval,corr6); // 70-90
1533
6f4a9952 1534 fIncMaxE = new TH2D("fIncMaxE","Inc",10,0,100,10,0,100);
1535 fOutputList->Add(fIncMaxE);
09022877 1536
c82f1b9d 1537 fIncReco = new TH2D("fIncReco","Inc",10,0,100,100,0,500);
5e38d973 1538 fOutputList->Add(fIncReco);
1539
c82f1b9d 1540 fPhoReco = new TH2D("fPhoReco","Pho",10,0,100,100,0,500);
1541 fOutputList->Add(fPhoReco);
1542
1543 fSamReco = new TH2D("fSamReco","Same",10,0,100,100,0,500);
1544 fOutputList->Add(fSamReco);
1545
1546 fIncRecoMaxE = new TH2D("fIncRecoMaxE","Inc",10,0,100,100,0,500);
1547 fOutputList->Add(fIncRecoMaxE);
1548
1549 fPhoRecoMaxE = new TH2D("fPhoRecoMaxE","Pho",10,0,100,100,0,500);
eabc1527 1550 fOutputList->Add(fPhoRecoMaxE);
5e38d973 1551
c82f1b9d 1552 fSamRecoMaxE = new TH2D("fSamRecoMaxE","Same",10,0,100,100,0,500);
eabc1527 1553 fOutputList->Add(fSamRecoMaxE);
5e38d973 1554
466d4f7a 1555 fPhoVertexReco_TPC = new TH2D("fPhoVertexReco_TPC","photon production Vertex mass selection TPC",40,0,20,250,0,50);
1556 fPhoVertexReco_TPC->Sumw2();
1557 fOutputList->Add(fPhoVertexReco_TPC);
1558
1559 fPhoVertexReco_TPC_Invmass = new TH2D("fPhoVertexReco_TPC_Invmass","photon production Vertex mass selection TPC Invmass",40,0,20,250,0,50);
1560 fPhoVertexReco_TPC_Invmass->Sumw2();
1561 fOutputList->Add(fPhoVertexReco_TPC_Invmass);
05ac9660 1562
1563 fPhoVertexReco_EMCal = new TH2D("fPhoVertexReco_EMCal","photon production Vertex mass selection",40,0,20,250,0,50);
7d495dc5 1564 fPhoVertexReco_EMCal->Sumw2();
05ac9660 1565 fOutputList->Add(fPhoVertexReco_EMCal);
1566
1567 fPhoVertexReco_Invmass = new TH2D("fPhoVertexReco_Invmass","photon production Vertex mass selection",40,0,20,250,0,50);
7d495dc5 1568 fPhoVertexReco_Invmass->Sumw2();
05ac9660 1569 fOutputList->Add(fPhoVertexReco_Invmass);
1570
1571 fPhoVertexReco_step0= new TH2D("fPhoVertexReco_step0","photon production Vertex mass selection",40,0,20,250,0,50);
d26226a1 1572 fPhoVertexReco_step0->Sumw2();
05ac9660 1573 fOutputList->Add(fPhoVertexReco_step0);
1574
1575 fPhoVertexReco_step1= new TH2D("fPhoVertexReco_step1","photon production Vertex mass selection",40,0,20,250,0,50);
d26226a1 1576 fPhoVertexReco_step1->Sumw2();
05ac9660 1577 fOutputList->Add(fPhoVertexReco_step1);
1578
f2fa5a00 1579 fMatchV0_0 = new TH1D("fMatchV0_0","V0 match",100,0,20);
1580 fOutputList->Add(fMatchV0_0);
1581
1582 fMatchV0_1 = new TH1D("fMatchV0_1","V0 match",100,0,20);
1583 fOutputList->Add(fMatchV0_1);
1584
1585 fMatchMC_0 = new TH1D("fMatchMC_0","MC match",100,0,20);
1586 fOutputList->Add(fMatchMC_0);
1587
1588 fMatchMC_1 = new TH1D("fMatchMC_1","MC match",100,0,20);
1589 fOutputList->Add(fMatchMC_1);
1590
36e9774f 1591 fpair = new TH2D("fpair","pair of associate",100,0,20,21,-10.5,10.5);
1592 fOutputList->Add(fpair);
1593
c5a4ddbf 1594 fFakeRejection0 = new TH2D("fFakeRejection0","TPC PID",2,-0.5,1.5,100,0,20);
1c9a06c8 1595 fFakeRejection0->Sumw2();
c5a4ddbf 1596 fOutputList->Add(fFakeRejection0);
1597
1598 fFakeRejection1 = new TH2D("fFakeRejection1","TPC PID + Tr match",2,-0.5,1.5,100,0,20);
1c9a06c8 1599 fFakeRejection1->Sumw2();
c5a4ddbf 1600 fOutputList->Add(fFakeRejection1);
1601
1602 fFakeRejection2 = new TH2D("fFakeRejection2","TPC PID + Tr match + E/p",2,-0.5,1.5,100,0,20);
1c9a06c8 1603 fFakeRejection2->Sumw2();
c5a4ddbf 1604 fOutputList->Add(fFakeRejection2);
1605
bfc7c23b 1606 PostData(1,fOutputList);
1607}
1608
1609//________________________________________________________________________
1610void AliAnalysisTaskHFECal::Terminate(Option_t *)
1611{
1612 // Info("Terminate");
1613 AliAnalysisTaskSE::Terminate();
1614}
1615
1616//________________________________________________________________________
1617Bool_t AliAnalysisTaskHFECal::ProcessCutStep(Int_t cutStep, AliVParticle *track)
1618{
1619 // Check single track cuts for a given cut step
1620 const Int_t kMCOffset = AliHFEcuts::kNcutStepsMCTrack;
1621 if(!fCFM->CheckParticleCuts(cutStep + kMCOffset, track)) return kFALSE;
1622 return kTRUE;
1623}
1624//_________________________________________
0f566acb 1625void AliAnalysisTaskHFECal::SelectPhotonicElectron(Int_t itrack, Double_t cent, AliESDtrack *track, Bool_t &fFlagPhotonic, Bool_t &fFlagConvinat, Double_t nSig, Double_t shower, Double_t ep, Double_t mce, Double_t w, Int_t ibgevent, Bool_t tagpi0, Bool_t tageta, Int_t iCal)
bfc7c23b 1626{
1627 //Identify non-heavy flavour electrons using Invariant mass method
1628
fb87d707 1629 fTrackCuts->SetAcceptKinkDaughters(kFALSE);
feffe705 1630 fTrackCuts->SetRequireTPCRefit(kTRUE);
85985bb0 1631 fTrackCuts->SetRequireITSRefit(kTRUE);
feffe705 1632 fTrackCuts->SetEtaRange(-0.9,0.9);
bfc7c23b 1633 fTrackCuts->SetMaxChi2PerClusterTPC(3.5);
feffe705 1634 fTrackCuts->SetMinNClustersTPC(90);
bfc7c23b 1635
c6bca3f2 1636 //const AliESDVertex *pVtx = fESD->GetPrimaryVertex();
6451f693 1637 Double_t bfield = fESD->GetMagneticField();
2f8ef315 1638 Double_t emass = 0.51*0.001; // (0.51 MeV)
6451f693 1639
93c189c5 1640 double ptEle = track->Pt(); //add
1641 if(ibgevent==0 && w > 0.0)
1642 {
1643 fpTCheck->Fill(ptEle,w);
1644 }
1645
bfc7c23b 1646 Bool_t flagPhotonicElec = kFALSE;
f09766dd 1647 Bool_t flagConvinatElec = kFALSE;
bfc7c23b 1648
8806ce6c 1649 int p1 = 0;
1650 if(mce==3)
1651 {
1652 Int_t label = TMath::Abs(track->GetLabel());
1653 TParticle* particle = stack->Particle(label);
1654 p1 = particle->GetFirstMother();
1655 }
1656
36e9774f 1657 int numULS = 0;
1658 int numLS = 0;
1659
5e65985c 1660 for(Int_t jTracks = 0; jTracks<fESD->GetNumberOfTracks(); jTracks++){
bfc7c23b 1661 AliESDtrack* trackAsso = fESD->GetTrack(jTracks);
1662 if (!trackAsso) {
1663 printf("ERROR: Could not receive track %d\n", jTracks);
1664 continue;
1665 }
5e65985c 1666 if(itrack==jTracks)continue;
d86a815c 1667 int jbgevent = 0;
5e65985c 1668
8806ce6c 1669 int p2 = 0;
1670 if(mce==3)
1671 {
1672 Int_t label2 = TMath::Abs(trackAsso->GetLabel());
1673 TParticle* particle2 = stack->Particle(label2);
d86a815c 1674 Bool_t MChijing_ass = fMC->IsFromBGEvent(label2);
1675 if(MChijing_ass)jbgevent =1;
8806ce6c 1676 if(particle2->GetFirstMother()>-1)
1677 p2 = particle2->GetFirstMother();
1678 }
1679
f09766dd 1680 Double_t dEdxAsso = -999., ptPrim=-999., ptAsso=-999., openingAngle = -999.;
bfc7c23b 1681 Double_t mass=999., width = -999;
1682 Bool_t fFlagLS=kFALSE, fFlagULS=kFALSE;
1683
93c189c5 1684 //ptPrim = track->Pt();
1685 ptPrim = ptEle;
f09766dd 1686
bfc7c23b 1687 dEdxAsso = trackAsso->GetTPCsignal();
1688 ptAsso = trackAsso->Pt();
1689 Int_t chargeAsso = trackAsso->Charge();
1690 Int_t charge = track->Charge();
1691
b12bc641 1692
6451f693 1693 if(ptAsso <fMimpTassCut) continue;
bfc7c23b 1694 if(!fTrackCuts->AcceptTrack(trackAsso)) continue;
6b2cee73 1695 double fTPCnSigmaAss = fPID->GetPIDResponse() ? fPID->GetPIDResponse()->NumberOfSigmasTPC(trackAsso, AliPID::kElectron) : 1000;
1696 //cout << "fTPCnSigmaAss = " << fTPCnSigmaAss << endl;
1697 //cout << "fTPCnSigmaAss Cut = " << fMimNsigassCut << endl;
1698 if(fTPCnSigmaAss <fMimNsigassCut || fTPCnSigmaAss>5) continue;
1699 //cout << "fTPCnSigmaAss a.f. cut = " << fTPCnSigmaAss << endl;
bfc7c23b 1700
1701 Int_t fPDGe1 = 11; Int_t fPDGe2 = 11;
1702 if(charge>0) fPDGe1 = -11;
1703 if(chargeAsso>0) fPDGe2 = -11;
b12bc641 1704
1705 //printf("chargeAsso = %d\n",chargeAsso);
1706 //printf("charge = %d\n",charge);
bfc7c23b 1707 if(charge == chargeAsso) fFlagLS = kTRUE;
1708 if(charge != chargeAsso) fFlagULS = kTRUE;
1709
b12bc641 1710 //printf("fFlagLS = %d\n",fFlagLS);
1711 //printf("fFlagULS = %d\n",fFlagULS);
cccab38f 1712 //printf("\n");
b12bc641 1713
6451f693 1714 AliKFParticle::SetField(bfield);
bfc7c23b 1715 AliKFParticle ge1(*track, fPDGe1);
1716 AliKFParticle ge2(*trackAsso, fPDGe2);
2f8ef315 1717
1718 if(fSetMassNonlinear)
1719 {
1720 ge1.SetNonlinearMassConstraint(emass);
1721 ge2.SetNonlinearMassConstraint(emass);
1722 }
1723
bfc7c23b 1724 AliKFParticle recg(ge1, ge2);
1725
2f8ef315 1726 /*
d19af75d 1727 // vertex
bfc7c23b 1728 AliKFVertex primV(*pVtx);
1729 primV += recg;
8079a103 1730 primV -= ge1;
1731 primV -= ge2;
bfc7c23b 1732 recg.SetProductionVertex(primV);
2f8ef315 1733 */
1734
d19af75d 1735 // check chi2
1736 if(recg.GetNDF()<1) continue;
a009053a 1737 Double_t chi2recg = recg.GetChi2()/recg.GetNDF();
1738 Double_t chi2cut = 3.0;
f9447b3b 1739
a009053a 1740 // mass.
1741 if(fSetMassConstraint)
1742 {
1743 recg.SetMassConstraint(0,0.0001);
1744 chi2cut = 30.0;
1745 }
f9447b3b 1746 recg.GetMass(mass,width);
1747
3c475a9d 1748 if(fSetMassWidthCut && width>1e10)continue;
1749
a009053a 1750 // angle
bfc7c23b 1751 openingAngle = ge1.GetAngle(ge2);
1752 if(fFlagLS) fOpeningAngleLS->Fill(openingAngle);
1753 if(fFlagULS) fOpeningAngleULS->Fill(openingAngle);
0f566acb 1754
c2801a73 1755 double ishower = 0;
1756 if(shower>0.0 && shower<0.3)ishower = 1;
1757
0f566acb 1758 if(iCal==1)
1759 {
1760 double phoinfo[9];
1761 phoinfo[0] = cent;
1762 phoinfo[1] = ptPrim;
1763 phoinfo[2] = mass;
1764 phoinfo[3] = nSig;
1765 phoinfo[4] = openingAngle;
1766 phoinfo[5] = ishower;
1767 phoinfo[6] = ep;
1768 phoinfo[7] = mce;
1769 phoinfo[8] = ptAsso;
1770
1771 if(fFlagLS) fInvmassLS->Fill(phoinfo);
1772 if(fFlagULS) fInvmassULS->Fill(phoinfo);
1773 if(fFlagLS && ibgevent==0 && jbgevent==0) fInvmassLSmc->Fill(phoinfo,w);
1774 if(fFlagULS && ibgevent==0 && jbgevent==0)
1775 {
1776 fInvmassULSmc->Fill(phoinfo,w);
1777 }
1778 //printf("fInvmassCut %f\n",fInvmassCut);
1779 //printf("openingAngle %f\n",fOpeningAngleCut);
1780 }
b12bc641 1781
f9447b3b 1782 // angle cut
ffa14dda 1783 if(openingAngle > fOpeningAngleCut) continue;
a009053a 1784 // chi2 cut
6ec70f8c 1785 //if(TMath::Sqrt(TMath::Abs(chi2recg))>chi2cut) continue;
1786 if(chi2recg>chi2cut) continue;
1787
1788 if(fFlagLS ) fInvmassLSreco->Fill(ptPrim,mass);
1789 if(fFlagULS) fInvmassULSreco->Fill(ptPrim,mass);
1790
36e9774f 1791 // check double count
1792 if(mass<fInvmassCut && fFlagULS)numULS++;
1793 if(mass<fInvmassCut && fFlagLS)numLS++;
1794
d86a815c 1795 // for real data
93c189c5 1796 //printf("mce =%f\n",mce);
0f566acb 1797
d86a815c 1798 if(mce<-0.5) // mce==-1. is real
0f566acb 1799 {
1800 //printf("Real data\n");
1801
1802 // count
1803 if(mass<fInvmassCut && fFlagULS && shower>0.0 && shower<0.3 && iCal==1)fPhoElecPtM20Mass->Fill(cent,ptPrim);
1804 if(mass<fInvmassCut && fFlagLS && shower>0.0 && shower<0.3 && iCal==1)fSameElecPtM20Mass->Fill(cent,ptPrim);
1805
1806 // flag
1807 if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec){
1808 flagPhotonicElec = kTRUE;
1809 }
1810 if(mass<fInvmassCut && fFlagLS && !flagConvinatElec){
1811 flagConvinatElec = kTRUE;
1812 }
1813 }
d86a815c 1814 // for MC data
1815 else
0f566acb 1816 {
1817 // count
1818 if(mass<fInvmassCut && fFlagULS && shower>0.0 && shower<0.3 && mce>2.1 && iCal==1)fPhoElecPtMCM20Mass->Fill(cent,ptPrim,w);
1819 if(mass<fInvmassCut && fFlagLS && shower>0.0 && shower<0.3 && mce>2.1 && iCal==1)fSameElecPtMCM20Mass->Fill(cent,ptPrim,w);
1820
1821 //printf("MC data\n");
1822
1823 if(w>0.0)
1824 {
1825 //cout << "tagpi0 = " << tagpi0 << " ; tageta = " << tageta << endl;
1826 if(iCal==1)
1827 {
1828 if(fFlagLS && ibgevent==0 && jbgevent==0 && tagpi0) fInvmassLSmc0->Fill(ptPrim,mass);
1829 if(fFlagULS && ibgevent==0 && jbgevent==0 && tagpi0) fInvmassULSmc0->Fill(ptPrim,mass);
1830 if(fFlagLS && ibgevent==0 && jbgevent==0 && tageta) fInvmassLSmc1->Fill(ptPrim,mass);
1831 if(fFlagULS && ibgevent==0 && jbgevent==0 && tageta) fInvmassULSmc1->Fill(ptPrim,mass);
1832 if(fFlagLS && ibgevent==0 && jbgevent==0 && (p1==p2) && tagpi0) fInvmassLSmc2->Fill(ptPrim,mass);
1833 if(fFlagULS && ibgevent==0 && jbgevent==0 && (p1==p2) && tagpi0) fInvmassULSmc2->Fill(ptPrim,mass);
1834 if(fFlagLS && ibgevent==0 && jbgevent==0 && (p1==p2) && tageta) fInvmassLSmc3->Fill(ptPrim,mass);
1835 if(fFlagULS && ibgevent==0 && jbgevent==0 && (p1==p2) && tageta) fInvmassULSmc3->Fill(ptPrim,mass);
1836 }
1837 }
1838
1839 if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec && (ibgevent==jbgevent)){
1840 //if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec && (p1==p2)){ //<--- only MC train (55,56) v5-03-68-AN , 69 & v5-05-70-AN (till 74)
1841 flagPhotonicElec = kTRUE;
1842 }
1843 if(mass<fInvmassCut && fFlagLS && !flagConvinatElec && (ibgevent==jbgevent)){
1844 flagConvinatElec = kTRUE;
1845 }
1846 }
d86a815c 1847
36e9774f 1848 } // end of associate loop
1849
6e1d7a8a 1850 if(numULS>0 || numLS>0)
1851 {
1852 int numPair = numULS-numLS;
0f566acb 1853 if(iCal==1)fpair->Fill(ptEle,numPair);
6e1d7a8a 1854 }
1855
0f566acb 1856 fFlagPhotonic = flagPhotonicElec;
1857 fFlagConvinat = flagConvinatElec;
6451f693 1858
1859}
1860
8806ce6c 1861//-------------------------------------------
46305ed6 1862
1863void AliAnalysisTaskHFECal::FindMother(TParticle* part, int &label, int &pid)
1864{
1865 //int label = 99999;
1866 //int pid = 99999;
1867
1868 if(part->GetFirstMother()>-1)
1869 {
1870 label = part->GetFirstMother();
1871 pid = stack->Particle(label)->GetPdgCode();
1872 }
70448e3b 1873 //cout << "Find Mother : label = " << label << " ; pid" << pid << endl;
46305ed6 1874}
1875
8806ce6c 1876double AliAnalysisTaskHFECal::GetMCweight(double mcPi0pT)
1877{
1878 double weight = 1.0;
1879
1880 if(mcPi0pT>0.0 && mcPi0pT<5.0)
1881 {
1882 weight = 0.323*mcPi0pT/(TMath::Exp(-1.6+0.767*mcPi0pT+0.0285*mcPi0pT*mcPi0pT));
1883 }
1884 else
1885 {
1886 weight = 115.0/(0.718*mcPi0pT*TMath::Power(mcPi0pT,3.65));
1887 }
1888 return weight;
1889}
fb87d707 1890
93c189c5 1891double AliAnalysisTaskHFECal::GetMCweightEta(double mcEtapT)
1892{
1893 double weight = 1.0;
1894
1895 weight = 223.3/TMath::Power((TMath::Exp(-0.17*mcEtapT-0.0322*mcEtapT*mcEtapT)+mcEtapT/1.69),5.65);
1896 return weight;
1897}
1898
1899
fb87d707 1900//_________________________________________
1901void AliAnalysisTaskHFECal::FindTriggerClusters()
1902{
bd6ee6dd 1903 //cout << "finding trigger patch" << endl;
fb87d707 1904 // constants
1905 const int nModuleCols = 2;
1906 const int nModuleRows = 5;
1907 const int nColsFeeModule = 48;
1908 const int nRowsFeeModule = 24;
1909 const int nColsFaltroModule = 24;
1910 const int nRowsFaltroModule = 12;
1911 //const int faltroWidthMax = 20;
1912
1913 // part 1, trigger extraction -------------------------------------
1914 Int_t globCol, globRow;
1915 //Int_t ntimes=0, nTrigChannel=0, nTrigChannelCut=0, trigInCut=0;
1916 Int_t ntimes=0, nTrigChannel=0, nTrigChannelCut=0;
1917
1918 //Int_t trigtimes[faltroWidthMax];
1919 Double_t cellTime[nColsFeeModule*nModuleCols][nRowsFeeModule*nModuleRows];
1920 Double_t cellEnergy[nColsFeeModule*nModuleCols][nRowsFeeModule*nModuleRows];
1921 //Double_t fTrigCutLow = 6;
1922 //Double_t fTrigCutHigh = 10;
1923 Double_t fTimeCutLow = 469e-09;
1924 Double_t fTimeCutHigh = 715e-09;
1925
1926 AliESDCaloTrigger * fCaloTrigger = fESD->GetCaloTrigger( "EMCAL" );
1927
1928 // erase trigger maps
1929 for(Int_t i = 0; i < nColsFaltroModule*nModuleCols; i++ )
1930 {
1931 for(Int_t j = 0; j < nRowsFaltroModule*nModuleRows; j++ )
1932 {
1933 ftriggersCut[i][j] = 0;
1934 ftriggers[i][j] = 0;
1935 ftriggersTime[i][j] = 0;
1936 }
1937 }
1938
1939 Int_t iglobCol=0, iglobRow=0;
1940 // go through triggers
1941 if( fCaloTrigger->GetEntries() > 0 )
1942 {
1943 // needs reset
1944 fCaloTrigger->Reset();
1945 while( fCaloTrigger->Next() )
1946 {
1947 fCaloTrigger->GetPosition( globCol, globRow );
1948 fCaloTrigger->GetNL0Times( ntimes );
1949 /*
1950 // no L0s
1951 if( ntimes < 1 ) continue;
1952 // get precise timings
1953 fCaloTrigger->GetL0Times( trigtimes );
1954 trigInCut = 0;
1955 for(Int_t i = 0; i < ntimes; i++ )
1956 {
1957 // save the first trigger time in channel
1958 if( i == 0 || triggersTime[globCol][globRow] > trigtimes[i] )
1959 triggersTime[globCol][globRow] = trigtimes[i];
1960 //printf("trigger times: %d\n",trigtimes[i]);
1961 // check if in cut
1962 if(trigtimes[i] > fTrigCutLow && trigtimes[i] < fTrigCutHigh )
1963 trigInCut = 1;
1964
1965 fTrigTimes->Fill(trigtimes[i]);
1966 }
1967 */
1968
1969 //L1 analysis from AliAnalysisTaskEMCALTriggerQA
1970 Int_t bit = 0;
1971 fCaloTrigger->GetTriggerBits(bit);
bd6ee6dd 1972 //cout << "bit = " << bit << endl;
fb87d707 1973
1974 Int_t ts = 0;
1975 fCaloTrigger->GetL1TimeSum(ts);
bd6ee6dd 1976 //cout << "ts = " << ts << endl;
fb87d707 1977 if (ts > 0)ftriggers[globCol][globRow] = 1;
1978 // number of triggered channels in event
1979 nTrigChannel++;
1980 // ... inside cut
1981 if(ts>0 && (bit >> 6 & 0x1))
1982 {
1983 iglobCol = globCol;
1984 iglobRow = globRow;
1985 nTrigChannelCut++;
bd6ee6dd 1986 //cout << "ts cut = " << ts << endl;
1987 //cout << "globCol = " << globCol << endl;
1988 //cout << "globRow = " << globRow << endl;
fb87d707 1989 ftriggersCut[globCol][globRow] = 1;
1990 }
1991
1992 } // calo trigger entries
1993 } // has calo trigger entries
1994
1995 // part 2 go through the clusters here -----------------------------------
bd6ee6dd 1996 //cout << " part 2 go through the clusters here ----------------------------------- " << endl;
fb87d707 1997 Int_t nCluster=0, nCell=0, iCell=0, gCell=0;
60d77596 1998 Short_t cellAddr, nSACell;
1999 Int_t mclabel;
fb87d707 2000 //Int_t nSACell, iSACell, mclabel;
2001 Int_t iSACell;
42c75692 2002 Double_t cellAmp=0, cellTimeT=0, clusterTime=0, efrac=0;
fb87d707 2003 Int_t nSupMod, nModule, nIphi, nIeta, iphi, ieta, gphi, geta, feta, fphi;
2004
2005 TRefArray *fCaloClusters = new TRefArray();
2006 fESD->GetEMCALClusters( fCaloClusters );
2007 nCluster = fCaloClusters->GetEntries();
2008
2009
2010 // save all cells times if there are clusters
2011 if( nCluster > 0 ){
2012 // erase time map
2013 for(Int_t i = 0; i < nColsFeeModule*nModuleCols; i++ ){
2014 for(Int_t j = 0; j < nRowsFeeModule*nModuleRows; j++ ){
2015 cellTime[i][j] = 0.;
2016 cellEnergy[i][j] = 0.;
2017 }
2018 }
2019
2020 // get cells
2021 AliESDCaloCells *fCaloCells = fESD->GetEMCALCells();
2022 //AliVCaloCells fCaloCells = fESD->GetEMCALCells();
2023 nSACell = fCaloCells->GetNumberOfCells();
2024 for(iSACell = 0; iSACell < nSACell; iSACell++ ){
2025 // get the cell info *fCal
2026 fCaloCells->GetCell( iSACell, cellAddr, cellAmp, cellTimeT , mclabel, efrac);
2027
2028 // get cell position
2029 fGeom->GetCellIndex( cellAddr, nSupMod, nModule, nIphi, nIeta );
2030 fGeom->GetCellPhiEtaIndexInSModule( nSupMod,nModule, nIphi, nIeta, iphi, ieta);
2031
2032 // convert co global phi eta
2033 gphi = iphi + nRowsFeeModule*(nSupMod/2);
2034 geta = ieta + nColsFeeModule*(nSupMod%2);
2035
2036 // save cell time and energy
2037 cellTime[geta][gphi] = cellTimeT;
2038 cellEnergy[geta][gphi] = cellAmp;
2039
2040 }
2041 }
2042
2043 Int_t nClusterTrig, nClusterTrigCut;
2044 UShort_t *cellAddrs;
2045 Double_t clsE=-999, clsEta=-999, clsPhi=-999;
2046 Float_t clsPos[3] = {0.,0.,0.};
2047
2048 for(Int_t icl=0; icl<fESD->GetNumberOfCaloClusters(); icl++)
2049 {
2050 AliESDCaloCluster *cluster = fESD->GetCaloCluster(icl);
2051 if(!cluster || !cluster->IsEMCAL()) continue;
2052
2053 // get cluster cells
2054 nCell = cluster->GetNCells();
2055
2056 // get cluster energy
2057 clsE = cluster->E();
2058
2059 // get cluster position
2060 cluster->GetPosition(clsPos);
2061 TVector3 clsPosVec(clsPos[0],clsPos[1],clsPos[2]);
2062 clsEta = clsPosVec.Eta();
2063 clsPhi = clsPosVec.Phi();
2064
2065 // get the cell addresses
2066 cellAddrs = cluster->GetCellsAbsId();
2067
2068 // check if the cluster contains cell, that was marked as triggered
2069 nClusterTrig = 0;
2070 nClusterTrigCut = 0;
2071
2072 // loop the cells to check, if cluser in acceptance
2073 // any cluster with a cell outside acceptance is not considered
2074 for( iCell = 0; iCell < nCell; iCell++ )
2075 {
42c75692 2076 // check hot cell
feffe705 2077 //if(clsE>6.0)fCellCheck->Fill(clsE,cellAddrs[iCell]);
42c75692 2078
fb87d707 2079 // get cell position
2080 fGeom->GetCellIndex( cellAddrs[iCell], nSupMod, nModule, nIphi, nIeta );
2081 fGeom->GetCellPhiEtaIndexInSModule( nSupMod,nModule, nIphi, nIeta, iphi, ieta);
2082
2083 // convert co global phi eta
2084 gphi = iphi + nRowsFeeModule*(nSupMod/2);
2085 geta = ieta + nColsFeeModule*(nSupMod%2);
2086
2087 if( cellTime[geta][gphi] > 0. ){
2088 clusterTime += cellTime[geta][gphi];
2089 gCell++;
2090 }
2091
2092 // get corresponding FALTRO
2093 fphi = gphi / 2;
2094 feta = geta / 2;
2095
bd6ee6dd 2096 //cout << "fphi = " << fphi << endl;
2097 //cout << "feta = " << feta << endl;
2098
fb87d707 2099 // try to match with a triggered
2100 if( ftriggers[feta][fphi]==1)
2101 { nClusterTrig++;
2102 }
2103 if( ftriggersCut[feta][fphi]==1)
2104 { nClusterTrigCut++;
2105 }
2106
bd6ee6dd 2107 //cout << "nClusterTrigCut : " << nClusterTrigCut << endl;
2108
fb87d707 2109 } // cells
2110
2111
2112 if( gCell > 0 )
2113 clusterTime = clusterTime / (Double_t)gCell;
2114 // fix the reconstriction code time 100ns jumps
2115 if( fESD->GetBunchCrossNumber() % 4 < 2 )
2116 clusterTime -= 0.0000001;
2117
feffe705 2118 //fClsETime->Fill(clsE,clusterTime);
2119 //fClsEBftTrigCut->Fill(clsE);
fb87d707 2120
2121 if(nClusterTrig>0){
feffe705 2122 //fClsETime1->Fill(clsE,clusterTime);
fb87d707 2123 }
2124
2125 if(nClusterTrig>0){
2126 cluster->SetChi2(1);
feffe705 2127 //fClsEAftTrigCut1->Fill(clsE);
fb87d707 2128 }
2129
2130 if(nClusterTrigCut>0){
2131 cluster->SetChi2(2);
feffe705 2132 //fClsEAftTrigCut2->Fill(clsE);
fb87d707 2133 }
2134
2135 if(nClusterTrigCut>0 && ( nCell > (1 + clsE / 3)))
2136 {
2137 cluster->SetChi2(3);
feffe705 2138 //fClsEAftTrigCut3->Fill(clsE);
fb87d707 2139 }
2140
2141 if(nClusterTrigCut>0 && (nCell > (1 + clsE / 3) )&&( clusterTime > fTimeCutLow && clusterTime < fTimeCutHigh ))
2142 {
2143 // cluster->SetChi2(4);
feffe705 2144 //fClsEAftTrigCut4->Fill(clsE);
fb87d707 2145 }
2146 if(nClusterTrigCut<1)
2147 {
2148 cluster->SetChi2(0);
2149
feffe705 2150 //fClsEAftTrigCut->Fill(clsE);
fb87d707 2151 }
2152
2153 } // clusters
2154}
2155
09022877 2156// <-------- only MC correction
8efacc22 2157double AliAnalysisTaskHFECal::MCEopMeanCorrection(double pTmc, float central)
8af58b0d 2158{
8efacc22 2159 TF1 *fcorr0 = new TF1("fcorr0","[0]*tanh([1]+[2]*x)");
2160 TF1 *fcorr1 = new TF1("fcorr1","[0]*tanh([1]+[2]*x)");
2161
8af58b0d 2162 double shift = 0.0;
8efacc22 2163
2164 if(central>0 && central<=10)
2165 {
2166 fcorr0->SetParameters(1.045,1.288,3.18e-01); //
2167 fcorr1->SetParameters(9.91e-01,3.466,2.344);
2168 }
09022877 2169 else if(central>10 && central<=20)
8efacc22 2170 {
2171 fcorr0->SetParameters(1.029,8.254e-01,4.07e-01);
2172 fcorr1->SetParameters(0.975,2.276,1.501e-01);
2173 }
2174 else if(central>20 && central<=30)
2175 {
2176 fcorr0->SetParameters(1.01,8.795e-01,3.904e-01);
2177 fcorr1->SetParameters(9.675e-01,1.654,2.583e-01);
2178 }
2179 else if(central>30 && central<=40)
2180 {
2181 fcorr0->SetParameters(1.00,1.466,2.305e-1);
2182 fcorr1->SetParameters(9.637e-01,1.473,2.754e-01);
2183 }
2184 else if(central>40 && central<=50)
2185 {
fd2126aa 2186 fcorr0->SetParameters(1.00,1.422,1.518e-01);
8efacc22 2187 fcorr1->SetParameters(9.59e-01,1.421,2.931e-01);
2188 }
2189
2190 else if(central>50 && central<=70)
2191 {
2192 fcorr0->SetParameters(0.989,2.495,2.167);
e1f0fb74 2193 fcorr1->SetParameters(0.961,1.734,1.438e-01);
8efacc22 2194 }
2195 else if(central>70 && central<=100)
2196 {
2197 fcorr0->SetParameters(0.981,-3.373,3.93327);
fd2126aa 2198 fcorr1->SetParameters(9.574e-01,1.698,1.58e-01);
8efacc22 2199 }
2200
2201
2202 shift = fcorr0->Eval(pTmc)-fcorr1->Eval(pTmc);
2203
89f41a30 2204 fcorr0->Delete();
2205 fcorr1->Delete();
2206
8af58b0d 2207 return shift;
2208}
fb87d707 2209
09022877 2210// <-------- only Data correction
2211double AliAnalysisTaskHFECal::NsigmaCorrection(double tmpeta, float central)
2212{
2213 int icent = 0;
2214
2215 if(central>=0 && central<10)
2216 {
2217 icent = 0;
2218 }
2219 else if(central>=10 && central<20)
2220 {
2221 icent = 1;
2222 }
2223 else if(central>=20 && central<30)
2224 {
2225 icent = 2;
2226 }
2227 else if(central>=30 && central<40)
2228 {
2229 icent = 3;
2230 }
2231 else if(central>=40 && central<50)
2232 {
2233 icent = 4;
2234 }
2235 else if(central>=50 && central<70)
2236 {
2237 icent = 5;
2238 }
2239 else
2240 {
2241 icent = 6;
2242 }
2243
2244 double shift = fnSigEtaCorr[icent]->Eval(tmpeta);
2245
2246 //cout << "eta correction"<< endl;
2247 //cout << "cent = "<< central<< endl;
2248 //cout << "icent = "<< icent << endl;
2249 //cout << "shift = "<< shift << endl;
2250
2251 return shift;
fb87d707 2252
09022877 2253}
fb87d707 2254
2255