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