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