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