]> 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
modify TPC cut parameter
[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 {
d113d7cd 647 double phoval[4];
44be9e1d 648 phoval[0] = cent;
649 phoval[1] = pt;
650 phoval[2] = fTPCnSigma;
d113d7cd 651 phoval[3] = iHijing;
44be9e1d 652
653 fIncpTMCpho->Fill(phoval);
654 if(fFlagPhotonicElec) fPhoElecPtMC->Fill(phoval);
655 if(fFlagConvinatElec) fSameElecPtMC->Fill(phoval);
e7d87aef 656
657 if(m20>0.0 && m20<0.3)
658 {
44be9e1d 659 fIncpTMCM20pho->Fill(phoval);
660 if(fFlagPhotonicElec) fPhoElecPtMCM20->Fill(phoval);
661 if(fFlagConvinatElec) fSameElecPtMCM20->Fill(phoval);
e4b0faf2 662 // pi0->g->e
663 if(mcWeight>-1)
664 {
665 fIncpTMCM20pho_pi0e->Fill(phoval,mcWeight);
666 if(fFlagPhotonicElec) fPhoElecPtMCM20_pi0e->Fill(phoval,mcWeight);
667 if(fFlagConvinatElec) fSameElecPtMCM20_pi0e->Fill(phoval,mcWeight);
668 }
e7d87aef 669 }
670 }
671 }
bfc7c23b 672 }
673 PostData(1, fOutputList);
674}
675//_________________________________________
676void AliAnalysisTaskHFECal::UserCreateOutputObjects()
677{
42c75692 678 //--- Check MC
679
e7d87aef 680 //Bool_t mcData = kFALSE;
42c75692 681 if(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler())
682 {
e7d87aef 683 fmcData = kTRUE;
42c75692 684 printf("+++++ MC Data available");
685 }
e7d87aef 686 if(fmcData)
42c75692 687 {
688 printf("++++++++= MC analysis \n");
689 }
690 else
691 {
692 printf("++++++++ real data analysis \n");
693 }
694
85985bb0 695 printf("+++++++ QA hist %d",fqahist);
696
fb87d707 697 //---- Geometry
698 fGeom = AliEMCALGeometry::GetInstance("EMCAL_COMPLETEV1");
699
bfc7c23b 700 //--------Initialize PID
42c75692 701 //fPID->SetHasMCData(kFALSE);
e7d87aef 702 fPID->SetHasMCData(fmcData);
bfc7c23b 703 if(!fPID->GetNumberOfPIDdetectors())
704 {
705 fPID->AddDetector("TPC", 0);
706 fPID->AddDetector("EMCAL", 1);
707 }
708
3db00c72 709 Double_t params[4];
78ff1095 710 const char *cutmodel;
3db00c72 711 cutmodel = "pol0";
712 params[0] = -1.0; //sigma min
9b3495ff 713 double maxnSig = 3.0;
714 if(fmcData)
715 {
716 params[0] = -5.0; //sigma min
717 maxnSig = 5.0;
718 }
719
720 for(Int_t a=0;a<11;a++)fPID->ConfigureTPCcentralityCut(a,cutmodel,params,maxnSig);
3db00c72 721
bfc7c23b 722 fPID->SortDetectors();
723 fPIDqa = new AliHFEpidQAmanager();
724 fPIDqa->Initialize(fPID);
a6df418c 725
726 //------- fcut --------------
727 fCuts = new AliHFEcuts();
728 fCuts->CreateStandardCuts();
729 //fCuts->SetMinNClustersTPC(100);
730 fCuts->SetMinNClustersTPC(90);
731 fCuts->SetMinRatioTPCclusters(0.6);
732 fCuts->SetTPCmodes(AliHFEextraCuts::kFound, AliHFEextraCuts::kFoundOverFindable);
733 //fCuts->SetMinNClustersITS(3);
734 fCuts->SetMinNClustersITS(2);
735 fCuts->SetCutITSpixel(AliHFEextraCuts::kAny);
736 fCuts->SetCheckITSLayerStatus(kFALSE);
737 fCuts->SetVertexRange(10.);
738 fCuts->SetTOFPIDStep(kFALSE);
739 fCuts->SetPtRange(2, 50);
740 fCuts->SetMaxImpactParam(3.,3.);
741
bfc7c23b 742 //--------Initialize correction Framework and Cuts
743 fCFM = new AliCFManager;
744 const Int_t kNcutSteps = AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kNcutStepsRecTrack + AliHFEcuts::kNcutStepsDETrack;
745 fCFM->SetNStepParticle(kNcutSteps);
746 for(Int_t istep = 0; istep < kNcutSteps; istep++)
747 fCFM->SetParticleCutsList(istep, NULL);
748
749 if(!fCuts){
750 AliWarning("Cuts not available. Default cuts will be used");
751 fCuts = new AliHFEcuts;
752 fCuts->CreateStandardCuts();
753 }
754 fCuts->Initialize(fCFM);
755
756 //---------Output Tlist
757 fOutputList = new TList();
758 fOutputList->SetOwner();
759 fOutputList->Add(fPIDqa->MakeList("PIDQA"));
760
761 fNoEvents = new TH1F("fNoEvents","",4,-0.5,3.5) ;
762 fOutputList->Add(fNoEvents);
763
02cd82a0 764 Int_t binsE[5] = {250, 100, 1000, 200, 10};
fb87d707 765 Double_t xminE[5] = {1.0, -1, 0.0, 0, -0.5};
766 Double_t xmaxE[5] = {3.5, 1, 100.0, 100, 9.5};
767 fEMCAccE = new THnSparseD("fEMCAccE","EMC acceptance & E;#eta;#phi;Energy;Centrality;trugCondition;",5,binsE,xminE,xmaxE);
85985bb0 768 if(fqahist==1)fOutputList->Add(fEMCAccE);
bfc7c23b 769
770 fTrkpt = new TH1F("fTrkpt","track pt",100,0,50);
771 fOutputList->Add(fTrkpt);
772
773 fTrackPtBefTrkCuts = new TH1F("fTrackPtBefTrkCuts","track pt before track cuts",100,0,50);
774 fOutputList->Add(fTrackPtBefTrkCuts);
775
776 fTrackPtAftTrkCuts = new TH1F("fTrackPtAftTrkCuts","track pt after track cuts",100,0,50);
777 fOutputList->Add(fTrackPtAftTrkCuts);
778
779 fTPCnsigma = new TH2F("fTPCnsigma", "TPC - n sigma",100,0,50,200,-10,10);
780 fOutputList->Add(fTPCnsigma);
781
782 fTrkEovPBef = new TH2F("fTrkEovPBef","track E/p before HFE pid",100,0,50,100,0,2);
783 fOutputList->Add(fTrkEovPBef);
784
785 fTrkEovPAft = new TH2F("fTrkEovPAft","track E/p after HFE pid",100,0,50,100,0,2);
786 fOutputList->Add(fTrkEovPAft);
787
42c75692 788 fdEdxBef = new TH2F("fdEdxBef","track dEdx vs p before HFE pid",100,0,50,200,-10,10);
bfc7c23b 789 fOutputList->Add(fdEdxBef);
790
42c75692 791 fdEdxAft = new TH2F("fdEdxAft","track dEdx vs p after HFE pid",100,0,50,200,-10,10);
bfc7c23b 792 fOutputList->Add(fdEdxAft);
793
02cd82a0 794 fIncpT = new TH2F("fIncpT","HFE pid electro vs. centrality",200,0,100,100,0,50);
bfc7c23b 795 fOutputList->Add(fIncpT);
796
02cd82a0 797 fIncpTM20 = new TH2F("fIncpTM20","HFE pid electro vs. centrality with M20",200,0,100,100,0,50);
fb87d707 798 fOutputList->Add(fIncpTM20);
4e01b68c 799
02cd82a0 800 Int_t nBinspho[9] = { 200, 100, 500, 12, 50, 4, 200, 8, 100};
b12bc641 801 Double_t minpho[9] = { 0., 0., 0., -2.5, 0, -0.5, 0,-1.5, 0};
802 Double_t maxpho[9] = {100., 50., 0.5, 3.5, 1, 3.5, 2, 6.5, 50};
f09766dd 803
b12bc641 804 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 805 fOutputList->Add(fInvmassLS);
806
b12bc641 807 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 808 fOutputList->Add(fInvmassULS);
809
810 fOpeningAngleLS = new TH1F("fOpeningAngleLS","Opening angle for LS pairs",100,0,1);
811 fOutputList->Add(fOpeningAngleLS);
812
813 fOpeningAngleULS = new TH1F("fOpeningAngleULS","Opening angle for ULS pairs",100,0,1);
814 fOutputList->Add(fOpeningAngleULS);
815
816 fPhotoElecPt = new TH1F("fPhotoElecPt", "photonic electron pt",100,0,50);
817 fOutputList->Add(fPhotoElecPt);
818
02cd82a0 819 fPhoElecPt = new TH2F("fPhoElecPt", "Pho-inclusive electron pt",200,0,100,100,0,50);
bfc7c23b 820 fOutputList->Add(fPhoElecPt);
821
02cd82a0 822 fPhoElecPtM20 = new TH2F("fPhoElecPtM20", "Pho-inclusive electron pt with M20",200,0,100,100,0,50);
fb87d707 823 fOutputList->Add(fPhoElecPtM20);
824
02cd82a0 825 fSameElecPt = new TH2F("fSameElecPt", "Same-inclusive electron pt",200,0,100,100,0,50);
f09766dd 826 fOutputList->Add(fSameElecPt);
827
02cd82a0 828 fSameElecPtM20 = new TH2F("fSameElecPtM20", "Same-inclusive electron pt with M20",200,0,100,100,0,50);
fb87d707 829 fOutputList->Add(fSameElecPtM20);
830
02cd82a0 831 fCent = new TH1F("fCent","Centrality",200,0,100) ;
bfc7c23b 832 fOutputList->Add(fCent);
833
834 // Make common binning
44be9e1d 835 const Double_t kMinP = 0.;
bfc7c23b 836 const Double_t kMaxP = 50.;
a6df418c 837 //const Double_t kTPCSigMim = 40.;
838 //const Double_t kTPCSigMax = 140.;
bfc7c23b 839
840 // 1st histogram: TPC dEdx with/without EMCAL (p, pT, TPC Signal, phi, eta, Sig, e/p, ,match, cell, M02, M20, Disp, Centrality, select)
a6df418c 841 Int_t nBins[16] = { 250, 10, 60, 20, 100, 300, 50, 40, 200, 200, 250, 200, 3, 5, 100, 8};
842 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 843 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};
844 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 845 if(fqahist==1)fOutputList->Add(fEleInfo);
bfc7c23b 846
fb87d707 847 //<--- Trigger info
feffe705 848 /*
fb87d707 849 fClsEBftTrigCut = new TH1F("fClsEBftTrigCut","cluster E before trigger selection",1000,0,100);
850 fOutputList->Add(fClsEBftTrigCut);
851
852 fClsEAftTrigCut = new TH1F("fClsEAftTrigCut","cluster E if cls has 0 trigcut channel",1000,0,100);
853 fOutputList->Add(fClsEAftTrigCut);
854
855 fClsEAftTrigCut1 = new TH1F("fClsEAftTrigCut1","cluster E if cls with trig channel",1000,0,100);
856 fOutputList->Add(fClsEAftTrigCut1);
857
858 fClsEAftTrigCut2 = new TH1F("fClsEAftTrigCut2","cluster E if cls with trigcut channel",1000,0,100);
859 fOutputList->Add(fClsEAftTrigCut2);
860
861 fClsEAftTrigCut3 = new TH1F("fClsEAftTrigCut3","cluster E if cls with trigcut channel + nCell>Ecorrect",1000,0,100);
862 fOutputList->Add(fClsEAftTrigCut3);
863
864 fClsEAftTrigCut4 = new TH1F("fClsEAftTrigCut4","cluster E if cls with trigcut channel + nCell>Ecorrect + cls time cut",1000,0,100);
865 fOutputList->Add(fClsEAftTrigCut4);
866
867 fClsETime = new TH2F("fClsETime", "Cls time vs E; E; time;",1000,0,100,1000,-0.0000002,0.0000009);
868 fOutputList->Add(fClsETime);
869
870 fClsETime1 = new TH2F("fClsETime1", "Cls time vs E if cls contains trigger channel; E; time;",1000,0,100,1000,-0.0000002,0.0000009);
871 fOutputList->Add(fClsETime1);
872
873 fTrigTimes = new TH1F("fTrigTimes", "Trigger time; time; N;",25,0,25);
874 fOutputList->Add(fTrigTimes);
875
42c75692 876 fCellCheck = new TH2F("fCellCheck", "Cell vs E; E GeV; Cell ID",10,6,26,12000,0,12000);
877 fOutputList->Add(fCellCheck);
feffe705 878 */
e7d87aef 879 //<---------- MC
880
02cd82a0 881 fInputHFEMC = new TH2F("fInputHFEMC","Input MC HFE pid electro vs. centrality",200,0,100,100,0,50);
e7d87aef 882 fOutputList->Add(fInputHFEMC);
883
02cd82a0 884 fInputAlle = new TH2F("fInputAlle","Input MC electro vs. centrality",200,0,100,100,0,50);
feffe705 885 fOutputList->Add(fInputAlle);
886
02cd82a0 887 fIncpTMChfe = new TH2F("fIncpTMChfe","MC HFE pid electro vs. centrality",200,0,100,100,0,50);
e7d87aef 888 fOutputList->Add(fIncpTMChfe);
889
02cd82a0 890 fIncpTMChfeAll = new TH2F("fIncpTMChfeAll","MC Alle pid electro vs. centrality",200,0,100,100,0,50);
3db00c72 891 fOutputList->Add(fIncpTMChfeAll);
892
02cd82a0 893 fIncpTMCM20hfe = new TH2F("fIncpTMCM20hfe","MC HFE pid electro vs. centrality with M20",200,0,100,100,0,50);
e7d87aef 894 fOutputList->Add(fIncpTMCM20hfe);
895
02cd82a0 896 fIncpTMCM20hfeAll = new TH2F("fIncpTMCM20hfeAll","MC Alle pid electro vs. centrality with M20",200,0,100,100,0,50);
3db00c72 897 fOutputList->Add(fIncpTMCM20hfeAll);
898
44be9e1d 899
d113d7cd 900 Int_t nBinspho2[4] = { 200, 100, 7, 3};
901 Double_t minpho2[4] = { 0., 0., -2.5, -0.5};
902 Double_t maxpho2[4] = {100., 50., 4.5, 2.5};
44be9e1d 903
d113d7cd 904 fIncpTMCpho = new THnSparseD("fIncpTMCpho","MC Pho pid electro vs. centrality",4,nBinspho2,minpho2,maxpho2);
e7d87aef 905 fOutputList->Add(fIncpTMCpho);
906
d113d7cd 907 fIncpTMCM20pho = new THnSparseD("fIncpTMCM20pho","MC Pho pid electro vs. centrality with M20",4,nBinspho2,minpho2,maxpho2);
e7d87aef 908 fOutputList->Add(fIncpTMCM20pho);
909
d113d7cd 910 fPhoElecPtMC = new THnSparseD("fPhoElecPtMC", "MC Pho-inclusive electron pt",4,nBinspho2,minpho2,maxpho2);
e7d87aef 911 fOutputList->Add(fPhoElecPtMC);
912
d113d7cd 913 fPhoElecPtMCM20 = new THnSparseD("fPhoElecPtMCM20", "MC Pho-inclusive electron pt with M20",4,nBinspho2,minpho2,maxpho2);
feffe705 914 fOutputList->Add(fPhoElecPtMCM20);
e7d87aef 915
d113d7cd 916 fSameElecPtMC = new THnSparseD("fSameElecPtMC", "MC Same-inclusive electron pt",4,nBinspho2,minpho2,maxpho2);
e7d87aef 917 fOutputList->Add(fSameElecPtMC);
918
d113d7cd 919 fSameElecPtMCM20 = new THnSparseD("fSameElecPtMCM20", "MC Same-inclusive electron pt with M20",4,nBinspho2,minpho2,maxpho2);
e7d87aef 920 fOutputList->Add(fSameElecPtMCM20);
921
e4b0faf2 922 fIncpTMCM20pho_pi0e = new THnSparseD("fIncpTMCM20pho_pi0e","MC Pho pi0->e pid electro vs. centrality with M20",4,nBinspho2,minpho2,maxpho2);
923 fOutputList->Add(fIncpTMCM20pho_pi0e);
924
925 fPhoElecPtMCM20_pi0e = new THnSparseD("fPhoElecPtMCM20_pi0e", "MC Pho-inclusive electron pt with M20 pi0->e",4,nBinspho2,minpho2,maxpho2);
926 fOutputList->Add(fPhoElecPtMCM20_pi0e);
927
928 fSameElecPtMCM20_pi0e = new THnSparseD("fSameElecPtMCM20_pi0e", "MC Same-inclusive electron pt pi0->e",4,nBinspho2,minpho2,maxpho2);
929 fOutputList->Add(fSameElecPtMCM20_pi0e);
930
a6df418c 931 CheckNclust = new TH1D("CheckNclust","cluster check",200,0,200);
932 fOutputList->Add(CheckNclust);
933
934 CheckNits = new TH1D("CheckNits","ITS cluster check",8,-0.5,7.5);
935 fOutputList->Add(CheckNits);
e7d87aef 936
d113d7cd 937 Hpi0pTcheck = new TH2D("Hpi0pTcheck","Pi0 pT from Hijing",100,0,50,3,-0.5,2.5);
938 fOutputList->Add(Hpi0pTcheck);
939
e4b0faf2 940 HphopTcheck = new TH2D("HphopTcheck","Pho pT from Hijing",100,0,50,3,-0.5,2.5);
941 fOutputList->Add(HphopTcheck);
942
50919258 943 fMomDtoE = new TH2D("fMomDtoE","D->E pT correlations;e p_{T} GeV/c;D p_{T} GeV/c",400,0,40,400,0,40);
944 fOutputList->Add(fMomDtoE);
d113d7cd 945
bfc7c23b 946 PostData(1,fOutputList);
947}
948
949//________________________________________________________________________
950void AliAnalysisTaskHFECal::Terminate(Option_t *)
951{
952 // Info("Terminate");
953 AliAnalysisTaskSE::Terminate();
954}
955
956//________________________________________________________________________
957Bool_t AliAnalysisTaskHFECal::ProcessCutStep(Int_t cutStep, AliVParticle *track)
958{
959 // Check single track cuts for a given cut step
960 const Int_t kMCOffset = AliHFEcuts::kNcutStepsMCTrack;
961 if(!fCFM->CheckParticleCuts(cutStep + kMCOffset, track)) return kFALSE;
962 return kTRUE;
963}
964//_________________________________________
f09766dd 965//void AliAnalysisTaskHFECal::SelectPhotonicElectron(Int_t itrack, Double_t cent, AliESDtrack *track, Bool_t &fFlagPhotonicElec)
c2801a73 966//void AliAnalysisTaskHFECal::SelectPhotonicElectron(Int_t itrack, Double_t cent, AliESDtrack *track, Bool_t &fFlagPhotonicElec, Bool_t &fFlagConvinatElec, Double_t nSig)
967void 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 968{
969 //Identify non-heavy flavour electrons using Invariant mass method
970
fb87d707 971 fTrackCuts->SetAcceptKinkDaughters(kFALSE);
feffe705 972 fTrackCuts->SetRequireTPCRefit(kTRUE);
85985bb0 973 fTrackCuts->SetRequireITSRefit(kTRUE);
feffe705 974 fTrackCuts->SetEtaRange(-0.9,0.9);
f09766dd 975 //fTrackCuts->SetRequireSigmaToVertex(kTRUE);
bfc7c23b 976 fTrackCuts->SetMaxChi2PerClusterTPC(3.5);
feffe705 977 fTrackCuts->SetMinNClustersTPC(90);
bfc7c23b 978
979 const AliESDVertex *pVtx = fESD->GetPrimaryVertex();
980
981 Bool_t flagPhotonicElec = kFALSE;
f09766dd 982 Bool_t flagConvinatElec = kFALSE;
bfc7c23b 983
5e65985c 984 //for(Int_t jTracks = itrack+1; jTracks<fESD->GetNumberOfTracks(); jTracks++){
985 for(Int_t jTracks = 0; jTracks<fESD->GetNumberOfTracks(); jTracks++){
bfc7c23b 986 AliESDtrack* trackAsso = fESD->GetTrack(jTracks);
987 if (!trackAsso) {
988 printf("ERROR: Could not receive track %d\n", jTracks);
989 continue;
990 }
5e65985c 991 if(itrack==jTracks)continue;
992
f09766dd 993 Double_t dEdxAsso = -999., ptPrim=-999., ptAsso=-999., openingAngle = -999.;
bfc7c23b 994 Double_t mass=999., width = -999;
995 Bool_t fFlagLS=kFALSE, fFlagULS=kFALSE;
996
f09766dd 997 ptPrim = track->Pt();
998
bfc7c23b 999 dEdxAsso = trackAsso->GetTPCsignal();
1000 ptAsso = trackAsso->Pt();
1001 Int_t chargeAsso = trackAsso->Charge();
1002 Int_t charge = track->Charge();
1003
b12bc641 1004
3db00c72 1005 if(ptAsso <0.5) continue;
bfc7c23b 1006 if(!fTrackCuts->AcceptTrack(trackAsso)) continue;
f09766dd 1007 if(dEdxAsso <65 || dEdxAsso>100) continue; //11a pass1
bfc7c23b 1008
1009 Int_t fPDGe1 = 11; Int_t fPDGe2 = 11;
1010 if(charge>0) fPDGe1 = -11;
1011 if(chargeAsso>0) fPDGe2 = -11;
b12bc641 1012
1013 //printf("chargeAsso = %d\n",chargeAsso);
1014 //printf("charge = %d\n",charge);
bfc7c23b 1015 if(charge == chargeAsso) fFlagLS = kTRUE;
1016 if(charge != chargeAsso) fFlagULS = kTRUE;
1017
b12bc641 1018 //printf("fFlagLS = %d\n",fFlagLS);
1019 //printf("fFlagULS = %d\n",fFlagULS);
1020 //printf("\n");
1021
bfc7c23b 1022 AliKFParticle ge1(*track, fPDGe1);
1023 AliKFParticle ge2(*trackAsso, fPDGe2);
1024 AliKFParticle recg(ge1, ge2);
1025
1026 if(recg.GetNDF()<1) continue;
1027 Double_t chi2recg = recg.GetChi2()/recg.GetNDF();
1028 if(TMath::Sqrt(TMath::Abs(chi2recg))>3.) continue;
1029
1030 AliKFVertex primV(*pVtx);
1031 primV += recg;
1032 recg.SetProductionVertex(primV);
1033
d113d7cd 1034 recg.SetMassConstraint(0,0.0001);
bfc7c23b 1035
1036 openingAngle = ge1.GetAngle(ge2);
1037 if(fFlagLS) fOpeningAngleLS->Fill(openingAngle);
1038 if(fFlagULS) fOpeningAngleULS->Fill(openingAngle);
1039
bfc7c23b 1040
1041 recg.GetMass(mass,width);
1042
c2801a73 1043 double ishower = 0;
1044 if(shower>0.0 && shower<0.3)ishower = 1;
1045
b12bc641 1046 double phoinfo[9];
f09766dd 1047 phoinfo[0] = cent;
1048 phoinfo[1] = ptPrim;
1049 phoinfo[2] = mass;
3db00c72 1050 phoinfo[3] = nSig;
4e01b68c 1051 phoinfo[4] = openingAngle;
c2801a73 1052 phoinfo[5] = ishower;
1053 phoinfo[6] = ep;
1054 phoinfo[7] = mce;
b12bc641 1055 phoinfo[8] = ptAsso;
f09766dd 1056
1057 if(fFlagLS) fInvmassLS->Fill(phoinfo);
1058 if(fFlagULS) fInvmassULS->Fill(phoinfo);
ffa14dda 1059
b12bc641 1060 //printf("fInvmassCut %f\n",fInvmassCut);
1061 //printf("openingAngle %f\n",fOpeningAngleCut);
1062
ffa14dda 1063 if(openingAngle > fOpeningAngleCut) continue;
1064
bfc7c23b 1065 if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec){
1066 flagPhotonicElec = kTRUE;
1067 }
f09766dd 1068 if(mass<fInvmassCut && fFlagLS && !flagConvinatElec){
1069 flagConvinatElec = kTRUE;
1070 }
bfc7c23b 1071
1072 }
1073 fFlagPhotonicElec = flagPhotonicElec;
f09766dd 1074 fFlagConvinatElec = flagConvinatElec;
bfc7c23b 1075
1076}
1077
fb87d707 1078
1079//_________________________________________
1080void AliAnalysisTaskHFECal::FindTriggerClusters()
1081{
1082 // constants
1083 const int nModuleCols = 2;
1084 const int nModuleRows = 5;
1085 const int nColsFeeModule = 48;
1086 const int nRowsFeeModule = 24;
1087 const int nColsFaltroModule = 24;
1088 const int nRowsFaltroModule = 12;
1089 //const int faltroWidthMax = 20;
1090
1091 // part 1, trigger extraction -------------------------------------
1092 Int_t globCol, globRow;
1093 //Int_t ntimes=0, nTrigChannel=0, nTrigChannelCut=0, trigInCut=0;
1094 Int_t ntimes=0, nTrigChannel=0, nTrigChannelCut=0;
1095
1096 //Int_t trigtimes[faltroWidthMax];
1097 Double_t cellTime[nColsFeeModule*nModuleCols][nRowsFeeModule*nModuleRows];
1098 Double_t cellEnergy[nColsFeeModule*nModuleCols][nRowsFeeModule*nModuleRows];
1099 //Double_t fTrigCutLow = 6;
1100 //Double_t fTrigCutHigh = 10;
1101 Double_t fTimeCutLow = 469e-09;
1102 Double_t fTimeCutHigh = 715e-09;
1103
1104 AliESDCaloTrigger * fCaloTrigger = fESD->GetCaloTrigger( "EMCAL" );
1105
1106 // erase trigger maps
1107 for(Int_t i = 0; i < nColsFaltroModule*nModuleCols; i++ )
1108 {
1109 for(Int_t j = 0; j < nRowsFaltroModule*nModuleRows; j++ )
1110 {
1111 ftriggersCut[i][j] = 0;
1112 ftriggers[i][j] = 0;
1113 ftriggersTime[i][j] = 0;
1114 }
1115 }
1116
1117 Int_t iglobCol=0, iglobRow=0;
1118 // go through triggers
1119 if( fCaloTrigger->GetEntries() > 0 )
1120 {
1121 // needs reset
1122 fCaloTrigger->Reset();
1123 while( fCaloTrigger->Next() )
1124 {
1125 fCaloTrigger->GetPosition( globCol, globRow );
1126 fCaloTrigger->GetNL0Times( ntimes );
1127 /*
1128 // no L0s
1129 if( ntimes < 1 ) continue;
1130 // get precise timings
1131 fCaloTrigger->GetL0Times( trigtimes );
1132 trigInCut = 0;
1133 for(Int_t i = 0; i < ntimes; i++ )
1134 {
1135 // save the first trigger time in channel
1136 if( i == 0 || triggersTime[globCol][globRow] > trigtimes[i] )
1137 triggersTime[globCol][globRow] = trigtimes[i];
1138 //printf("trigger times: %d\n",trigtimes[i]);
1139 // check if in cut
1140 if(trigtimes[i] > fTrigCutLow && trigtimes[i] < fTrigCutHigh )
1141 trigInCut = 1;
1142
1143 fTrigTimes->Fill(trigtimes[i]);
1144 }
1145 */
1146
1147 //L1 analysis from AliAnalysisTaskEMCALTriggerQA
1148 Int_t bit = 0;
1149 fCaloTrigger->GetTriggerBits(bit);
1150
1151 Int_t ts = 0;
1152 fCaloTrigger->GetL1TimeSum(ts);
1153 if (ts > 0)ftriggers[globCol][globRow] = 1;
1154 // number of triggered channels in event
1155 nTrigChannel++;
1156 // ... inside cut
1157 if(ts>0 && (bit >> 6 & 0x1))
1158 {
1159 iglobCol = globCol;
1160 iglobRow = globRow;
1161 nTrigChannelCut++;
1162 ftriggersCut[globCol][globRow] = 1;
1163 }
1164
1165 } // calo trigger entries
1166 } // has calo trigger entries
1167
1168 // part 2 go through the clusters here -----------------------------------
1169 Int_t nCluster=0, nCell=0, iCell=0, gCell=0;
1170 Short_t cellAddr, nSACell, mclabel;
1171 //Int_t nSACell, iSACell, mclabel;
1172 Int_t iSACell;
42c75692 1173 Double_t cellAmp=0, cellTimeT=0, clusterTime=0, efrac=0;
fb87d707 1174 Int_t nSupMod, nModule, nIphi, nIeta, iphi, ieta, gphi, geta, feta, fphi;
1175
1176 TRefArray *fCaloClusters = new TRefArray();
1177 fESD->GetEMCALClusters( fCaloClusters );
1178 nCluster = fCaloClusters->GetEntries();
1179
1180
1181 // save all cells times if there are clusters
1182 if( nCluster > 0 ){
1183 // erase time map
1184 for(Int_t i = 0; i < nColsFeeModule*nModuleCols; i++ ){
1185 for(Int_t j = 0; j < nRowsFeeModule*nModuleRows; j++ ){
1186 cellTime[i][j] = 0.;
1187 cellEnergy[i][j] = 0.;
1188 }
1189 }
1190
1191 // get cells
1192 AliESDCaloCells *fCaloCells = fESD->GetEMCALCells();
1193 //AliVCaloCells fCaloCells = fESD->GetEMCALCells();
1194 nSACell = fCaloCells->GetNumberOfCells();
1195 for(iSACell = 0; iSACell < nSACell; iSACell++ ){
1196 // get the cell info *fCal
1197 fCaloCells->GetCell( iSACell, cellAddr, cellAmp, cellTimeT , mclabel, efrac);
1198
1199 // get cell position
1200 fGeom->GetCellIndex( cellAddr, nSupMod, nModule, nIphi, nIeta );
1201 fGeom->GetCellPhiEtaIndexInSModule( nSupMod,nModule, nIphi, nIeta, iphi, ieta);
1202
1203 // convert co global phi eta
1204 gphi = iphi + nRowsFeeModule*(nSupMod/2);
1205 geta = ieta + nColsFeeModule*(nSupMod%2);
1206
1207 // save cell time and energy
1208 cellTime[geta][gphi] = cellTimeT;
1209 cellEnergy[geta][gphi] = cellAmp;
1210
1211 }
1212 }
1213
1214 Int_t nClusterTrig, nClusterTrigCut;
1215 UShort_t *cellAddrs;
1216 Double_t clsE=-999, clsEta=-999, clsPhi=-999;
1217 Float_t clsPos[3] = {0.,0.,0.};
1218
1219 for(Int_t icl=0; icl<fESD->GetNumberOfCaloClusters(); icl++)
1220 {
1221 AliESDCaloCluster *cluster = fESD->GetCaloCluster(icl);
1222 if(!cluster || !cluster->IsEMCAL()) continue;
1223
1224 // get cluster cells
1225 nCell = cluster->GetNCells();
1226
1227 // get cluster energy
1228 clsE = cluster->E();
1229
1230 // get cluster position
1231 cluster->GetPosition(clsPos);
1232 TVector3 clsPosVec(clsPos[0],clsPos[1],clsPos[2]);
1233 clsEta = clsPosVec.Eta();
1234 clsPhi = clsPosVec.Phi();
1235
1236 // get the cell addresses
1237 cellAddrs = cluster->GetCellsAbsId();
1238
1239 // check if the cluster contains cell, that was marked as triggered
1240 nClusterTrig = 0;
1241 nClusterTrigCut = 0;
1242
1243 // loop the cells to check, if cluser in acceptance
1244 // any cluster with a cell outside acceptance is not considered
1245 for( iCell = 0; iCell < nCell; iCell++ )
1246 {
42c75692 1247 // check hot cell
feffe705 1248 //if(clsE>6.0)fCellCheck->Fill(clsE,cellAddrs[iCell]);
42c75692 1249
fb87d707 1250 // get cell position
1251 fGeom->GetCellIndex( cellAddrs[iCell], nSupMod, nModule, nIphi, nIeta );
1252 fGeom->GetCellPhiEtaIndexInSModule( nSupMod,nModule, nIphi, nIeta, iphi, ieta);
1253
1254 // convert co global phi eta
1255 gphi = iphi + nRowsFeeModule*(nSupMod/2);
1256 geta = ieta + nColsFeeModule*(nSupMod%2);
1257
1258 if( cellTime[geta][gphi] > 0. ){
1259 clusterTime += cellTime[geta][gphi];
1260 gCell++;
1261 }
1262
1263 // get corresponding FALTRO
1264 fphi = gphi / 2;
1265 feta = geta / 2;
1266
1267 // try to match with a triggered
1268 if( ftriggers[feta][fphi]==1)
1269 { nClusterTrig++;
1270 }
1271 if( ftriggersCut[feta][fphi]==1)
1272 { nClusterTrigCut++;
1273 }
1274
1275 } // cells
1276
1277
1278 if( gCell > 0 )
1279 clusterTime = clusterTime / (Double_t)gCell;
1280 // fix the reconstriction code time 100ns jumps
1281 if( fESD->GetBunchCrossNumber() % 4 < 2 )
1282 clusterTime -= 0.0000001;
1283
feffe705 1284 //fClsETime->Fill(clsE,clusterTime);
1285 //fClsEBftTrigCut->Fill(clsE);
fb87d707 1286
1287 if(nClusterTrig>0){
feffe705 1288 //fClsETime1->Fill(clsE,clusterTime);
fb87d707 1289 }
1290
1291 if(nClusterTrig>0){
1292 cluster->SetChi2(1);
feffe705 1293 //fClsEAftTrigCut1->Fill(clsE);
fb87d707 1294 }
1295
1296 if(nClusterTrigCut>0){
1297 cluster->SetChi2(2);
feffe705 1298 //fClsEAftTrigCut2->Fill(clsE);
fb87d707 1299 }
1300
1301 if(nClusterTrigCut>0 && ( nCell > (1 + clsE / 3)))
1302 {
1303 cluster->SetChi2(3);
feffe705 1304 //fClsEAftTrigCut3->Fill(clsE);
fb87d707 1305 }
1306
1307 if(nClusterTrigCut>0 && (nCell > (1 + clsE / 3) )&&( clusterTime > fTimeCutLow && clusterTime < fTimeCutHigh ))
1308 {
1309 // cluster->SetChi2(4);
feffe705 1310 //fClsEAftTrigCut4->Fill(clsE);
fb87d707 1311 }
1312 if(nClusterTrigCut<1)
1313 {
1314 cluster->SetChi2(0);
1315
feffe705 1316 //fClsEAftTrigCut->Fill(clsE);
fb87d707 1317 }
1318
1319 } // clusters
1320}
1321
1322
1323
1324
1325
1326
1327