git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
- improve protection against data in impossible padrows (stack 2)
[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)
fb87d707 89 ,fGeom(0)
bfc7c23b 90 ,fOutputList(0)
85985bb0 91 ,fqahist(1)
bfc7c23b 92 ,fTrackCuts(0)
93 ,fCuts(0)
94 ,fIdentifiedAsOutInz(kFALSE)
95 ,fPassTheEventCut(kFALSE)
96 ,fRejectKinkMother(kFALSE)
e7d87aef 97 ,fmcData(kFALSE)
bfc7c23b 98 ,fVz(0.0)
99 ,fCFM(0)
100 ,fPID(0)
101 ,fPIDqa(0)
102 ,fOpeningAngleCut(0.1)
103 ,fInvmassCut(0.01)
104 ,fNoEvents(0)
105 ,fEMCAccE(0)
106 ,fTrkpt(0)
107 ,fTrkEovPBef(0)
108 ,fTrkEovPAft(0)
109 ,fdEdxBef(0)
110 ,fdEdxAft(0)
111 ,fIncpT(0)
fb87d707 112 ,fIncpTM20(0)
bfc7c23b 113 ,fInvmassLS(0)
114 ,fInvmassULS(0)
115 ,fOpeningAngleLS(0)
116 ,fOpeningAngleULS(0)
117 ,fPhotoElecPt(0)
118 ,fPhoElecPt(0)
fb87d707 119 ,fPhoElecPtM20(0)
f09766dd 120 ,fSameElecPt(0)
fb87d707 121 ,fSameElecPtM20(0)
bfc7c23b 122 ,fTrackPtBefTrkCuts(0)
123 ,fTrackPtAftTrkCuts(0)
124 ,fTPCnsigma(0)
125 ,fCent(0)
126 ,fEleInfo(0)
feffe705 127 /*
fb87d707 128 ,fClsEBftTrigCut(0)
129 ,fClsEAftTrigCut(0)
130 ,fClsEAftTrigCut1(0)
131 ,fClsEAftTrigCut2(0)
132 ,fClsEAftTrigCut3(0)
133 ,fClsEAftTrigCut4(0)
134 ,fClsETime(0)
135 ,fClsETime1(0)
136 ,fTrigTimes(0)
42c75692 137 ,fCellCheck(0)
feffe705 138 */
e7d87aef 139 ,fInputHFEMC(0)
feffe705 140 ,fInputAlle(0)
e7d87aef 141 ,fIncpTMChfe(0)
3db00c72 142 ,fIncpTMChfeAll(0)
e7d87aef 143 ,fIncpTMCM20hfe(0)
3db00c72 144 ,fIncpTMCM20hfeAll(0)
e7d87aef 145 ,fIncpTMCpho(0)
146 ,fIncpTMCM20pho(0)
147 ,fPhoElecPtMC(0)
148 ,fPhoElecPtMCM20(0)
149 ,fSameElecPtMC(0)
150 ,fSameElecPtMCM20(0)
e4b0faf2 151 ,fIncpTMCM20pho_pi0e(0)
152 ,fPhoElecPtMCM20_pi0e(0)
153 ,fSameElecPtMCM20_pi0e(0)
a6df418c 154 ,CheckNclust(0)
155 ,CheckNits(0)
d113d7cd 156 ,Hpi0pTcheck(0)
e4b0faf2 157 ,HphopTcheck(0)
50919258 158 ,fMomDtoE(0)
bfc7c23b 159{
160 //Named constructor
161
162 fPID = new AliHFEpid("hfePid");
163 fTrackCuts = new AliESDtrackCuts();
164
165 // Define input and output slots here
166 // Input slot #0 works with a TChain
167 DefineInput(0, TChain::Class());
168 // Output slot #0 id reserved by the base class for AOD
169 // Output slot #1 writes into a TH1 container
170 // DefineOutput(1, TH1I::Class());
171 DefineOutput(1, TList::Class());
172 // DefineOutput(3, TTree::Class());
173}
174
175//________________________________________________________________________
176AliAnalysisTaskHFECal::AliAnalysisTaskHFECal()
177 : AliAnalysisTaskSE("DefaultAnalysis_AliAnalysisTaskHFECal")
178 ,fESD(0)
e7d87aef 179 ,fMC(0)
fb87d707 180 ,fGeom(0)
bfc7c23b 181 ,fOutputList(0)
85985bb0 182 ,fqahist(1)
bfc7c23b 183 ,fTrackCuts(0)
184 ,fCuts(0)
185 ,fIdentifiedAsOutInz(kFALSE)
186 ,fPassTheEventCut(kFALSE)
187 ,fRejectKinkMother(kFALSE)
e7d87aef 188 ,fmcData(kFALSE)
bfc7c23b 189 ,fVz(0.0)
190 ,fCFM(0)
191 ,fPID(0)
192 ,fPIDqa(0)
193 ,fOpeningAngleCut(0.1)
194 ,fInvmassCut(0.01)
195 ,fNoEvents(0)
196 ,fEMCAccE(0)
197 ,fTrkpt(0)
198 ,fTrkEovPBef(0)
199 ,fTrkEovPAft(0)
200 ,fdEdxBef(0)
201 ,fdEdxAft(0)
202 ,fIncpT(0)
fb87d707 203 ,fIncpTM20(0)
bfc7c23b 204 ,fInvmassLS(0)
205 ,fInvmassULS(0)
206 ,fOpeningAngleLS(0)
207 ,fOpeningAngleULS(0)
208 ,fPhotoElecPt(0)
209 ,fPhoElecPt(0)
fb87d707 210 ,fPhoElecPtM20(0)
f09766dd 211 ,fSameElecPt(0)
fb87d707 212 ,fSameElecPtM20(0)
bfc7c23b 213 ,fTrackPtBefTrkCuts(0)
214 ,fTrackPtAftTrkCuts(0)
215 ,fTPCnsigma(0)
216 ,fCent(0)
217 ,fEleInfo(0)
feffe705 218 /*
fb87d707 219 ,fClsEBftTrigCut(0)
220 ,fClsEAftTrigCut(0)
221 ,fClsEAftTrigCut1(0)
222 ,fClsEAftTrigCut2(0)
223 ,fClsEAftTrigCut3(0)
224 ,fClsEAftTrigCut4(0)
225 ,fClsETime(0)
226 ,fClsETime1(0)
227 ,fTrigTimes(0)
42c75692 228 ,fCellCheck(0)
feffe705 229 */
e7d87aef 230 ,fInputHFEMC(0)
feffe705 231 ,fInputAlle(0)
e7d87aef 232 ,fIncpTMChfe(0)
3db00c72 233 ,fIncpTMChfeAll(0)
e7d87aef 234 ,fIncpTMCM20hfe(0)
3db00c72 235 ,fIncpTMCM20hfeAll(0)
e7d87aef 236 ,fIncpTMCpho(0)
237 ,fIncpTMCM20pho(0)
238 ,fPhoElecPtMC(0)
239 ,fPhoElecPtMCM20(0)
240 ,fSameElecPtMC(0)
241 ,fSameElecPtMCM20(0)
e4b0faf2 242 ,fIncpTMCM20pho_pi0e(0)
243 ,fPhoElecPtMCM20_pi0e(0)
244 ,fSameElecPtMCM20_pi0e(0)
a6df418c 245 ,CheckNclust(0)
246 ,CheckNits(0)
d113d7cd 247 ,Hpi0pTcheck(0)
e4b0faf2 248 ,HphopTcheck(0)
50919258 249 ,fMomDtoE(0)
bfc7c23b 250{
251 //Default constructor
252 fPID = new AliHFEpid("hfePid");
253
254 fTrackCuts = new AliESDtrackCuts();
255
256 // Constructor
257 // Define input and output slots here
258 // Input slot #0 works with a TChain
259 DefineInput(0, TChain::Class());
260 // Output slot #0 id reserved by the base class for AOD
261 // Output slot #1 writes into a TH1 container
262 // DefineOutput(1, TH1I::Class());
263 DefineOutput(1, TList::Class());
264 //DefineOutput(3, TTree::Class());
265}
266//_________________________________________
267
268AliAnalysisTaskHFECal::~AliAnalysisTaskHFECal()
269{
270 //Destructor
271
272 delete fOutputList;
fb87d707 273 delete fGeom;
bfc7c23b 274 delete fPID;
a6df418c 275 delete fCuts;
bfc7c23b 276 delete fCFM;
277 delete fPIDqa;
278 delete fTrackCuts;
279}
280//_________________________________________
281
282void AliAnalysisTaskHFECal::UserExec(Option_t*)
283{
284 //Main loop
285 //Called for each event
286
287 // create pointer to event
288 fESD = dynamic_cast<AliESDEvent*>(InputEvent());
289 if (!fESD) {
290 printf("ERROR: fESD not available\n");
291 return;
292 }
293
294 if(!fCuts){
295 AliError("HFE cuts not available");
296 return;
297 }
298
299 if(!fPID->IsInitialized()){
300 // Initialize PID with the given run number
301 AliWarning("PID not initialised, get from Run no");
302 fPID->InitializePID(fESD->GetRunNumber());
303 }
304
e7d87aef 305 if(fmcData)fMC = MCEvent();
306 AliStack* stack = NULL;
307 if(fmcData && fMC)stack = fMC->Stack();
308
309 Float_t cent = -1.;
310 AliCentrality *centrality = fESD->GetCentrality();
311 cent = centrality->GetCentralityPercentile("V0M");
312
313 //---- fill MC track info
314 if(fmcData && fMC)
315 {
316 Int_t nParticles = stack->GetNtrack();
317 for (Int_t iParticle = 0; iParticle < nParticles; iParticle++) {
318 TParticle* particle = stack->Particle(iParticle);
319 int fPDG = particle->GetPdgCode();
320 double mcZvertex = fMC->GetPrimaryVertex()->GetZ();
321 double pTMC = particle->Pt();
feffe705 322 double proR = particle->R();
9b3495ff 323 double etaMC = particle->Eta();
324 if(fabs(etaMC)>0.6)continue;
e7d87aef 325 Bool_t mcInDtoE= kFALSE;
326 Bool_t mcInBtoE= kFALSE;
327
d113d7cd 328 Bool_t MChijing = fMC->IsFromBGEvent(iParticle);
fd426741 329 //if(!MChijing)printf("not MC hijing");
d113d7cd 330 int iHijing = 1;
331 if(!MChijing)iHijing = 0;
332 if(fPDG==111)Hpi0pTcheck->Fill(pTMC,iHijing);
333
e4b0faf2 334 if(particle->GetFirstMother()>-1 && fPDG==22)
335 {
336 int parentPID = stack->Particle(particle->GetFirstMother())->GetPdgCode();
337 if(parentPID==111 || parentPID==221)HphopTcheck->Fill(pTMC,iHijing); // pi0->g & eta->g
338 }
339
e7d87aef 340 if(particle->GetFirstMother()>-1 && fabs(fPDG)==11)
341 {
342 int parentPID = stack->Particle(particle->GetFirstMother())->GetPdgCode();
343 if((fabs(parentPID)==411 || fabs(parentPID)==413 || fabs(parentPID)==421 || fabs(parentPID)==423 || fabs(parentPID)==431)&& fabs(fPDG)==11)mcInDtoE = kTRUE;
344 if((fabs(parentPID)==511 || fabs(parentPID)==513 || fabs(parentPID)==521 || fabs(parentPID)==523 || fabs(parentPID)==531)&& fabs(fPDG)==11)mcInBtoE = kTRUE;
345 if((mcInBtoE || mcInDtoE) && fabs(mcZvertex)<10.0)fInputHFEMC->Fill(cent,pTMC);
346 }
347
feffe705 348 if(proR<7.0 && fabs(fPDG)==11)fInputAlle->Fill(cent,pTMC);
349
e7d87aef 350 }
351 }
352
bfc7c23b 353 fNoEvents->Fill(0);
354
355 Int_t fNOtrks = fESD->GetNumberOfTracks();
356 const AliESDVertex *pVtx = fESD->GetPrimaryVertex();
357
358 Double_t pVtxZ = -999;
359 pVtxZ = pVtx->GetZ();
360
361 if(TMath::Abs(pVtxZ)>10) return;
362 fNoEvents->Fill(1);
363
364 if(fNOtrks<2) return;
365 fNoEvents->Fill(2);
366
367 AliPIDResponse *pidResponse = fInputHandler->GetPIDResponse();
368 if(!pidResponse){
369 AliDebug(1, "Using default PID Response");
370 pidResponse = AliHFEtools::GetDefaultPID(kFALSE, fInputEvent->IsA() == AliAODEvent::Class());
371 }
372
373 fPID->SetPIDResponse(pidResponse);
374
375 fCFM->SetRecEventInfo(fESD);
376
e7d87aef 377 //Float_t cent = -1.;
378 //AliCentrality *centrality = fESD->GetCentrality();
379 //cent = centrality->GetCentralityPercentile("V0M");
bfc7c23b 380 fCent->Fill(cent);
381
b12bc641 382 //if(cent>90.) return;
bfc7c23b 383
384 // Calorimeter info.
385
42c75692 386 FindTriggerClusters();
387
bfc7c23b 388 // make EMCAL array
389 for(Int_t iCluster=0; iCluster<fESD->GetNumberOfCaloClusters(); iCluster++)
390 {
391 AliESDCaloCluster *clust = fESD->GetCaloCluster(iCluster);
392 if(clust->IsEMCAL())
393 {
394 double clustE = clust->E();
395 float emcx[3]; // cluster pos
396 clust->GetPosition(emcx);
397 TVector3 clustpos(emcx[0],emcx[1],emcx[2]);
398 double emcphi = clustpos.Phi();
399 double emceta = clustpos.Eta();
fb87d707 400 double calInfo[5];
401 calInfo[0] = emcphi; calInfo[1] = emceta; calInfo[2] = clustE; calInfo[3] = cent; calInfo[4] = clust->Chi2();
85985bb0 402 //if(clustE>1.5)fEMCAccE->Fill(calInfo);
44be9e1d 403 //if(fqahist==1 && clustE>1.5)fEMCAccE->Fill(calInfo);
bfc7c23b 404 }
405 }
406
407 // Track loop
408 for (Int_t iTracks = 0; iTracks < fESD->GetNumberOfTracks(); iTracks++) {
409 AliESDtrack* track = fESD->GetTrack(iTracks);
410 if (!track) {
411 printf("ERROR: Could not receive track %d\n", iTracks);
412 continue;
413 }
e7d87aef 414
415 Bool_t mcPho = kFALSE;
416 Bool_t mcDtoE= kFALSE;
417 Bool_t mcBtoE= kFALSE;
418 double mcele = -1.;
44be9e1d 419 double mcpT = 0.0;
b12bc641 420 double mcMompT = 0.0;
e4b0faf2 421 //double mcGrandMompT = 0.0;
422 double mcWeight = -10.0;
d113d7cd 423
424 int iHijing = 1;
425
e7d87aef 426 if(fmcData && fMC && stack)
427 {
feffe705 428 Int_t label = TMath::Abs(track->GetLabel());
d113d7cd 429
430 Bool_t MChijing = fMC->IsFromBGEvent(label);
431 if(!MChijing)iHijing = 0;
432
e7d87aef 433 TParticle* particle = stack->Particle(label);
434 int mcpid = particle->GetPdgCode();
44be9e1d 435 mcpT = particle->Pt();
9b3495ff 436 //printf("MCpid = %d",mcpid);
e7d87aef 437 if(particle->GetFirstMother()>-1)
438 {
e4b0faf2 439 int parentlabel = particle->GetFirstMother();
e7d87aef 440 int parentPID = stack->Particle(particle->GetFirstMother())->GetPdgCode();
b12bc641 441 mcMompT = stack->Particle(particle->GetFirstMother())->Pt();
e7d87aef 442 if((parentPID==22 || parentPID==111 || parentPID==221)&& fabs(mcpid)==11)mcPho = kTRUE;
443 if((fabs(parentPID)==411 || fabs(parentPID)==413 || fabs(parentPID)==421 || fabs(parentPID)==423 || fabs(parentPID)==431)&& fabs(mcpid)==11)mcDtoE = kTRUE;
444 if((fabs(parentPID)==511 || fabs(parentPID)==513 || fabs(parentPID)==521 || fabs(parentPID)==523 || fabs(parentPID)==531)&& fabs(mcpid)==11)mcBtoE = kTRUE;
e4b0faf2 445
50919258 446 // make D->e pT correlation
447 if(mcDtoE)fMomDtoE->Fill(mcpT,mcMompT);
448
449 // pi->e (Dalitz)
450 if(parentPID==111 && fabs(mcpid)==11 && mcMompT>0.0)
e4b0faf2 451 {
452 if(mcMompT>0.0 && mcMompT<5.0)
453 {
454 mcWeight = 0.323*mcMompT/(TMath::Exp(-1.6+0.767*mcMompT+0.0285*mcMompT*mcMompT));
455 }
456 else
457 {
458 mcWeight = 115.0/(0.718*mcMompT*TMath::Power(mcMompT,3.65));
459 }
460 }
461
462 // access grand parent pi0->g->e
463 TParticle* particle_parent = stack->Particle(parentlabel); // get parent pointer
464 if(particle_parent->GetFirstMother()>-1 && parentPID==22 && fabs(mcpid)==11) // get grand parent g->e
465 {
466 int grand_parentPID = stack->Particle(particle_parent->GetFirstMother())->GetPdgCode();
467 if(grand_parentPID==111)
468 {
469 //mcGrandMompT = stack->Particle(particle_parent->GetFirstMother())->Pt();
470 double pTtmp = stack->Particle(particle_parent->GetFirstMother())->Pt();
471 //if(mcGrandMompT>0.0 && mcGrandMompT<5.0)
50919258 472 if(pTtmp>0.0)
473 {
474 if(pTtmp<5.0)
475 {
476 mcWeight = 0.323*pTtmp/(TMath::Exp(-1.6+0.767*pTtmp+0.0285*pTtmp*pTtmp));
477 }
478 else
e4b0faf2 479 {
480 mcWeight = 115.0/(0.718*pTtmp*TMath::Power(pTtmp,3.65));
481 }
50919258 482 }
e4b0faf2 483 }
484 }
e7d87aef 485 }
9b3495ff 486 if(fabs(mcpid)==11 && mcDtoE)mcele= 1.;
487 if(fabs(mcpid)==11 && mcBtoE)mcele= 2.;
488 if(fabs(mcpid)==11 && mcPho)mcele= 3.;
e7d87aef 489 }
490
50919258 491 //printf("weight = %f\n",mcWeight);
492
9b3495ff 493 if(TMath::Abs(track->Eta())>0.6) continue;
f09766dd 494 if(TMath::Abs(track->Pt()<2.0)) continue;
bfc7c23b 495
496 fTrackPtBefTrkCuts->Fill(track->Pt());
497 // RecKine: ITSTPC cuts
498 if(!ProcessCutStep(AliHFEcuts::kStepRecKineITSTPC, track)) continue;
499
500 //RecKink
501 if(fRejectKinkMother) { // Quick and dirty fix to reject both kink mothers and daughters
502 if(track->GetKinkIndex(0) != 0) continue;
503 }
504
505 // RecPrim
506 if(!ProcessCutStep(AliHFEcuts::kStepRecPrim, track)) continue;
507
508 // HFEcuts: ITS layers cuts
509 if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsITS, track)) continue;
510
f09766dd 511 // HFE cuts: TPC PID cleanup
512 if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTPC, track)) continue;
513
a6df418c 514 int nTPCcl = track->GetTPCNcls();
6e92c428 515 int nTPCclF = track->GetTPCNclsF();
a6df418c 516 int nITS = track->GetNcls(0);
bfc7c23b 517
518 fTrackPtAftTrkCuts->Fill(track->Pt());
519
520 Double_t mom = -999., eop=-999., pt = -999., dEdx=-999., fTPCnSigma=-10, phi=-999., eta=-999.;
feffe705 521 pt = track->Pt();
522 if(pt<2.0)continue;
bfc7c23b 523
524 // Track extrapolation
525
bfc7c23b 526 Int_t charge = track->Charge();
527 fTrkpt->Fill(pt);
528 mom = track->P();
529 phi = track->Phi();
530 eta = track->Eta();
531 dEdx = track->GetTPCsignal();
532 fTPCnSigma = fPID->GetPIDResponse() ? fPID->GetPIDResponse()->NumberOfSigmasTPC(track, AliPID::kElectron) : 1000;
533
534 double ncells = -1.0;
535 double m20 = -1.0;
536 double m02 = -1.0;
537 double disp = -1.0;
538 double rmatch = -1.0;
539 double nmatch = -1.0;
f09766dd 540 double oppstatus = 0.0;
bfc7c23b 541
f09766dd 542 Bool_t fFlagPhotonicElec = kFALSE;
543 Bool_t fFlagConvinatElec = kFALSE;
bfc7c23b 544
545 Int_t clsId = track->GetEMCALcluster();
546 if (clsId>0){
547 AliESDCaloCluster *clust = fESD->GetCaloCluster(clsId);
548 if(clust && clust->IsEMCAL()){
549
550 double clustE = clust->E();
551 eop = clustE/fabs(mom);
552 //double clustT = clust->GetTOF();
553 ncells = clust->GetNCells();
554 m02 = clust->GetM02();
555 m20 = clust->GetM20();
556 disp = clust->GetDispersion();
557 double delphi = clust->GetTrackDx();
558 double deleta = clust->GetTrackDz();
559 rmatch = sqrt(pow(delphi,2)+pow(deleta,2));
560 nmatch = clust->GetNTracksMatched();
561
c2801a73 562 if(fTPCnSigma>-1.5 && fTPCnSigma<3.0)
563 {
564 SelectPhotonicElectron(iTracks,cent,track,fFlagPhotonicElec,fFlagConvinatElec,fTPCnSigma,m20,eop,mcele);
565 }
566 if(fFlagPhotonicElec)oppstatus = 1.0;
567 if(fFlagConvinatElec)oppstatus = 2.0;
568 if(fFlagPhotonicElec && fFlagConvinatElec)oppstatus = 3.0;
569
feffe705 570 double valdedx[16];
a6df418c 571 valdedx[0] = pt; valdedx[1] = nITS; valdedx[2] = phi; valdedx[3] = eta; valdedx[4] = fTPCnSigma;
6e92c428 572 valdedx[5] = eop; valdedx[6] = rmatch; valdedx[7] = ncells, valdedx[8] = nTPCclF; valdedx[9] = m20; valdedx[10] = mcpT;
a6df418c 573 valdedx[11] = cent; valdedx[12] = charge; valdedx[13] = oppstatus; valdedx[14] = nTPCcl;
feffe705 574 valdedx[15] = mcele;
85985bb0 575 if(fqahist==1)fEleInfo->Fill(valdedx);
bfc7c23b 576
577
578 }
579 }
a6df418c 580
581 if(nITS<2.5)continue;
582 if(nTPCcl<100)continue;
583
584 CheckNclust->Fill(nTPCcl);
585 CheckNits->Fill(nITS);
586
42c75692 587 fdEdxBef->Fill(mom,fTPCnSigma);
bfc7c23b 588 fTPCnsigma->Fill(mom,fTPCnSigma);
f09766dd 589 if(fTPCnSigma >= -1.0 && fTPCnSigma <= 3)fTrkEovPBef->Fill(pt,eop);
bfc7c23b 590
bfc7c23b 591 Int_t pidpassed = 1;
592
593
594 //--- track accepted
595 AliHFEpidObject hfetrack;
596 hfetrack.SetAnalysisType(AliHFEpidObject::kESDanalysis);
597 hfetrack.SetRecTrack(track);
3db00c72 598 double binct = 10.5;
599 if((0.0< cent) && (cent<5.0)) binct = 0.5;
600 if((5.0< cent) && (cent<10.0)) binct = 1.5;
601 if((10.0< cent) && (cent<20.0)) binct = 2.5;
602 if((20.0< cent) && (cent<30.0)) binct = 3.5;
603 if((30.0< cent) && (cent<40.0)) binct = 4.5;
604 if((40.0< cent) && (cent<50.0)) binct = 5.5;
605 if((50.0< cent) && (cent<60.0)) binct = 6.5;
606 if((60.0< cent) && (cent<70.0)) binct = 7.5;
607 if((70.0< cent) && (cent<80.0)) binct = 8.5;
608 if((80.0< cent) && (cent<90.0)) binct = 9.5;
609 if((90.0< cent) && (cent<100.0)) binct = 10.5;
610
611 hfetrack.SetCentrality((int)binct); //added
bfc7c23b 612 hfetrack.SetPbPb();
613 if(!fPID->IsSelected(&hfetrack, NULL, "", fPIDqa)) pidpassed = 0;
614
615 if(pidpassed==0) continue;
616
617 fTrkEovPAft->Fill(pt,eop);
42c75692 618 fdEdxAft->Fill(mom,fTPCnSigma);
bfc7c23b 619
e7d87aef 620 // Fill real data
fb87d707 621 fIncpT->Fill(cent,pt);
bfc7c23b 622 if(fFlagPhotonicElec) fPhoElecPt->Fill(cent,pt);
f09766dd 623 if(fFlagConvinatElec) fSameElecPt->Fill(cent,pt);
fb87d707 624
625 if(m20>0.0 && m20<0.3)
626 {
627 fIncpTM20->Fill(cent,pt);
628 if(fFlagPhotonicElec) fPhoElecPtM20->Fill(cent,pt);
629 if(fFlagConvinatElec) fSameElecPtM20->Fill(cent,pt);
630 }
feffe705 631
e7d87aef 632 // MC
e4b0faf2 633 if(mcele>0) // select MC electrons
e7d87aef 634 {
3db00c72 635
636 fIncpTMChfeAll->Fill(cent,pt);
637 if(m20>0.0 && m20<0.3)fIncpTMCM20hfeAll->Fill(cent,pt);
638
e4b0faf2 639 if(mcBtoE || mcDtoE) // select B->e & D->e
e7d87aef 640 {
641 fIncpTMChfe->Fill(cent,pt);
642 if(m20>0.0 && m20<0.3)fIncpTMCM20hfe->Fill(cent,pt);
643 }
e4b0faf2 644
645 if(mcPho) // select photonic electrons
e7d87aef 646 {
987053ce 647 double phoval[5];
44be9e1d 648 phoval[0] = cent;
649 phoval[1] = pt;
650 phoval[2] = fTPCnSigma;
d113d7cd 651 phoval[3] = iHijing;
987053ce 652 phoval[4] = mcMompT;
44be9e1d 653
654 fIncpTMCpho->Fill(phoval);
655 if(fFlagPhotonicElec) fPhoElecPtMC->Fill(phoval);
656 if(fFlagConvinatElec) fSameElecPtMC->Fill(phoval);
e7d87aef 657
658 if(m20>0.0 && m20<0.3)
659 {
44be9e1d 660 fIncpTMCM20pho->Fill(phoval);
661 if(fFlagPhotonicElec) fPhoElecPtMCM20->Fill(phoval);
662 if(fFlagConvinatElec) fSameElecPtMCM20->Fill(phoval);
e4b0faf2 663 // pi0->g->e
664 if(mcWeight>-1)
665 {
987053ce 666 if(iHijing==1)mcWeight = 1.0;
e4b0faf2 667 fIncpTMCM20pho_pi0e->Fill(phoval,mcWeight);
668 if(fFlagPhotonicElec) fPhoElecPtMCM20_pi0e->Fill(phoval,mcWeight);
669 if(fFlagConvinatElec) fSameElecPtMCM20_pi0e->Fill(phoval,mcWeight);
670 }
e7d87aef 671 }
672 }
673 }
bfc7c23b 674 }
675 PostData(1, fOutputList);
676}
677//_________________________________________
678void AliAnalysisTaskHFECal::UserCreateOutputObjects()
679{
42c75692 680 //--- Check MC
681
e7d87aef 682 //Bool_t mcData = kFALSE;
42c75692 683 if(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler())
684 {
e7d87aef 685 fmcData = kTRUE;
42c75692 686 printf("+++++ MC Data available");
687 }
e7d87aef 688 if(fmcData)
42c75692 689 {
690 printf("++++++++= MC analysis \n");
691 }
692 else
693 {
694 printf("++++++++ real data analysis \n");
695 }
696
85985bb0 697 printf("+++++++ QA hist %d",fqahist);
698
fb87d707 699 //---- Geometry
700 fGeom = AliEMCALGeometry::GetInstance("EMCAL_COMPLETEV1");
701
bfc7c23b 702 //--------Initialize PID
42c75692 703 //fPID->SetHasMCData(kFALSE);
e7d87aef 704 fPID->SetHasMCData(fmcData);
bfc7c23b 705 if(!fPID->GetNumberOfPIDdetectors())
706 {
707 fPID->AddDetector("TPC", 0);
708 fPID->AddDetector("EMCAL", 1);
709 }
710
3db00c72 711 Double_t params[4];
78ff1095 712 const char *cutmodel;
3db00c72 713 cutmodel = "pol0";
714 params[0] = -1.0; //sigma min
9b3495ff 715 double maxnSig = 3.0;
716 if(fmcData)
717 {
718 params[0] = -5.0; //sigma min
719 maxnSig = 5.0;
720 }
721
722 for(Int_t a=0;a<11;a++)fPID->ConfigureTPCcentralityCut(a,cutmodel,params,maxnSig);
3db00c72 723
bfc7c23b 724 fPID->SortDetectors();
725 fPIDqa = new AliHFEpidQAmanager();
726 fPIDqa->Initialize(fPID);
a6df418c 727
728 //------- fcut --------------
729 fCuts = new AliHFEcuts();
730 fCuts->CreateStandardCuts();
731 //fCuts->SetMinNClustersTPC(100);
732 fCuts->SetMinNClustersTPC(90);
733 fCuts->SetMinRatioTPCclusters(0.6);
734 fCuts->SetTPCmodes(AliHFEextraCuts::kFound, AliHFEextraCuts::kFoundOverFindable);
735 //fCuts->SetMinNClustersITS(3);
736 fCuts->SetMinNClustersITS(2);
737 fCuts->SetCutITSpixel(AliHFEextraCuts::kAny);
738 fCuts->SetCheckITSLayerStatus(kFALSE);
739 fCuts->SetVertexRange(10.);
740 fCuts->SetTOFPIDStep(kFALSE);
741 fCuts->SetPtRange(2, 50);
742 fCuts->SetMaxImpactParam(3.,3.);
743
bfc7c23b 744 //--------Initialize correction Framework and Cuts
745 fCFM = new AliCFManager;
746 const Int_t kNcutSteps = AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kNcutStepsRecTrack + AliHFEcuts::kNcutStepsDETrack;
747 fCFM->SetNStepParticle(kNcutSteps);
748 for(Int_t istep = 0; istep < kNcutSteps; istep++)
749 fCFM->SetParticleCutsList(istep, NULL);
750
751 if(!fCuts){
752 AliWarning("Cuts not available. Default cuts will be used");
753 fCuts = new AliHFEcuts;
754 fCuts->CreateStandardCuts();
755 }
756 fCuts->Initialize(fCFM);
757
758 //---------Output Tlist
759 fOutputList = new TList();
760 fOutputList->SetOwner();
761 fOutputList->Add(fPIDqa->MakeList("PIDQA"));
762
763 fNoEvents = new TH1F("fNoEvents","",4,-0.5,3.5) ;
764 fOutputList->Add(fNoEvents);
765
02cd82a0 766 Int_t binsE[5] = {250, 100, 1000, 200, 10};
fb87d707 767 Double_t xminE[5] = {1.0, -1, 0.0, 0, -0.5};
768 Double_t xmaxE[5] = {3.5, 1, 100.0, 100, 9.5};
769 fEMCAccE = new THnSparseD("fEMCAccE","EMC acceptance & E;#eta;#phi;Energy;Centrality;trugCondition;",5,binsE,xminE,xmaxE);
85985bb0 770 if(fqahist==1)fOutputList->Add(fEMCAccE);
bfc7c23b 771
772 fTrkpt = new TH1F("fTrkpt","track pt",100,0,50);
773 fOutputList->Add(fTrkpt);
774
775 fTrackPtBefTrkCuts = new TH1F("fTrackPtBefTrkCuts","track pt before track cuts",100,0,50);
776 fOutputList->Add(fTrackPtBefTrkCuts);
777
778 fTrackPtAftTrkCuts = new TH1F("fTrackPtAftTrkCuts","track pt after track cuts",100,0,50);
779 fOutputList->Add(fTrackPtAftTrkCuts);
780
781 fTPCnsigma = new TH2F("fTPCnsigma", "TPC - n sigma",100,0,50,200,-10,10);
782 fOutputList->Add(fTPCnsigma);
783
784 fTrkEovPBef = new TH2F("fTrkEovPBef","track E/p before HFE pid",100,0,50,100,0,2);
785 fOutputList->Add(fTrkEovPBef);
786
787 fTrkEovPAft = new TH2F("fTrkEovPAft","track E/p after HFE pid",100,0,50,100,0,2);
788 fOutputList->Add(fTrkEovPAft);
789
42c75692 790 fdEdxBef = new TH2F("fdEdxBef","track dEdx vs p before HFE pid",100,0,50,200,-10,10);
bfc7c23b 791 fOutputList->Add(fdEdxBef);
792
42c75692 793 fdEdxAft = new TH2F("fdEdxAft","track dEdx vs p after HFE pid",100,0,50,200,-10,10);
bfc7c23b 794 fOutputList->Add(fdEdxAft);
795
02cd82a0 796 fIncpT = new TH2F("fIncpT","HFE pid electro vs. centrality",200,0,100,100,0,50);
bfc7c23b 797 fOutputList->Add(fIncpT);
798
02cd82a0 799 fIncpTM20 = new TH2F("fIncpTM20","HFE pid electro vs. centrality with M20",200,0,100,100,0,50);
fb87d707 800 fOutputList->Add(fIncpTM20);
4e01b68c 801
02cd82a0 802 Int_t nBinspho[9] = { 200, 100, 500, 12, 50, 4, 200, 8, 100};
b12bc641 803 Double_t minpho[9] = { 0., 0., 0., -2.5, 0, -0.5, 0,-1.5, 0};
804 Double_t maxpho[9] = {100., 50., 0.5, 3.5, 1, 3.5, 2, 6.5, 50};
f09766dd 805
b12bc641 806 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);
bfc7c23b 807 fOutputList->Add(fInvmassLS);
808
b12bc641 809 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);
bfc7c23b 810 fOutputList->Add(fInvmassULS);
811
812 fOpeningAngleLS = new TH1F("fOpeningAngleLS","Opening angle for LS pairs",100,0,1);
813 fOutputList->Add(fOpeningAngleLS);
814
815 fOpeningAngleULS = new TH1F("fOpeningAngleULS","Opening angle for ULS pairs",100,0,1);
816 fOutputList->Add(fOpeningAngleULS);
817
818 fPhotoElecPt = new TH1F("fPhotoElecPt", "photonic electron pt",100,0,50);
819 fOutputList->Add(fPhotoElecPt);
820
02cd82a0 821 fPhoElecPt = new TH2F("fPhoElecPt", "Pho-inclusive electron pt",200,0,100,100,0,50);
bfc7c23b 822 fOutputList->Add(fPhoElecPt);
823
02cd82a0 824 fPhoElecPtM20 = new TH2F("fPhoElecPtM20", "Pho-inclusive electron pt with M20",200,0,100,100,0,50);
fb87d707 825 fOutputList->Add(fPhoElecPtM20);
826
02cd82a0 827 fSameElecPt = new TH2F("fSameElecPt", "Same-inclusive electron pt",200,0,100,100,0,50);
f09766dd 828 fOutputList->Add(fSameElecPt);
829
02cd82a0 830 fSameElecPtM20 = new TH2F("fSameElecPtM20", "Same-inclusive electron pt with M20",200,0,100,100,0,50);
fb87d707 831 fOutputList->Add(fSameElecPtM20);
832
02cd82a0 833 fCent = new TH1F("fCent","Centrality",200,0,100) ;
bfc7c23b 834 fOutputList->Add(fCent);
835
836 // Make common binning
44be9e1d 837 const Double_t kMinP = 0.;
bfc7c23b 838 const Double_t kMaxP = 50.;
a6df418c 839 //const Double_t kTPCSigMim = 40.;
840 //const Double_t kTPCSigMax = 140.;
bfc7c23b 841
842 // 1st histogram: TPC dEdx with/without EMCAL (p, pT, TPC Signal, phi, eta, Sig, e/p, ,match, cell, M02, M20, Disp, Centrality, select)
a6df418c 843 Int_t nBins[16] = { 250, 10, 60, 20, 100, 300, 50, 40, 200, 200, 250, 200, 3, 5, 100, 8};
844 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 845 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};
846 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 847 if(fqahist==1)fOutputList->Add(fEleInfo);
bfc7c23b 848
fb87d707 849 //<--- Trigger info
feffe705 850 /*
fb87d707 851 fClsEBftTrigCut = new TH1F("fClsEBftTrigCut","cluster E before trigger selection",1000,0,100);
852 fOutputList->Add(fClsEBftTrigCut);
853
854 fClsEAftTrigCut = new TH1F("fClsEAftTrigCut","cluster E if cls has 0 trigcut channel",1000,0,100);
855 fOutputList->Add(fClsEAftTrigCut);
856
857 fClsEAftTrigCut1 = new TH1F("fClsEAftTrigCut1","cluster E if cls with trig channel",1000,0,100);
858 fOutputList->Add(fClsEAftTrigCut1);
859
860 fClsEAftTrigCut2 = new TH1F("fClsEAftTrigCut2","cluster E if cls with trigcut channel",1000,0,100);
861 fOutputList->Add(fClsEAftTrigCut2);
862
863 fClsEAftTrigCut3 = new TH1F("fClsEAftTrigCut3","cluster E if cls with trigcut channel + nCell>Ecorrect",1000,0,100);
864 fOutputList->Add(fClsEAftTrigCut3);
865
866 fClsEAftTrigCut4 = new TH1F("fClsEAftTrigCut4","cluster E if cls with trigcut channel + nCell>Ecorrect + cls time cut",1000,0,100);
867 fOutputList->Add(fClsEAftTrigCut4);
868
869 fClsETime = new TH2F("fClsETime", "Cls time vs E; E; time;",1000,0,100,1000,-0.0000002,0.0000009);
870 fOutputList->Add(fClsETime);
871
872 fClsETime1 = new TH2F("fClsETime1", "Cls time vs E if cls contains trigger channel; E; time;",1000,0,100,1000,-0.0000002,0.0000009);
873 fOutputList->Add(fClsETime1);
874
875 fTrigTimes = new TH1F("fTrigTimes", "Trigger time; time; N;",25,0,25);
876 fOutputList->Add(fTrigTimes);
877
42c75692 878 fCellCheck = new TH2F("fCellCheck", "Cell vs E; E GeV; Cell ID",10,6,26,12000,0,12000);
879 fOutputList->Add(fCellCheck);
feffe705 880 */
e7d87aef 881 //<---------- MC
882
02cd82a0 883 fInputHFEMC = new TH2F("fInputHFEMC","Input MC HFE pid electro vs. centrality",200,0,100,100,0,50);
e7d87aef 884 fOutputList->Add(fInputHFEMC);
885
02cd82a0 886 fInputAlle = new TH2F("fInputAlle","Input MC electro vs. centrality",200,0,100,100,0,50);
feffe705 887 fOutputList->Add(fInputAlle);
888
02cd82a0 889 fIncpTMChfe = new TH2F("fIncpTMChfe","MC HFE pid electro vs. centrality",200,0,100,100,0,50);
e7d87aef 890 fOutputList->Add(fIncpTMChfe);
891
02cd82a0 892 fIncpTMChfeAll = new TH2F("fIncpTMChfeAll","MC Alle pid electro vs. centrality",200,0,100,100,0,50);
3db00c72 893 fOutputList->Add(fIncpTMChfeAll);
894
02cd82a0 895 fIncpTMCM20hfe = new TH2F("fIncpTMCM20hfe","MC HFE pid electro vs. centrality with M20",200,0,100,100,0,50);
e7d87aef 896 fOutputList->Add(fIncpTMCM20hfe);
897
02cd82a0 898 fIncpTMCM20hfeAll = new TH2F("fIncpTMCM20hfeAll","MC Alle pid electro vs. centrality with M20",200,0,100,100,0,50);
3db00c72 899 fOutputList->Add(fIncpTMCM20hfeAll);
900
44be9e1d 901
987053ce 902 Int_t nBinspho2[5] = { 200, 100, 7, 3, 100};
903 Double_t minpho2[5] = { 0., 0., -2.5, -0.5, 0.};
904 Double_t maxpho2[5] = {100., 50., 4.5, 2.5, 50.};
44be9e1d 905
987053ce 906 fIncpTMCpho = new THnSparseD("fIncpTMCpho","MC Pho pid electro vs. centrality",5,nBinspho2,minpho2,maxpho2);
e7d87aef 907 fOutputList->Add(fIncpTMCpho);
908
987053ce 909 fIncpTMCM20pho = new THnSparseD("fIncpTMCM20pho","MC Pho pid electro vs. centrality with M20",5,nBinspho2,minpho2,maxpho2);
e7d87aef 910 fOutputList->Add(fIncpTMCM20pho);
911
987053ce 912 fPhoElecPtMC = new THnSparseD("fPhoElecPtMC", "MC Pho-inclusive electron pt",5,nBinspho2,minpho2,maxpho2);
e7d87aef 913 fOutputList->Add(fPhoElecPtMC);
914
987053ce 915 fPhoElecPtMCM20 = new THnSparseD("fPhoElecPtMCM20", "MC Pho-inclusive electron pt with M20",5,nBinspho2,minpho2,maxpho2);
feffe705 916 fOutputList->Add(fPhoElecPtMCM20);
e7d87aef 917
987053ce 918 fSameElecPtMC = new THnSparseD("fSameElecPtMC", "MC Same-inclusive electron pt",5,nBinspho2,minpho2,maxpho2);
e7d87aef 919 fOutputList->Add(fSameElecPtMC);
920
987053ce 921 fSameElecPtMCM20 = new THnSparseD("fSameElecPtMCM20", "MC Same-inclusive electron pt with M20",5,nBinspho2,minpho2,maxpho2);
e7d87aef 922 fOutputList->Add(fSameElecPtMCM20);
923
987053ce 924 fIncpTMCM20pho_pi0e = new THnSparseD("fIncpTMCM20pho_pi0e","MC Pho pi0->e pid electro vs. centrality with M20",5,nBinspho2,minpho2,maxpho2);
925 fIncpTMCM20pho_pi0e->Sumw2();
e4b0faf2 926 fOutputList->Add(fIncpTMCM20pho_pi0e);
927
987053ce 928 fPhoElecPtMCM20_pi0e = new THnSparseD("fPhoElecPtMCM20_pi0e", "MC Pho-inclusive electron pt with M20 pi0->e",5,nBinspho2,minpho2,maxpho2);
929 fPhoElecPtMCM20_pi0e->Sumw2();
e4b0faf2 930 fOutputList->Add(fPhoElecPtMCM20_pi0e);
931
987053ce 932 fSameElecPtMCM20_pi0e = new THnSparseD("fSameElecPtMCM20_pi0e", "MC Same-inclusive electron pt pi0->e",5,nBinspho2,minpho2,maxpho2);
933 fSameElecPtMCM20_pi0e->Sumw2();
e4b0faf2 934 fOutputList->Add(fSameElecPtMCM20_pi0e);
935
a6df418c 936 CheckNclust = new TH1D("CheckNclust","cluster check",200,0,200);
937 fOutputList->Add(CheckNclust);
938
939 CheckNits = new TH1D("CheckNits","ITS cluster check",8,-0.5,7.5);
940 fOutputList->Add(CheckNits);
e7d87aef 941
d113d7cd 942 Hpi0pTcheck = new TH2D("Hpi0pTcheck","Pi0 pT from Hijing",100,0,50,3,-0.5,2.5);
943 fOutputList->Add(Hpi0pTcheck);
944
e4b0faf2 945 HphopTcheck = new TH2D("HphopTcheck","Pho pT from Hijing",100,0,50,3,-0.5,2.5);
946 fOutputList->Add(HphopTcheck);
947
50919258 948 fMomDtoE = new TH2D("fMomDtoE","D->E pT correlations;e p_{T} GeV/c;D p_{T} GeV/c",400,0,40,400,0,40);
949 fOutputList->Add(fMomDtoE);
d113d7cd 950
bfc7c23b 951 PostData(1,fOutputList);
952}
953
954//________________________________________________________________________
955void AliAnalysisTaskHFECal::Terminate(Option_t *)
956{
957 // Info("Terminate");
958 AliAnalysisTaskSE::Terminate();
959}
960
961//________________________________________________________________________
962Bool_t AliAnalysisTaskHFECal::ProcessCutStep(Int_t cutStep, AliVParticle *track)
963{
964 // Check single track cuts for a given cut step
965 const Int_t kMCOffset = AliHFEcuts::kNcutStepsMCTrack;
966 if(!fCFM->CheckParticleCuts(cutStep + kMCOffset, track)) return kFALSE;
967 return kTRUE;
968}
969//_________________________________________
f09766dd 970//void AliAnalysisTaskHFECal::SelectPhotonicElectron(Int_t itrack, Double_t cent, AliESDtrack *track, Bool_t &fFlagPhotonicElec)
c2801a73 971//void AliAnalysisTaskHFECal::SelectPhotonicElectron(Int_t itrack, Double_t cent, AliESDtrack *track, Bool_t &fFlagPhotonicElec, Bool_t &fFlagConvinatElec, Double_t nSig)
972void 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)
bfc7c23b 973{
974 //Identify non-heavy flavour electrons using Invariant mass method
975
fb87d707 976 fTrackCuts->SetAcceptKinkDaughters(kFALSE);
feffe705 977 fTrackCuts->SetRequireTPCRefit(kTRUE);
85985bb0 978 fTrackCuts->SetRequireITSRefit(kTRUE);
feffe705 979 fTrackCuts->SetEtaRange(-0.9,0.9);
f09766dd 980 //fTrackCuts->SetRequireSigmaToVertex(kTRUE);
bfc7c23b 981 fTrackCuts->SetMaxChi2PerClusterTPC(3.5);
feffe705 982 fTrackCuts->SetMinNClustersTPC(90);
bfc7c23b 983
984 const AliESDVertex *pVtx = fESD->GetPrimaryVertex();
985
986 Bool_t flagPhotonicElec = kFALSE;
f09766dd 987 Bool_t flagConvinatElec = kFALSE;
bfc7c23b 988
5e65985c 989 //for(Int_t jTracks = itrack+1; jTracks<fESD->GetNumberOfTracks(); jTracks++){
990 for(Int_t jTracks = 0; jTracks<fESD->GetNumberOfTracks(); jTracks++){
bfc7c23b 991 AliESDtrack* trackAsso = fESD->GetTrack(jTracks);
992 if (!trackAsso) {
993 printf("ERROR: Could not receive track %d\n", jTracks);
994 continue;
995 }
5e65985c 996 if(itrack==jTracks)continue;
997
f09766dd 998 Double_t dEdxAsso = -999., ptPrim=-999., ptAsso=-999., openingAngle = -999.;
bfc7c23b 999 Double_t mass=999., width = -999;
1000 Bool_t fFlagLS=kFALSE, fFlagULS=kFALSE;
1001
f09766dd 1002 ptPrim = track->Pt();
1003
bfc7c23b 1004 dEdxAsso = trackAsso->GetTPCsignal();
1005 ptAsso = trackAsso->Pt();
1006 Int_t chargeAsso = trackAsso->Charge();
1007 Int_t charge = track->Charge();
1008
b12bc641 1009
3db00c72 1010 if(ptAsso <0.5) continue;
bfc7c23b 1011 if(!fTrackCuts->AcceptTrack(trackAsso)) continue;
f09766dd 1012 if(dEdxAsso <65 || dEdxAsso>100) continue; //11a pass1
bfc7c23b 1013
1014 Int_t fPDGe1 = 11; Int_t fPDGe2 = 11;
1015 if(charge>0) fPDGe1 = -11;
1016 if(chargeAsso>0) fPDGe2 = -11;
b12bc641 1017
1018 //printf("chargeAsso = %d\n",chargeAsso);
1019 //printf("charge = %d\n",charge);
bfc7c23b 1020 if(charge == chargeAsso) fFlagLS = kTRUE;
1021 if(charge != chargeAsso) fFlagULS = kTRUE;
1022
b12bc641 1023 //printf("fFlagLS = %d\n",fFlagLS);
1024 //printf("fFlagULS = %d\n",fFlagULS);
1025 //printf("\n");
1026
bfc7c23b 1027 AliKFParticle ge1(*track, fPDGe1);
1028 AliKFParticle ge2(*trackAsso, fPDGe2);
1029 AliKFParticle recg(ge1, ge2);
1030
1031 if(recg.GetNDF()<1) continue;
1032 Double_t chi2recg = recg.GetChi2()/recg.GetNDF();
1033 if(TMath::Sqrt(TMath::Abs(chi2recg))>3.) continue;
1034
1035 AliKFVertex primV(*pVtx);
1036 primV += recg;
1037 recg.SetProductionVertex(primV);
1038
d113d7cd 1039 recg.SetMassConstraint(0,0.0001);
bfc7c23b 1040
1041 openingAngle = ge1.GetAngle(ge2);
1042 if(fFlagLS) fOpeningAngleLS->Fill(openingAngle);
1043 if(fFlagULS) fOpeningAngleULS->Fill(openingAngle);
1044
bfc7c23b 1045
1046 recg.GetMass(mass,width);
1047
c2801a73 1048 double ishower = 0;
1049 if(shower>0.0 && shower<0.3)ishower = 1;
1050
b12bc641 1051 double phoinfo[9];
f09766dd 1052 phoinfo[0] = cent;
1053 phoinfo[1] = ptPrim;
1054 phoinfo[2] = mass;
3db00c72 1055 phoinfo[3] = nSig;
4e01b68c 1056 phoinfo[4] = openingAngle;
c2801a73 1057 phoinfo[5] = ishower;
1058 phoinfo[6] = ep;
1059 phoinfo[7] = mce;
b12bc641 1060 phoinfo[8] = ptAsso;
f09766dd 1061
1062 if(fFlagLS) fInvmassLS->Fill(phoinfo);
1063 if(fFlagULS) fInvmassULS->Fill(phoinfo);
ffa14dda 1064
b12bc641 1065 //printf("fInvmassCut %f\n",fInvmassCut);
1066 //printf("openingAngle %f\n",fOpeningAngleCut);
1067
ffa14dda 1068 if(openingAngle > fOpeningAngleCut) continue;
1069
bfc7c23b 1070 if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec){
1071 flagPhotonicElec = kTRUE;
1072 }
f09766dd 1073 if(mass<fInvmassCut && fFlagLS && !flagConvinatElec){
1074 flagConvinatElec = kTRUE;
1075 }
bfc7c23b 1076
1077 }
1078 fFlagPhotonicElec = flagPhotonicElec;
f09766dd 1079 fFlagConvinatElec = flagConvinatElec;
bfc7c23b 1080
1081}
1082
fb87d707 1083
1084//_________________________________________
1085void AliAnalysisTaskHFECal::FindTriggerClusters()
1086{
1087 // constants
1088 const int nModuleCols = 2;
1089 const int nModuleRows = 5;
1090 const int nColsFeeModule = 48;
1091 const int nRowsFeeModule = 24;
1092 const int nColsFaltroModule = 24;
1093 const int nRowsFaltroModule = 12;
1094 //const int faltroWidthMax = 20;
1095
1096 // part 1, trigger extraction -------------------------------------
1097 Int_t globCol, globRow;
1098 //Int_t ntimes=0, nTrigChannel=0, nTrigChannelCut=0, trigInCut=0;
1099 Int_t ntimes=0, nTrigChannel=0, nTrigChannelCut=0;
1100
1101 //Int_t trigtimes[faltroWidthMax];
1102 Double_t cellTime[nColsFeeModule*nModuleCols][nRowsFeeModule*nModuleRows];
1103 Double_t cellEnergy[nColsFeeModule*nModuleCols][nRowsFeeModule*nModuleRows];
1104 //Double_t fTrigCutLow = 6;
1105 //Double_t fTrigCutHigh = 10;
1106 Double_t fTimeCutLow = 469e-09;
1107 Double_t fTimeCutHigh = 715e-09;
1108
1109 AliESDCaloTrigger * fCaloTrigger = fESD->GetCaloTrigger( "EMCAL" );
1110
1111 // erase trigger maps
1112 for(Int_t i = 0; i < nColsFaltroModule*nModuleCols; i++ )
1113 {
1114 for(Int_t j = 0; j < nRowsFaltroModule*nModuleRows; j++ )
1115 {
1116 ftriggersCut[i][j] = 0;
1117 ftriggers[i][j] = 0;
1118 ftriggersTime[i][j] = 0;
1119 }
1120 }
1121
1122 Int_t iglobCol=0, iglobRow=0;
1123 // go through triggers
1124 if( fCaloTrigger->GetEntries() > 0 )
1125 {
1126 // needs reset
1127 fCaloTrigger->Reset();
1128 while( fCaloTrigger->Next() )
1129 {
1130 fCaloTrigger->GetPosition( globCol, globRow );
1131 fCaloTrigger->GetNL0Times( ntimes );
1132 /*
1133 // no L0s
1134 if( ntimes < 1 ) continue;
1135 // get precise timings
1136 fCaloTrigger->GetL0Times( trigtimes );
1137 trigInCut = 0;
1138 for(Int_t i = 0; i < ntimes; i++ )
1139 {
1140 // save the first trigger time in channel
1141 if( i == 0 || triggersTime[globCol][globRow] > trigtimes[i] )
1142 triggersTime[globCol][globRow] = trigtimes[i];
1143 //printf("trigger times: %d\n",trigtimes[i]);
1144 // check if in cut
1145 if(trigtimes[i] > fTrigCutLow && trigtimes[i] < fTrigCutHigh )
1146 trigInCut = 1;
1147
1148 fTrigTimes->Fill(trigtimes[i]);
1149 }
1150 */
1151
1152 //L1 analysis from AliAnalysisTaskEMCALTriggerQA
1153 Int_t bit = 0;
1154 fCaloTrigger->GetTriggerBits(bit);
1155
1156 Int_t ts = 0;
1157 fCaloTrigger->GetL1TimeSum(ts);
1158 if (ts > 0)ftriggers[globCol][globRow] = 1;
1159 // number of triggered channels in event
1160 nTrigChannel++;
1161 // ... inside cut
1162 if(ts>0 && (bit >> 6 & 0x1))
1163 {
1164 iglobCol = globCol;
1165 iglobRow = globRow;
1166 nTrigChannelCut++;
1167 ftriggersCut[globCol][globRow] = 1;
1168 }
1169
1170 } // calo trigger entries
1171 } // has calo trigger entries
1172
1173 // part 2 go through the clusters here -----------------------------------
1174 Int_t nCluster=0, nCell=0, iCell=0, gCell=0;
1175 Short_t cellAddr, nSACell, mclabel;
1176 //Int_t nSACell, iSACell, mclabel;
1177 Int_t iSACell;
42c75692 1178 Double_t cellAmp=0, cellTimeT=0, clusterTime=0, efrac=0;
fb87d707 1179 Int_t nSupMod, nModule, nIphi, nIeta, iphi, ieta, gphi, geta, feta, fphi;
1180
1181 TRefArray *fCaloClusters = new TRefArray();
1182 fESD->GetEMCALClusters( fCaloClusters );
1183 nCluster = fCaloClusters->GetEntries();
1184
1185
1186 // save all cells times if there are clusters
1187 if( nCluster > 0 ){
1188 // erase time map
1189 for(Int_t i = 0; i < nColsFeeModule*nModuleCols; i++ ){
1190 for(Int_t j = 0; j < nRowsFeeModule*nModuleRows; j++ ){
1191 cellTime[i][j] = 0.;
1192 cellEnergy[i][j] = 0.;
1193 }
1194 }
1195
1196 // get cells
1197 AliESDCaloCells *fCaloCells = fESD->GetEMCALCells();
1198 //AliVCaloCells fCaloCells = fESD->GetEMCALCells();
1199 nSACell = fCaloCells->GetNumberOfCells();
1200 for(iSACell = 0; iSACell < nSACell; iSACell++ ){
1201 // get the cell info *fCal
1202 fCaloCells->GetCell( iSACell, cellAddr, cellAmp, cellTimeT , mclabel, efrac);
1203
1204 // get cell position
1205 fGeom->GetCellIndex( cellAddr, nSupMod, nModule, nIphi, nIeta );
1206 fGeom->GetCellPhiEtaIndexInSModule( nSupMod,nModule, nIphi, nIeta, iphi, ieta);
1207
1208 // convert co global phi eta
1209 gphi = iphi + nRowsFeeModule*(nSupMod/2);
1210 geta = ieta + nColsFeeModule*(nSupMod%2);
1211
1212 // save cell time and energy
1213 cellTime[geta][gphi] = cellTimeT;
1214 cellEnergy[geta][gphi] = cellAmp;
1215
1216 }
1217 }
1218
1219 Int_t nClusterTrig, nClusterTrigCut;
1220 UShort_t *cellAddrs;
1221 Double_t clsE=-999, clsEta=-999, clsPhi=-999;
1222 Float_t clsPos[3] = {0.,0.,0.};
1223
1224 for(Int_t icl=0; icl<fESD->GetNumberOfCaloClusters(); icl++)
1225 {
1226 AliESDCaloCluster *cluster = fESD->GetCaloCluster(icl);
1227 if(!cluster || !cluster->IsEMCAL()) continue;
1228
1229 // get cluster cells
1230 nCell = cluster->GetNCells();
1231
1232 // get cluster energy
1233 clsE = cluster->E();
1234
1235 // get cluster position
1236 cluster->GetPosition(clsPos);
1237 TVector3 clsPosVec(clsPos[0],clsPos[1],clsPos[2]);
1238 clsEta = clsPosVec.Eta();
1239 clsPhi = clsPosVec.Phi();
1240
1241 // get the cell addresses
1242 cellAddrs = cluster->GetCellsAbsId();
1243
1244 // check if the cluster contains cell, that was marked as triggered
1245 nClusterTrig = 0;
1246 nClusterTrigCut = 0;
1247
1248 // loop the cells to check, if cluser in acceptance
1249 // any cluster with a cell outside acceptance is not considered
1250 for( iCell = 0; iCell < nCell; iCell++ )
1251 {
42c75692 1252 // check hot cell
feffe705 1253 //if(clsE>6.0)fCellCheck->Fill(clsE,cellAddrs[iCell]);
42c75692 1254
fb87d707 1255 // get cell position
1256 fGeom->GetCellIndex( cellAddrs[iCell], nSupMod, nModule, nIphi, nIeta );
1257 fGeom->GetCellPhiEtaIndexInSModule( nSupMod,nModule, nIphi, nIeta, iphi, ieta);
1258
1259 // convert co global phi eta
1260 gphi = iphi + nRowsFeeModule*(nSupMod/2);
1261 geta = ieta + nColsFeeModule*(nSupMod%2);
1262
1263 if( cellTime[geta][gphi] > 0. ){
1264 clusterTime += cellTime[geta][gphi];
1265 gCell++;
1266 }
1267
1268 // get corresponding FALTRO
1269 fphi = gphi / 2;
1270 feta = geta / 2;
1271
1272 // try to match with a triggered
1273 if( ftriggers[feta][fphi]==1)
1274 { nClusterTrig++;
1275 }
1276 if( ftriggersCut[feta][fphi]==1)
1277 { nClusterTrigCut++;
1278 }
1279
1280 } // cells
1281
1282
1283 if( gCell > 0 )
1284 clusterTime = clusterTime / (Double_t)gCell;
1285 // fix the reconstriction code time 100ns jumps
1286 if( fESD->GetBunchCrossNumber() % 4 < 2 )
1287 clusterTime -= 0.0000001;
1288
feffe705 1289 //fClsETime->Fill(clsE,clusterTime);
1290 //fClsEBftTrigCut->Fill(clsE);
fb87d707 1291
1292 if(nClusterTrig>0){
feffe705 1293 //fClsETime1->Fill(clsE,clusterTime);
fb87d707 1294 }
1295
1296 if(nClusterTrig>0){
1297 cluster->SetChi2(1);
feffe705 1298 //fClsEAftTrigCut1->Fill(clsE);
fb87d707 1299 }
1300
1301 if(nClusterTrigCut>0){
1302 cluster->SetChi2(2);
feffe705 1303 //fClsEAftTrigCut2->Fill(clsE);
fb87d707 1304 }
1305
1306 if(nClusterTrigCut>0 && ( nCell > (1 + clsE / 3)))
1307 {
1308 cluster->SetChi2(3);
feffe705 1309 //fClsEAftTrigCut3->Fill(clsE);
fb87d707 1310 }
1311
1312 if(nClusterTrigCut>0 && (nCell > (1 + clsE / 3) )&&( clusterTime > fTimeCutLow && clusterTime < fTimeCutHigh ))
1313 {
1314 // cluster->SetChi2(4);
feffe705 1315 //fClsEAftTrigCut4->Fill(clsE);
fb87d707 1316 }
1317 if(nClusterTrigCut<1)
1318 {
1319 cluster->SetChi2(0);
1320
feffe705 1321 //fClsEAftTrigCut->Fill(clsE);
fb87d707 1322 }
1323
1324 } // clusters
1325}
1326
1327
1328
1329
1330
1331
1332