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