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