]> 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
mixed events modified
[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"
30
31#include "AliAnalysisTask.h"
32#include "AliAnalysisManager.h"
33
34#include "AliESDEvent.h"
35#include "AliESDHandler.h"
36#include "AliAODEvent.h"
37#include "AliAODHandler.h"
38
e7d87aef 39#include "AliMCEventHandler.h"
40#include "AliMCEvent.h"
41#include "AliMCParticle.h"
42
bfc7c23b 43#include "AliAnalysisTaskHFECal.h"
44#include "TGeoGlobalMagField.h"
45#include "AliLog.h"
46#include "AliAnalysisTaskSE.h"
47#include "TRefArray.h"
48#include "TVector.h"
49#include "AliESDInputHandler.h"
50#include "AliESDpid.h"
51#include "AliESDtrackCuts.h"
52#include "AliPhysicsSelection.h"
53#include "AliESDCaloCluster.h"
54#include "AliAODCaloCluster.h"
55#include "AliEMCALRecoUtils.h"
56#include "AliEMCALGeometry.h"
57#include "AliGeomManager.h"
58#include "stdio.h"
59#include "TGeoManager.h"
60#include "iostream"
61#include "fstream"
62
63#include "AliEMCALTrack.h"
64#include "AliMagF.h"
65
66#include "AliKFParticle.h"
67#include "AliKFVertex.h"
68
69#include "AliPID.h"
70#include "AliPIDResponse.h"
71#include "AliHFEcontainer.h"
72#include "AliHFEcuts.h"
73#include "AliHFEpid.h"
74#include "AliHFEpidBase.h"
75#include "AliHFEpidQAmanager.h"
76#include "AliHFEtools.h"
77#include "AliCFContainer.h"
78#include "AliCFManager.h"
79
e7d87aef 80#include "AliStack.h"
81
bfc7c23b 82#include "AliCentrality.h"
83
84ClassImp(AliAnalysisTaskHFECal)
85//________________________________________________________________________
86AliAnalysisTaskHFECal::AliAnalysisTaskHFECal(const char *name)
87 : AliAnalysisTaskSE(name)
88 ,fESD(0)
e7d87aef 89 ,fMC(0)
46305ed6 90 ,stack(0)
fb87d707 91 ,fGeom(0)
bfc7c23b 92 ,fOutputList(0)
5e0e45b3 93 ,fqahist(1)
bfc7c23b 94 ,fTrackCuts(0)
95 ,fCuts(0)
96 ,fIdentifiedAsOutInz(kFALSE)
97 ,fPassTheEventCut(kFALSE)
98 ,fRejectKinkMother(kFALSE)
e7d87aef 99 ,fmcData(kFALSE)
bfc7c23b 100 ,fVz(0.0)
101 ,fCFM(0)
102 ,fPID(0)
103 ,fPIDqa(0)
104 ,fOpeningAngleCut(0.1)
105 ,fInvmassCut(0.01)
106 ,fNoEvents(0)
107 ,fEMCAccE(0)
f4e0d2d5 108 ,hEMCAccE(0)
bfc7c23b 109 ,fTrkpt(0)
110 ,fTrkEovPBef(0)
111 ,fTrkEovPAft(0)
112 ,fdEdxBef(0)
113 ,fdEdxAft(0)
114 ,fIncpT(0)
fb87d707 115 ,fIncpTM20(0)
bfc7c23b 116 ,fInvmassLS(0)
117 ,fInvmassULS(0)
8806ce6c 118 ,fInvmassLSmc(0)
119 ,fInvmassULSmc(0)
93c189c5 120 ,fInvmassLSmc0(0)
121 ,fInvmassLSmc1(0)
122 ,fInvmassLSmc2(0)
123 ,fInvmassLSmc3(0)
124 ,fInvmassULSmc0(0)
125 ,fInvmassULSmc1(0)
126 ,fInvmassULSmc2(0)
127 ,fInvmassULSmc3(0)
bfc7c23b 128 ,fOpeningAngleLS(0)
129 ,fOpeningAngleULS(0)
130 ,fPhotoElecPt(0)
131 ,fPhoElecPt(0)
fb87d707 132 ,fPhoElecPtM20(0)
f09766dd 133 ,fSameElecPt(0)
fb87d707 134 ,fSameElecPtM20(0)
bfc7c23b 135 ,fTrackPtBefTrkCuts(0)
136 ,fTrackPtAftTrkCuts(0)
137 ,fTPCnsigma(0)
138 ,fCent(0)
139 ,fEleInfo(0)
feffe705 140 /*
fb87d707 141 ,fClsEBftTrigCut(0)
142 ,fClsEAftTrigCut(0)
143 ,fClsEAftTrigCut1(0)
144 ,fClsEAftTrigCut2(0)
145 ,fClsEAftTrigCut3(0)
146 ,fClsEAftTrigCut4(0)
147 ,fClsETime(0)
148 ,fClsETime1(0)
149 ,fTrigTimes(0)
42c75692 150 ,fCellCheck(0)
feffe705 151 */
e7d87aef 152 ,fInputHFEMC(0)
feffe705 153 ,fInputAlle(0)
e7d87aef 154 ,fIncpTMChfe(0)
3db00c72 155 ,fIncpTMChfeAll(0)
e7d87aef 156 ,fIncpTMCM20hfe(0)
3db00c72 157 ,fIncpTMCM20hfeAll(0)
60544aea 158 ,fIncpTMCM20hfeCheck(0)
bd6ee6dd 159 ,fInputHFEMC_weight(0)
160 ,fIncpTMCM20hfeCheck_weight(0)
e7d87aef 161 ,fIncpTMCpho(0)
162 ,fIncpTMCM20pho(0)
163 ,fPhoElecPtMC(0)
164 ,fPhoElecPtMCM20(0)
165 ,fSameElecPtMC(0)
166 ,fSameElecPtMCM20(0)
e4b0faf2 167 ,fIncpTMCM20pho_pi0e(0)
168 ,fPhoElecPtMCM20_pi0e(0)
169 ,fSameElecPtMCM20_pi0e(0)
93c189c5 170 ,fIncpTMCM20pho_eta(0)
171 ,fPhoElecPtMCM20_eta(0)
172 ,fSameElecPtMCM20_eta(0)
697ecf6b 173 ,fIncpTMCpho_pi0e_TPC(0)
174 ,fPhoElecPtMC_pi0e_TPC(0)
175 ,fSameElecPtMC_pi0e_TPC(0)
a6df418c 176 ,CheckNclust(0)
177 ,CheckNits(0)
d113d7cd 178 ,Hpi0pTcheck(0)
19325033 179 ,HETApTcheck(0)
bd6ee6dd 180 ,HphopTcheck(0)
5e0e45b3 181 ,HDpTcheck(0)
182 ,HBpTcheck(0)
93c189c5 183 ,fpTCheck(0)
50919258 184 ,fMomDtoE(0)
70448e3b 185 ,fLabelCheck(0)
186 ,fgeoFake(0)
59d998de 187 ,ftimingEle(0)
bfc7c23b 188{
189 //Named constructor
190
191 fPID = new AliHFEpid("hfePid");
192 fTrackCuts = new AliESDtrackCuts();
193
194 // Define input and output slots here
195 // Input slot #0 works with a TChain
196 DefineInput(0, TChain::Class());
197 // Output slot #0 id reserved by the base class for AOD
198 // Output slot #1 writes into a TH1 container
199 // DefineOutput(1, TH1I::Class());
200 DefineOutput(1, TList::Class());
201 // DefineOutput(3, TTree::Class());
202}
203
204//________________________________________________________________________
205AliAnalysisTaskHFECal::AliAnalysisTaskHFECal()
206 : AliAnalysisTaskSE("DefaultAnalysis_AliAnalysisTaskHFECal")
207 ,fESD(0)
e7d87aef 208 ,fMC(0)
46305ed6 209 ,stack(0)
fb87d707 210 ,fGeom(0)
bfc7c23b 211 ,fOutputList(0)
5e0e45b3 212 ,fqahist(1)
bfc7c23b 213 ,fTrackCuts(0)
214 ,fCuts(0)
215 ,fIdentifiedAsOutInz(kFALSE)
216 ,fPassTheEventCut(kFALSE)
217 ,fRejectKinkMother(kFALSE)
e7d87aef 218 ,fmcData(kFALSE)
bfc7c23b 219 ,fVz(0.0)
220 ,fCFM(0)
221 ,fPID(0)
222 ,fPIDqa(0)
223 ,fOpeningAngleCut(0.1)
224 ,fInvmassCut(0.01)
225 ,fNoEvents(0)
226 ,fEMCAccE(0)
f4e0d2d5 227 ,hEMCAccE(0)
bfc7c23b 228 ,fTrkpt(0)
229 ,fTrkEovPBef(0)
230 ,fTrkEovPAft(0)
231 ,fdEdxBef(0)
232 ,fdEdxAft(0)
233 ,fIncpT(0)
fb87d707 234 ,fIncpTM20(0)
bfc7c23b 235 ,fInvmassLS(0)
236 ,fInvmassULS(0)
8806ce6c 237 ,fInvmassLSmc(0)
238 ,fInvmassULSmc(0)
93c189c5 239 ,fInvmassLSmc0(0)
240 ,fInvmassLSmc1(0)
241 ,fInvmassLSmc2(0)
242 ,fInvmassLSmc3(0)
243 ,fInvmassULSmc0(0)
244 ,fInvmassULSmc1(0)
245 ,fInvmassULSmc2(0)
246 ,fInvmassULSmc3(0)
bfc7c23b 247 ,fOpeningAngleLS(0)
248 ,fOpeningAngleULS(0)
249 ,fPhotoElecPt(0)
250 ,fPhoElecPt(0)
fb87d707 251 ,fPhoElecPtM20(0)
f09766dd 252 ,fSameElecPt(0)
fb87d707 253 ,fSameElecPtM20(0)
bfc7c23b 254 ,fTrackPtBefTrkCuts(0)
255 ,fTrackPtAftTrkCuts(0)
256 ,fTPCnsigma(0)
257 ,fCent(0)
258 ,fEleInfo(0)
feffe705 259 /*
fb87d707 260 ,fClsEBftTrigCut(0)
261 ,fClsEAftTrigCut(0)
262 ,fClsEAftTrigCut1(0)
263 ,fClsEAftTrigCut2(0)
264 ,fClsEAftTrigCut3(0)
265 ,fClsEAftTrigCut4(0)
266 ,fClsETime(0)
267 ,fClsETime1(0)
268 ,fTrigTimes(0)
42c75692 269 ,fCellCheck(0)
feffe705 270 */
e7d87aef 271 ,fInputHFEMC(0)
feffe705 272 ,fInputAlle(0)
e7d87aef 273 ,fIncpTMChfe(0)
3db00c72 274 ,fIncpTMChfeAll(0)
e7d87aef 275 ,fIncpTMCM20hfe(0)
3db00c72 276 ,fIncpTMCM20hfeAll(0)
60544aea 277 ,fIncpTMCM20hfeCheck(0)
bd6ee6dd 278 ,fInputHFEMC_weight(0)
279 ,fIncpTMCM20hfeCheck_weight(0)
e7d87aef 280 ,fIncpTMCpho(0)
281 ,fIncpTMCM20pho(0)
282 ,fPhoElecPtMC(0)
283 ,fPhoElecPtMCM20(0)
284 ,fSameElecPtMC(0)
285 ,fSameElecPtMCM20(0)
e4b0faf2 286 ,fIncpTMCM20pho_pi0e(0)
287 ,fPhoElecPtMCM20_pi0e(0)
288 ,fSameElecPtMCM20_pi0e(0)
93c189c5 289 ,fIncpTMCM20pho_eta(0)
290 ,fPhoElecPtMCM20_eta(0)
291 ,fSameElecPtMCM20_eta(0)
697ecf6b 292 ,fIncpTMCpho_pi0e_TPC(0)
293 ,fPhoElecPtMC_pi0e_TPC(0)
294 ,fSameElecPtMC_pi0e_TPC(0)
a6df418c 295 ,CheckNclust(0)
296 ,CheckNits(0)
d113d7cd 297 ,Hpi0pTcheck(0)
19325033 298 ,HETApTcheck(0)
e4b0faf2 299 ,HphopTcheck(0)
5e0e45b3 300 ,HDpTcheck(0)
301 ,HBpTcheck(0)
93c189c5 302 ,fpTCheck(0)
50919258 303 ,fMomDtoE(0)
70448e3b 304 ,fLabelCheck(0)
305 ,fgeoFake(0)
59d998de 306 ,ftimingEle(0)
bfc7c23b 307{
308 //Default constructor
309 fPID = new AliHFEpid("hfePid");
310
311 fTrackCuts = new AliESDtrackCuts();
312
313 // Constructor
314 // Define input and output slots here
315 // Input slot #0 works with a TChain
316 DefineInput(0, TChain::Class());
317 // Output slot #0 id reserved by the base class for AOD
318 // Output slot #1 writes into a TH1 container
319 // DefineOutput(1, TH1I::Class());
320 DefineOutput(1, TList::Class());
321 //DefineOutput(3, TTree::Class());
322}
323//_________________________________________
324
325AliAnalysisTaskHFECal::~AliAnalysisTaskHFECal()
326{
327 //Destructor
328
329 delete fOutputList;
fb87d707 330 delete fGeom;
bfc7c23b 331 delete fPID;
a6df418c 332 delete fCuts;
bfc7c23b 333 delete fCFM;
334 delete fPIDqa;
335 delete fTrackCuts;
336}
337//_________________________________________
338
339void AliAnalysisTaskHFECal::UserExec(Option_t*)
340{
341 //Main loop
342 //Called for each event
343
344 // create pointer to event
345 fESD = dynamic_cast<AliESDEvent*>(InputEvent());
346 if (!fESD) {
347 printf("ERROR: fESD not available\n");
348 return;
349 }
350
351 if(!fCuts){
352 AliError("HFE cuts not available");
353 return;
354 }
355
356 if(!fPID->IsInitialized()){
357 // Initialize PID with the given run number
358 AliWarning("PID not initialised, get from Run no");
359 fPID->InitializePID(fESD->GetRunNumber());
360 }
361
e7d87aef 362 if(fmcData)fMC = MCEvent();
46305ed6 363 //AliStack* stack = NULL;
e7d87aef 364 if(fmcData && fMC)stack = fMC->Stack();
365
366 Float_t cent = -1.;
367 AliCentrality *centrality = fESD->GetCentrality();
368 cent = centrality->GetCentralityPercentile("V0M");
369
370 //---- fill MC track info
371 if(fmcData && fMC)
372 {
373 Int_t nParticles = stack->GetNtrack();
374 for (Int_t iParticle = 0; iParticle < nParticles; iParticle++) {
375 TParticle* particle = stack->Particle(iParticle);
376 int fPDG = particle->GetPdgCode();
377 double mcZvertex = fMC->GetPrimaryVertex()->GetZ();
378 double pTMC = particle->Pt();
feffe705 379 double proR = particle->R();
9b3495ff 380 double etaMC = particle->Eta();
381 if(fabs(etaMC)>0.6)continue;
e7d87aef 382 Bool_t mcInDtoE= kFALSE;
383 Bool_t mcInBtoE= kFALSE;
384
d113d7cd 385 Bool_t MChijing = fMC->IsFromBGEvent(iParticle);
fd426741 386 //if(!MChijing)printf("not MC hijing");
d113d7cd 387 int iHijing = 1;
388 if(!MChijing)iHijing = 0;
68d718e7 389 double mcphoinfo[5];
390 mcphoinfo[0] = cent;
391 mcphoinfo[1] = pTMC;
392 mcphoinfo[2] = iHijing;
393 //if(fPDG==111)Hpi0pTcheck->Fill(pTMC,iHijing);
394 //if(fPDG==221)HETApTcheck->Fill(pTMC,iHijing);
395 if(fPDG==111)Hpi0pTcheck->Fill(mcphoinfo);
396 if(fPDG==221)HETApTcheck->Fill(mcphoinfo);
bd6ee6dd 397 if(fabs(fPDG)==411 || fabs(fPDG)==413 || fabs(fPDG)==421 || fabs(fPDG)==423 || fabs(fPDG)==431)HDpTcheck->Fill(pTMC,iHijing);
398 if(fabs(fPDG)==511 || fabs(fPDG)==513 || fabs(fPDG)==521 || fabs(fPDG)==523 || fabs(fPDG)==531)HBpTcheck->Fill(pTMC,iHijing);
d113d7cd 399
e4b0faf2 400 if(particle->GetFirstMother()>-1 && fPDG==22)
401 {
402 int parentPID = stack->Particle(particle->GetFirstMother())->GetPdgCode();
403 if(parentPID==111 || parentPID==221)HphopTcheck->Fill(pTMC,iHijing); // pi0->g & eta->g
404 }
405
e7d87aef 406 if(particle->GetFirstMother()>-1 && fabs(fPDG)==11)
407 {
408 int parentPID = stack->Particle(particle->GetFirstMother())->GetPdgCode();
bd6ee6dd 409 double pTMCparent = stack->Particle(particle->GetFirstMother())->Pt();
e7d87aef 410 if((fabs(parentPID)==411 || fabs(parentPID)==413 || fabs(parentPID)==421 || fabs(parentPID)==423 || fabs(parentPID)==431)&& fabs(fPDG)==11)mcInDtoE = kTRUE;
411 if((fabs(parentPID)==511 || fabs(parentPID)==513 || fabs(parentPID)==521 || fabs(parentPID)==523 || fabs(parentPID)==531)&& fabs(fPDG)==11)mcInBtoE = kTRUE;
bd6ee6dd 412 if((mcInBtoE || mcInDtoE) && fabs(mcZvertex)<10.0)
413 {
414 fInputHFEMC->Fill(cent,pTMC);
415 double mcinfo[5];
416 mcinfo[0] = cent;
417 mcinfo[1] = pTMC;
418 mcinfo[2] = 0.0;
419 mcinfo[3] = iHijing;
420 mcinfo[4] = pTMCparent;
421 fInputHFEMC_weight->Fill(mcinfo);
422 }
e7d87aef 423 }
424
5e0e45b3 425
feffe705 426 if(proR<7.0 && fabs(fPDG)==11)fInputAlle->Fill(cent,pTMC);
427
e7d87aef 428 }
429 }
430
bfc7c23b 431 fNoEvents->Fill(0);
432
433 Int_t fNOtrks = fESD->GetNumberOfTracks();
434 const AliESDVertex *pVtx = fESD->GetPrimaryVertex();
435
436 Double_t pVtxZ = -999;
437 pVtxZ = pVtx->GetZ();
438
439 if(TMath::Abs(pVtxZ)>10) return;
440 fNoEvents->Fill(1);
441
f4e0d2d5 442 if(fNOtrks<1) return;
bfc7c23b 443 fNoEvents->Fill(2);
444
445 AliPIDResponse *pidResponse = fInputHandler->GetPIDResponse();
446 if(!pidResponse){
447 AliDebug(1, "Using default PID Response");
448 pidResponse = AliHFEtools::GetDefaultPID(kFALSE, fInputEvent->IsA() == AliAODEvent::Class());
449 }
450
451 fPID->SetPIDResponse(pidResponse);
452
453 fCFM->SetRecEventInfo(fESD);
454
e7d87aef 455 //Float_t cent = -1.;
456 //AliCentrality *centrality = fESD->GetCentrality();
457 //cent = centrality->GetCentralityPercentile("V0M");
bfc7c23b 458 fCent->Fill(cent);
459
b12bc641 460 //if(cent>90.) return;
bfc7c23b 461
462 // Calorimeter info.
463
42c75692 464 FindTriggerClusters();
465
bfc7c23b 466 // make EMCAL array
467 for(Int_t iCluster=0; iCluster<fESD->GetNumberOfCaloClusters(); iCluster++)
468 {
469 AliESDCaloCluster *clust = fESD->GetCaloCluster(iCluster);
68d718e7 470 if(clust && clust->IsEMCAL())
bfc7c23b 471 {
472 double clustE = clust->E();
473 float emcx[3]; // cluster pos
474 clust->GetPosition(emcx);
475 TVector3 clustpos(emcx[0],emcx[1],emcx[2]);
476 double emcphi = clustpos.Phi();
477 double emceta = clustpos.Eta();
fb87d707 478 double calInfo[5];
479 calInfo[0] = emcphi; calInfo[1] = emceta; calInfo[2] = clustE; calInfo[3] = cent; calInfo[4] = clust->Chi2();
f4e0d2d5 480 //fEMCAccE->Fill(calInfo);
481 //if(clustE>3.0)fEMCAccE->Fill(calInfo);
44be9e1d 482 //if(fqahist==1 && clustE>1.5)fEMCAccE->Fill(calInfo);
f4e0d2d5 483 hEMCAccE->Fill(cent,clustE);
bfc7c23b 484 }
485 }
486
487 // Track loop
488 for (Int_t iTracks = 0; iTracks < fESD->GetNumberOfTracks(); iTracks++) {
489 AliESDtrack* track = fESD->GetTrack(iTracks);
490 if (!track) {
491 printf("ERROR: Could not receive track %d\n", iTracks);
492 continue;
493 }
e7d87aef 494
46305ed6 495 int parentlabel = 99999;
496 int parentPID = 99999;
497 int grand_parentlabel = 99999;
498 int grand_parentPID = 99999;
e7d87aef 499 Bool_t mcPho = kFALSE;
500 Bool_t mcDtoE= kFALSE;
501 Bool_t mcBtoE= kFALSE;
93c189c5 502 Bool_t mcOrgPi0 = kFALSE;
503 Bool_t mcOrgEta = kFALSE;
e7d87aef 504 double mcele = -1.;
44be9e1d 505 double mcpT = 0.0;
b12bc641 506 double mcMompT = 0.0;
e4b0faf2 507 //double mcGrandMompT = 0.0;
508 double mcWeight = -10.0;
d113d7cd 509
510 int iHijing = 1;
70448e3b 511 int mcLabel = -1;
d113d7cd 512
e7d87aef 513 if(fmcData && fMC && stack)
514 {
feffe705 515 Int_t label = TMath::Abs(track->GetLabel());
70448e3b 516 mcLabel = track->GetLabel();
bd6ee6dd 517
8efacc22 518 if(mcLabel>-1)
bd6ee6dd 519 {
520
521 Bool_t MChijing = fMC->IsFromBGEvent(label);
522 if(!MChijing)iHijing = 0;
523
524 TParticle* particle = stack->Particle(label);
525 int mcpid = particle->GetPdgCode();
526 mcpT = particle->Pt();
527 //printf("MCpid = %d",mcpid);
528 if(particle->GetFirstMother()>-1)
529 {
530 //int parentlabel = particle->GetFirstMother();
531 //int parentPID = stack->Particle(particle->GetFirstMother())->GetPdgCode();
532 //mcMompT = stack->Particle(particle->GetFirstMother())->Pt();
533 FindMother(particle, parentlabel, parentPID);
534 mcMompT = stack->Particle(parentlabel)->Pt();
535 if((parentPID==22 || parentPID==111 || parentPID==221)&& fabs(mcpid)==11)mcPho = kTRUE;
536 if((fabs(parentPID)==411 || fabs(parentPID)==413 || fabs(parentPID)==421 || fabs(parentPID)==423 || fabs(parentPID)==431)&& fabs(mcpid)==11)mcDtoE = kTRUE;
537 if((fabs(parentPID)==511 || fabs(parentPID)==513 || fabs(parentPID)==521 || fabs(parentPID)==523 || fabs(parentPID)==531)&& fabs(mcpid)==11)mcBtoE = kTRUE;
538
539 // make D->e pT correlation
540 if(mcDtoE)fMomDtoE->Fill(mcpT,mcMompT);
541
542 //cout << "check PID = " << parentPID << endl;
543 //cout << "check pho = " << mcPho << endl;
544 //cout << "check D or B = " << mcDtoE << endl;
545 // pi->e (Dalitz)
546 if(parentPID==111 && fabs(mcpid)==11 && mcMompT>0.0)
547 {
548 //cout << "find pi0->e " << endl;
549 mcOrgPi0 = kTRUE;
550 mcWeight = GetMCweight(mcMompT);
551 }
552 // eta->e (Dalitz)
553 if(parentPID==221 && fabs(mcpid)==11 && mcMompT>0.0)
554 {
555 //cout << "find Eta->e " << endl;
556 mcOrgEta = kTRUE;
557 mcWeight = GetMCweightEta(mcMompT);
558 }
559
560 // access grand parent
561 TParticle* particle_parent = stack->Particle(parentlabel); // get parent pointer
562 //if(particle_parent->GetFirstMother()>-1 && parentPID==22 && fabs(mcpid)==11) // get grand parent g->e
563 if(particle_parent->GetFirstMother()>-1 && (parentPID==22 || parentPID==111) && fabs(mcpid)==11) // get grand parent g->e
564 {
565 //int grand_parentPID = stack->Particle(particle_parent->GetFirstMother())->GetPdgCode();
566 //double pTtmp = stack->Particle(particle_parent->GetFirstMother())->Pt();
567 FindMother(particle_parent, grand_parentlabel, grand_parentPID);
568 double mcGrandpT = stack->Particle(grand_parentlabel)->Pt();
569 if(grand_parentPID==111 && mcGrandpT>0.0)
570 {
571 // check eta->pi0 decay !
572 int grand2_parentlabel = 99999; int grand2_parentPID = 99999;
573 TParticle* particle_grand = stack->Particle(grand_parentlabel); // get parent pointer
574 FindMother(particle_grand, grand2_parentlabel, grand2_parentPID);
575 if(grand2_parentPID==221)
576 {
577 //cout << "find Eta->e " << endl;
578 double mcGrandpT2 = stack->Particle(grand2_parentlabel)->Pt();
579 mcOrgEta = kTRUE;
580 mcWeight = GetMCweight(mcGrandpT2);
581 mcMompT = mcGrandpT2;
582 }
583 else
584 {
585 //cout << "find pi0->e " << endl;
586 mcOrgPi0 = kTRUE;
587 mcWeight = GetMCweight(mcGrandpT);
588 mcMompT = mcGrandpT;
589 }
590 }
591
592 if(grand_parentPID==221 && mcGrandpT>0.0)
593 {
594 //cout << "find Eta->e " << endl;
595 mcOrgEta = kTRUE;
596 mcOrgPi0 = kFALSE;
597 mcWeight = GetMCweightEta(mcGrandpT);
598 mcMompT = mcGrandpT;
599 }
600 }
601 }
602
cc1605fe 603 //cout << "===================="<<endl;
604 //cout << "mcDtoE : " << mcDtoE << endl;
605 //cout << "mcBtoE : " << mcBtoE << endl;
606 //cout << "mcPho : " << mcPho << endl;
607
608 if(fabs(mcpid)==11)mcele= 0.;
609 //cout << "check e: " << mcele << endl;
bd6ee6dd 610 if(fabs(mcpid)==11 && mcDtoE)mcele= 1.;
cc1605fe 611 //cout << "check D->e: " << mcele << endl;
bd6ee6dd 612 if(fabs(mcpid)==11 && mcBtoE)mcele= 2.;
cc1605fe 613 //cout << "check B->e: " << mcele << endl;
bd6ee6dd 614 if(fabs(mcpid)==11 && mcPho)mcele= 3.;
cc1605fe 615 //cout << "check Pho->e: " << mcele << endl;
bd6ee6dd 616
fd2126aa 617 cout << "check PID " << endl;
618 if(fabs(mcpid)!=11)
619 {
620 cout << "!= 11" << endl;
621 cout << mcpid << endl;
622 }
623 if(mcele==-1)
624 {
625 cout << "mcele==-1" << endl;
626 cout << mcele << endl;
627 cout << mcpid << endl;
628 }
629
bd6ee6dd 630 } // end of mcLabel>-1
631
632 } // end of MC info.
e7d87aef 633
2b4460a6 634 //cout << "Pi0 = " << mcOrgPi0 << " ; Eta = " << mcOrgEta << endl;
50919258 635 //printf("weight = %f\n",mcWeight);
636
9b3495ff 637 if(TMath::Abs(track->Eta())>0.6) continue;
f09766dd 638 if(TMath::Abs(track->Pt()<2.0)) continue;
bfc7c23b 639
640 fTrackPtBefTrkCuts->Fill(track->Pt());
641 // RecKine: ITSTPC cuts
642 if(!ProcessCutStep(AliHFEcuts::kStepRecKineITSTPC, track)) continue;
643
644 //RecKink
645 if(fRejectKinkMother) { // Quick and dirty fix to reject both kink mothers and daughters
646 if(track->GetKinkIndex(0) != 0) continue;
647 }
648
649 // RecPrim
650 if(!ProcessCutStep(AliHFEcuts::kStepRecPrim, track)) continue;
651
652 // HFEcuts: ITS layers cuts
653 if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsITS, track)) continue;
654
f09766dd 655 // HFE cuts: TPC PID cleanup
656 if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTPC, track)) continue;
657
a6df418c 658 int nTPCcl = track->GetTPCNcls();
96167a04 659 //int nTPCclF = track->GetTPCNclsF(); // warnings
a6df418c 660 int nITS = track->GetNcls(0);
bfc7c23b 661
662 fTrackPtAftTrkCuts->Fill(track->Pt());
663
664 Double_t mom = -999., eop=-999., pt = -999., dEdx=-999., fTPCnSigma=-10, phi=-999., eta=-999.;
feffe705 665 pt = track->Pt();
666 if(pt<2.0)continue;
bfc7c23b 667
668 // Track extrapolation
669
bfc7c23b 670 Int_t charge = track->Charge();
671 fTrkpt->Fill(pt);
672 mom = track->P();
673 phi = track->Phi();
674 eta = track->Eta();
675 dEdx = track->GetTPCsignal();
676 fTPCnSigma = fPID->GetPIDResponse() ? fPID->GetPIDResponse()->NumberOfSigmasTPC(track, AliPID::kElectron) : 1000;
8af58b0d 677
bfc7c23b 678 double ncells = -1.0;
679 double m20 = -1.0;
680 double m02 = -1.0;
681 double disp = -1.0;
682 double rmatch = -1.0;
683 double nmatch = -1.0;
f09766dd 684 double oppstatus = 0.0;
59d998de 685 double emctof = 0.0;
bfc7c23b 686
f09766dd 687 Bool_t fFlagPhotonicElec = kFALSE;
688 Bool_t fFlagConvinatElec = kFALSE;
bfc7c23b 689
690 Int_t clsId = track->GetEMCALcluster();
691 if (clsId>0){
692 AliESDCaloCluster *clust = fESD->GetCaloCluster(clsId);
693 if(clust && clust->IsEMCAL()){
694
695 double clustE = clust->E();
696 eop = clustE/fabs(mom);
8efacc22 697 cout << "eop org = "<< eop << endl;
698 if(mcLabel>-1.0)
699 {
700 double mceopcorr = MCEopMeanCorrection(pt,cent);
701 eop += mceopcorr;
702 }
703 cout << "eop corr = " << eop << endl;
704
bfc7c23b 705 //double clustT = clust->GetTOF();
706 ncells = clust->GetNCells();
707 m02 = clust->GetM02();
708 m20 = clust->GetM20();
709 disp = clust->GetDispersion();
710 double delphi = clust->GetTrackDx();
711 double deleta = clust->GetTrackDz();
712 rmatch = sqrt(pow(delphi,2)+pow(deleta,2));
713 nmatch = clust->GetNTracksMatched();
59d998de 714 emctof = clust->GetTOF();
715 //cout << "emctof = " << emctof << endl;
bfc7c23b 716
c2801a73 717 if(fTPCnSigma>-1.5 && fTPCnSigma<3.0)
718 {
2b4460a6 719 SelectPhotonicElectron(iTracks,cent,track,fFlagPhotonicElec,fFlagConvinatElec,fTPCnSigma,m20,eop,mcele,mcWeight,iHijing,mcOrgPi0,mcOrgEta);
c2801a73 720 }
721 if(fFlagPhotonicElec)oppstatus = 1.0;
722 if(fFlagConvinatElec)oppstatus = 2.0;
723 if(fFlagPhotonicElec && fFlagConvinatElec)oppstatus = 3.0;
724
feffe705 725 double valdedx[16];
a6df418c 726 valdedx[0] = pt; valdedx[1] = nITS; valdedx[2] = phi; valdedx[3] = eta; valdedx[4] = fTPCnSigma;
615ffe11 727 //valdedx[5] = eop; valdedx[6] = rmatch; valdedx[7] = ncells, valdedx[8] = nTPCclF; valdedx[9] = m20; valdedx[10] = mcpT;
728 valdedx[5] = eop; valdedx[6] = rmatch; valdedx[7] = ncells, valdedx[8] = nmatch; valdedx[9] = m20; valdedx[10] = mcpT;
a6df418c 729 valdedx[11] = cent; valdedx[12] = charge; valdedx[13] = oppstatus; valdedx[14] = nTPCcl;
feffe705 730 valdedx[15] = mcele;
85985bb0 731 if(fqahist==1)fEleInfo->Fill(valdedx);
bfc7c23b 732
733
734 }
735 }
a6df418c 736
737 if(nITS<2.5)continue;
738 if(nTPCcl<100)continue;
739
740 CheckNclust->Fill(nTPCcl);
741 CheckNits->Fill(nITS);
742
42c75692 743 fdEdxBef->Fill(mom,fTPCnSigma);
bfc7c23b 744 fTPCnsigma->Fill(mom,fTPCnSigma);
f09766dd 745 if(fTPCnSigma >= -1.0 && fTPCnSigma <= 3)fTrkEovPBef->Fill(pt,eop);
bfc7c23b 746
bfc7c23b 747 Int_t pidpassed = 1;
697ecf6b 748
749 // check reco eff. with TPC
750
751 double phoval[5];
752 phoval[0] = cent;
753 phoval[1] = pt;
754 phoval[2] = fTPCnSigma;
755 phoval[3] = iHijing;
756 phoval[4] = mcMompT;
757
758 if((fTPCnSigma >= -1.0 && fTPCnSigma <= 3) && mcele>0 && mcPho && mcOrgPi0)
759 {
760 if(iHijing==1)mcWeight = 1.0;
761 fIncpTMCpho_pi0e_TPC->Fill(phoval,mcWeight);
762 if(fFlagPhotonicElec) fPhoElecPtMC_pi0e_TPC->Fill(phoval,mcWeight);
763 if(fFlagConvinatElec) fSameElecPtMC_pi0e_TPC->Fill(phoval,mcWeight);
764 }
bfc7c23b 765
766 //--- track accepted
767 AliHFEpidObject hfetrack;
768 hfetrack.SetAnalysisType(AliHFEpidObject::kESDanalysis);
769 hfetrack.SetRecTrack(track);
3db00c72 770 double binct = 10.5;
771 if((0.0< cent) && (cent<5.0)) binct = 0.5;
772 if((5.0< cent) && (cent<10.0)) binct = 1.5;
773 if((10.0< cent) && (cent<20.0)) binct = 2.5;
774 if((20.0< cent) && (cent<30.0)) binct = 3.5;
775 if((30.0< cent) && (cent<40.0)) binct = 4.5;
776 if((40.0< cent) && (cent<50.0)) binct = 5.5;
777 if((50.0< cent) && (cent<60.0)) binct = 6.5;
778 if((60.0< cent) && (cent<70.0)) binct = 7.5;
779 if((70.0< cent) && (cent<80.0)) binct = 8.5;
780 if((80.0< cent) && (cent<90.0)) binct = 9.5;
781 if((90.0< cent) && (cent<100.0)) binct = 10.5;
782
783 hfetrack.SetCentrality((int)binct); //added
bfc7c23b 784 hfetrack.SetPbPb();
785 if(!fPID->IsSelected(&hfetrack, NULL, "", fPIDqa)) pidpassed = 0;
786
787 if(pidpassed==0) continue;
788
789 fTrkEovPAft->Fill(pt,eop);
42c75692 790 fdEdxAft->Fill(mom,fTPCnSigma);
bfc7c23b 791
e7d87aef 792 // Fill real data
fb87d707 793 fIncpT->Fill(cent,pt);
bfc7c23b 794 if(fFlagPhotonicElec) fPhoElecPt->Fill(cent,pt);
f09766dd 795 if(fFlagConvinatElec) fSameElecPt->Fill(cent,pt);
fb87d707 796
797 if(m20>0.0 && m20<0.3)
798 {
59d998de 799 fIncpTM20->Fill(cent,pt);
800 ftimingEle->Fill(pt,emctof);
fb87d707 801 if(fFlagPhotonicElec) fPhoElecPtM20->Fill(cent,pt);
802 if(fFlagConvinatElec) fSameElecPtM20->Fill(cent,pt);
803 }
feffe705 804
e7d87aef 805 // MC
70448e3b 806 // check label for electron candidiates
807
808 int idlabel = 1;
809 if(mcLabel==0)idlabel = 0;
810 fLabelCheck->Fill(pt,idlabel);
811 if(mcLabel==0)fgeoFake->Fill(phi,eta);
812
e4b0faf2 813 if(mcele>0) // select MC electrons
e7d87aef 814 {
3db00c72 815
816 fIncpTMChfeAll->Fill(cent,pt);
817 if(m20>0.0 && m20<0.3)fIncpTMCM20hfeAll->Fill(cent,pt);
818
e4b0faf2 819 if(mcBtoE || mcDtoE) // select B->e & D->e
e7d87aef 820 {
821 fIncpTMChfe->Fill(cent,pt);
bd6ee6dd 822 //if(m20>0.0 && m20<0.3)fIncpTMCM20hfe->Fill(cent,pt);
823 //if(m20>0.0 && m20<0.3)fIncpTMCM20hfeCheck->Fill(cent,mcpT);
824 if(m20>0.0 && m20<0.3)
825 {
68d718e7 826 cout << "MC label = " << mcLabel << endl;
bd6ee6dd 827 fIncpTMCM20hfe->Fill(cent,pt);
828 fIncpTMCM20hfeCheck->Fill(cent,mcpT);
829 fIncpTMCM20hfeCheck_weight->Fill(phoval);
830 }
e7d87aef 831 }
e4b0faf2 832
833 if(mcPho) // select photonic electrons
e7d87aef 834 {
44be9e1d 835
836 fIncpTMCpho->Fill(phoval);
837 if(fFlagPhotonicElec) fPhoElecPtMC->Fill(phoval);
838 if(fFlagConvinatElec) fSameElecPtMC->Fill(phoval);
e7d87aef 839
840 if(m20>0.0 && m20<0.3)
841 {
44be9e1d 842 fIncpTMCM20pho->Fill(phoval);
843 if(fFlagPhotonicElec) fPhoElecPtMCM20->Fill(phoval);
844 if(fFlagConvinatElec) fSameElecPtMCM20->Fill(phoval);
e4b0faf2 845 // pi0->g->e
846 if(mcWeight>-1)
847 {
987053ce 848 if(iHijing==1)mcWeight = 1.0;
93c189c5 849 if(mcOrgPi0)
850 {
697ecf6b 851 fIncpTMCM20pho_pi0e->Fill(phoval,mcWeight);
852 if(fFlagPhotonicElec) fPhoElecPtMCM20_pi0e->Fill(phoval,mcWeight);
853 if(fFlagConvinatElec) fSameElecPtMCM20_pi0e->Fill(phoval,mcWeight);
854 //fIncpTMCM20pho_pi0e->Fill(phoval); // v5-04-02-AN & v5-04-06-AN
855 //if(fFlagPhotonicElec) fPhoElecPtMCM20_pi0e->Fill(phoval);
856 //if(fFlagConvinatElec) fSameElecPtMCM20_pi0e->Fill(phoval);
93c189c5 857 }
858 if(mcOrgEta)
859 {
697ecf6b 860 fIncpTMCM20pho_eta->Fill(phoval,mcWeight);
861 if(fFlagPhotonicElec) fPhoElecPtMCM20_eta->Fill(phoval,mcWeight);
862 if(fFlagConvinatElec) fSameElecPtMCM20_eta->Fill(phoval,mcWeight);
863 //fIncpTMCM20pho_eta->Fill(phoval);
864 //if(fFlagPhotonicElec) fPhoElecPtMCM20_eta->Fill(phoval);
865 //if(fFlagConvinatElec) fSameElecPtMCM20_eta->Fill(phoval);
93c189c5 866 }
867 // --- eta
e4b0faf2 868 }
e7d87aef 869 }
870 }
871 }
bfc7c23b 872 }
873 PostData(1, fOutputList);
874}
875//_________________________________________
876void AliAnalysisTaskHFECal::UserCreateOutputObjects()
877{
42c75692 878 //--- Check MC
879
e7d87aef 880 //Bool_t mcData = kFALSE;
42c75692 881 if(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler())
882 {
e7d87aef 883 fmcData = kTRUE;
42c75692 884 printf("+++++ MC Data available");
885 }
e7d87aef 886 if(fmcData)
42c75692 887 {
888 printf("++++++++= MC analysis \n");
889 }
890 else
891 {
892 printf("++++++++ real data analysis \n");
893 }
894
85985bb0 895 printf("+++++++ QA hist %d",fqahist);
896
fb87d707 897 //---- Geometry
898 fGeom = AliEMCALGeometry::GetInstance("EMCAL_COMPLETEV1");
899
bfc7c23b 900 //--------Initialize PID
42c75692 901 //fPID->SetHasMCData(kFALSE);
e7d87aef 902 fPID->SetHasMCData(fmcData);
bfc7c23b 903 if(!fPID->GetNumberOfPIDdetectors())
904 {
905 fPID->AddDetector("TPC", 0);
906 fPID->AddDetector("EMCAL", 1);
907 }
908
3db00c72 909 Double_t params[4];
78ff1095 910 const char *cutmodel;
3db00c72 911 cutmodel = "pol0";
912 params[0] = -1.0; //sigma min
9b3495ff 913 double maxnSig = 3.0;
914 if(fmcData)
915 {
916 params[0] = -5.0; //sigma min
917 maxnSig = 5.0;
918 }
919
920 for(Int_t a=0;a<11;a++)fPID->ConfigureTPCcentralityCut(a,cutmodel,params,maxnSig);
3db00c72 921
bfc7c23b 922 fPID->SortDetectors();
923 fPIDqa = new AliHFEpidQAmanager();
924 fPIDqa->Initialize(fPID);
a6df418c 925
926 //------- fcut --------------
927 fCuts = new AliHFEcuts();
928 fCuts->CreateStandardCuts();
929 //fCuts->SetMinNClustersTPC(100);
930 fCuts->SetMinNClustersTPC(90);
931 fCuts->SetMinRatioTPCclusters(0.6);
932 fCuts->SetTPCmodes(AliHFEextraCuts::kFound, AliHFEextraCuts::kFoundOverFindable);
933 //fCuts->SetMinNClustersITS(3);
934 fCuts->SetMinNClustersITS(2);
935 fCuts->SetCutITSpixel(AliHFEextraCuts::kAny);
936 fCuts->SetCheckITSLayerStatus(kFALSE);
937 fCuts->SetVertexRange(10.);
938 fCuts->SetTOFPIDStep(kFALSE);
939 fCuts->SetPtRange(2, 50);
940 fCuts->SetMaxImpactParam(3.,3.);
941
bfc7c23b 942 //--------Initialize correction Framework and Cuts
943 fCFM = new AliCFManager;
944 const Int_t kNcutSteps = AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kNcutStepsRecTrack + AliHFEcuts::kNcutStepsDETrack;
945 fCFM->SetNStepParticle(kNcutSteps);
946 for(Int_t istep = 0; istep < kNcutSteps; istep++)
947 fCFM->SetParticleCutsList(istep, NULL);
948
949 if(!fCuts){
950 AliWarning("Cuts not available. Default cuts will be used");
951 fCuts = new AliHFEcuts;
952 fCuts->CreateStandardCuts();
953 }
954 fCuts->Initialize(fCFM);
955
956 //---------Output Tlist
957 fOutputList = new TList();
958 fOutputList->SetOwner();
959 fOutputList->Add(fPIDqa->MakeList("PIDQA"));
960
961 fNoEvents = new TH1F("fNoEvents","",4,-0.5,3.5) ;
962 fOutputList->Add(fNoEvents);
963
02cd82a0 964 Int_t binsE[5] = {250, 100, 1000, 200, 10};
fb87d707 965 Double_t xminE[5] = {1.0, -1, 0.0, 0, -0.5};
966 Double_t xmaxE[5] = {3.5, 1, 100.0, 100, 9.5};
967 fEMCAccE = new THnSparseD("fEMCAccE","EMC acceptance & E;#eta;#phi;Energy;Centrality;trugCondition;",5,binsE,xminE,xmaxE);
85985bb0 968 if(fqahist==1)fOutputList->Add(fEMCAccE);
bfc7c23b 969
f4e0d2d5 970 hEMCAccE = new TH2F("hEMCAccE","Cluster Energy",200,0,100,100,0,20);
971 fOutputList->Add(hEMCAccE);
972
bfc7c23b 973 fTrkpt = new TH1F("fTrkpt","track pt",100,0,50);
974 fOutputList->Add(fTrkpt);
975
976 fTrackPtBefTrkCuts = new TH1F("fTrackPtBefTrkCuts","track pt before track cuts",100,0,50);
977 fOutputList->Add(fTrackPtBefTrkCuts);
978
979 fTrackPtAftTrkCuts = new TH1F("fTrackPtAftTrkCuts","track pt after track cuts",100,0,50);
980 fOutputList->Add(fTrackPtAftTrkCuts);
981
982 fTPCnsigma = new TH2F("fTPCnsigma", "TPC - n sigma",100,0,50,200,-10,10);
983 fOutputList->Add(fTPCnsigma);
984
985 fTrkEovPBef = new TH2F("fTrkEovPBef","track E/p before HFE pid",100,0,50,100,0,2);
986 fOutputList->Add(fTrkEovPBef);
987
988 fTrkEovPAft = new TH2F("fTrkEovPAft","track E/p after HFE pid",100,0,50,100,0,2);
989 fOutputList->Add(fTrkEovPAft);
990
42c75692 991 fdEdxBef = new TH2F("fdEdxBef","track dEdx vs p before HFE pid",100,0,50,200,-10,10);
bfc7c23b 992 fOutputList->Add(fdEdxBef);
993
42c75692 994 fdEdxAft = new TH2F("fdEdxAft","track dEdx vs p after HFE pid",100,0,50,200,-10,10);
bfc7c23b 995 fOutputList->Add(fdEdxAft);
996
02cd82a0 997 fIncpT = new TH2F("fIncpT","HFE pid electro vs. centrality",200,0,100,100,0,50);
bfc7c23b 998 fOutputList->Add(fIncpT);
999
02cd82a0 1000 fIncpTM20 = new TH2F("fIncpTM20","HFE pid electro vs. centrality with M20",200,0,100,100,0,50);
fb87d707 1001 fOutputList->Add(fIncpTM20);
4e01b68c 1002
02cd82a0 1003 Int_t nBinspho[9] = { 200, 100, 500, 12, 50, 4, 200, 8, 100};
b12bc641 1004 Double_t minpho[9] = { 0., 0., 0., -2.5, 0, -0.5, 0,-1.5, 0};
1005 Double_t maxpho[9] = {100., 50., 0.5, 3.5, 1, 3.5, 2, 6.5, 50};
f09766dd 1006
b12bc641 1007 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 1008 if(fqahist==1)fOutputList->Add(fInvmassLS);
bfc7c23b 1009
b12bc641 1010 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 1011 if(fqahist==1)fOutputList->Add(fInvmassULS);
bfc7c23b 1012
8806ce6c 1013 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 1014 if(fqahist==1)fOutputList->Add(fInvmassLSmc);
8806ce6c 1015
1016 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 1017 if(fqahist==1)fOutputList->Add(fInvmassULSmc);
8806ce6c 1018
93c189c5 1019 fInvmassLSmc0 = new TH2D("fInvmassLSmc0", "Inv mass of LS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2)",20,0,20,500,0,0.5 );
1020 fInvmassLSmc0->Sumw2();
1021 fOutputList->Add(fInvmassLSmc0);
1022
1023 fInvmassLSmc1 = new TH2D("fInvmassLSmc1", "Inv mass of LS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2)",20,0,20,500,0,0.5 );
1024 fInvmassLSmc1->Sumw2();
1025 fOutputList->Add(fInvmassLSmc1);
1026
1027 fInvmassLSmc2 = new TH2D("fInvmassLSmc2", "Inv mass of LS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2)",20,0,20,500,0,0.5 );
1028 fInvmassLSmc2->Sumw2();
1029 fOutputList->Add(fInvmassLSmc2);
1030
1031 fInvmassLSmc3 = new TH2D("fInvmassLSmc3", "Inv mass of LS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2)",20,0,20,500,0,0.5 );
1032 fInvmassLSmc3->Sumw2();
1033 fOutputList->Add(fInvmassLSmc3);
1034
1035 fInvmassULSmc0 = new TH2D("fInvmassULSmc0", "Inv mass of ULS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2)",20,0,20,500,0,0.5 );
1036 fInvmassULSmc0->Sumw2();
1037 fOutputList->Add(fInvmassULSmc0);
1038
1039 fInvmassULSmc1 = new TH2D("fInvmassULSmc1", "Inv mass of ULS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2)",20,0,20,500,0,0.5 );
1040 fInvmassULSmc1->Sumw2();
1041 fOutputList->Add(fInvmassULSmc1);
1042
1043 fInvmassULSmc2 = new TH2D("fInvmassULSmc2", "Inv mass of ULS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2)",20,0,20,500,0,0.5 );
1044 fInvmassULSmc2->Sumw2();
1045 fOutputList->Add(fInvmassULSmc2);
1046
1047 fInvmassULSmc3 = new TH2D("fInvmassULSmc3", "Inv mass of ULS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2)",20,0,20,500,0,0.5 );
1048 fInvmassULSmc3->Sumw2();
93c189c5 1049 fOutputList->Add(fInvmassULSmc3);
1050
bfc7c23b 1051 fOpeningAngleLS = new TH1F("fOpeningAngleLS","Opening angle for LS pairs",100,0,1);
1052 fOutputList->Add(fOpeningAngleLS);
1053
1054 fOpeningAngleULS = new TH1F("fOpeningAngleULS","Opening angle for ULS pairs",100,0,1);
1055 fOutputList->Add(fOpeningAngleULS);
1056
1057 fPhotoElecPt = new TH1F("fPhotoElecPt", "photonic electron pt",100,0,50);
1058 fOutputList->Add(fPhotoElecPt);
1059
02cd82a0 1060 fPhoElecPt = new TH2F("fPhoElecPt", "Pho-inclusive electron pt",200,0,100,100,0,50);
bfc7c23b 1061 fOutputList->Add(fPhoElecPt);
1062
02cd82a0 1063 fPhoElecPtM20 = new TH2F("fPhoElecPtM20", "Pho-inclusive electron pt with M20",200,0,100,100,0,50);
fb87d707 1064 fOutputList->Add(fPhoElecPtM20);
1065
02cd82a0 1066 fSameElecPt = new TH2F("fSameElecPt", "Same-inclusive electron pt",200,0,100,100,0,50);
f09766dd 1067 fOutputList->Add(fSameElecPt);
1068
02cd82a0 1069 fSameElecPtM20 = new TH2F("fSameElecPtM20", "Same-inclusive electron pt with M20",200,0,100,100,0,50);
fb87d707 1070 fOutputList->Add(fSameElecPtM20);
1071
02cd82a0 1072 fCent = new TH1F("fCent","Centrality",200,0,100) ;
bfc7c23b 1073 fOutputList->Add(fCent);
1074
1075 // Make common binning
44be9e1d 1076 const Double_t kMinP = 0.;
bfc7c23b 1077 const Double_t kMaxP = 50.;
a6df418c 1078 //const Double_t kTPCSigMim = 40.;
1079 //const Double_t kTPCSigMax = 140.;
bfc7c23b 1080
1081 // 1st histogram: TPC dEdx with/without EMCAL (p, pT, TPC Signal, phi, eta, Sig, e/p, ,match, cell, M02, M20, Disp, Centrality, select)
a6df418c 1082 Int_t nBins[16] = { 250, 10, 60, 20, 100, 300, 50, 40, 200, 200, 250, 200, 3, 5, 100, 8};
1083 Double_t min[16] = {kMinP, -0.5, 1.0, -1.0, -6.0, 0, 0, 0, 0.0, 0.0, 0.0, 0, -1.5, -0.5, 80, -1.5};
6e92c428 1084 Double_t max[16] = {kMaxP, 9.5, 4.0, 1.0, 4.0, 3.0, 0.1, 40, 200, 2.0, 50.0, 100, 1.5, 4.5, 180, 6.5};
1085 fEleInfo = new THnSparseD("fEleInfo", "Electron Info; pT [GeV/c]; TPC signal;phi;eta;nSig; E/p;Rmatch;Ncell;clsF;M20;mcpT;Centrality;charge;opp;same;trigCond;MCele", 16, nBins, min, max);
85985bb0 1086 if(fqahist==1)fOutputList->Add(fEleInfo);
bfc7c23b 1087
fb87d707 1088 //<--- Trigger info
feffe705 1089 /*
fb87d707 1090 fClsEBftTrigCut = new TH1F("fClsEBftTrigCut","cluster E before trigger selection",1000,0,100);
1091 fOutputList->Add(fClsEBftTrigCut);
1092
1093 fClsEAftTrigCut = new TH1F("fClsEAftTrigCut","cluster E if cls has 0 trigcut channel",1000,0,100);
1094 fOutputList->Add(fClsEAftTrigCut);
1095
1096 fClsEAftTrigCut1 = new TH1F("fClsEAftTrigCut1","cluster E if cls with trig channel",1000,0,100);
1097 fOutputList->Add(fClsEAftTrigCut1);
1098
1099 fClsEAftTrigCut2 = new TH1F("fClsEAftTrigCut2","cluster E if cls with trigcut channel",1000,0,100);
1100 fOutputList->Add(fClsEAftTrigCut2);
1101
1102 fClsEAftTrigCut3 = new TH1F("fClsEAftTrigCut3","cluster E if cls with trigcut channel + nCell>Ecorrect",1000,0,100);
1103 fOutputList->Add(fClsEAftTrigCut3);
1104
1105 fClsEAftTrigCut4 = new TH1F("fClsEAftTrigCut4","cluster E if cls with trigcut channel + nCell>Ecorrect + cls time cut",1000,0,100);
1106 fOutputList->Add(fClsEAftTrigCut4);
1107
1108 fClsETime = new TH2F("fClsETime", "Cls time vs E; E; time;",1000,0,100,1000,-0.0000002,0.0000009);
1109 fOutputList->Add(fClsETime);
1110
1111 fClsETime1 = new TH2F("fClsETime1", "Cls time vs E if cls contains trigger channel; E; time;",1000,0,100,1000,-0.0000002,0.0000009);
1112 fOutputList->Add(fClsETime1);
1113
1114 fTrigTimes = new TH1F("fTrigTimes", "Trigger time; time; N;",25,0,25);
1115 fOutputList->Add(fTrigTimes);
1116
42c75692 1117 fCellCheck = new TH2F("fCellCheck", "Cell vs E; E GeV; Cell ID",10,6,26,12000,0,12000);
1118 fOutputList->Add(fCellCheck);
feffe705 1119 */
e7d87aef 1120 //<---------- MC
1121
02cd82a0 1122 fInputHFEMC = new TH2F("fInputHFEMC","Input MC HFE pid electro vs. centrality",200,0,100,100,0,50);
e7d87aef 1123 fOutputList->Add(fInputHFEMC);
1124
02cd82a0 1125 fInputAlle = new TH2F("fInputAlle","Input MC electro vs. centrality",200,0,100,100,0,50);
feffe705 1126 fOutputList->Add(fInputAlle);
1127
02cd82a0 1128 fIncpTMChfe = new TH2F("fIncpTMChfe","MC HFE pid electro vs. centrality",200,0,100,100,0,50);
e7d87aef 1129 fOutputList->Add(fIncpTMChfe);
1130
02cd82a0 1131 fIncpTMChfeAll = new TH2F("fIncpTMChfeAll","MC Alle pid electro vs. centrality",200,0,100,100,0,50);
3db00c72 1132 fOutputList->Add(fIncpTMChfeAll);
1133
02cd82a0 1134 fIncpTMCM20hfe = new TH2F("fIncpTMCM20hfe","MC HFE pid electro vs. centrality with M20",200,0,100,100,0,50);
e7d87aef 1135 fOutputList->Add(fIncpTMCM20hfe);
1136
02cd82a0 1137 fIncpTMCM20hfeAll = new TH2F("fIncpTMCM20hfeAll","MC Alle pid electro vs. centrality with M20",200,0,100,100,0,50);
3db00c72 1138 fOutputList->Add(fIncpTMCM20hfeAll);
1139
60544aea 1140 fIncpTMCM20hfeCheck = new TH2F("fIncpTMCM20hfeCheck","MC HFE pid electro vs. centrality with M20 Check",200,0,100,100,0,50);
1141 fOutputList->Add(fIncpTMCM20hfeCheck);
44be9e1d 1142
2b4460a6 1143 Int_t nBinspho2[5] = { 200, 100, 7, 3, 700};
987053ce 1144 Double_t minpho2[5] = { 0., 0., -2.5, -0.5, 0.};
2b4460a6 1145 Double_t maxpho2[5] = {100., 50., 4.5, 2.5, 70.};
bd6ee6dd 1146
1147 fInputHFEMC_weight = new THnSparseD("fInputHFEMC_weight", "MC HFE electron pt",5,nBinspho2,minpho2,maxpho2);
1148 fOutputList->Add(fInputHFEMC_weight);
1149
1150 fIncpTMCM20hfeCheck_weight = new THnSparseD("fIncpTMCM20hfeCheck_weight", "HFE electron pt with M20",5,nBinspho2,minpho2,maxpho2);
1151 fOutputList->Add(fIncpTMCM20hfeCheck_weight);
1152
987053ce 1153 fIncpTMCpho = new THnSparseD("fIncpTMCpho","MC Pho pid electro vs. centrality",5,nBinspho2,minpho2,maxpho2);
e7d87aef 1154 fOutputList->Add(fIncpTMCpho);
1155
987053ce 1156 fIncpTMCM20pho = new THnSparseD("fIncpTMCM20pho","MC Pho pid electro vs. centrality with M20",5,nBinspho2,minpho2,maxpho2);
e7d87aef 1157 fOutputList->Add(fIncpTMCM20pho);
1158
987053ce 1159 fPhoElecPtMC = new THnSparseD("fPhoElecPtMC", "MC Pho-inclusive electron pt",5,nBinspho2,minpho2,maxpho2);
e7d87aef 1160 fOutputList->Add(fPhoElecPtMC);
1161
987053ce 1162 fPhoElecPtMCM20 = new THnSparseD("fPhoElecPtMCM20", "MC Pho-inclusive electron pt with M20",5,nBinspho2,minpho2,maxpho2);
feffe705 1163 fOutputList->Add(fPhoElecPtMCM20);
e7d87aef 1164
987053ce 1165 fSameElecPtMC = new THnSparseD("fSameElecPtMC", "MC Same-inclusive electron pt",5,nBinspho2,minpho2,maxpho2);
e7d87aef 1166 fOutputList->Add(fSameElecPtMC);
1167
987053ce 1168 fSameElecPtMCM20 = new THnSparseD("fSameElecPtMCM20", "MC Same-inclusive electron pt with M20",5,nBinspho2,minpho2,maxpho2);
e7d87aef 1169 fOutputList->Add(fSameElecPtMCM20);
1170
987053ce 1171 fIncpTMCM20pho_pi0e = new THnSparseD("fIncpTMCM20pho_pi0e","MC Pho pi0->e pid electro vs. centrality with M20",5,nBinspho2,minpho2,maxpho2);
1172 fIncpTMCM20pho_pi0e->Sumw2();
e4b0faf2 1173 fOutputList->Add(fIncpTMCM20pho_pi0e);
1174
987053ce 1175 fPhoElecPtMCM20_pi0e = new THnSparseD("fPhoElecPtMCM20_pi0e", "MC Pho-inclusive electron pt with M20 pi0->e",5,nBinspho2,minpho2,maxpho2);
1176 fPhoElecPtMCM20_pi0e->Sumw2();
e4b0faf2 1177 fOutputList->Add(fPhoElecPtMCM20_pi0e);
1178
987053ce 1179 fSameElecPtMCM20_pi0e = new THnSparseD("fSameElecPtMCM20_pi0e", "MC Same-inclusive electron pt pi0->e",5,nBinspho2,minpho2,maxpho2);
1180 fSameElecPtMCM20_pi0e->Sumw2();
e4b0faf2 1181 fOutputList->Add(fSameElecPtMCM20_pi0e);
697ecf6b 1182 //
93c189c5 1183 fIncpTMCM20pho_eta = new THnSparseD("fIncpTMCM20pho_eta","MC Pho pi0->e pid electro vs. centrality with M20",5,nBinspho2,minpho2,maxpho2);
1184 fIncpTMCM20pho_eta->Sumw2();
1185 fOutputList->Add(fIncpTMCM20pho_eta);
1186
1187 fPhoElecPtMCM20_eta = new THnSparseD("fPhoElecPtMCM20_eta", "MC Pho-inclusive electron pt with M20 pi0->e",5,nBinspho2,minpho2,maxpho2);
1188 fPhoElecPtMCM20_eta->Sumw2();
1189 fOutputList->Add(fPhoElecPtMCM20_eta);
1190
1191 fSameElecPtMCM20_eta = new THnSparseD("fSameElecPtMCM20_eta", "MC Same-inclusive electron pt pi0->e",5,nBinspho2,minpho2,maxpho2);
1192 fSameElecPtMCM20_eta->Sumw2();
1193 fOutputList->Add(fSameElecPtMCM20_eta);
697ecf6b 1194 // ------------
1195 fIncpTMCpho_pi0e_TPC = new THnSparseD("fIncpTMCpho_pi0e_TPC","MC Pho pi0->e pid electro vs. centrality with M20",5,nBinspho2,minpho2,maxpho2);
1196 fIncpTMCpho_pi0e_TPC->Sumw2();
1197 fOutputList->Add(fIncpTMCpho_pi0e_TPC);
1198
1199 fPhoElecPtMC_pi0e_TPC = new THnSparseD("fPhoElecPtMC_pi0e_TPC", "MC Pho-inclusive electron pt with pi0->e",5,nBinspho2,minpho2,maxpho2);
1200 fPhoElecPtMC_pi0e_TPC->Sumw2();
1201 fOutputList->Add(fPhoElecPtMC_pi0e_TPC);
1202
1203 fSameElecPtMC_pi0e_TPC = new THnSparseD("fSameElecPtMC_pi0e_TPC", "MC Same-inclusive electron pt pi0->e",5,nBinspho2,minpho2,maxpho2);
1204 fSameElecPtMC_pi0e_TPC->Sumw2();
1205 fOutputList->Add(fSameElecPtMC_pi0e_TPC);
1206 //-------------
e4b0faf2 1207
a6df418c 1208 CheckNclust = new TH1D("CheckNclust","cluster check",200,0,200);
1209 fOutputList->Add(CheckNclust);
1210
1211 CheckNits = new TH1D("CheckNits","ITS cluster check",8,-0.5,7.5);
1212 fOutputList->Add(CheckNits);
68d718e7 1213 /*
d113d7cd 1214 Hpi0pTcheck = new TH2D("Hpi0pTcheck","Pi0 pT from Hijing",100,0,50,3,-0.5,2.5);
1215 fOutputList->Add(Hpi0pTcheck);
1216
19325033 1217 HETApTcheck = new TH2D("HETApTcheck","Eta pT from Hijing",100,0,50,3,-0.5,2.5);
1218 fOutputList->Add(HETApTcheck);
68d718e7 1219 */
1220
1221 Int_t nBinspho3[3] = { 200, 100, 3};
1222 Double_t minpho3[3] = { 0., 0., -0.5};
1223 Double_t maxpho3[3] = {100., 50., 2.5};
19325033 1224
68d718e7 1225 Hpi0pTcheck = new THnSparseD("Hpi0pTcheck","Pi0 pT from Hijing",3,nBinspho3,minpho3,maxpho3);
1226 fOutputList->Add(Hpi0pTcheck);
1227
1228 HETApTcheck = new THnSparseD("HETApTcheck","Eta pT from Hijing",3,nBinspho3,minpho3,maxpho3);
1229 fOutputList->Add(HETApTcheck);
1230 //--
e4b0faf2 1231 HphopTcheck = new TH2D("HphopTcheck","Pho pT from Hijing",100,0,50,3,-0.5,2.5);
1232 fOutputList->Add(HphopTcheck);
5e0e45b3 1233 //
1234 HDpTcheck = new TH2D("HDpTcheck","D pT from Hijing",100,0,50,3,-0.5,2.5);
1235 fOutputList->Add(HDpTcheck);
1236
1237 HBpTcheck = new TH2D("HBpTcheck","B pT from Hijing",100,0,50,3,-0.5,2.5);
1238 fOutputList->Add(HBpTcheck);
1239 //
e4b0faf2 1240
93c189c5 1241 fpTCheck = new TH1D("fpTCheck","pT check",500,0,50);
1242 fOutputList->Add(fpTCheck);
1243
50919258 1244 fMomDtoE = new TH2D("fMomDtoE","D->E pT correlations;e p_{T} GeV/c;D p_{T} GeV/c",400,0,40,400,0,40);
1245 fOutputList->Add(fMomDtoE);
d113d7cd 1246
70448e3b 1247 fLabelCheck = new TH2D("fLabelCheck","MC label",50,0,50,5,-1.5,3.5);
1248 fOutputList->Add(fLabelCheck);
1249
1250 fgeoFake = new TH2D("fgeoFake","Label==0 eta and phi",628,0,6.28,200,-1,1);
1251 fOutputList->Add(fgeoFake);
1252
59d998de 1253 ftimingEle = new TH2D("ftimingEle","electron TOF",100,0,20,100,1e-7,1e-6);
1254 fOutputList->Add(ftimingEle);
1255
bfc7c23b 1256 PostData(1,fOutputList);
1257}
1258
1259//________________________________________________________________________
1260void AliAnalysisTaskHFECal::Terminate(Option_t *)
1261{
1262 // Info("Terminate");
1263 AliAnalysisTaskSE::Terminate();
1264}
1265
1266//________________________________________________________________________
1267Bool_t AliAnalysisTaskHFECal::ProcessCutStep(Int_t cutStep, AliVParticle *track)
1268{
1269 // Check single track cuts for a given cut step
1270 const Int_t kMCOffset = AliHFEcuts::kNcutStepsMCTrack;
1271 if(!fCFM->CheckParticleCuts(cutStep + kMCOffset, track)) return kFALSE;
1272 return kTRUE;
1273}
1274//_________________________________________
f09766dd 1275//void AliAnalysisTaskHFECal::SelectPhotonicElectron(Int_t itrack, Double_t cent, AliESDtrack *track, Bool_t &fFlagPhotonicElec)
c2801a73 1276//void AliAnalysisTaskHFECal::SelectPhotonicElectron(Int_t itrack, Double_t cent, AliESDtrack *track, Bool_t &fFlagPhotonicElec, Bool_t &fFlagConvinatElec, Double_t nSig)
2b4460a6 1277void 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 1278{
1279 //Identify non-heavy flavour electrons using Invariant mass method
1280
fb87d707 1281 fTrackCuts->SetAcceptKinkDaughters(kFALSE);
feffe705 1282 fTrackCuts->SetRequireTPCRefit(kTRUE);
85985bb0 1283 fTrackCuts->SetRequireITSRefit(kTRUE);
feffe705 1284 fTrackCuts->SetEtaRange(-0.9,0.9);
f09766dd 1285 //fTrackCuts->SetRequireSigmaToVertex(kTRUE);
bfc7c23b 1286 fTrackCuts->SetMaxChi2PerClusterTPC(3.5);
feffe705 1287 fTrackCuts->SetMinNClustersTPC(90);
bfc7c23b 1288
1289 const AliESDVertex *pVtx = fESD->GetPrimaryVertex();
1290
93c189c5 1291 double ptEle = track->Pt(); //add
1292 if(ibgevent==0 && w > 0.0)
1293 {
1294 fpTCheck->Fill(ptEle,w);
1295 }
1296
bfc7c23b 1297 Bool_t flagPhotonicElec = kFALSE;
f09766dd 1298 Bool_t flagConvinatElec = kFALSE;
bfc7c23b 1299
8806ce6c 1300 int p1 = 0;
1301 if(mce==3)
1302 {
1303 Int_t label = TMath::Abs(track->GetLabel());
1304 TParticle* particle = stack->Particle(label);
1305 p1 = particle->GetFirstMother();
1306 }
1307
5e65985c 1308 //for(Int_t jTracks = itrack+1; jTracks<fESD->GetNumberOfTracks(); jTracks++){
1309 for(Int_t jTracks = 0; jTracks<fESD->GetNumberOfTracks(); jTracks++){
bfc7c23b 1310 AliESDtrack* trackAsso = fESD->GetTrack(jTracks);
1311 if (!trackAsso) {
1312 printf("ERROR: Could not receive track %d\n", jTracks);
1313 continue;
1314 }
5e65985c 1315 if(itrack==jTracks)continue;
d86a815c 1316 int jbgevent = 0;
5e65985c 1317
8806ce6c 1318 int p2 = 0;
1319 if(mce==3)
1320 {
1321 Int_t label2 = TMath::Abs(trackAsso->GetLabel());
1322 TParticle* particle2 = stack->Particle(label2);
d86a815c 1323 Bool_t MChijing_ass = fMC->IsFromBGEvent(label2);
1324 if(MChijing_ass)jbgevent =1;
8806ce6c 1325 if(particle2->GetFirstMother()>-1)
1326 p2 = particle2->GetFirstMother();
1327 }
1328
f09766dd 1329 Double_t dEdxAsso = -999., ptPrim=-999., ptAsso=-999., openingAngle = -999.;
bfc7c23b 1330 Double_t mass=999., width = -999;
1331 Bool_t fFlagLS=kFALSE, fFlagULS=kFALSE;
1332
93c189c5 1333 //ptPrim = track->Pt();
1334 ptPrim = ptEle;
f09766dd 1335
bfc7c23b 1336 dEdxAsso = trackAsso->GetTPCsignal();
1337 ptAsso = trackAsso->Pt();
1338 Int_t chargeAsso = trackAsso->Charge();
1339 Int_t charge = track->Charge();
1340
b12bc641 1341
3db00c72 1342 if(ptAsso <0.5) continue;
bfc7c23b 1343 if(!fTrackCuts->AcceptTrack(trackAsso)) continue;
f09766dd 1344 if(dEdxAsso <65 || dEdxAsso>100) continue; //11a pass1
bfc7c23b 1345
1346 Int_t fPDGe1 = 11; Int_t fPDGe2 = 11;
1347 if(charge>0) fPDGe1 = -11;
1348 if(chargeAsso>0) fPDGe2 = -11;
b12bc641 1349
1350 //printf("chargeAsso = %d\n",chargeAsso);
1351 //printf("charge = %d\n",charge);
bfc7c23b 1352 if(charge == chargeAsso) fFlagLS = kTRUE;
1353 if(charge != chargeAsso) fFlagULS = kTRUE;
1354
b12bc641 1355 //printf("fFlagLS = %d\n",fFlagLS);
1356 //printf("fFlagULS = %d\n",fFlagULS);
1357 //printf("\n");
1358
bfc7c23b 1359 AliKFParticle ge1(*track, fPDGe1);
1360 AliKFParticle ge2(*trackAsso, fPDGe2);
1361 AliKFParticle recg(ge1, ge2);
1362
1363 if(recg.GetNDF()<1) continue;
1364 Double_t chi2recg = recg.GetChi2()/recg.GetNDF();
1365 if(TMath::Sqrt(TMath::Abs(chi2recg))>3.) continue;
1366
d86a815c 1367 // for v5-03-70-AN comment
bfc7c23b 1368 AliKFVertex primV(*pVtx);
1369 primV += recg;
1370 recg.SetProductionVertex(primV);
1371
d113d7cd 1372 recg.SetMassConstraint(0,0.0001);
d86a815c 1373
bfc7c23b 1374 openingAngle = ge1.GetAngle(ge2);
1375 if(fFlagLS) fOpeningAngleLS->Fill(openingAngle);
1376 if(fFlagULS) fOpeningAngleULS->Fill(openingAngle);
1377
bfc7c23b 1378
1379 recg.GetMass(mass,width);
1380
c2801a73 1381 double ishower = 0;
1382 if(shower>0.0 && shower<0.3)ishower = 1;
1383
b12bc641 1384 double phoinfo[9];
f09766dd 1385 phoinfo[0] = cent;
1386 phoinfo[1] = ptPrim;
1387 phoinfo[2] = mass;
3db00c72 1388 phoinfo[3] = nSig;
4e01b68c 1389 phoinfo[4] = openingAngle;
c2801a73 1390 phoinfo[5] = ishower;
1391 phoinfo[6] = ep;
1392 phoinfo[7] = mce;
b12bc641 1393 phoinfo[8] = ptAsso;
f09766dd 1394
1395 if(fFlagLS) fInvmassLS->Fill(phoinfo);
1396 if(fFlagULS) fInvmassULS->Fill(phoinfo);
d86a815c 1397 if(fFlagLS && ibgevent==0 && jbgevent==0) fInvmassLSmc->Fill(phoinfo,w);
93c189c5 1398 if(fFlagULS && ibgevent==0 && jbgevent==0)
1399 {
1400 fInvmassULSmc->Fill(phoinfo,w);
1401 }
b12bc641 1402 //printf("fInvmassCut %f\n",fInvmassCut);
1403 //printf("openingAngle %f\n",fOpeningAngleCut);
1404
ffa14dda 1405 if(openingAngle > fOpeningAngleCut) continue;
bfc7c23b 1406
d86a815c 1407 // for real data
93c189c5 1408 //printf("mce =%f\n",mce);
d86a815c 1409 if(mce<-0.5) // mce==-1. is real
1410 {
93c189c5 1411 //printf("Real data\n");
d86a815c 1412 if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec){
1413 //if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec && (p1==p2)){ <--- only MC train (55,56) v5-03-68-AN & 69 for check
1414 flagPhotonicElec = kTRUE;
1415 }
1416 if(mass<fInvmassCut && fFlagLS && !flagConvinatElec){
1417 flagConvinatElec = kTRUE;
1418 }
1419 }
1420 // for MC data
1421 else
1422 {
93c189c5 1423 //printf("MC data\n");
1424
1425 if(w>0.0)
1426 {
2b4460a6 1427 //cout << "tagpi0 = " << tagpi0 << " ; tageta = " << tageta << endl;
1428 if(fFlagLS && ibgevent==0 && jbgevent==0 && tagpi0) fInvmassLSmc0->Fill(ptPrim,mass,w);
1429 if(fFlagULS && ibgevent==0 && jbgevent==0 && tagpi0) fInvmassULSmc0->Fill(ptPrim,mass,w);
1430 if(fFlagLS && ibgevent==0 && jbgevent==0 && tageta) fInvmassLSmc1->Fill(ptPrim,mass,w);
1431 if(fFlagULS && ibgevent==0 && jbgevent==0 && tageta) fInvmassULSmc1->Fill(ptPrim,mass,w);
1432 if(fFlagLS && ibgevent==0 && jbgevent==0 && (p1==p2) && tagpi0) fInvmassLSmc2->Fill(ptPrim,mass,w);
1433 if(fFlagULS && ibgevent==0 && jbgevent==0 && (p1==p2) && tagpi0) fInvmassULSmc2->Fill(ptPrim,mass,w);
1434 if(fFlagLS && ibgevent==0 && jbgevent==0 && (p1==p2) && tageta) fInvmassLSmc3->Fill(ptPrim,mass,w);
1435 if(fFlagULS && ibgevent==0 && jbgevent==0 && (p1==p2) && tageta) fInvmassULSmc3->Fill(ptPrim,mass,w);
93c189c5 1436 }
d86a815c 1437 if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec && (ibgevent==jbgevent)){
1438 //if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec && (p1==p2)){ <--- only MC train (55,56) v5-03-68-AN & 69 for check
1439 flagPhotonicElec = kTRUE;
1440 }
1441 if(mass<fInvmassCut && fFlagLS && !flagConvinatElec && (ibgevent==jbgevent)){
1442 flagConvinatElec = kTRUE;
1443 }
1444 }
1445
bfc7c23b 1446 }
1447 fFlagPhotonicElec = flagPhotonicElec;
f09766dd 1448 fFlagConvinatElec = flagConvinatElec;
bfc7c23b 1449
1450}
8806ce6c 1451//-------------------------------------------
46305ed6 1452
1453void AliAnalysisTaskHFECal::FindMother(TParticle* part, int &label, int &pid)
1454{
1455 //int label = 99999;
1456 //int pid = 99999;
1457
1458 if(part->GetFirstMother()>-1)
1459 {
1460 label = part->GetFirstMother();
1461 pid = stack->Particle(label)->GetPdgCode();
1462 }
70448e3b 1463 //cout << "Find Mother : label = " << label << " ; pid" << pid << endl;
46305ed6 1464}
1465
8806ce6c 1466double AliAnalysisTaskHFECal::GetMCweight(double mcPi0pT)
1467{
1468 double weight = 1.0;
1469
1470 if(mcPi0pT>0.0 && mcPi0pT<5.0)
1471 {
1472 weight = 0.323*mcPi0pT/(TMath::Exp(-1.6+0.767*mcPi0pT+0.0285*mcPi0pT*mcPi0pT));
1473 }
1474 else
1475 {
1476 weight = 115.0/(0.718*mcPi0pT*TMath::Power(mcPi0pT,3.65));
1477 }
1478 return weight;
1479}
fb87d707 1480
93c189c5 1481double AliAnalysisTaskHFECal::GetMCweightEta(double mcEtapT)
1482{
1483 double weight = 1.0;
1484
1485 weight = 223.3/TMath::Power((TMath::Exp(-0.17*mcEtapT-0.0322*mcEtapT*mcEtapT)+mcEtapT/1.69),5.65);
1486 return weight;
1487}
1488
1489
fb87d707 1490//_________________________________________
1491void AliAnalysisTaskHFECal::FindTriggerClusters()
1492{
bd6ee6dd 1493 //cout << "finding trigger patch" << endl;
fb87d707 1494 // constants
1495 const int nModuleCols = 2;
1496 const int nModuleRows = 5;
1497 const int nColsFeeModule = 48;
1498 const int nRowsFeeModule = 24;
1499 const int nColsFaltroModule = 24;
1500 const int nRowsFaltroModule = 12;
1501 //const int faltroWidthMax = 20;
1502
1503 // part 1, trigger extraction -------------------------------------
1504 Int_t globCol, globRow;
1505 //Int_t ntimes=0, nTrigChannel=0, nTrigChannelCut=0, trigInCut=0;
1506 Int_t ntimes=0, nTrigChannel=0, nTrigChannelCut=0;
1507
1508 //Int_t trigtimes[faltroWidthMax];
1509 Double_t cellTime[nColsFeeModule*nModuleCols][nRowsFeeModule*nModuleRows];
1510 Double_t cellEnergy[nColsFeeModule*nModuleCols][nRowsFeeModule*nModuleRows];
1511 //Double_t fTrigCutLow = 6;
1512 //Double_t fTrigCutHigh = 10;
1513 Double_t fTimeCutLow = 469e-09;
1514 Double_t fTimeCutHigh = 715e-09;
1515
1516 AliESDCaloTrigger * fCaloTrigger = fESD->GetCaloTrigger( "EMCAL" );
1517
1518 // erase trigger maps
1519 for(Int_t i = 0; i < nColsFaltroModule*nModuleCols; i++ )
1520 {
1521 for(Int_t j = 0; j < nRowsFaltroModule*nModuleRows; j++ )
1522 {
1523 ftriggersCut[i][j] = 0;
1524 ftriggers[i][j] = 0;
1525 ftriggersTime[i][j] = 0;
1526 }
1527 }
1528
1529 Int_t iglobCol=0, iglobRow=0;
1530 // go through triggers
1531 if( fCaloTrigger->GetEntries() > 0 )
1532 {
1533 // needs reset
1534 fCaloTrigger->Reset();
1535 while( fCaloTrigger->Next() )
1536 {
1537 fCaloTrigger->GetPosition( globCol, globRow );
1538 fCaloTrigger->GetNL0Times( ntimes );
1539 /*
1540 // no L0s
1541 if( ntimes < 1 ) continue;
1542 // get precise timings
1543 fCaloTrigger->GetL0Times( trigtimes );
1544 trigInCut = 0;
1545 for(Int_t i = 0; i < ntimes; i++ )
1546 {
1547 // save the first trigger time in channel
1548 if( i == 0 || triggersTime[globCol][globRow] > trigtimes[i] )
1549 triggersTime[globCol][globRow] = trigtimes[i];
1550 //printf("trigger times: %d\n",trigtimes[i]);
1551 // check if in cut
1552 if(trigtimes[i] > fTrigCutLow && trigtimes[i] < fTrigCutHigh )
1553 trigInCut = 1;
1554
1555 fTrigTimes->Fill(trigtimes[i]);
1556 }
1557 */
1558
1559 //L1 analysis from AliAnalysisTaskEMCALTriggerQA
1560 Int_t bit = 0;
1561 fCaloTrigger->GetTriggerBits(bit);
bd6ee6dd 1562 //cout << "bit = " << bit << endl;
fb87d707 1563
1564 Int_t ts = 0;
1565 fCaloTrigger->GetL1TimeSum(ts);
bd6ee6dd 1566 //cout << "ts = " << ts << endl;
fb87d707 1567 if (ts > 0)ftriggers[globCol][globRow] = 1;
1568 // number of triggered channels in event
1569 nTrigChannel++;
1570 // ... inside cut
1571 if(ts>0 && (bit >> 6 & 0x1))
1572 {
1573 iglobCol = globCol;
1574 iglobRow = globRow;
1575 nTrigChannelCut++;
bd6ee6dd 1576 //cout << "ts cut = " << ts << endl;
1577 //cout << "globCol = " << globCol << endl;
1578 //cout << "globRow = " << globRow << endl;
fb87d707 1579 ftriggersCut[globCol][globRow] = 1;
1580 }
1581
1582 } // calo trigger entries
1583 } // has calo trigger entries
1584
1585 // part 2 go through the clusters here -----------------------------------
bd6ee6dd 1586 //cout << " part 2 go through the clusters here ----------------------------------- " << endl;
fb87d707 1587 Int_t nCluster=0, nCell=0, iCell=0, gCell=0;
60d77596 1588 Short_t cellAddr, nSACell;
1589 Int_t mclabel;
fb87d707 1590 //Int_t nSACell, iSACell, mclabel;
1591 Int_t iSACell;
42c75692 1592 Double_t cellAmp=0, cellTimeT=0, clusterTime=0, efrac=0;
fb87d707 1593 Int_t nSupMod, nModule, nIphi, nIeta, iphi, ieta, gphi, geta, feta, fphi;
1594
1595 TRefArray *fCaloClusters = new TRefArray();
1596 fESD->GetEMCALClusters( fCaloClusters );
1597 nCluster = fCaloClusters->GetEntries();
1598
1599
1600 // save all cells times if there are clusters
1601 if( nCluster > 0 ){
1602 // erase time map
1603 for(Int_t i = 0; i < nColsFeeModule*nModuleCols; i++ ){
1604 for(Int_t j = 0; j < nRowsFeeModule*nModuleRows; j++ ){
1605 cellTime[i][j] = 0.;
1606 cellEnergy[i][j] = 0.;
1607 }
1608 }
1609
1610 // get cells
1611 AliESDCaloCells *fCaloCells = fESD->GetEMCALCells();
1612 //AliVCaloCells fCaloCells = fESD->GetEMCALCells();
1613 nSACell = fCaloCells->GetNumberOfCells();
1614 for(iSACell = 0; iSACell < nSACell; iSACell++ ){
1615 // get the cell info *fCal
1616 fCaloCells->GetCell( iSACell, cellAddr, cellAmp, cellTimeT , mclabel, efrac);
1617
1618 // get cell position
1619 fGeom->GetCellIndex( cellAddr, nSupMod, nModule, nIphi, nIeta );
1620 fGeom->GetCellPhiEtaIndexInSModule( nSupMod,nModule, nIphi, nIeta, iphi, ieta);
1621
1622 // convert co global phi eta
1623 gphi = iphi + nRowsFeeModule*(nSupMod/2);
1624 geta = ieta + nColsFeeModule*(nSupMod%2);
1625
1626 // save cell time and energy
1627 cellTime[geta][gphi] = cellTimeT;
1628 cellEnergy[geta][gphi] = cellAmp;
1629
1630 }
1631 }
1632
1633 Int_t nClusterTrig, nClusterTrigCut;
1634 UShort_t *cellAddrs;
1635 Double_t clsE=-999, clsEta=-999, clsPhi=-999;
1636 Float_t clsPos[3] = {0.,0.,0.};
1637
1638 for(Int_t icl=0; icl<fESD->GetNumberOfCaloClusters(); icl++)
1639 {
1640 AliESDCaloCluster *cluster = fESD->GetCaloCluster(icl);
1641 if(!cluster || !cluster->IsEMCAL()) continue;
1642
1643 // get cluster cells
1644 nCell = cluster->GetNCells();
1645
1646 // get cluster energy
1647 clsE = cluster->E();
1648
1649 // get cluster position
1650 cluster->GetPosition(clsPos);
1651 TVector3 clsPosVec(clsPos[0],clsPos[1],clsPos[2]);
1652 clsEta = clsPosVec.Eta();
1653 clsPhi = clsPosVec.Phi();
1654
1655 // get the cell addresses
1656 cellAddrs = cluster->GetCellsAbsId();
1657
1658 // check if the cluster contains cell, that was marked as triggered
1659 nClusterTrig = 0;
1660 nClusterTrigCut = 0;
1661
1662 // loop the cells to check, if cluser in acceptance
1663 // any cluster with a cell outside acceptance is not considered
1664 for( iCell = 0; iCell < nCell; iCell++ )
1665 {
42c75692 1666 // check hot cell
feffe705 1667 //if(clsE>6.0)fCellCheck->Fill(clsE,cellAddrs[iCell]);
42c75692 1668
fb87d707 1669 // get cell position
1670 fGeom->GetCellIndex( cellAddrs[iCell], nSupMod, nModule, nIphi, nIeta );
1671 fGeom->GetCellPhiEtaIndexInSModule( nSupMod,nModule, nIphi, nIeta, iphi, ieta);
1672
1673 // convert co global phi eta
1674 gphi = iphi + nRowsFeeModule*(nSupMod/2);
1675 geta = ieta + nColsFeeModule*(nSupMod%2);
1676
1677 if( cellTime[geta][gphi] > 0. ){
1678 clusterTime += cellTime[geta][gphi];
1679 gCell++;
1680 }
1681
1682 // get corresponding FALTRO
1683 fphi = gphi / 2;
1684 feta = geta / 2;
1685
bd6ee6dd 1686 //cout << "fphi = " << fphi << endl;
1687 //cout << "feta = " << feta << endl;
1688
fb87d707 1689 // try to match with a triggered
1690 if( ftriggers[feta][fphi]==1)
1691 { nClusterTrig++;
1692 }
1693 if( ftriggersCut[feta][fphi]==1)
1694 { nClusterTrigCut++;
1695 }
1696
bd6ee6dd 1697 //cout << "nClusterTrigCut : " << nClusterTrigCut << endl;
1698
fb87d707 1699 } // cells
1700
1701
1702 if( gCell > 0 )
1703 clusterTime = clusterTime / (Double_t)gCell;
1704 // fix the reconstriction code time 100ns jumps
1705 if( fESD->GetBunchCrossNumber() % 4 < 2 )
1706 clusterTime -= 0.0000001;
1707
feffe705 1708 //fClsETime->Fill(clsE,clusterTime);
1709 //fClsEBftTrigCut->Fill(clsE);
fb87d707 1710
1711 if(nClusterTrig>0){
feffe705 1712 //fClsETime1->Fill(clsE,clusterTime);
fb87d707 1713 }
1714
1715 if(nClusterTrig>0){
1716 cluster->SetChi2(1);
feffe705 1717 //fClsEAftTrigCut1->Fill(clsE);
fb87d707 1718 }
1719
1720 if(nClusterTrigCut>0){
1721 cluster->SetChi2(2);
feffe705 1722 //fClsEAftTrigCut2->Fill(clsE);
fb87d707 1723 }
1724
1725 if(nClusterTrigCut>0 && ( nCell > (1 + clsE / 3)))
1726 {
1727 cluster->SetChi2(3);
feffe705 1728 //fClsEAftTrigCut3->Fill(clsE);
fb87d707 1729 }
1730
1731 if(nClusterTrigCut>0 && (nCell > (1 + clsE / 3) )&&( clusterTime > fTimeCutLow && clusterTime < fTimeCutHigh ))
1732 {
1733 // cluster->SetChi2(4);
feffe705 1734 //fClsEAftTrigCut4->Fill(clsE);
fb87d707 1735 }
1736 if(nClusterTrigCut<1)
1737 {
1738 cluster->SetChi2(0);
1739
feffe705 1740 //fClsEAftTrigCut->Fill(clsE);
fb87d707 1741 }
1742
1743 } // clusters
1744}
1745
1746
8efacc22 1747double AliAnalysisTaskHFECal::MCEopMeanCorrection(double pTmc, float central)
8af58b0d 1748{
8efacc22 1749 TF1 *fcorr0 = new TF1("fcorr0","[0]*tanh([1]+[2]*x)");
1750 TF1 *fcorr1 = new TF1("fcorr1","[0]*tanh([1]+[2]*x)");
1751
8af58b0d 1752 double shift = 0.0;
8efacc22 1753
1754 if(central>0 && central<=10)
1755 {
1756 fcorr0->SetParameters(1.045,1.288,3.18e-01); //
1757 fcorr1->SetParameters(9.91e-01,3.466,2.344);
1758 }
1759
1760 if(central>10 && central<=20)
1761 {
1762 fcorr0->SetParameters(1.029,8.254e-01,4.07e-01);
1763 fcorr1->SetParameters(0.975,2.276,1.501e-01);
1764 }
1765 else if(central>20 && central<=30)
1766 {
1767 fcorr0->SetParameters(1.01,8.795e-01,3.904e-01);
1768 fcorr1->SetParameters(9.675e-01,1.654,2.583e-01);
1769 }
1770 else if(central>30 && central<=40)
1771 {
1772 fcorr0->SetParameters(1.00,1.466,2.305e-1);
1773 fcorr1->SetParameters(9.637e-01,1.473,2.754e-01);
1774 }
1775 else if(central>40 && central<=50)
1776 {
fd2126aa 1777 fcorr0->SetParameters(1.00,1.422,1.518e-01);
8efacc22 1778 fcorr1->SetParameters(9.59e-01,1.421,2.931e-01);
1779 }
1780
1781 else if(central>50 && central<=70)
1782 {
1783 fcorr0->SetParameters(0.989,2.495,2.167);
1784 fcorr1->SetParameters(0.961e-1,1.734,1.438e-01);
1785 }
1786 else if(central>70 && central<=100)
1787 {
1788 fcorr0->SetParameters(0.981,-3.373,3.93327);
fd2126aa 1789 fcorr1->SetParameters(9.574e-01,1.698,1.58e-01);
8efacc22 1790 }
1791
1792
1793 shift = fcorr0->Eval(pTmc)-fcorr1->Eval(pTmc);
1794
8af58b0d 1795 return shift;
1796}
fb87d707 1797
1798
1799
1800