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