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