]> 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
[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)
05ac9660 216 ,fPhoVertexReco_HFE(0)
217 ,fPhoVertexReco_EMCal(0)
218 ,fPhoVertexReco_Invmass(0)
219 ,fPhoVertexReco_step0(0)
220 ,fPhoVertexReco_step1(0)
f2fa5a00 221 ,fMatchV0_0(0)
222 ,fMatchV0_1(0)
223 ,fMatchMC_0(0)
224 ,fMatchMC_1(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)
05ac9660 370 ,fPhoVertexReco_HFE(0)
371 ,fPhoVertexReco_EMCal(0)
372 ,fPhoVertexReco_Invmass(0)
373 ,fPhoVertexReco_step0(0)
374 ,fPhoVertexReco_step1(0)
f2fa5a00 375 ,fMatchV0_0(0)
376 ,fMatchV0_1(0)
377 ,fMatchMC_0(0)
378 ,fMatchMC_1(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;
fbe402d7 712 //if(TMath::Abs(track->Pt()<2.5)) continue;
713 if(TMath::Abs(track->Pt()<0.1)) continue;
bfc7c23b 714
05ac9660 715 int nITS = track->GetNcls(0);
716
bfc7c23b 717 fTrackPtBefTrkCuts->Fill(track->Pt());
f2111233 718
05ac9660 719 /*
f2fa5a00 720 UChar_t itsPixel = track->GetITSClusterMap();
05ac9660 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
bfc7c23b 729 // RecKine: ITSTPC cuts
730 if(!ProcessCutStep(AliHFEcuts::kStepRecKineITSTPC, track)) continue;
731
732 //RecKink
733 if(fRejectKinkMother) { // Quick and dirty fix to reject both kink mothers and daughters
734 if(track->GetKinkIndex(0) != 0) continue;
735 }
736
737 // RecPrim
738 if(!ProcessCutStep(AliHFEcuts::kStepRecPrim, track)) continue;
739
740 // HFEcuts: ITS layers cuts
fbe402d7 741 if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsITS, track)) continue;
d26226a1 742 if(mcPho && iHijing==0)fPhoVertexReco_step0->Fill(track->Pt(),conv_proR); // check MC vertex
743 if(mcPho && iHijing==1)fPhoVertexReco_step1->Fill(track->Pt(),conv_proR); // check MC vertex
bfc7c23b 744
f09766dd 745 // HFE cuts: TPC PID cleanup
746 if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTPC, track)) continue;
747
a6df418c 748 int nTPCcl = track->GetTPCNcls();
96167a04 749 //int nTPCclF = track->GetTPCNclsF(); // warnings
05ac9660 750 //int nITS = track->GetNcls(0);
751
f2fa5a00 752 if(mcPho && iHijing==1)fPhoVertexReco_HFE->Fill(track->Pt(),conv_proR); // check MC vertex
05ac9660 753
bfc7c23b 754 fTrackPtAftTrkCuts->Fill(track->Pt());
755
756 Double_t mom = -999., eop=-999., pt = -999., dEdx=-999., fTPCnSigma=-10, phi=-999., eta=-999.;
feffe705 757 pt = track->Pt();
d26226a1 758 if(pt<2.5)continue;
fbe402d7 759 if(pt<0.1)continue;
bfc7c23b 760
6f4a9952 761 //Int_t charge = track->Charge();
bfc7c23b 762 fTrkpt->Fill(pt);
763 mom = track->P();
764 phi = track->Phi();
765 eta = track->Eta();
ae88597c 766 float dca_xy;
767 float dca_z;
768 track->GetImpactParameters(dca_xy,dca_z);
769
bfc7c23b 770 dEdx = track->GetTPCsignal();
771 fTPCnSigma = fPID->GetPIDResponse() ? fPID->GetPIDResponse()->NumberOfSigmasTPC(track, AliPID::kElectron) : 1000;
8af58b0d 772
09022877 773 //cout << "nSigma correctoon-----" << endl;
774 //cout << "org = " << fTPCnSigma << endl;
9eae9ea5 775 /*
09022877 776 if(!fmcData) // nsigma eta correction
777 {
778 double nSigexpCorr = NsigmaCorrection(eta,cent);
779 fTPCnSigma -= nSigexpCorr;
780 }
9eae9ea5 781 */
09022877 782 //cout << "correction = " << fTPCnSigma << endl;
783
89f41a30 784 //--- track cluster match
785
bfc7c23b 786 double ncells = -1.0;
787 double m20 = -1.0;
788 double m02 = -1.0;
789 double disp = -1.0;
790 double rmatch = -1.0;
791 double nmatch = -1.0;
89f41a30 792 //double oppstatus = 0.0;
59d998de 793 double emctof = 0.0;
6f4a9952 794 Bool_t MaxEmatch = kFALSE;
bfc7c23b 795
bfc7c23b 796 Int_t clsId = track->GetEMCALcluster();
797 if (clsId>0){
798 AliESDCaloCluster *clust = fESD->GetCaloCluster(clsId);
799 if(clust && clust->IsEMCAL()){
800
801 double clustE = clust->E();
6f4a9952 802 if(clustE==maxE)MaxEmatch = kTRUE;
bfc7c23b 803 eop = clustE/fabs(mom);
8efacc22 804
bfc7c23b 805 //double clustT = clust->GetTOF();
806 ncells = clust->GetNCells();
807 m02 = clust->GetM02();
808 m20 = clust->GetM20();
809 disp = clust->GetDispersion();
810 double delphi = clust->GetTrackDx();
811 double deleta = clust->GetTrackDz();
812 rmatch = sqrt(pow(delphi,2)+pow(deleta,2));
813 nmatch = clust->GetNTracksMatched();
59d998de 814 emctof = clust->GetTOF();
815 //cout << "emctof = " << emctof << endl;
89f41a30 816 cout << "eop org = "<< eop << endl;
817 double eoporg = eop;
818 if(fmcData)
9eae9ea5 819 {
820 double mceopcorr = MCEopMeanCorrection(pt,cent);
821 eop += mceopcorr;
822 }
89f41a30 823 cout << "eop corr = " << eop << endl;
9eae9ea5 824
feffe705 825 double valdedx[16];
a6df418c 826 valdedx[0] = pt; valdedx[1] = nITS; valdedx[2] = phi; valdedx[3] = eta; valdedx[4] = fTPCnSigma;
05ac9660 827 //valdedx[5] = eop; valdedx[6] = rmatch; valdedx[7] = ncells, valdedx[8] = nmatch; valdedx[9] = m20; valdedx[10] = mcpT;
828 valdedx[5] = eop; valdedx[6] = rmatch; valdedx[7] = dca_xy, valdedx[8] = dca_z; valdedx[9] = m20; valdedx[10] = mcpT;
89f41a30 829 valdedx[11] = cent; valdedx[12] = dEdx; valdedx[13] = eoporg; valdedx[14] = nTPCcl;
feffe705 830 valdedx[15] = mcele;
89f41a30 831 fEleInfo->Fill(valdedx);
bfc7c23b 832
833 }
834 }
e1f0fb74 835
836 //Get Cal info PID response
9eae9ea5 837 /*
e1f0fb74 838 double eop2;
839 double ss[4];
840 Double_t nSigmaEop = fPID->GetPIDResponse()->NumberOfSigmasEMCAL(track,AliPID::kElectron,eop2,ss);
841 if(fTPCnSigma>-1.5 && fTPCnSigma<3.0 && nITS>2.5 && nTPCcl>100)
842 {
843 double valEop[3];
844 valEop[0] = cent;
845 valEop[1] = pt;
846 valEop[2] = nSigmaEop;
847 fElenSigma->Fill(valEop);
848 }
9eae9ea5 849 */
89f41a30 850 // --- tarck cut & e ID
e1f0fb74 851
a6df418c 852 if(nITS<2.5)continue;
853 if(nTPCcl<100)continue;
854
855 CheckNclust->Fill(nTPCcl);
856 CheckNits->Fill(nITS);
ae88597c 857 CheckDCA->Fill(dca_xy,dca_z);
f2111233 858 // check production vertex of photons
f2111233 859
a6df418c 860
42c75692 861 fdEdxBef->Fill(mom,fTPCnSigma);
bfc7c23b 862 fTPCnsigma->Fill(mom,fTPCnSigma);
f09766dd 863 if(fTPCnSigma >= -1.0 && fTPCnSigma <= 3)fTrkEovPBef->Fill(pt,eop);
bfc7c23b 864
89f41a30 865 //--- track accepted by HFE
bfc7c23b 866
89f41a30 867 Int_t pidpassed = 1;
bfc7c23b 868 AliHFEpidObject hfetrack;
869 hfetrack.SetAnalysisType(AliHFEpidObject::kESDanalysis);
870 hfetrack.SetRecTrack(track);
3db00c72 871 double binct = 10.5;
872 if((0.0< cent) && (cent<5.0)) binct = 0.5;
873 if((5.0< cent) && (cent<10.0)) binct = 1.5;
874 if((10.0< cent) && (cent<20.0)) binct = 2.5;
875 if((20.0< cent) && (cent<30.0)) binct = 3.5;
876 if((30.0< cent) && (cent<40.0)) binct = 4.5;
877 if((40.0< cent) && (cent<50.0)) binct = 5.5;
878 if((50.0< cent) && (cent<60.0)) binct = 6.5;
879 if((60.0< cent) && (cent<70.0)) binct = 7.5;
880 if((70.0< cent) && (cent<80.0)) binct = 8.5;
881 if((80.0< cent) && (cent<90.0)) binct = 9.5;
882 if((90.0< cent) && (cent<100.0)) binct = 10.5;
883
884 hfetrack.SetCentrality((int)binct); //added
bfc7c23b 885 hfetrack.SetPbPb();
886 if(!fPID->IsSelected(&hfetrack, NULL, "", fPIDqa)) pidpassed = 0;
887
f2fa5a00 888 if(pidpassed==0) continue; // nSigma rejection
89f41a30 889
890 //--- photonic ID
891
892 Bool_t fFlagPhotonicElec = kFALSE;
893 Bool_t fFlagConvinatElec = kFALSE;
894
895 if(fSetKFpart)
896 {
897 SelectPhotonicElectron(iTracks,cent,track,fFlagPhotonicElec,fFlagConvinatElec,fTPCnSigma,m20,eop,mcele,mcWeight,iHijing,mcOrgPi0,mcOrgEta);
898 }
899 else
900 {
901 SelectPhotonicElectron2(iTracks,cent,track,fFlagPhotonicElec,fFlagConvinatElec,fTPCnSigma,m20,eop,mcele,mcWeight,iHijing,mcOrgPi0,mcOrgEta);
902 }
903
904 //--- check reco eff. with TPC
905 double phoval[5];
906 phoval[0] = cent;
907 phoval[1] = pt;
908 phoval[2] = fTPCnSigma;
909 phoval[3] = iHijing;
910 phoval[4] = mcMompT;
911
2f8ef315 912 //if((fTPCnSigma >= -5.0 && fTPCnSigma <= 5) && mcele>-1 && mcPho && mcOrgPi0)
913 if((fTPCnSigma >= -5.0 && fTPCnSigma <= 5) && (mcOrgPi0 || mcOrgEta))
89f41a30 914 {
915 if(iHijing==1)mcWeight = 1.0;
2f8ef315 916 if(mcOrgPi0)
917 {
918 fIncpTMCpho_pi0e_TPC->Fill(phoval,mcWeight);
919 if(fFlagPhotonicElec) fPhoElecPtMC_pi0e_TPC->Fill(phoval,mcWeight);
920 if(fFlagConvinatElec) fSameElecPtMC_pi0e_TPC->Fill(phoval,mcWeight);
921 }
922 if(mcOrgEta)
923 {
924 fIncpTMCpho_eta_TPC->Fill(phoval,mcWeight);
925 if(fFlagPhotonicElec) fPhoElecPtMC_eta_TPC->Fill(phoval,mcWeight);
926 if(fFlagConvinatElec) fSameElecPtMC_eta_TPC->Fill(phoval,mcWeight);
927 }
89f41a30 928 }
929
f2fa5a00 930 //--- matchin check
931 double emcphimim = 1.396;
932 double emcphimax = 3.14;
933 if(phi>emcphimim && phi<emcphimax)
934 {
d26226a1 935 if(fFlagPhotonicElec)fMatchV0_0->Fill(pt); // data
936 if(mcele>2.1)fMatchMC_0->Fill(pt,mcWeight); // MC
937
938 if(eop>=0.0) // have a match
f2fa5a00 939 {
940 if(fFlagPhotonicElec)fMatchV0_1->Fill(pt);
d26226a1 941 if(mcele>2.1)fMatchMC_1->Fill(pt,mcWeight);
f2fa5a00 942 }
943 }
944
89f41a30 945 //+++++++ E/p cut ++++++++++++++++
946
947 if(eop<0.9 || eop>1.3)continue;
948
bfc7c23b 949 fTrkEovPAft->Fill(pt,eop);
42c75692 950 fdEdxAft->Fill(mom,fTPCnSigma);
bfc7c23b 951
e7d87aef 952 // Fill real data
fb87d707 953 fIncpT->Fill(cent,pt);
bfc7c23b 954 if(fFlagPhotonicElec) fPhoElecPt->Fill(cent,pt);
f09766dd 955 if(fFlagConvinatElec) fSameElecPt->Fill(cent,pt);
fb87d707 956
957 if(m20>0.0 && m20<0.3)
958 {
59d998de 959 fIncpTM20->Fill(cent,pt);
960 ftimingEle->Fill(pt,emctof);
fb87d707 961 if(fFlagPhotonicElec) fPhoElecPtM20->Fill(cent,pt);
962 if(fFlagConvinatElec) fSameElecPtM20->Fill(cent,pt);
ebef1da9 963 if(!fFlagPhotonicElec) fSemiElecPtM20->Fill(cent,pt);
fb87d707 964 }
feffe705 965
5e38d973 966
e7d87aef 967 // MC
70448e3b 968 // check label for electron candidiates
969
970 int idlabel = 1;
971 if(mcLabel==0)idlabel = 0;
972 fLabelCheck->Fill(pt,idlabel);
973 if(mcLabel==0)fgeoFake->Fill(phi,eta);
974
d19af75d 975 if(mcLabel<0 && m20>0.0 && m20<0.3 && fTPCnSigma>-1 && fTPCnSigma<3)
976 {
977 fFakeTrk0->Fill(cent,pt);
978 }
979
980 if(mcele>-1) // select MC electrons
e7d87aef 981 {
3db00c72 982
983 fIncpTMChfeAll->Fill(cent,pt);
984 if(m20>0.0 && m20<0.3)fIncpTMCM20hfeAll->Fill(cent,pt);
d19af75d 985 if(m20>0.0 && m20<0.3 && fTPCnSigma>-1 && fTPCnSigma<3)fFakeTrk1->Fill(cent,pt);
3db00c72 986
e4b0faf2 987 if(mcBtoE || mcDtoE) // select B->e & D->e
e7d87aef 988 {
989 fIncpTMChfe->Fill(cent,pt);
bd6ee6dd 990 if(m20>0.0 && m20<0.3)
991 {
09022877 992 //cout << "MC label = " << mcLabel << endl;
bd6ee6dd 993 fIncpTMCM20hfe->Fill(cent,pt);
994 fIncpTMCM20hfeCheck->Fill(cent,mcpT);
995 fIncpTMCM20hfeCheck_weight->Fill(phoval);
996 }
e7d87aef 997 }
e4b0faf2 998
999 if(mcPho) // select photonic electrons
e7d87aef 1000 {
44be9e1d 1001
1002 fIncpTMCpho->Fill(phoval);
1003 if(fFlagPhotonicElec) fPhoElecPtMC->Fill(phoval);
1004 if(fFlagConvinatElec) fSameElecPtMC->Fill(phoval);
e7d87aef 1005
1006 if(m20>0.0 && m20<0.3)
1007 {
44be9e1d 1008 fIncpTMCM20pho->Fill(phoval);
1009 if(fFlagPhotonicElec) fPhoElecPtMCM20->Fill(phoval);
1010 if(fFlagConvinatElec) fSameElecPtMCM20->Fill(phoval);
e4b0faf2 1011 // pi0->g->e
1012 if(mcWeight>-1)
1013 {
987053ce 1014 if(iHijing==1)mcWeight = 1.0;
93c189c5 1015 if(mcOrgPi0)
1016 {
697ecf6b 1017 fIncpTMCM20pho_pi0e->Fill(phoval,mcWeight);
1018 if(fFlagPhotonicElec) fPhoElecPtMCM20_pi0e->Fill(phoval,mcWeight);
1019 if(fFlagConvinatElec) fSameElecPtMCM20_pi0e->Fill(phoval,mcWeight);
05ac9660 1020
1021 // check production vertex
1022 fPhoVertexReco_EMCal->Fill(track->Pt(),conv_proR);
1023 if(fFlagPhotonicElec)fPhoVertexReco_Invmass->Fill(track->Pt(),conv_proR);
1024
1025 }
bd50bf24 1026 // --- eta
93c189c5 1027 if(mcOrgEta)
1028 {
697ecf6b 1029 fIncpTMCM20pho_eta->Fill(phoval,mcWeight);
1030 if(fFlagPhotonicElec) fPhoElecPtMCM20_eta->Fill(phoval,mcWeight);
1031 if(fFlagConvinatElec) fSameElecPtMCM20_eta->Fill(phoval,mcWeight);
93c189c5 1032 }
bd50bf24 1033 // check production vertex
e4b0faf2 1034 }
e7d87aef 1035 }
1036 }
1037 }
bfc7c23b 1038 }
1039 PostData(1, fOutputList);
1040}
1041//_________________________________________
1042void AliAnalysisTaskHFECal::UserCreateOutputObjects()
1043{
42c75692 1044 //--- Check MC
1045
e7d87aef 1046 //Bool_t mcData = kFALSE;
42c75692 1047 if(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler())
1048 {
e7d87aef 1049 fmcData = kTRUE;
42c75692 1050 printf("+++++ MC Data available");
1051 }
e7d87aef 1052 if(fmcData)
42c75692 1053 {
1054 printf("++++++++= MC analysis \n");
1055 }
1056 else
1057 {
1058 printf("++++++++ real data analysis \n");
1059 }
1060
85985bb0 1061 printf("+++++++ QA hist %d",fqahist);
1062
fb87d707 1063 //---- Geometry
1064 fGeom = AliEMCALGeometry::GetInstance("EMCAL_COMPLETEV1");
1065
bfc7c23b 1066 //--------Initialize PID
42c75692 1067 //fPID->SetHasMCData(kFALSE);
e7d87aef 1068 fPID->SetHasMCData(fmcData);
bfc7c23b 1069 if(!fPID->GetNumberOfPIDdetectors())
1070 {
1071 fPID->AddDetector("TPC", 0);
89f41a30 1072 //fPID->AddDetector("EMCAL", 1); <--- apply PID selection in task
bfc7c23b 1073 }
1074
3db00c72 1075 Double_t params[4];
78ff1095 1076 const char *cutmodel;
3db00c72 1077 cutmodel = "pol0";
1078 params[0] = -1.0; //sigma min
9b3495ff 1079 double maxnSig = 3.0;
1080 if(fmcData)
1081 {
1082 params[0] = -5.0; //sigma min
1083 maxnSig = 5.0;
1084 }
1085
1086 for(Int_t a=0;a<11;a++)fPID->ConfigureTPCcentralityCut(a,cutmodel,params,maxnSig);
3db00c72 1087
bfc7c23b 1088 fPID->SortDetectors();
1089 fPIDqa = new AliHFEpidQAmanager();
1090 fPIDqa->Initialize(fPID);
a6df418c 1091
1092 //------- fcut --------------
1093 fCuts = new AliHFEcuts();
1094 fCuts->CreateStandardCuts();
1095 //fCuts->SetMinNClustersTPC(100);
1096 fCuts->SetMinNClustersTPC(90);
1097 fCuts->SetMinRatioTPCclusters(0.6);
1098 fCuts->SetTPCmodes(AliHFEextraCuts::kFound, AliHFEextraCuts::kFoundOverFindable);
1099 //fCuts->SetMinNClustersITS(3);
1100 fCuts->SetMinNClustersITS(2);
05ac9660 1101 fCuts->SetProductionVertex(0,50,0,50);
a6df418c 1102 fCuts->SetCutITSpixel(AliHFEextraCuts::kAny);
1103 fCuts->SetCheckITSLayerStatus(kFALSE);
1104 fCuts->SetVertexRange(10.);
1105 fCuts->SetTOFPIDStep(kFALSE);
fbe402d7 1106 //fCuts->SetPtRange(2, 50);
1107 fCuts->SetPtRange(0.1, 50);
a6df418c 1108 fCuts->SetMaxImpactParam(3.,3.);
1109
bfc7c23b 1110 //--------Initialize correction Framework and Cuts
1111 fCFM = new AliCFManager;
1112 const Int_t kNcutSteps = AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kNcutStepsRecTrack + AliHFEcuts::kNcutStepsDETrack;
1113 fCFM->SetNStepParticle(kNcutSteps);
1114 for(Int_t istep = 0; istep < kNcutSteps; istep++)
1115 fCFM->SetParticleCutsList(istep, NULL);
1116
1117 if(!fCuts){
1118 AliWarning("Cuts not available. Default cuts will be used");
1119 fCuts = new AliHFEcuts;
1120 fCuts->CreateStandardCuts();
1121 }
1122 fCuts->Initialize(fCFM);
1123
1124 //---------Output Tlist
1125 fOutputList = new TList();
1126 fOutputList->SetOwner();
1127 fOutputList->Add(fPIDqa->MakeList("PIDQA"));
1128
1129 fNoEvents = new TH1F("fNoEvents","",4,-0.5,3.5) ;
1130 fOutputList->Add(fNoEvents);
1131
02cd82a0 1132 Int_t binsE[5] = {250, 100, 1000, 200, 10};
fb87d707 1133 Double_t xminE[5] = {1.0, -1, 0.0, 0, -0.5};
1134 Double_t xmaxE[5] = {3.5, 1, 100.0, 100, 9.5};
1135 fEMCAccE = new THnSparseD("fEMCAccE","EMC acceptance & E;#eta;#phi;Energy;Centrality;trugCondition;",5,binsE,xminE,xmaxE);
c82f1b9d 1136 //if(fqahist==1)fOutputList->Add(fEMCAccE);
bfc7c23b 1137
f4e0d2d5 1138 hEMCAccE = new TH2F("hEMCAccE","Cluster Energy",200,0,100,100,0,20);
1139 fOutputList->Add(hEMCAccE);
1140
bfc7c23b 1141 fTrkpt = new TH1F("fTrkpt","track pt",100,0,50);
1142 fOutputList->Add(fTrkpt);
1143
1144 fTrackPtBefTrkCuts = new TH1F("fTrackPtBefTrkCuts","track pt before track cuts",100,0,50);
1145 fOutputList->Add(fTrackPtBefTrkCuts);
1146
1147 fTrackPtAftTrkCuts = new TH1F("fTrackPtAftTrkCuts","track pt after track cuts",100,0,50);
1148 fOutputList->Add(fTrackPtAftTrkCuts);
1149
1150 fTPCnsigma = new TH2F("fTPCnsigma", "TPC - n sigma",100,0,50,200,-10,10);
1151 fOutputList->Add(fTPCnsigma);
1152
1153 fTrkEovPBef = new TH2F("fTrkEovPBef","track E/p before HFE pid",100,0,50,100,0,2);
1154 fOutputList->Add(fTrkEovPBef);
1155
1156 fTrkEovPAft = new TH2F("fTrkEovPAft","track E/p after HFE pid",100,0,50,100,0,2);
1157 fOutputList->Add(fTrkEovPAft);
1158
42c75692 1159 fdEdxBef = new TH2F("fdEdxBef","track dEdx vs p before HFE pid",100,0,50,200,-10,10);
bfc7c23b 1160 fOutputList->Add(fdEdxBef);
1161
42c75692 1162 fdEdxAft = new TH2F("fdEdxAft","track dEdx vs p after HFE pid",100,0,50,200,-10,10);
bfc7c23b 1163 fOutputList->Add(fdEdxAft);
1164
02cd82a0 1165 fIncpT = new TH2F("fIncpT","HFE pid electro vs. centrality",200,0,100,100,0,50);
bfc7c23b 1166 fOutputList->Add(fIncpT);
1167
02cd82a0 1168 fIncpTM20 = new TH2F("fIncpTM20","HFE pid electro vs. centrality with M20",200,0,100,100,0,50);
fb87d707 1169 fOutputList->Add(fIncpTM20);
4e01b68c 1170
5e38d973 1171 Int_t nBinspho[9] = { 10, 30, 600, 60, 50, 4, 40, 8, 30};
1172 Double_t minpho[9] = { 0., 0., -0.1, 40, 0, -0.5, 0,-1.5, 0};
1173 Double_t maxpho[9] = {100., 30., 0.5, 100, 1, 3.5, 2, 6.5, 30};
f09766dd 1174
b12bc641 1175 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 1176 if(fqahist==1)fOutputList->Add(fInvmassLS);
bfc7c23b 1177
b12bc641 1178 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 1179 if(fqahist==1)fOutputList->Add(fInvmassULS);
bfc7c23b 1180
8806ce6c 1181 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 1182 if(fqahist==1)fOutputList->Add(fInvmassLSmc);
8806ce6c 1183
1184 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 1185 if(fqahist==1)fOutputList->Add(fInvmassULSmc);
8806ce6c 1186
d89866fd 1187 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 1188 fInvmassLSreco->Sumw2();
1189 fOutputList->Add(fInvmassLSreco);
1190
d89866fd 1191 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 1192 fInvmassULSreco->Sumw2();
1193 fOutputList->Add(fInvmassULSreco);
1194
d89866fd 1195 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 1196 fInvmassLSmc0->Sumw2();
1197 fOutputList->Add(fInvmassLSmc0);
1198
d89866fd 1199 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 1200 fInvmassLSmc1->Sumw2();
1201 fOutputList->Add(fInvmassLSmc1);
1202
d89866fd 1203 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 1204 fInvmassLSmc2->Sumw2();
1205 fOutputList->Add(fInvmassLSmc2);
1206
d89866fd 1207 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 1208 fInvmassLSmc3->Sumw2();
1209 fOutputList->Add(fInvmassLSmc3);
1210
d89866fd 1211 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 1212 fInvmassULSmc0->Sumw2();
1213 fOutputList->Add(fInvmassULSmc0);
1214
d89866fd 1215 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 1216 fInvmassULSmc1->Sumw2();
1217 fOutputList->Add(fInvmassULSmc1);
1218
d89866fd 1219 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 1220 fInvmassULSmc2->Sumw2();
1221 fOutputList->Add(fInvmassULSmc2);
1222
d89866fd 1223 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 1224 fInvmassULSmc3->Sumw2();
93c189c5 1225 fOutputList->Add(fInvmassULSmc3);
1226
bfc7c23b 1227 fOpeningAngleLS = new TH1F("fOpeningAngleLS","Opening angle for LS pairs",100,0,1);
1228 fOutputList->Add(fOpeningAngleLS);
1229
1230 fOpeningAngleULS = new TH1F("fOpeningAngleULS","Opening angle for ULS pairs",100,0,1);
1231 fOutputList->Add(fOpeningAngleULS);
1232
1233 fPhotoElecPt = new TH1F("fPhotoElecPt", "photonic electron pt",100,0,50);
1234 fOutputList->Add(fPhotoElecPt);
1235
02cd82a0 1236 fPhoElecPt = new TH2F("fPhoElecPt", "Pho-inclusive electron pt",200,0,100,100,0,50);
bfc7c23b 1237 fOutputList->Add(fPhoElecPt);
1238
02cd82a0 1239 fPhoElecPtM20 = new TH2F("fPhoElecPtM20", "Pho-inclusive electron pt with M20",200,0,100,100,0,50);
fb87d707 1240 fOutputList->Add(fPhoElecPtM20);
1241
02cd82a0 1242 fSameElecPt = new TH2F("fSameElecPt", "Same-inclusive electron pt",200,0,100,100,0,50);
f09766dd 1243 fOutputList->Add(fSameElecPt);
1244
02cd82a0 1245 fSameElecPtM20 = new TH2F("fSameElecPtM20", "Same-inclusive electron pt with M20",200,0,100,100,0,50);
fb87d707 1246 fOutputList->Add(fSameElecPtM20);
1247
ebef1da9 1248 fSemiElecPtM20 = new TH2F("fSemiElecPtM20", "Semi-inclusive electron pt with M20",200,0,100,100,0,50);
1249 fOutputList->Add(fSemiElecPtM20);
1250
02cd82a0 1251 fCent = new TH1F("fCent","Centrality",200,0,100) ;
bfc7c23b 1252 fOutputList->Add(fCent);
1253
1254 // Make common binning
44be9e1d 1255 const Double_t kMinP = 0.;
d19af75d 1256 const Double_t kMaxP = 20.;
bfc7c23b 1257
05ac9660 1258 //+++ 1st histogram: TPC dEdx with/without EMCAL (p, pT, TPC Signal, phi, eta, Sig, e/p, ,match, cell, M02, M20, Disp, Centrality, select)
1259 // 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)
1260 Int_t nBins[16] = { 100, 7, 60, 20, 90, 100, 25, 60, 60, 100, 40, 10, 250, 100, 100, 8};
1261 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};
1262 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};
1263 //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);
1264 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 1265 //if(fqahist==1)fOutputList->Add(fEleInfo);
bfc7c23b 1266
e1f0fb74 1267 // Make common binning
1268 Int_t nBinsEop[3] = { 10, 50, 100};
1269 Double_t minEop[3] = { 0, 0, -5};
1270 Double_t maxEop[3] = {100, 50, 5};
1271 fElenSigma= new THnSparseD("fElenSigma", "Electron nSigma; cent; pT [GeV/c]; nSigma", 3, nBinsEop, minEop, maxEop);
1272 fOutputList->Add(fElenSigma);
1273
1274
fb87d707 1275 //<--- Trigger info
feffe705 1276 /*
fb87d707 1277 fClsEBftTrigCut = new TH1F("fClsEBftTrigCut","cluster E before trigger selection",1000,0,100);
1278 fOutputList->Add(fClsEBftTrigCut);
1279
1280 fClsEAftTrigCut = new TH1F("fClsEAftTrigCut","cluster E if cls has 0 trigcut channel",1000,0,100);
1281 fOutputList->Add(fClsEAftTrigCut);
1282
1283 fClsEAftTrigCut1 = new TH1F("fClsEAftTrigCut1","cluster E if cls with trig channel",1000,0,100);
1284 fOutputList->Add(fClsEAftTrigCut1);
1285
1286 fClsEAftTrigCut2 = new TH1F("fClsEAftTrigCut2","cluster E if cls with trigcut channel",1000,0,100);
1287 fOutputList->Add(fClsEAftTrigCut2);
1288
1289 fClsEAftTrigCut3 = new TH1F("fClsEAftTrigCut3","cluster E if cls with trigcut channel + nCell>Ecorrect",1000,0,100);
1290 fOutputList->Add(fClsEAftTrigCut3);
1291
1292 fClsEAftTrigCut4 = new TH1F("fClsEAftTrigCut4","cluster E if cls with trigcut channel + nCell>Ecorrect + cls time cut",1000,0,100);
1293 fOutputList->Add(fClsEAftTrigCut4);
1294
1295 fClsETime = new TH2F("fClsETime", "Cls time vs E; E; time;",1000,0,100,1000,-0.0000002,0.0000009);
1296 fOutputList->Add(fClsETime);
1297
1298 fClsETime1 = new TH2F("fClsETime1", "Cls time vs E if cls contains trigger channel; E; time;",1000,0,100,1000,-0.0000002,0.0000009);
1299 fOutputList->Add(fClsETime1);
1300
1301 fTrigTimes = new TH1F("fTrigTimes", "Trigger time; time; N;",25,0,25);
1302 fOutputList->Add(fTrigTimes);
1303
42c75692 1304 fCellCheck = new TH2F("fCellCheck", "Cell vs E; E GeV; Cell ID",10,6,26,12000,0,12000);
1305 fOutputList->Add(fCellCheck);
feffe705 1306 */
e7d87aef 1307 //<---------- MC
1308
02cd82a0 1309 fInputHFEMC = new TH2F("fInputHFEMC","Input MC HFE pid electro vs. centrality",200,0,100,100,0,50);
e7d87aef 1310 fOutputList->Add(fInputHFEMC);
1311
02cd82a0 1312 fInputAlle = new TH2F("fInputAlle","Input MC electro vs. centrality",200,0,100,100,0,50);
feffe705 1313 fOutputList->Add(fInputAlle);
1314
02cd82a0 1315 fIncpTMChfe = new TH2F("fIncpTMChfe","MC HFE pid electro vs. centrality",200,0,100,100,0,50);
e7d87aef 1316 fOutputList->Add(fIncpTMChfe);
1317
02cd82a0 1318 fIncpTMChfeAll = new TH2F("fIncpTMChfeAll","MC Alle pid electro vs. centrality",200,0,100,100,0,50);
3db00c72 1319 fOutputList->Add(fIncpTMChfeAll);
1320
02cd82a0 1321 fIncpTMCM20hfe = new TH2F("fIncpTMCM20hfe","MC HFE pid electro vs. centrality with M20",200,0,100,100,0,50);
e7d87aef 1322 fOutputList->Add(fIncpTMCM20hfe);
1323
02cd82a0 1324 fIncpTMCM20hfeAll = new TH2F("fIncpTMCM20hfeAll","MC Alle pid electro vs. centrality with M20",200,0,100,100,0,50);
3db00c72 1325 fOutputList->Add(fIncpTMCM20hfeAll);
1326
60544aea 1327 fIncpTMCM20hfeCheck = new TH2F("fIncpTMCM20hfeCheck","MC HFE pid electro vs. centrality with M20 Check",200,0,100,100,0,50);
1328 fOutputList->Add(fIncpTMCM20hfeCheck);
44be9e1d 1329
2b4460a6 1330 Int_t nBinspho2[5] = { 200, 100, 7, 3, 700};
987053ce 1331 Double_t minpho2[5] = { 0., 0., -2.5, -0.5, 0.};
2b4460a6 1332 Double_t maxpho2[5] = {100., 50., 4.5, 2.5, 70.};
bd6ee6dd 1333
1334 fInputHFEMC_weight = new THnSparseD("fInputHFEMC_weight", "MC HFE electron pt",5,nBinspho2,minpho2,maxpho2);
1335 fOutputList->Add(fInputHFEMC_weight);
1336
1337 fIncpTMCM20hfeCheck_weight = new THnSparseD("fIncpTMCM20hfeCheck_weight", "HFE electron pt with M20",5,nBinspho2,minpho2,maxpho2);
1338 fOutputList->Add(fIncpTMCM20hfeCheck_weight);
1339
987053ce 1340 fIncpTMCpho = new THnSparseD("fIncpTMCpho","MC Pho pid electro vs. centrality",5,nBinspho2,minpho2,maxpho2);
e7d87aef 1341 fOutputList->Add(fIncpTMCpho);
1342
987053ce 1343 fIncpTMCM20pho = new THnSparseD("fIncpTMCM20pho","MC Pho pid electro vs. centrality with M20",5,nBinspho2,minpho2,maxpho2);
e7d87aef 1344 fOutputList->Add(fIncpTMCM20pho);
1345
987053ce 1346 fPhoElecPtMC = new THnSparseD("fPhoElecPtMC", "MC Pho-inclusive electron pt",5,nBinspho2,minpho2,maxpho2);
e7d87aef 1347 fOutputList->Add(fPhoElecPtMC);
1348
987053ce 1349 fPhoElecPtMCM20 = new THnSparseD("fPhoElecPtMCM20", "MC Pho-inclusive electron pt with M20",5,nBinspho2,minpho2,maxpho2);
feffe705 1350 fOutputList->Add(fPhoElecPtMCM20);
e7d87aef 1351
987053ce 1352 fSameElecPtMC = new THnSparseD("fSameElecPtMC", "MC Same-inclusive electron pt",5,nBinspho2,minpho2,maxpho2);
e7d87aef 1353 fOutputList->Add(fSameElecPtMC);
1354
987053ce 1355 fSameElecPtMCM20 = new THnSparseD("fSameElecPtMCM20", "MC Same-inclusive electron pt with M20",5,nBinspho2,minpho2,maxpho2);
e7d87aef 1356 fOutputList->Add(fSameElecPtMCM20);
1357
987053ce 1358 fIncpTMCM20pho_pi0e = new THnSparseD("fIncpTMCM20pho_pi0e","MC Pho pi0->e pid electro vs. centrality with M20",5,nBinspho2,minpho2,maxpho2);
1359 fIncpTMCM20pho_pi0e->Sumw2();
e4b0faf2 1360 fOutputList->Add(fIncpTMCM20pho_pi0e);
1361
987053ce 1362 fPhoElecPtMCM20_pi0e = new THnSparseD("fPhoElecPtMCM20_pi0e", "MC Pho-inclusive electron pt with M20 pi0->e",5,nBinspho2,minpho2,maxpho2);
1363 fPhoElecPtMCM20_pi0e->Sumw2();
e4b0faf2 1364 fOutputList->Add(fPhoElecPtMCM20_pi0e);
1365
987053ce 1366 fSameElecPtMCM20_pi0e = new THnSparseD("fSameElecPtMCM20_pi0e", "MC Same-inclusive electron pt pi0->e",5,nBinspho2,minpho2,maxpho2);
1367 fSameElecPtMCM20_pi0e->Sumw2();
e4b0faf2 1368 fOutputList->Add(fSameElecPtMCM20_pi0e);
697ecf6b 1369 //
93c189c5 1370 fIncpTMCM20pho_eta = new THnSparseD("fIncpTMCM20pho_eta","MC Pho pi0->e pid electro vs. centrality with M20",5,nBinspho2,minpho2,maxpho2);
1371 fIncpTMCM20pho_eta->Sumw2();
1372 fOutputList->Add(fIncpTMCM20pho_eta);
1373
1374 fPhoElecPtMCM20_eta = new THnSparseD("fPhoElecPtMCM20_eta", "MC Pho-inclusive electron pt with M20 pi0->e",5,nBinspho2,minpho2,maxpho2);
1375 fPhoElecPtMCM20_eta->Sumw2();
1376 fOutputList->Add(fPhoElecPtMCM20_eta);
1377
1378 fSameElecPtMCM20_eta = new THnSparseD("fSameElecPtMCM20_eta", "MC Same-inclusive electron pt pi0->e",5,nBinspho2,minpho2,maxpho2);
1379 fSameElecPtMCM20_eta->Sumw2();
1380 fOutputList->Add(fSameElecPtMCM20_eta);
697ecf6b 1381 // ------------
1382 fIncpTMCpho_pi0e_TPC = new THnSparseD("fIncpTMCpho_pi0e_TPC","MC Pho pi0->e pid electro vs. centrality with M20",5,nBinspho2,minpho2,maxpho2);
1383 fIncpTMCpho_pi0e_TPC->Sumw2();
1384 fOutputList->Add(fIncpTMCpho_pi0e_TPC);
1385
1386 fPhoElecPtMC_pi0e_TPC = new THnSparseD("fPhoElecPtMC_pi0e_TPC", "MC Pho-inclusive electron pt with pi0->e",5,nBinspho2,minpho2,maxpho2);
1387 fPhoElecPtMC_pi0e_TPC->Sumw2();
1388 fOutputList->Add(fPhoElecPtMC_pi0e_TPC);
1389
1390 fSameElecPtMC_pi0e_TPC = new THnSparseD("fSameElecPtMC_pi0e_TPC", "MC Same-inclusive electron pt pi0->e",5,nBinspho2,minpho2,maxpho2);
1391 fSameElecPtMC_pi0e_TPC->Sumw2();
1392 fOutputList->Add(fSameElecPtMC_pi0e_TPC);
2f8ef315 1393
1394 fIncpTMCpho_eta_TPC = new THnSparseD("fIncpTMCpho_eta_TPC","MC Pho pi0->e pid electro vs. centrality with M20",5,nBinspho2,minpho2,maxpho2);
1395 fIncpTMCpho_eta_TPC->Sumw2();
1396 fOutputList->Add(fIncpTMCpho_eta_TPC);
1397
1398 fPhoElecPtMC_eta_TPC = new THnSparseD("fPhoElecPtMC_eta_TPC", "MC Pho-inclusive electron pt with pi0->e",5,nBinspho2,minpho2,maxpho2);
1399 fPhoElecPtMC_eta_TPC->Sumw2();
1400 fOutputList->Add(fPhoElecPtMC_eta_TPC);
1401
1402 fSameElecPtMC_eta_TPC = new THnSparseD("fSameElecPtMC_eta_TPC", "MC Same-inclusive electron pt pi0->e",5,nBinspho2,minpho2,maxpho2);
1403 fSameElecPtMC_eta_TPC->Sumw2();
1404 fOutputList->Add(fSameElecPtMC_eta_TPC);
1405
1406
697ecf6b 1407 //-------------
e4b0faf2 1408
a6df418c 1409 CheckNclust = new TH1D("CheckNclust","cluster check",200,0,200);
1410 fOutputList->Add(CheckNclust);
1411
1412 CheckNits = new TH1D("CheckNits","ITS cluster check",8,-0.5,7.5);
1413 fOutputList->Add(CheckNits);
ae88597c 1414
1415 CheckDCA = new TH2D("CheckDCA","DCA check",200,-5,5,200,-5,5);
1416 fOutputList->Add(CheckDCA);
68d718e7 1417 /*
d113d7cd 1418 Hpi0pTcheck = new TH2D("Hpi0pTcheck","Pi0 pT from Hijing",100,0,50,3,-0.5,2.5);
1419 fOutputList->Add(Hpi0pTcheck);
1420
19325033 1421 HETApTcheck = new TH2D("HETApTcheck","Eta pT from Hijing",100,0,50,3,-0.5,2.5);
1422 fOutputList->Add(HETApTcheck);
68d718e7 1423 */
1424
1425 Int_t nBinspho3[3] = { 200, 100, 3};
1426 Double_t minpho3[3] = { 0., 0., -0.5};
1427 Double_t maxpho3[3] = {100., 50., 2.5};
19325033 1428
68d718e7 1429 Hpi0pTcheck = new THnSparseD("Hpi0pTcheck","Pi0 pT from Hijing",3,nBinspho3,minpho3,maxpho3);
1430 fOutputList->Add(Hpi0pTcheck);
1431
1432 HETApTcheck = new THnSparseD("HETApTcheck","Eta pT from Hijing",3,nBinspho3,minpho3,maxpho3);
1433 fOutputList->Add(HETApTcheck);
1434 //--
e4b0faf2 1435 HphopTcheck = new TH2D("HphopTcheck","Pho pT from Hijing",100,0,50,3,-0.5,2.5);
1436 fOutputList->Add(HphopTcheck);
5e0e45b3 1437 //
1438 HDpTcheck = new TH2D("HDpTcheck","D pT from Hijing",100,0,50,3,-0.5,2.5);
1439 fOutputList->Add(HDpTcheck);
1440
1441 HBpTcheck = new TH2D("HBpTcheck","B pT from Hijing",100,0,50,3,-0.5,2.5);
1442 fOutputList->Add(HBpTcheck);
1443 //
e4b0faf2 1444
93c189c5 1445 fpTCheck = new TH1D("fpTCheck","pT check",500,0,50);
1446 fOutputList->Add(fpTCheck);
1447
50919258 1448 fMomDtoE = new TH2D("fMomDtoE","D->E pT correlations;e p_{T} GeV/c;D p_{T} GeV/c",400,0,40,400,0,40);
1449 fOutputList->Add(fMomDtoE);
d113d7cd 1450
70448e3b 1451 fLabelCheck = new TH2D("fLabelCheck","MC label",50,0,50,5,-1.5,3.5);
1452 fOutputList->Add(fLabelCheck);
1453
1454 fgeoFake = new TH2D("fgeoFake","Label==0 eta and phi",628,0,6.28,200,-1,1);
1455 fOutputList->Add(fgeoFake);
1456
d19af75d 1457 fFakeTrk0 = new TH2D("fFakeTrk0","fake trakcs",10,0,100,20,0,20);
1458 fOutputList->Add(fFakeTrk0);
1459
1460 fFakeTrk1 = new TH2D("fFakeTrk1","true all e a.f. eID",10,0,100,20,0,20);
1461 fOutputList->Add(fFakeTrk1);
1462
59d998de 1463 ftimingEle = new TH2D("ftimingEle","electron TOF",100,0,20,100,1e-7,1e-6);
1464 fOutputList->Add(ftimingEle);
1465
09022877 1466 // eta correction
1467 // note: parameters 01/31new.TPCnSigmaEtaDep
1468 // 70-90 delta_eta = 0.2
1469
1470 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};
1471 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)
1472 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)
1473 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)
1474 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)
1475 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)
1476 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)
1477 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)
1478
1479 fnSigEtaCorr[0] = new TGraphErrors(12,etaval,corr0); // 0-10
1480 fnSigEtaCorr[1] = new TGraphErrors(12,etaval,corr1); // 10-20
1481 fnSigEtaCorr[2] = new TGraphErrors(12,etaval,corr2); // 20-30
1482 fnSigEtaCorr[3] = new TGraphErrors(12,etaval,corr3); // 30-40
1483 fnSigEtaCorr[4] = new TGraphErrors(12,etaval,corr4); // 40-50
1484 fnSigEtaCorr[5] = new TGraphErrors(12,etaval,corr5); // 50-70
1485 fnSigEtaCorr[6] = new TGraphErrors(12,etaval,corr6); // 70-90
1486
6f4a9952 1487 fIncMaxE = new TH2D("fIncMaxE","Inc",10,0,100,10,0,100);
1488 fOutputList->Add(fIncMaxE);
09022877 1489
c82f1b9d 1490 fIncReco = new TH2D("fIncReco","Inc",10,0,100,100,0,500);
5e38d973 1491 fOutputList->Add(fIncReco);
1492
c82f1b9d 1493 fPhoReco = new TH2D("fPhoReco","Pho",10,0,100,100,0,500);
1494 fOutputList->Add(fPhoReco);
1495
1496 fSamReco = new TH2D("fSamReco","Same",10,0,100,100,0,500);
1497 fOutputList->Add(fSamReco);
1498
1499 fIncRecoMaxE = new TH2D("fIncRecoMaxE","Inc",10,0,100,100,0,500);
1500 fOutputList->Add(fIncRecoMaxE);
1501
1502 fPhoRecoMaxE = new TH2D("fPhoRecoMaxE","Pho",10,0,100,100,0,500);
eabc1527 1503 fOutputList->Add(fPhoRecoMaxE);
5e38d973 1504
c82f1b9d 1505 fSamRecoMaxE = new TH2D("fSamRecoMaxE","Same",10,0,100,100,0,500);
eabc1527 1506 fOutputList->Add(fSamRecoMaxE);
5e38d973 1507
05ac9660 1508 fPhoVertexReco_HFE = new TH2D("fPhoVertexReco_HFE","photon production Vertex mass selection",40,0,20,250,0,50);
1509 fOutputList->Add(fPhoVertexReco_HFE);
1510
1511 fPhoVertexReco_EMCal = new TH2D("fPhoVertexReco_EMCal","photon production Vertex mass selection",40,0,20,250,0,50);
1512 fOutputList->Add(fPhoVertexReco_EMCal);
1513
1514 fPhoVertexReco_Invmass = new TH2D("fPhoVertexReco_Invmass","photon production Vertex mass selection",40,0,20,250,0,50);
1515 fOutputList->Add(fPhoVertexReco_Invmass);
1516
1517 fPhoVertexReco_step0= new TH2D("fPhoVertexReco_step0","photon production Vertex mass selection",40,0,20,250,0,50);
d26226a1 1518 fPhoVertexReco_step0->Sumw2();
05ac9660 1519 fOutputList->Add(fPhoVertexReco_step0);
1520
1521 fPhoVertexReco_step1= new TH2D("fPhoVertexReco_step1","photon production Vertex mass selection",40,0,20,250,0,50);
d26226a1 1522 fPhoVertexReco_step1->Sumw2();
05ac9660 1523 fOutputList->Add(fPhoVertexReco_step1);
1524
f2fa5a00 1525 fMatchV0_0 = new TH1D("fMatchV0_0","V0 match",100,0,20);
1526 fOutputList->Add(fMatchV0_0);
1527
1528 fMatchV0_1 = new TH1D("fMatchV0_1","V0 match",100,0,20);
1529 fOutputList->Add(fMatchV0_1);
1530
1531 fMatchMC_0 = new TH1D("fMatchMC_0","MC match",100,0,20);
1532 fOutputList->Add(fMatchMC_0);
1533
1534 fMatchMC_1 = new TH1D("fMatchMC_1","MC match",100,0,20);
1535 fOutputList->Add(fMatchMC_1);
1536
bfc7c23b 1537 PostData(1,fOutputList);
1538}
1539
1540//________________________________________________________________________
1541void AliAnalysisTaskHFECal::Terminate(Option_t *)
1542{
1543 // Info("Terminate");
1544 AliAnalysisTaskSE::Terminate();
1545}
1546
1547//________________________________________________________________________
1548Bool_t AliAnalysisTaskHFECal::ProcessCutStep(Int_t cutStep, AliVParticle *track)
1549{
1550 // Check single track cuts for a given cut step
1551 const Int_t kMCOffset = AliHFEcuts::kNcutStepsMCTrack;
1552 if(!fCFM->CheckParticleCuts(cutStep + kMCOffset, track)) return kFALSE;
1553 return kTRUE;
1554}
1555//_________________________________________
2b4460a6 1556void 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 1557{
1558 //Identify non-heavy flavour electrons using Invariant mass method
1559
fb87d707 1560 fTrackCuts->SetAcceptKinkDaughters(kFALSE);
feffe705 1561 fTrackCuts->SetRequireTPCRefit(kTRUE);
85985bb0 1562 fTrackCuts->SetRequireITSRefit(kTRUE);
feffe705 1563 fTrackCuts->SetEtaRange(-0.9,0.9);
bfc7c23b 1564 fTrackCuts->SetMaxChi2PerClusterTPC(3.5);
feffe705 1565 fTrackCuts->SetMinNClustersTPC(90);
bfc7c23b 1566
c6bca3f2 1567 //const AliESDVertex *pVtx = fESD->GetPrimaryVertex();
6451f693 1568 Double_t bfield = fESD->GetMagneticField();
2f8ef315 1569 Double_t emass = 0.51*0.001; // (0.51 MeV)
6451f693 1570
93c189c5 1571 double ptEle = track->Pt(); //add
1572 if(ibgevent==0 && w > 0.0)
1573 {
1574 fpTCheck->Fill(ptEle,w);
1575 }
1576
bfc7c23b 1577 Bool_t flagPhotonicElec = kFALSE;
f09766dd 1578 Bool_t flagConvinatElec = kFALSE;
bfc7c23b 1579
8806ce6c 1580 int p1 = 0;
1581 if(mce==3)
1582 {
1583 Int_t label = TMath::Abs(track->GetLabel());
1584 TParticle* particle = stack->Particle(label);
1585 p1 = particle->GetFirstMother();
1586 }
1587
5e65985c 1588 //for(Int_t jTracks = itrack+1; jTracks<fESD->GetNumberOfTracks(); jTracks++){
1589 for(Int_t jTracks = 0; jTracks<fESD->GetNumberOfTracks(); jTracks++){
bfc7c23b 1590 AliESDtrack* trackAsso = fESD->GetTrack(jTracks);
1591 if (!trackAsso) {
1592 printf("ERROR: Could not receive track %d\n", jTracks);
1593 continue;
1594 }
5e65985c 1595 if(itrack==jTracks)continue;
d86a815c 1596 int jbgevent = 0;
5e65985c 1597
8806ce6c 1598 int p2 = 0;
1599 if(mce==3)
1600 {
1601 Int_t label2 = TMath::Abs(trackAsso->GetLabel());
1602 TParticle* particle2 = stack->Particle(label2);
d86a815c 1603 Bool_t MChijing_ass = fMC->IsFromBGEvent(label2);
1604 if(MChijing_ass)jbgevent =1;
8806ce6c 1605 if(particle2->GetFirstMother()>-1)
1606 p2 = particle2->GetFirstMother();
1607 }
1608
f09766dd 1609 Double_t dEdxAsso = -999., ptPrim=-999., ptAsso=-999., openingAngle = -999.;
bfc7c23b 1610 Double_t mass=999., width = -999;
1611 Bool_t fFlagLS=kFALSE, fFlagULS=kFALSE;
1612
93c189c5 1613 //ptPrim = track->Pt();
1614 ptPrim = ptEle;
f09766dd 1615
bfc7c23b 1616 dEdxAsso = trackAsso->GetTPCsignal();
1617 ptAsso = trackAsso->Pt();
1618 Int_t chargeAsso = trackAsso->Charge();
1619 Int_t charge = track->Charge();
1620
b12bc641 1621
6451f693 1622 if(ptAsso <fMimpTassCut) continue;
bfc7c23b 1623 if(!fTrackCuts->AcceptTrack(trackAsso)) continue;
6b2cee73 1624 //if(dEdxAsso <65 || dEdxAsso>100) continue;
1625 double fTPCnSigmaAss = fPID->GetPIDResponse() ? fPID->GetPIDResponse()->NumberOfSigmasTPC(trackAsso, AliPID::kElectron) : 1000;
1626 //cout << "fTPCnSigmaAss = " << fTPCnSigmaAss << endl;
1627 //cout << "fTPCnSigmaAss Cut = " << fMimNsigassCut << endl;
1628 if(fTPCnSigmaAss <fMimNsigassCut || fTPCnSigmaAss>5) continue;
1629 //cout << "fTPCnSigmaAss a.f. cut = " << fTPCnSigmaAss << endl;
bfc7c23b 1630
1631 Int_t fPDGe1 = 11; Int_t fPDGe2 = 11;
1632 if(charge>0) fPDGe1 = -11;
1633 if(chargeAsso>0) fPDGe2 = -11;
b12bc641 1634
1635 //printf("chargeAsso = %d\n",chargeAsso);
1636 //printf("charge = %d\n",charge);
bfc7c23b 1637 if(charge == chargeAsso) fFlagLS = kTRUE;
1638 if(charge != chargeAsso) fFlagULS = kTRUE;
1639
b12bc641 1640 //printf("fFlagLS = %d\n",fFlagLS);
1641 //printf("fFlagULS = %d\n",fFlagULS);
cccab38f 1642 //printf("\n");
b12bc641 1643
6451f693 1644 AliKFParticle::SetField(bfield);
bfc7c23b 1645 AliKFParticle ge1(*track, fPDGe1);
1646 AliKFParticle ge2(*trackAsso, fPDGe2);
2f8ef315 1647
1648 if(fSetMassNonlinear)
1649 {
1650 ge1.SetNonlinearMassConstraint(emass);
1651 ge2.SetNonlinearMassConstraint(emass);
1652 }
1653
bfc7c23b 1654 AliKFParticle recg(ge1, ge2);
1655
2f8ef315 1656 /*
d19af75d 1657 // vertex
bfc7c23b 1658 AliKFVertex primV(*pVtx);
1659 primV += recg;
8079a103 1660 primV -= ge1;
1661 primV -= ge2;
bfc7c23b 1662 recg.SetProductionVertex(primV);
2f8ef315 1663 */
1664
d19af75d 1665 // check chi2
1666 if(recg.GetNDF()<1) continue;
a009053a 1667 Double_t chi2recg = recg.GetChi2()/recg.GetNDF();
1668 Double_t chi2cut = 3.0;
f9447b3b 1669
a009053a 1670 // mass.
1671 if(fSetMassConstraint)
1672 {
1673 recg.SetMassConstraint(0,0.0001);
1674 chi2cut = 30.0;
1675 }
f9447b3b 1676 recg.GetMass(mass,width);
1677
3c475a9d 1678 if(fSetMassWidthCut && width>1e10)continue;
1679
a009053a 1680 // angle
bfc7c23b 1681 openingAngle = ge1.GetAngle(ge2);
1682 if(fFlagLS) fOpeningAngleLS->Fill(openingAngle);
1683 if(fFlagULS) fOpeningAngleULS->Fill(openingAngle);
1684
c2801a73 1685 double ishower = 0;
1686 if(shower>0.0 && shower<0.3)ishower = 1;
1687
b12bc641 1688 double phoinfo[9];
f09766dd 1689 phoinfo[0] = cent;
1690 phoinfo[1] = ptPrim;
1691 phoinfo[2] = mass;
6f4a9952 1692 phoinfo[3] = nSig;
1693 //phoinfo[3] = dEdxAsso;
4e01b68c 1694 phoinfo[4] = openingAngle;
c2801a73 1695 phoinfo[5] = ishower;
1696 phoinfo[6] = ep;
1697 phoinfo[7] = mce;
b12bc641 1698 phoinfo[8] = ptAsso;
f09766dd 1699
1700 if(fFlagLS) fInvmassLS->Fill(phoinfo);
1701 if(fFlagULS) fInvmassULS->Fill(phoinfo);
d86a815c 1702 if(fFlagLS && ibgevent==0 && jbgevent==0) fInvmassLSmc->Fill(phoinfo,w);
93c189c5 1703 if(fFlagULS && ibgevent==0 && jbgevent==0)
1704 {
1705 fInvmassULSmc->Fill(phoinfo,w);
1706 }
b12bc641 1707 //printf("fInvmassCut %f\n",fInvmassCut);
1708 //printf("openingAngle %f\n",fOpeningAngleCut);
1709
f9447b3b 1710 // angle cut
ffa14dda 1711 if(openingAngle > fOpeningAngleCut) continue;
a009053a 1712 // chi2 cut
6ec70f8c 1713 //if(TMath::Sqrt(TMath::Abs(chi2recg))>chi2cut) continue;
1714 if(chi2recg>chi2cut) continue;
1715
1716 if(fFlagLS ) fInvmassLSreco->Fill(ptPrim,mass);
1717 if(fFlagULS) fInvmassULSreco->Fill(ptPrim,mass);
1718
d86a815c 1719 // for real data
93c189c5 1720 //printf("mce =%f\n",mce);
d86a815c 1721 if(mce<-0.5) // mce==-1. is real
1722 {
93c189c5 1723 //printf("Real data\n");
d86a815c 1724 if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec){
d86a815c 1725 flagPhotonicElec = kTRUE;
1726 }
1727 if(mass<fInvmassCut && fFlagLS && !flagConvinatElec){
1728 flagConvinatElec = kTRUE;
1729 }
1730 }
1731 // for MC data
1732 else
1733 {
93c189c5 1734 //printf("MC data\n");
1735
1736 if(w>0.0)
1737 {
2b4460a6 1738 //cout << "tagpi0 = " << tagpi0 << " ; tageta = " << tageta << endl;
6ec70f8c 1739 if(fFlagLS && ibgevent==0 && jbgevent==0 && tagpi0) fInvmassLSmc0->Fill(ptPrim,mass);
1740 if(fFlagULS && ibgevent==0 && jbgevent==0 && tagpi0) fInvmassULSmc0->Fill(ptPrim,mass);
1741 if(fFlagLS && ibgevent==0 && jbgevent==0 && tageta) fInvmassLSmc1->Fill(ptPrim,mass);
1742 if(fFlagULS && ibgevent==0 && jbgevent==0 && tageta) fInvmassULSmc1->Fill(ptPrim,mass);
1743 if(fFlagLS && ibgevent==0 && jbgevent==0 && (p1==p2) && tagpi0) fInvmassLSmc2->Fill(ptPrim,mass);
1744 if(fFlagULS && ibgevent==0 && jbgevent==0 && (p1==p2) && tagpi0) fInvmassULSmc2->Fill(ptPrim,mass);
1745 if(fFlagLS && ibgevent==0 && jbgevent==0 && (p1==p2) && tageta) fInvmassLSmc3->Fill(ptPrim,mass);
1746 if(fFlagULS && ibgevent==0 && jbgevent==0 && (p1==p2) && tageta) fInvmassULSmc3->Fill(ptPrim,mass);
93c189c5 1747 }
6ec70f8c 1748
05ac9660 1749 if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec && (ibgevent==jbgevent)){
1750 //if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec && (p1==p2)){ //<--- only MC train (55,56) v5-03-68-AN , 69 & v5-05-70-AN (till 74)
d86a815c 1751 flagPhotonicElec = kTRUE;
1752 }
1753 if(mass<fInvmassCut && fFlagLS && !flagConvinatElec && (ibgevent==jbgevent)){
1754 flagConvinatElec = kTRUE;
1755 }
1756 }
1757
bfc7c23b 1758 }
1759 fFlagPhotonicElec = flagPhotonicElec;
f09766dd 1760 fFlagConvinatElec = flagConvinatElec;
bfc7c23b 1761
1762}
6451f693 1763
1764//_________________________________________
1765void 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)
1766{
1767 //Identify non-heavy flavour electrons using Invariant mass method
1768
1769 fTrackCuts->SetAcceptKinkDaughters(kFALSE);
1770 fTrackCuts->SetRequireTPCRefit(kTRUE);
1771 fTrackCuts->SetRequireITSRefit(kTRUE);
1772 fTrackCuts->SetEtaRange(-0.9,0.9);
1773 //fTrackCuts->SetRequireSigmaToVertex(kTRUE);
1774 fTrackCuts->SetMaxChi2PerClusterTPC(3.5);
1775 fTrackCuts->SetMinNClustersTPC(90);
1776
6451f693 1777 Double_t eMass = TDatabasePDG::Instance()->GetParticle(11)->Mass(); //Electron mass in GeV
1778 Double_t bfield = fESD->GetMagneticField();
1779
1780 double ptEle = track->Pt(); //add
1781 if(ibgevent==0 && w > 0.0)
1782 {
1783 fpTCheck->Fill(ptEle,w);
1784 }
1785
1786 Bool_t flagPhotonicElec = kFALSE;
1787 Bool_t flagConvinatElec = kFALSE;
1788
1789 int p1 = 0;
1790 if(mce==3)
1791 {
1792 Int_t label = TMath::Abs(track->GetLabel());
1793 TParticle* particle = stack->Particle(label);
1794 p1 = particle->GetFirstMother();
1795 }
1796
1797 //for(Int_t jTracks = itrack+1; jTracks<fESD->GetNumberOfTracks(); jTracks++){
1798 for(Int_t jTracks = 0; jTracks<fESD->GetNumberOfTracks(); jTracks++){
1799 AliESDtrack* trackAsso = fESD->GetTrack(jTracks);
1800 if (!trackAsso) {
1801 printf("ERROR: Could not receive track %d\n", jTracks);
1802 continue;
1803 }
1804 if(itrack==jTracks)continue;
1805 int jbgevent = 0;
1806
1807 int p2 = 0;
1808 if(mce==3)
1809 {
1810 Int_t label2 = TMath::Abs(trackAsso->GetLabel());
1811 TParticle* particle2 = stack->Particle(label2);
1812 Bool_t MChijing_ass = fMC->IsFromBGEvent(label2);
1813 if(MChijing_ass)jbgevent =1;
1814 if(particle2->GetFirstMother()>-1)
1815 p2 = particle2->GetFirstMother();
1816 }
1817
1818 Double_t dEdxAsso = -999., ptPrim=-999., ptAsso=-999., openingAngle = -999.;
1819 //Double_t mass=999., width = -999;
1820 Double_t mass=999.;
1821 Bool_t fFlagLS=kFALSE, fFlagULS=kFALSE;
1822
1823 //ptPrim = track->Pt();
1824 ptPrim = ptEle;
1825
1826 dEdxAsso = trackAsso->GetTPCsignal();
1827 ptAsso = trackAsso->Pt();
1828 Int_t chargeAsso = trackAsso->Charge();
1829 Int_t charge = track->Charge();
1830
1831
1832 if(ptAsso <0.5) continue;
1833 if(!fTrackCuts->AcceptTrack(trackAsso)) continue;
1834 if(dEdxAsso <65 || dEdxAsso>100) continue; //11a pass1
1835
1836 Int_t fPDGe1 = 11; Int_t fPDGe2 = 11;
1837 if(charge>0) fPDGe1 = -11;
1838 if(chargeAsso>0) fPDGe2 = -11;
1839
1840 //printf("chargeAsso = %d\n",chargeAsso);
1841 //printf("charge = %d\n",charge);
1842 if(charge == chargeAsso) fFlagLS = kTRUE;
1843 if(charge != chargeAsso) fFlagULS = kTRUE;
1844
1845 // mass cal 0
1846 double xt1 = 0.0, xt2 = 0.0;
1847 Double_t p1at[3] = {0,0,0};
1848 Double_t p2at[3] = {0,0,0};
1849 //double dca12 = trackAsso->GetDCA(track,bfield,xt2,xt1); //DCA track1-track2
1850 double kHasdcaT1 = track->GetPxPyPzAt(xt1,bfield,p1at); //Track1
1851 double kHasdcaT2 = trackAsso->GetPxPyPzAt(xt2,bfield,p2at); //Track2
1852 if(!kHasdcaT1 || !kHasdcaT2) AliWarning("It could be a problem in the extrapolation");
1853 //cout << "dca = " << dca12 << endl;
1854 // 3D
1855 TLorentzVector electron1;
1856 TLorentzVector electron2;
1857 TLorentzVector mother;
1858
1859 electron1.SetXYZM(p1at[0], p1at[1], p1at[2], eMass);
1860 electron2.SetXYZM(p2at[0], p2at[1], p2at[2], eMass);
1861 openingAngle = TVector2::Phi_0_2pi(electron1.Angle(electron2.Vect()));
4ade0749 1862
6451f693 1863 mother = electron1 + electron2;
4ade0749 1864 //double invmassAtDCA = mother.M();
6451f693 1865
1866 // 2D
1867 TLorentzVector electron1_2D;
1868 TLorentzVector electron2_2D;
1869 TLorentzVector mother_2D;
1870
1871 double pT1at = sqrt(pow(p1at[0],2)+pow(p1at[1],2));
1872 double pT2at = sqrt(pow(p2at[0],2)+pow(p2at[1],2));
1873
1874 electron1_2D.SetXYZM(pT1at, 0.0, p1at[2], eMass);
1875 electron2_2D.SetXYZM(pT2at, 0.0, p2at[2], eMass);
1876
1877 mother_2D = electron1_2D + electron2_2D;
1878 double invmassAtDCA_2D = mother_2D.M();
1879 mass = invmassAtDCA_2D;
1880
1881 // angle
1882 if(fFlagLS) fOpeningAngleLS->Fill(openingAngle);
1883 if(fFlagULS) fOpeningAngleULS->Fill(openingAngle);
1884
1885 double ishower = 0;
1886 if(shower>0.0 && shower<0.3)ishower = 1;
1887
1888 double phoinfo[9];
1889 phoinfo[0] = cent;
1890 phoinfo[1] = ptPrim;
1891 phoinfo[2] = mass;
1892 phoinfo[3] = nSig;
1893 //phoinfo[3] = dEdxAsso;
1894 phoinfo[4] = openingAngle;
1895 phoinfo[5] = ishower;
1896 phoinfo[6] = ep;
1897 phoinfo[7] = mce;
1898 phoinfo[8] = ptAsso;
1899
1900 if(fFlagLS) fInvmassLS->Fill(phoinfo);
1901 if(fFlagULS) fInvmassULS->Fill(phoinfo);
1902 if(fFlagLS && ibgevent==0 && jbgevent==0) fInvmassLSmc->Fill(phoinfo,w);
1903 if(fFlagULS && ibgevent==0 && jbgevent==0)
1904 {
1905 fInvmassULSmc->Fill(phoinfo,w);
1906 }
1907 //printf("fInvmassCut %f\n",fInvmassCut);
1908 //printf("openingAngle %f\n",fOpeningAngleCut);
1909
1910 // angle cut
1911 if(openingAngle > fOpeningAngleCut) continue;
1912 // chi2 cut
1913 //if(TMath::Sqrt(TMath::Abs(chi2recg))>chi2cut) continue;
1914 //if(chi2recg>chi2cut) continue;
1915
1916 if(fFlagLS ) fInvmassLSreco->Fill(ptPrim,mass);
1917 if(fFlagULS) fInvmassULSreco->Fill(ptPrim,mass);
1918
1919 // for real data
1920 //printf("mce =%f\n",mce);
1921 if(mce<-0.5) // mce==-1. is real
1922 {
1923 //printf("Real data\n");
1924 if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec){
1925 flagPhotonicElec = kTRUE;
1926 }
1927 if(mass<fInvmassCut && fFlagLS && !flagConvinatElec){
1928 flagConvinatElec = kTRUE;
1929 }
1930 }
1931 // for MC data
1932 else
1933 {
1934 //printf("MC data\n");
1935
1936 if(w>0.0)
1937 {
1938 //cout << "tagpi0 = " << tagpi0 << " ; tageta = " << tageta << endl;
1939 if(fFlagLS && ibgevent==0 && jbgevent==0 && tagpi0) fInvmassLSmc0->Fill(ptPrim,mass);
1940 if(fFlagULS && ibgevent==0 && jbgevent==0 && tagpi0) fInvmassULSmc0->Fill(ptPrim,mass);
1941 if(fFlagLS && ibgevent==0 && jbgevent==0 && tageta) fInvmassLSmc1->Fill(ptPrim,mass);
1942 if(fFlagULS && ibgevent==0 && jbgevent==0 && tageta) fInvmassULSmc1->Fill(ptPrim,mass);
1943 if(fFlagLS && ibgevent==0 && jbgevent==0 && (p1==p2) && tagpi0) fInvmassLSmc2->Fill(ptPrim,mass);
1944 if(fFlagULS && ibgevent==0 && jbgevent==0 && (p1==p2) && tagpi0) fInvmassULSmc2->Fill(ptPrim,mass);
1945 if(fFlagLS && ibgevent==0 && jbgevent==0 && (p1==p2) && tageta) fInvmassLSmc3->Fill(ptPrim,mass);
1946 if(fFlagULS && ibgevent==0 && jbgevent==0 && (p1==p2) && tageta) fInvmassULSmc3->Fill(ptPrim,mass);
1947 }
1948
1949 if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec && (ibgevent==jbgevent)){
1950 //if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec && (p1==p2)){ <--- only MC train (55,56) v5-03-68-AN & 69 for check
1951 flagPhotonicElec = kTRUE;
1952 }
1953 if(mass<fInvmassCut && fFlagLS && !flagConvinatElec && (ibgevent==jbgevent)){
1954 flagConvinatElec = kTRUE;
1955 }
1956 }
1957
1958 }
1959 fFlagPhotonicElec = flagPhotonicElec;
1960 fFlagConvinatElec = flagConvinatElec;
1961
1962}
1963
8806ce6c 1964//-------------------------------------------
46305ed6 1965
1966void AliAnalysisTaskHFECal::FindMother(TParticle* part, int &label, int &pid)
1967{
1968 //int label = 99999;
1969 //int pid = 99999;
1970
1971 if(part->GetFirstMother()>-1)
1972 {
1973 label = part->GetFirstMother();
1974 pid = stack->Particle(label)->GetPdgCode();
1975 }
70448e3b 1976 //cout << "Find Mother : label = " << label << " ; pid" << pid << endl;
46305ed6 1977}
1978
8806ce6c 1979double AliAnalysisTaskHFECal::GetMCweight(double mcPi0pT)
1980{
1981 double weight = 1.0;
1982
1983 if(mcPi0pT>0.0 && mcPi0pT<5.0)
1984 {
1985 weight = 0.323*mcPi0pT/(TMath::Exp(-1.6+0.767*mcPi0pT+0.0285*mcPi0pT*mcPi0pT));
1986 }
1987 else
1988 {
1989 weight = 115.0/(0.718*mcPi0pT*TMath::Power(mcPi0pT,3.65));
1990 }
1991 return weight;
1992}
fb87d707 1993
93c189c5 1994double AliAnalysisTaskHFECal::GetMCweightEta(double mcEtapT)
1995{
1996 double weight = 1.0;
1997
1998 weight = 223.3/TMath::Power((TMath::Exp(-0.17*mcEtapT-0.0322*mcEtapT*mcEtapT)+mcEtapT/1.69),5.65);
1999 return weight;
2000}
2001
2002
fb87d707 2003//_________________________________________
2004void AliAnalysisTaskHFECal::FindTriggerClusters()
2005{
bd6ee6dd 2006 //cout << "finding trigger patch" << endl;
fb87d707 2007 // constants
2008 const int nModuleCols = 2;
2009 const int nModuleRows = 5;
2010 const int nColsFeeModule = 48;
2011 const int nRowsFeeModule = 24;
2012 const int nColsFaltroModule = 24;
2013 const int nRowsFaltroModule = 12;
2014 //const int faltroWidthMax = 20;
2015
2016 // part 1, trigger extraction -------------------------------------
2017 Int_t globCol, globRow;
2018 //Int_t ntimes=0, nTrigChannel=0, nTrigChannelCut=0, trigInCut=0;
2019 Int_t ntimes=0, nTrigChannel=0, nTrigChannelCut=0;
2020
2021 //Int_t trigtimes[faltroWidthMax];
2022 Double_t cellTime[nColsFeeModule*nModuleCols][nRowsFeeModule*nModuleRows];
2023 Double_t cellEnergy[nColsFeeModule*nModuleCols][nRowsFeeModule*nModuleRows];
2024 //Double_t fTrigCutLow = 6;
2025 //Double_t fTrigCutHigh = 10;
2026 Double_t fTimeCutLow = 469e-09;
2027 Double_t fTimeCutHigh = 715e-09;
2028
2029 AliESDCaloTrigger * fCaloTrigger = fESD->GetCaloTrigger( "EMCAL" );
2030
2031 // erase trigger maps
2032 for(Int_t i = 0; i < nColsFaltroModule*nModuleCols; i++ )
2033 {
2034 for(Int_t j = 0; j < nRowsFaltroModule*nModuleRows; j++ )
2035 {
2036 ftriggersCut[i][j] = 0;
2037 ftriggers[i][j] = 0;
2038 ftriggersTime[i][j] = 0;
2039 }
2040 }
2041
2042 Int_t iglobCol=0, iglobRow=0;
2043 // go through triggers
2044 if( fCaloTrigger->GetEntries() > 0 )
2045 {
2046 // needs reset
2047 fCaloTrigger->Reset();
2048 while( fCaloTrigger->Next() )
2049 {
2050 fCaloTrigger->GetPosition( globCol, globRow );
2051 fCaloTrigger->GetNL0Times( ntimes );
2052 /*
2053 // no L0s
2054 if( ntimes < 1 ) continue;
2055 // get precise timings
2056 fCaloTrigger->GetL0Times( trigtimes );
2057 trigInCut = 0;
2058 for(Int_t i = 0; i < ntimes; i++ )
2059 {
2060 // save the first trigger time in channel
2061 if( i == 0 || triggersTime[globCol][globRow] > trigtimes[i] )
2062 triggersTime[globCol][globRow] = trigtimes[i];
2063 //printf("trigger times: %d\n",trigtimes[i]);
2064 // check if in cut
2065 if(trigtimes[i] > fTrigCutLow && trigtimes[i] < fTrigCutHigh )
2066 trigInCut = 1;
2067
2068 fTrigTimes->Fill(trigtimes[i]);
2069 }
2070 */
2071
2072 //L1 analysis from AliAnalysisTaskEMCALTriggerQA
2073 Int_t bit = 0;
2074 fCaloTrigger->GetTriggerBits(bit);
bd6ee6dd 2075 //cout << "bit = " << bit << endl;
fb87d707 2076
2077 Int_t ts = 0;
2078 fCaloTrigger->GetL1TimeSum(ts);
bd6ee6dd 2079 //cout << "ts = " << ts << endl;
fb87d707 2080 if (ts > 0)ftriggers[globCol][globRow] = 1;
2081 // number of triggered channels in event
2082 nTrigChannel++;
2083 // ... inside cut
2084 if(ts>0 && (bit >> 6 & 0x1))
2085 {
2086 iglobCol = globCol;
2087 iglobRow = globRow;
2088 nTrigChannelCut++;
bd6ee6dd 2089 //cout << "ts cut = " << ts << endl;
2090 //cout << "globCol = " << globCol << endl;
2091 //cout << "globRow = " << globRow << endl;
fb87d707 2092 ftriggersCut[globCol][globRow] = 1;
2093 }
2094
2095 } // calo trigger entries
2096 } // has calo trigger entries
2097
2098 // part 2 go through the clusters here -----------------------------------
bd6ee6dd 2099 //cout << " part 2 go through the clusters here ----------------------------------- " << endl;
fb87d707 2100 Int_t nCluster=0, nCell=0, iCell=0, gCell=0;
60d77596 2101 Short_t cellAddr, nSACell;
2102 Int_t mclabel;
fb87d707 2103 //Int_t nSACell, iSACell, mclabel;
2104 Int_t iSACell;
42c75692 2105 Double_t cellAmp=0, cellTimeT=0, clusterTime=0, efrac=0;
fb87d707 2106 Int_t nSupMod, nModule, nIphi, nIeta, iphi, ieta, gphi, geta, feta, fphi;
2107
2108 TRefArray *fCaloClusters = new TRefArray();
2109 fESD->GetEMCALClusters( fCaloClusters );
2110 nCluster = fCaloClusters->GetEntries();
2111
2112
2113 // save all cells times if there are clusters
2114 if( nCluster > 0 ){
2115 // erase time map
2116 for(Int_t i = 0; i < nColsFeeModule*nModuleCols; i++ ){
2117 for(Int_t j = 0; j < nRowsFeeModule*nModuleRows; j++ ){
2118 cellTime[i][j] = 0.;
2119 cellEnergy[i][j] = 0.;
2120 }
2121 }
2122
2123 // get cells
2124 AliESDCaloCells *fCaloCells = fESD->GetEMCALCells();
2125 //AliVCaloCells fCaloCells = fESD->GetEMCALCells();
2126 nSACell = fCaloCells->GetNumberOfCells();
2127 for(iSACell = 0; iSACell < nSACell; iSACell++ ){
2128 // get the cell info *fCal
2129 fCaloCells->GetCell( iSACell, cellAddr, cellAmp, cellTimeT , mclabel, efrac);
2130
2131 // get cell position
2132 fGeom->GetCellIndex( cellAddr, nSupMod, nModule, nIphi, nIeta );
2133 fGeom->GetCellPhiEtaIndexInSModule( nSupMod,nModule, nIphi, nIeta, iphi, ieta);
2134
2135 // convert co global phi eta
2136 gphi = iphi + nRowsFeeModule*(nSupMod/2);
2137 geta = ieta + nColsFeeModule*(nSupMod%2);
2138
2139 // save cell time and energy
2140 cellTime[geta][gphi] = cellTimeT;
2141 cellEnergy[geta][gphi] = cellAmp;
2142
2143 }
2144 }
2145
2146 Int_t nClusterTrig, nClusterTrigCut;
2147 UShort_t *cellAddrs;
2148 Double_t clsE=-999, clsEta=-999, clsPhi=-999;
2149 Float_t clsPos[3] = {0.,0.,0.};
2150
2151 for(Int_t icl=0; icl<fESD->GetNumberOfCaloClusters(); icl++)
2152 {
2153 AliESDCaloCluster *cluster = fESD->GetCaloCluster(icl);
2154 if(!cluster || !cluster->IsEMCAL()) continue;
2155
2156 // get cluster cells
2157 nCell = cluster->GetNCells();
2158
2159 // get cluster energy
2160 clsE = cluster->E();
2161
2162 // get cluster position
2163 cluster->GetPosition(clsPos);
2164 TVector3 clsPosVec(clsPos[0],clsPos[1],clsPos[2]);
2165 clsEta = clsPosVec.Eta();
2166 clsPhi = clsPosVec.Phi();
2167
2168 // get the cell addresses
2169 cellAddrs = cluster->GetCellsAbsId();
2170
2171 // check if the cluster contains cell, that was marked as triggered
2172 nClusterTrig = 0;
2173 nClusterTrigCut = 0;
2174
2175 // loop the cells to check, if cluser in acceptance
2176 // any cluster with a cell outside acceptance is not considered
2177 for( iCell = 0; iCell < nCell; iCell++ )
2178 {
42c75692 2179 // check hot cell
feffe705 2180 //if(clsE>6.0)fCellCheck->Fill(clsE,cellAddrs[iCell]);
42c75692 2181
fb87d707 2182 // get cell position
2183 fGeom->GetCellIndex( cellAddrs[iCell], nSupMod, nModule, nIphi, nIeta );
2184 fGeom->GetCellPhiEtaIndexInSModule( nSupMod,nModule, nIphi, nIeta, iphi, ieta);
2185
2186 // convert co global phi eta
2187 gphi = iphi + nRowsFeeModule*(nSupMod/2);
2188 geta = ieta + nColsFeeModule*(nSupMod%2);
2189
2190 if( cellTime[geta][gphi] > 0. ){
2191 clusterTime += cellTime[geta][gphi];
2192 gCell++;
2193 }
2194
2195 // get corresponding FALTRO
2196 fphi = gphi / 2;
2197 feta = geta / 2;
2198
bd6ee6dd 2199 //cout << "fphi = " << fphi << endl;
2200 //cout << "feta = " << feta << endl;
2201
fb87d707 2202 // try to match with a triggered
2203 if( ftriggers[feta][fphi]==1)
2204 { nClusterTrig++;
2205 }
2206 if( ftriggersCut[feta][fphi]==1)
2207 { nClusterTrigCut++;
2208 }
2209
bd6ee6dd 2210 //cout << "nClusterTrigCut : " << nClusterTrigCut << endl;
2211
fb87d707 2212 } // cells
2213
2214
2215 if( gCell > 0 )
2216 clusterTime = clusterTime / (Double_t)gCell;
2217 // fix the reconstriction code time 100ns jumps
2218 if( fESD->GetBunchCrossNumber() % 4 < 2 )
2219 clusterTime -= 0.0000001;
2220
feffe705 2221 //fClsETime->Fill(clsE,clusterTime);
2222 //fClsEBftTrigCut->Fill(clsE);
fb87d707 2223
2224 if(nClusterTrig>0){
feffe705 2225 //fClsETime1->Fill(clsE,clusterTime);
fb87d707 2226 }
2227
2228 if(nClusterTrig>0){
2229 cluster->SetChi2(1);
feffe705 2230 //fClsEAftTrigCut1->Fill(clsE);
fb87d707 2231 }
2232
2233 if(nClusterTrigCut>0){
2234 cluster->SetChi2(2);
feffe705 2235 //fClsEAftTrigCut2->Fill(clsE);
fb87d707 2236 }
2237
2238 if(nClusterTrigCut>0 && ( nCell > (1 + clsE / 3)))
2239 {
2240 cluster->SetChi2(3);
feffe705 2241 //fClsEAftTrigCut3->Fill(clsE);
fb87d707 2242 }
2243
2244 if(nClusterTrigCut>0 && (nCell > (1 + clsE / 3) )&&( clusterTime > fTimeCutLow && clusterTime < fTimeCutHigh ))
2245 {
2246 // cluster->SetChi2(4);
feffe705 2247 //fClsEAftTrigCut4->Fill(clsE);
fb87d707 2248 }
2249 if(nClusterTrigCut<1)
2250 {
2251 cluster->SetChi2(0);
2252
feffe705 2253 //fClsEAftTrigCut->Fill(clsE);
fb87d707 2254 }
2255
2256 } // clusters
2257}
2258
09022877 2259// <-------- only MC correction
8efacc22 2260double AliAnalysisTaskHFECal::MCEopMeanCorrection(double pTmc, float central)
8af58b0d 2261{
8efacc22 2262 TF1 *fcorr0 = new TF1("fcorr0","[0]*tanh([1]+[2]*x)");
2263 TF1 *fcorr1 = new TF1("fcorr1","[0]*tanh([1]+[2]*x)");
2264
8af58b0d 2265 double shift = 0.0;
8efacc22 2266
2267 if(central>0 && central<=10)
2268 {
2269 fcorr0->SetParameters(1.045,1.288,3.18e-01); //
2270 fcorr1->SetParameters(9.91e-01,3.466,2.344);
2271 }
09022877 2272 else if(central>10 && central<=20)
8efacc22 2273 {
2274 fcorr0->SetParameters(1.029,8.254e-01,4.07e-01);
2275 fcorr1->SetParameters(0.975,2.276,1.501e-01);
2276 }
2277 else if(central>20 && central<=30)
2278 {
2279 fcorr0->SetParameters(1.01,8.795e-01,3.904e-01);
2280 fcorr1->SetParameters(9.675e-01,1.654,2.583e-01);
2281 }
2282 else if(central>30 && central<=40)
2283 {
2284 fcorr0->SetParameters(1.00,1.466,2.305e-1);
2285 fcorr1->SetParameters(9.637e-01,1.473,2.754e-01);
2286 }
2287 else if(central>40 && central<=50)
2288 {
fd2126aa 2289 fcorr0->SetParameters(1.00,1.422,1.518e-01);
8efacc22 2290 fcorr1->SetParameters(9.59e-01,1.421,2.931e-01);
2291 }
2292
2293 else if(central>50 && central<=70)
2294 {
2295 fcorr0->SetParameters(0.989,2.495,2.167);
e1f0fb74 2296 fcorr1->SetParameters(0.961,1.734,1.438e-01);
8efacc22 2297 }
2298 else if(central>70 && central<=100)
2299 {
2300 fcorr0->SetParameters(0.981,-3.373,3.93327);
fd2126aa 2301 fcorr1->SetParameters(9.574e-01,1.698,1.58e-01);
8efacc22 2302 }
2303
2304
2305 shift = fcorr0->Eval(pTmc)-fcorr1->Eval(pTmc);
2306
89f41a30 2307 fcorr0->Delete();
2308 fcorr1->Delete();
2309
8af58b0d 2310 return shift;
2311}
fb87d707 2312
09022877 2313// <-------- only Data correction
2314double AliAnalysisTaskHFECal::NsigmaCorrection(double tmpeta, float central)
2315{
2316 int icent = 0;
2317
2318 if(central>=0 && central<10)
2319 {
2320 icent = 0;
2321 }
2322 else if(central>=10 && central<20)
2323 {
2324 icent = 1;
2325 }
2326 else if(central>=20 && central<30)
2327 {
2328 icent = 2;
2329 }
2330 else if(central>=30 && central<40)
2331 {
2332 icent = 3;
2333 }
2334 else if(central>=40 && central<50)
2335 {
2336 icent = 4;
2337 }
2338 else if(central>=50 && central<70)
2339 {
2340 icent = 5;
2341 }
2342 else
2343 {
2344 icent = 6;
2345 }
2346
2347 double shift = fnSigEtaCorr[icent]->Eval(tmpeta);
2348
2349 //cout << "eta correction"<< endl;
2350 //cout << "cent = "<< central<< endl;
2351 //cout << "icent = "<< icent << endl;
2352 //cout << "shift = "<< shift << endl;
2353
2354 return shift;
fb87d707 2355
09022877 2356}
fb87d707 2357
2358