updates from Claude
[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
38#include "AliAnalysisTaskHFECal.h"
39#include "TGeoGlobalMagField.h"
40#include "AliLog.h"
41#include "AliAnalysisTaskSE.h"
42#include "TRefArray.h"
43#include "TVector.h"
44#include "AliESDInputHandler.h"
45#include "AliESDpid.h"
46#include "AliESDtrackCuts.h"
47#include "AliPhysicsSelection.h"
48#include "AliESDCaloCluster.h"
49#include "AliAODCaloCluster.h"
50#include "AliEMCALRecoUtils.h"
51#include "AliEMCALGeometry.h"
52#include "AliGeomManager.h"
53#include "stdio.h"
54#include "TGeoManager.h"
55#include "iostream"
56#include "fstream"
57
58#include "AliEMCALTrack.h"
59#include "AliMagF.h"
60
61#include "AliKFParticle.h"
62#include "AliKFVertex.h"
63
64#include "AliPID.h"
65#include "AliPIDResponse.h"
66#include "AliHFEcontainer.h"
67#include "AliHFEcuts.h"
68#include "AliHFEpid.h"
69#include "AliHFEpidBase.h"
70#include "AliHFEpidQAmanager.h"
71#include "AliHFEtools.h"
72#include "AliCFContainer.h"
73#include "AliCFManager.h"
74
75#include "AliCentrality.h"
76
77ClassImp(AliAnalysisTaskHFECal)
78//________________________________________________________________________
79AliAnalysisTaskHFECal::AliAnalysisTaskHFECal(const char *name)
80 : AliAnalysisTaskSE(name)
81 ,fESD(0)
82 ,fOutputList(0)
83 ,fTrackCuts(0)
84 ,fCuts(0)
85 ,fIdentifiedAsOutInz(kFALSE)
86 ,fPassTheEventCut(kFALSE)
87 ,fRejectKinkMother(kFALSE)
88 ,fVz(0.0)
89 ,fCFM(0)
90 ,fPID(0)
91 ,fPIDqa(0)
92 ,fOpeningAngleCut(0.1)
93 ,fInvmassCut(0.01)
94 ,fNoEvents(0)
95 ,fEMCAccE(0)
96 ,fTrkpt(0)
97 ,fTrkEovPBef(0)
98 ,fTrkEovPAft(0)
99 ,fdEdxBef(0)
100 ,fdEdxAft(0)
101 ,fIncpT(0)
102 ,fInvmassLS(0)
103 ,fInvmassULS(0)
104 ,fOpeningAngleLS(0)
105 ,fOpeningAngleULS(0)
106 ,fPhotoElecPt(0)
107 ,fPhoElecPt(0)
108 ,fTrackPtBefTrkCuts(0)
109 ,fTrackPtAftTrkCuts(0)
110 ,fTPCnsigma(0)
111 ,fCent(0)
112 ,fEleInfo(0)
113{
114 //Named constructor
115
116 fPID = new AliHFEpid("hfePid");
117 fTrackCuts = new AliESDtrackCuts();
118
119 // Define input and output slots here
120 // Input slot #0 works with a TChain
121 DefineInput(0, TChain::Class());
122 // Output slot #0 id reserved by the base class for AOD
123 // Output slot #1 writes into a TH1 container
124 // DefineOutput(1, TH1I::Class());
125 DefineOutput(1, TList::Class());
126 // DefineOutput(3, TTree::Class());
127}
128
129//________________________________________________________________________
130AliAnalysisTaskHFECal::AliAnalysisTaskHFECal()
131 : AliAnalysisTaskSE("DefaultAnalysis_AliAnalysisTaskHFECal")
132 ,fESD(0)
133 ,fOutputList(0)
134 ,fTrackCuts(0)
135 ,fCuts(0)
136 ,fIdentifiedAsOutInz(kFALSE)
137 ,fPassTheEventCut(kFALSE)
138 ,fRejectKinkMother(kFALSE)
139 ,fVz(0.0)
140 ,fCFM(0)
141 ,fPID(0)
142 ,fPIDqa(0)
143 ,fOpeningAngleCut(0.1)
144 ,fInvmassCut(0.01)
145 ,fNoEvents(0)
146 ,fEMCAccE(0)
147 ,fTrkpt(0)
148 ,fTrkEovPBef(0)
149 ,fTrkEovPAft(0)
150 ,fdEdxBef(0)
151 ,fdEdxAft(0)
152 ,fIncpT(0)
153 ,fInvmassLS(0)
154 ,fInvmassULS(0)
155 ,fOpeningAngleLS(0)
156 ,fOpeningAngleULS(0)
157 ,fPhotoElecPt(0)
158 ,fPhoElecPt(0)
159 ,fTrackPtBefTrkCuts(0)
160 ,fTrackPtAftTrkCuts(0)
161 ,fTPCnsigma(0)
162 ,fCent(0)
163 ,fEleInfo(0)
164{
165 //Default constructor
166 fPID = new AliHFEpid("hfePid");
167
168 fTrackCuts = new AliESDtrackCuts();
169
170 // Constructor
171 // Define input and output slots here
172 // Input slot #0 works with a TChain
173 DefineInput(0, TChain::Class());
174 // Output slot #0 id reserved by the base class for AOD
175 // Output slot #1 writes into a TH1 container
176 // DefineOutput(1, TH1I::Class());
177 DefineOutput(1, TList::Class());
178 //DefineOutput(3, TTree::Class());
179}
180//_________________________________________
181
182AliAnalysisTaskHFECal::~AliAnalysisTaskHFECal()
183{
184 //Destructor
185
186 delete fOutputList;
187 delete fPID;
188 delete fCFM;
189 delete fPIDqa;
190 delete fTrackCuts;
191}
192//_________________________________________
193
194void AliAnalysisTaskHFECal::UserExec(Option_t*)
195{
196 //Main loop
197 //Called for each event
198
199 // create pointer to event
200 fESD = dynamic_cast<AliESDEvent*>(InputEvent());
201 if (!fESD) {
202 printf("ERROR: fESD not available\n");
203 return;
204 }
205
206 if(!fCuts){
207 AliError("HFE cuts not available");
208 return;
209 }
210
211 if(!fPID->IsInitialized()){
212 // Initialize PID with the given run number
213 AliWarning("PID not initialised, get from Run no");
214 fPID->InitializePID(fESD->GetRunNumber());
215 }
216
217 fNoEvents->Fill(0);
218
219 Int_t fNOtrks = fESD->GetNumberOfTracks();
220 const AliESDVertex *pVtx = fESD->GetPrimaryVertex();
221
222 Double_t pVtxZ = -999;
223 pVtxZ = pVtx->GetZ();
224
225 if(TMath::Abs(pVtxZ)>10) return;
226 fNoEvents->Fill(1);
227
228 if(fNOtrks<2) return;
229 fNoEvents->Fill(2);
230
231 AliPIDResponse *pidResponse = fInputHandler->GetPIDResponse();
232 if(!pidResponse){
233 AliDebug(1, "Using default PID Response");
234 pidResponse = AliHFEtools::GetDefaultPID(kFALSE, fInputEvent->IsA() == AliAODEvent::Class());
235 }
236
237 fPID->SetPIDResponse(pidResponse);
238
239 fCFM->SetRecEventInfo(fESD);
240
241 Float_t cent = -1.;
242 AliCentrality *centrality = fESD->GetCentrality();
243 cent = centrality->GetCentralityPercentile("V0M");
244 fCent->Fill(cent);
245
246 if(cent>90.) return;
247
248 // Calorimeter info.
249
250 // make EMCAL array
251 for(Int_t iCluster=0; iCluster<fESD->GetNumberOfCaloClusters(); iCluster++)
252 {
253 AliESDCaloCluster *clust = fESD->GetCaloCluster(iCluster);
254 if(clust->IsEMCAL())
255 {
256 double clustE = clust->E();
257 float emcx[3]; // cluster pos
258 clust->GetPosition(emcx);
259 TVector3 clustpos(emcx[0],emcx[1],emcx[2]);
260 double emcphi = clustpos.Phi();
261 double emceta = clustpos.Eta();
5e65985c 262 double calInfo[4];
bfc7c23b 263 calInfo[0] = emcphi; calInfo[1] = emceta; calInfo[2] = clustE; calInfo[3] = cent;
264 fEMCAccE->Fill(calInfo);
265 }
266 }
267
268 // Track loop
269 for (Int_t iTracks = 0; iTracks < fESD->GetNumberOfTracks(); iTracks++) {
270 AliESDtrack* track = fESD->GetTrack(iTracks);
271 if (!track) {
272 printf("ERROR: Could not receive track %d\n", iTracks);
273 continue;
274 }
275
276 if(TMath::Abs(track->Eta())>0.7) continue;
277
278 fTrackPtBefTrkCuts->Fill(track->Pt());
279 // RecKine: ITSTPC cuts
280 if(!ProcessCutStep(AliHFEcuts::kStepRecKineITSTPC, track)) continue;
281
282 //RecKink
283 if(fRejectKinkMother) { // Quick and dirty fix to reject both kink mothers and daughters
284 if(track->GetKinkIndex(0) != 0) continue;
285 }
286
287 // RecPrim
288 if(!ProcessCutStep(AliHFEcuts::kStepRecPrim, track)) continue;
289
290 // HFEcuts: ITS layers cuts
291 if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsITS, track)) continue;
292
293
294 fTrackPtAftTrkCuts->Fill(track->Pt());
295
296 Double_t mom = -999., eop=-999., pt = -999., dEdx=-999., fTPCnSigma=-10, phi=-999., eta=-999.;
297
298 // Track extrapolation
299
300 pt = track->Pt();
301 if(pt<2.0)continue;
302 Int_t charge = track->Charge();
303 fTrkpt->Fill(pt);
304 mom = track->P();
305 phi = track->Phi();
306 eta = track->Eta();
307 dEdx = track->GetTPCsignal();
308 fTPCnSigma = fPID->GetPIDResponse() ? fPID->GetPIDResponse()->NumberOfSigmasTPC(track, AliPID::kElectron) : 1000;
309
310 double ncells = -1.0;
311 double m20 = -1.0;
312 double m02 = -1.0;
313 double disp = -1.0;
314 double rmatch = -1.0;
315 double nmatch = -1.0;
316
317
318 Int_t clsId = track->GetEMCALcluster();
319 if (clsId>0){
320 AliESDCaloCluster *clust = fESD->GetCaloCluster(clsId);
321 if(clust && clust->IsEMCAL()){
322
323 double clustE = clust->E();
324 eop = clustE/fabs(mom);
325 //double clustT = clust->GetTOF();
326 ncells = clust->GetNCells();
327 m02 = clust->GetM02();
328 m20 = clust->GetM20();
329 disp = clust->GetDispersion();
330 double delphi = clust->GetTrackDx();
331 double deleta = clust->GetTrackDz();
332 rmatch = sqrt(pow(delphi,2)+pow(deleta,2));
333 nmatch = clust->GetNTracksMatched();
334
335 double valdedx[14];
336 valdedx[0] = mom; valdedx[1] = pt; valdedx[2] = dEdx; valdedx[3] = phi; valdedx[4] = eta; valdedx[5] = fTPCnSigma;
337 valdedx[6] = eop; valdedx[7] = rmatch; valdedx[8] = ncells, valdedx[9] = m02; valdedx[10] = m20; valdedx[11] = disp;
338 valdedx[12] = cent; valdedx[13] = charge;
339 fEleInfo->Fill(valdedx);
340
341
342 }
343 }
344
345 fdEdxBef->Fill(mom,dEdx);
346 fTPCnsigma->Fill(mom,fTPCnSigma);
347
348
349 // HFE cuts: TPC PID cleanup
350 if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTPC, track)) continue;
351
352 if(fTPCnSigma >= 1.0 && fTPCnSigma <= 3)fTrkEovPBef->Fill(pt,eop);
353 Int_t pidpassed = 1;
354
355
356 //--- track accepted
357 AliHFEpidObject hfetrack;
358 hfetrack.SetAnalysisType(AliHFEpidObject::kESDanalysis);
359 hfetrack.SetRecTrack(track);
360 hfetrack.SetPbPb();
361 if(!fPID->IsSelected(&hfetrack, NULL, "", fPIDqa)) pidpassed = 0;
362
363 if(pidpassed==0) continue;
364
365 fTrkEovPAft->Fill(pt,eop);
366 fdEdxAft->Fill(mom,dEdx);
367 fIncpT->Fill(cent,pt);
368
369 Bool_t fFlagPhotonicElec = kFALSE;
370 SelectPhotonicElectron(iTracks,cent,track,fFlagPhotonicElec);
371
372 if(fFlagPhotonicElec) fPhoElecPt->Fill(cent,pt);
373 }
374 PostData(1, fOutputList);
375}
376//_________________________________________
377void AliAnalysisTaskHFECal::UserCreateOutputObjects()
378{
379 //--------Initialize PID
380 fPID->SetHasMCData(kFALSE);
381 if(!fPID->GetNumberOfPIDdetectors())
382 {
383 fPID->AddDetector("TPC", 0);
384 fPID->AddDetector("EMCAL", 1);
385 }
386
387 fPID->SortDetectors();
388 fPIDqa = new AliHFEpidQAmanager();
389 fPIDqa->Initialize(fPID);
390
391 //--------Initialize correction Framework and Cuts
392 fCFM = new AliCFManager;
393 const Int_t kNcutSteps = AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kNcutStepsRecTrack + AliHFEcuts::kNcutStepsDETrack;
394 fCFM->SetNStepParticle(kNcutSteps);
395 for(Int_t istep = 0; istep < kNcutSteps; istep++)
396 fCFM->SetParticleCutsList(istep, NULL);
397
398 if(!fCuts){
399 AliWarning("Cuts not available. Default cuts will be used");
400 fCuts = new AliHFEcuts;
401 fCuts->CreateStandardCuts();
402 }
403 fCuts->Initialize(fCFM);
404
405 //---------Output Tlist
406 fOutputList = new TList();
407 fOutputList->SetOwner();
408 fOutputList->Add(fPIDqa->MakeList("PIDQA"));
409
410 fNoEvents = new TH1F("fNoEvents","",4,-0.5,3.5) ;
411 fOutputList->Add(fNoEvents);
412
413 Int_t binsE[4] = {250, 100, 1000, 100};
414 Double_t xminE[4] = {1.0, -1, 0.0, 0};
415 Double_t xmaxE[4] = {3.5, 1, 100.0, 100};
416 fEMCAccE = new THnSparseD("fEMCAccE","EMC acceptance & E;#eta;#phi;Energy;Centrality",4,binsE,xminE,xmaxE);
417 fOutputList->Add(fEMCAccE);
418
419 fTrkpt = new TH1F("fTrkpt","track pt",100,0,50);
420 fOutputList->Add(fTrkpt);
421
422 fTrackPtBefTrkCuts = new TH1F("fTrackPtBefTrkCuts","track pt before track cuts",100,0,50);
423 fOutputList->Add(fTrackPtBefTrkCuts);
424
425 fTrackPtAftTrkCuts = new TH1F("fTrackPtAftTrkCuts","track pt after track cuts",100,0,50);
426 fOutputList->Add(fTrackPtAftTrkCuts);
427
428 fTPCnsigma = new TH2F("fTPCnsigma", "TPC - n sigma",100,0,50,200,-10,10);
429 fOutputList->Add(fTPCnsigma);
430
431 fTrkEovPBef = new TH2F("fTrkEovPBef","track E/p before HFE pid",100,0,50,100,0,2);
432 fOutputList->Add(fTrkEovPBef);
433
434 fTrkEovPAft = new TH2F("fTrkEovPAft","track E/p after HFE pid",100,0,50,100,0,2);
435 fOutputList->Add(fTrkEovPAft);
436
437 fdEdxBef = new TH2F("fdEdxBef","track dEdx vs p before HFE pid",100,0,50,150,0,150);
438 fOutputList->Add(fdEdxBef);
439
440 fdEdxAft = new TH2F("fdEdxAft","track dEdx vs p after HFE pid",100,0,50,150,0,150);
441 fOutputList->Add(fdEdxAft);
442
443 fIncpT = new TH2F("fIncpT","HFE pid electro vs. centrality",100,0,100,100,0,50);
444 fOutputList->Add(fIncpT);
445
446 fInvmassLS = new TH2F("fInvmassLS", "Inv mass of LS (e,e); mass(GeV/c^2); counts;", 100,0,100,500,0,0.5);
447 fOutputList->Add(fInvmassLS);
448
449 fInvmassULS = new TH2F("fInvmassULS", "Inv mass of ULS (e,e); mass(GeV/c^2); counts;", 100,0,100,500,0,0.5);
450 fOutputList->Add(fInvmassULS);
451
452 fOpeningAngleLS = new TH1F("fOpeningAngleLS","Opening angle for LS pairs",100,0,1);
453 fOutputList->Add(fOpeningAngleLS);
454
455 fOpeningAngleULS = new TH1F("fOpeningAngleULS","Opening angle for ULS pairs",100,0,1);
456 fOutputList->Add(fOpeningAngleULS);
457
458 fPhotoElecPt = new TH1F("fPhotoElecPt", "photonic electron pt",100,0,50);
459 fOutputList->Add(fPhotoElecPt);
460
461 fPhoElecPt = new TH2F("fPhoElecPt", "Semi-inclusive electron pt",100,0,100,100,0,50);
462 fOutputList->Add(fPhoElecPt);
463
464 fCent = new TH1F("fCent","Centrality",100,0,100) ;
465 fOutputList->Add(fCent);
466
467 // Make common binning
468 const Double_t kMinP = 0.;
469 const Double_t kMaxP = 50.;
470 const Double_t kTPCSigMim = 40.;
471 const Double_t kTPCSigMax = 140.;
472
473 // 1st histogram: TPC dEdx with/without EMCAL (p, pT, TPC Signal, phi, eta, Sig, e/p, ,match, cell, M02, M20, Disp, Centrality, select)
474 Int_t nBins[14] = { 500, 500, 200, 60, 20, 600, 300, 100, 40, 200, 200, 200, 100, 3};
475 Double_t min[14] = {kMinP, kMinP, kTPCSigMim, 1.0, -1.0, -8.0, 0, 0, 0, 0.0, 0.0, 0.0, 0, -1.5};
476 Double_t max[14] = {kMaxP, kMaxP, kTPCSigMax, 4.0, 1.0, 4.0, 3.0, 0.1, 40, 2.0, 2.0, 2.0, 100, 1.5};
477 fEleInfo = new THnSparseD("fEleInfo", "p [GeV/c]; pT [GeV/c]; TPC signal;phi;eta;nSig; E/p;Rmatch;Ncell;M02;M20;Disp; Centrality; charge", 14, nBins, min, max);
478 fOutputList->Add(fEleInfo);
479
480//_________________________________________________________
481
482 PostData(1,fOutputList);
483}
484
485//________________________________________________________________________
486void AliAnalysisTaskHFECal::Terminate(Option_t *)
487{
488 // Info("Terminate");
489 AliAnalysisTaskSE::Terminate();
490}
491
492//________________________________________________________________________
493Bool_t AliAnalysisTaskHFECal::ProcessCutStep(Int_t cutStep, AliVParticle *track)
494{
495 // Check single track cuts for a given cut step
496 const Int_t kMCOffset = AliHFEcuts::kNcutStepsMCTrack;
497 if(!fCFM->CheckParticleCuts(cutStep + kMCOffset, track)) return kFALSE;
498 return kTRUE;
499}
500//_________________________________________
501void AliAnalysisTaskHFECal::SelectPhotonicElectron(Int_t itrack, Double_t cent, AliESDtrack *track, Bool_t &fFlagPhotonicElec)
502{
503 //Identify non-heavy flavour electrons using Invariant mass method
504
505 fTrackCuts->SetAcceptKinkDaughters(kFALSE);
506 fTrackCuts->SetRequireTPCRefit(kTRUE);
507 fTrackCuts->SetEtaRange(-0.7,0.7);
508 fTrackCuts->SetRequireSigmaToVertex(kTRUE);
509 fTrackCuts->SetMaxChi2PerClusterTPC(3.5);
510 fTrackCuts->SetMinNClustersTPC(100);
511
512 const AliESDVertex *pVtx = fESD->GetPrimaryVertex();
513
514 Bool_t flagPhotonicElec = kFALSE;
515
5e65985c 516 //for(Int_t jTracks = itrack+1; jTracks<fESD->GetNumberOfTracks(); jTracks++){
517 for(Int_t jTracks = 0; jTracks<fESD->GetNumberOfTracks(); jTracks++){
bfc7c23b 518 AliESDtrack* trackAsso = fESD->GetTrack(jTracks);
519 if (!trackAsso) {
520 printf("ERROR: Could not receive track %d\n", jTracks);
521 continue;
522 }
5e65985c 523 if(itrack==jTracks)continue;
524
bfc7c23b 525 Double_t dEdxAsso = -999., ptAsso=-999., openingAngle = -999.;
526 Double_t mass=999., width = -999;
527 Bool_t fFlagLS=kFALSE, fFlagULS=kFALSE;
528
529 dEdxAsso = trackAsso->GetTPCsignal();
530 ptAsso = trackAsso->Pt();
531 Int_t chargeAsso = trackAsso->Charge();
532 Int_t charge = track->Charge();
533
5e65985c 534 //if(ptAsso <0.3) continue;
535 if(ptAsso <0.5) continue;
bfc7c23b 536 if(!fTrackCuts->AcceptTrack(trackAsso)) continue;
537 if(dEdxAsso <70 || dEdxAsso>100) continue; //11a pass1
538
539 Int_t fPDGe1 = 11; Int_t fPDGe2 = 11;
540 if(charge>0) fPDGe1 = -11;
541 if(chargeAsso>0) fPDGe2 = -11;
542
543 if(charge == chargeAsso) fFlagLS = kTRUE;
544 if(charge != chargeAsso) fFlagULS = kTRUE;
545
546 AliKFParticle ge1(*track, fPDGe1);
547 AliKFParticle ge2(*trackAsso, fPDGe2);
548 AliKFParticle recg(ge1, ge2);
549
550 if(recg.GetNDF()<1) continue;
551 Double_t chi2recg = recg.GetChi2()/recg.GetNDF();
552 if(TMath::Sqrt(TMath::Abs(chi2recg))>3.) continue;
553
554 AliKFVertex primV(*pVtx);
555 primV += recg;
556 recg.SetProductionVertex(primV);
557
558 recg.SetMassConstraint(0,0.0001);
559
560 openingAngle = ge1.GetAngle(ge2);
561 if(fFlagLS) fOpeningAngleLS->Fill(openingAngle);
562 if(fFlagULS) fOpeningAngleULS->Fill(openingAngle);
563
564 if(openingAngle > fOpeningAngleCut) continue;
565
566 recg.GetMass(mass,width);
567
568 if(fFlagLS) fInvmassLS->Fill(cent,mass);
569 if(fFlagULS) fInvmassULS->Fill(cent,mass);
570
571 if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec){
572 flagPhotonicElec = kTRUE;
573 }
574
575 }
576 fFlagPhotonicElec = flagPhotonicElec;
577
578}
579