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