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