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