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