git://git.uio.no / u/mrichter/AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Updates
[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)
fb87d707 82 ,fGeom(0)
bfc7c23b 83 ,fOutputList(0)
84 ,fTrackCuts(0)
85 ,fCuts(0)
86 ,fIdentifiedAsOutInz(kFALSE)
87 ,fPassTheEventCut(kFALSE)
88 ,fRejectKinkMother(kFALSE)
89 ,fVz(0.0)
90 ,fCFM(0)
91 ,fPID(0)
92 ,fPIDqa(0)
93 ,fOpeningAngleCut(0.1)
94 ,fInvmassCut(0.01)
95 ,fNoEvents(0)
96 ,fEMCAccE(0)
97 ,fTrkpt(0)
98 ,fTrkEovPBef(0)
99 ,fTrkEovPAft(0)
100 ,fdEdxBef(0)
101 ,fdEdxAft(0)
102 ,fIncpT(0)
fb87d707 103 ,fIncpTM20(0)
bfc7c23b 104 ,fInvmassLS(0)
105 ,fInvmassULS(0)
106 ,fOpeningAngleLS(0)
107 ,fOpeningAngleULS(0)
108 ,fPhotoElecPt(0)
109 ,fPhoElecPt(0)
fb87d707 110 ,fPhoElecPtM20(0)
f09766dd 111 ,fSameElecPt(0)
fb87d707 112 ,fSameElecPtM20(0)
bfc7c23b 113 ,fTrackPtBefTrkCuts(0)
114 ,fTrackPtAftTrkCuts(0)
115 ,fTPCnsigma(0)
116 ,fCent(0)
117 ,fEleInfo(0)
fb87d707 118 ,fClsEBftTrigCut(0)
119 ,fClsEAftTrigCut(0)
120 ,fClsEAftTrigCut1(0)
121 ,fClsEAftTrigCut2(0)
122 ,fClsEAftTrigCut3(0)
123 ,fClsEAftTrigCut4(0)
124 ,fClsETime(0)
125 ,fClsETime1(0)
126 ,fTrigTimes(0)
bfc7c23b 127{
128 //Named constructor
129
130 fPID = new AliHFEpid("hfePid");
131 fTrackCuts = new AliESDtrackCuts();
132
133 // Define input and output slots here
134 // Input slot #0 works with a TChain
135 DefineInput(0, TChain::Class());
136 // Output slot #0 id reserved by the base class for AOD
137 // Output slot #1 writes into a TH1 container
138 // DefineOutput(1, TH1I::Class());
139 DefineOutput(1, TList::Class());
140 // DefineOutput(3, TTree::Class());
141}
142
143//________________________________________________________________________
144AliAnalysisTaskHFECal::AliAnalysisTaskHFECal()
145 : AliAnalysisTaskSE("DefaultAnalysis_AliAnalysisTaskHFECal")
146 ,fESD(0)
fb87d707 147 ,fGeom(0)
bfc7c23b 148 ,fOutputList(0)
149 ,fTrackCuts(0)
150 ,fCuts(0)
151 ,fIdentifiedAsOutInz(kFALSE)
152 ,fPassTheEventCut(kFALSE)
153 ,fRejectKinkMother(kFALSE)
154 ,fVz(0.0)
155 ,fCFM(0)
156 ,fPID(0)
157 ,fPIDqa(0)
158 ,fOpeningAngleCut(0.1)
159 ,fInvmassCut(0.01)
160 ,fNoEvents(0)
161 ,fEMCAccE(0)
162 ,fTrkpt(0)
163 ,fTrkEovPBef(0)
164 ,fTrkEovPAft(0)
165 ,fdEdxBef(0)
166 ,fdEdxAft(0)
167 ,fIncpT(0)
fb87d707 168 ,fIncpTM20(0)
bfc7c23b 169 ,fInvmassLS(0)
170 ,fInvmassULS(0)
171 ,fOpeningAngleLS(0)
172 ,fOpeningAngleULS(0)
173 ,fPhotoElecPt(0)
174 ,fPhoElecPt(0)
fb87d707 175 ,fPhoElecPtM20(0)
f09766dd 176 ,fSameElecPt(0)
fb87d707 177 ,fSameElecPtM20(0)
bfc7c23b 178 ,fTrackPtBefTrkCuts(0)
179 ,fTrackPtAftTrkCuts(0)
180 ,fTPCnsigma(0)
181 ,fCent(0)
182 ,fEleInfo(0)
fb87d707 183 ,fClsEBftTrigCut(0)
184 ,fClsEAftTrigCut(0)
185 ,fClsEAftTrigCut1(0)
186 ,fClsEAftTrigCut2(0)
187 ,fClsEAftTrigCut3(0)
188 ,fClsEAftTrigCut4(0)
189 ,fClsETime(0)
190 ,fClsETime1(0)
191 ,fTrigTimes(0)
bfc7c23b 192{
193 //Default constructor
194 fPID = new AliHFEpid("hfePid");
195
196 fTrackCuts = new AliESDtrackCuts();
197
198 // Constructor
199 // Define input and output slots here
200 // Input slot #0 works with a TChain
201 DefineInput(0, TChain::Class());
202 // Output slot #0 id reserved by the base class for AOD
203 // Output slot #1 writes into a TH1 container
204 // DefineOutput(1, TH1I::Class());
205 DefineOutput(1, TList::Class());
206 //DefineOutput(3, TTree::Class());
207}
208//_________________________________________
209
210AliAnalysisTaskHFECal::~AliAnalysisTaskHFECal()
211{
212 //Destructor
213
214 delete fOutputList;
fb87d707 215 delete fGeom;
bfc7c23b 216 delete fPID;
217 delete fCFM;
218 delete fPIDqa;
219 delete fTrackCuts;
220}
221//_________________________________________
222
223void AliAnalysisTaskHFECal::UserExec(Option_t*)
224{
225 //Main loop
226 //Called for each event
227
228 // create pointer to event
229 fESD = dynamic_cast<AliESDEvent*>(InputEvent());
230 if (!fESD) {
231 printf("ERROR: fESD not available\n");
232 return;
233 }
234
235 if(!fCuts){
236 AliError("HFE cuts not available");
237 return;
238 }
239
240 if(!fPID->IsInitialized()){
241 // Initialize PID with the given run number
242 AliWarning("PID not initialised, get from Run no");
243 fPID->InitializePID(fESD->GetRunNumber());
244 }
245
246 fNoEvents->Fill(0);
247
248 Int_t fNOtrks = fESD->GetNumberOfTracks();
249 const AliESDVertex *pVtx = fESD->GetPrimaryVertex();
250
251 Double_t pVtxZ = -999;
252 pVtxZ = pVtx->GetZ();
253
254 if(TMath::Abs(pVtxZ)>10) return;
255 fNoEvents->Fill(1);
256
257 if(fNOtrks<2) return;
258 fNoEvents->Fill(2);
259
260 AliPIDResponse *pidResponse = fInputHandler->GetPIDResponse();
261 if(!pidResponse){
262 AliDebug(1, "Using default PID Response");
263 pidResponse = AliHFEtools::GetDefaultPID(kFALSE, fInputEvent->IsA() == AliAODEvent::Class());
264 }
265
266 fPID->SetPIDResponse(pidResponse);
267
268 fCFM->SetRecEventInfo(fESD);
269
270 Float_t cent = -1.;
271 AliCentrality *centrality = fESD->GetCentrality();
272 cent = centrality->GetCentralityPercentile("V0M");
273 fCent->Fill(cent);
274
275 if(cent>90.) return;
276
277 // Calorimeter info.
278
279 // make EMCAL array
280 for(Int_t iCluster=0; iCluster<fESD->GetNumberOfCaloClusters(); iCluster++)
281 {
282 AliESDCaloCluster *clust = fESD->GetCaloCluster(iCluster);
283 if(clust->IsEMCAL())
284 {
285 double clustE = clust->E();
286 float emcx[3]; // cluster pos
287 clust->GetPosition(emcx);
288 TVector3 clustpos(emcx[0],emcx[1],emcx[2]);
289 double emcphi = clustpos.Phi();
290 double emceta = clustpos.Eta();
fb87d707 291 double calInfo[5];
292 calInfo[0] = emcphi; calInfo[1] = emceta; calInfo[2] = clustE; calInfo[3] = cent; calInfo[4] = clust->Chi2();
f09766dd 293 if(clustE>1.5)fEMCAccE->Fill(calInfo);
bfc7c23b 294 }
295 }
296
297 // Track loop
298 for (Int_t iTracks = 0; iTracks < fESD->GetNumberOfTracks(); iTracks++) {
299 AliESDtrack* track = fESD->GetTrack(iTracks);
300 if (!track) {
301 printf("ERROR: Could not receive track %d\n", iTracks);
302 continue;
303 }
304
305 if(TMath::Abs(track->Eta())>0.7) continue;
f09766dd 306 if(TMath::Abs(track->Pt()<2.0)) continue;
bfc7c23b 307
308 fTrackPtBefTrkCuts->Fill(track->Pt());
309 // RecKine: ITSTPC cuts
310 if(!ProcessCutStep(AliHFEcuts::kStepRecKineITSTPC, track)) continue;
311
312 //RecKink
313 if(fRejectKinkMother) { // Quick and dirty fix to reject both kink mothers and daughters
314 if(track->GetKinkIndex(0) != 0) continue;
315 }
316
317 // RecPrim
318 if(!ProcessCutStep(AliHFEcuts::kStepRecPrim, track)) continue;
319
320 // HFEcuts: ITS layers cuts
321 if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsITS, track)) continue;
322
f09766dd 323 // HFE cuts: TPC PID cleanup
324 if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTPC, track)) continue;
325
bfc7c23b 326
327 fTrackPtAftTrkCuts->Fill(track->Pt());
328
329 Double_t mom = -999., eop=-999., pt = -999., dEdx=-999., fTPCnSigma=-10, phi=-999., eta=-999.;
fb87d707 330 //pt = track->Pt();
331 //if(pt<2.0)continue;
bfc7c23b 332
333 // Track extrapolation
334
bfc7c23b 335 Int_t charge = track->Charge();
336 fTrkpt->Fill(pt);
337 mom = track->P();
338 phi = track->Phi();
339 eta = track->Eta();
340 dEdx = track->GetTPCsignal();
341 fTPCnSigma = fPID->GetPIDResponse() ? fPID->GetPIDResponse()->NumberOfSigmasTPC(track, AliPID::kElectron) : 1000;
342
343 double ncells = -1.0;
344 double m20 = -1.0;
345 double m02 = -1.0;
346 double disp = -1.0;
347 double rmatch = -1.0;
348 double nmatch = -1.0;
f09766dd 349 double oppstatus = 0.0;
350 double samestatus = 0.0;
bfc7c23b 351
f09766dd 352 Bool_t fFlagPhotonicElec = kFALSE;
353 Bool_t fFlagConvinatElec = kFALSE;
354 SelectPhotonicElectron(iTracks,cent,track,fFlagPhotonicElec,fFlagConvinatElec);
355 if(fFlagPhotonicElec)oppstatus = 1.0;
356 if(fFlagConvinatElec)samestatus = 1.0;
bfc7c23b 357
358 Int_t clsId = track->GetEMCALcluster();
359 if (clsId>0){
360 AliESDCaloCluster *clust = fESD->GetCaloCluster(clsId);
361 if(clust && clust->IsEMCAL()){
362
363 double clustE = clust->E();
364 eop = clustE/fabs(mom);
365 //double clustT = clust->GetTOF();
366 ncells = clust->GetNCells();
367 m02 = clust->GetM02();
368 m20 = clust->GetM20();
369 disp = clust->GetDispersion();
370 double delphi = clust->GetTrackDx();
371 double deleta = clust->GetTrackDz();
372 rmatch = sqrt(pow(delphi,2)+pow(deleta,2));
373 nmatch = clust->GetNTracksMatched();
374
f09766dd 375 double valdedx[16];
fb87d707 376 valdedx[0] = pt; valdedx[1] = dEdx; valdedx[2] = phi; valdedx[3] = eta; valdedx[4] = fTPCnSigma;
377 valdedx[5] = eop; valdedx[6] = rmatch; valdedx[7] = ncells, valdedx[8] = m02; valdedx[9] = m20; valdedx[10] = disp;
378 valdedx[11] = cent; valdedx[12] = charge; valdedx[13] = oppstatus; valdedx[14] = samestatus; valdedx[15] = clust->Chi2();
bfc7c23b 379 fEleInfo->Fill(valdedx);
380
381
382 }
383 }
384
385 fdEdxBef->Fill(mom,dEdx);
386 fTPCnsigma->Fill(mom,fTPCnSigma);
f09766dd 387 if(fTPCnSigma >= -1.0 && fTPCnSigma <= 3)fTrkEovPBef->Fill(pt,eop);
bfc7c23b 388
bfc7c23b 389 Int_t pidpassed = 1;
390
391
392 //--- track accepted
393 AliHFEpidObject hfetrack;
394 hfetrack.SetAnalysisType(AliHFEpidObject::kESDanalysis);
395 hfetrack.SetRecTrack(track);
fb87d707 396 int centf = (int)cent;
397 hfetrack.SetCentrality(centf); //added
bfc7c23b 398 hfetrack.SetPbPb();
399 if(!fPID->IsSelected(&hfetrack, NULL, "", fPIDqa)) pidpassed = 0;
400
401 if(pidpassed==0) continue;
402
403 fTrkEovPAft->Fill(pt,eop);
404 fdEdxAft->Fill(mom,dEdx);
bfc7c23b 405
fb87d707 406 fIncpT->Fill(cent,pt);
bfc7c23b 407 if(fFlagPhotonicElec) fPhoElecPt->Fill(cent,pt);
f09766dd 408 if(fFlagConvinatElec) fSameElecPt->Fill(cent,pt);
fb87d707 409
410 if(m20>0.0 && m20<0.3)
411 {
412 fIncpTM20->Fill(cent,pt);
413 if(fFlagPhotonicElec) fPhoElecPtM20->Fill(cent,pt);
414 if(fFlagConvinatElec) fSameElecPtM20->Fill(cent,pt);
415 }
bfc7c23b 416 }
417 PostData(1, fOutputList);
418}
419//_________________________________________
420void AliAnalysisTaskHFECal::UserCreateOutputObjects()
421{
fb87d707 422 //---- Geometry
423 fGeom = AliEMCALGeometry::GetInstance("EMCAL_COMPLETEV1");
424
bfc7c23b 425 //--------Initialize PID
426 fPID->SetHasMCData(kFALSE);
427 if(!fPID->GetNumberOfPIDdetectors())
428 {
429 fPID->AddDetector("TPC", 0);
430 fPID->AddDetector("EMCAL", 1);
431 }
432
433 fPID->SortDetectors();
434 fPIDqa = new AliHFEpidQAmanager();
435 fPIDqa->Initialize(fPID);
436
437 //--------Initialize correction Framework and Cuts
438 fCFM = new AliCFManager;
439 const Int_t kNcutSteps = AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kNcutStepsRecTrack + AliHFEcuts::kNcutStepsDETrack;
440 fCFM->SetNStepParticle(kNcutSteps);
441 for(Int_t istep = 0; istep < kNcutSteps; istep++)
442 fCFM->SetParticleCutsList(istep, NULL);
443
444 if(!fCuts){
445 AliWarning("Cuts not available. Default cuts will be used");
446 fCuts = new AliHFEcuts;
447 fCuts->CreateStandardCuts();
448 }
449 fCuts->Initialize(fCFM);
450
451 //---------Output Tlist
452 fOutputList = new TList();
453 fOutputList->SetOwner();
454 fOutputList->Add(fPIDqa->MakeList("PIDQA"));
455
456 fNoEvents = new TH1F("fNoEvents","",4,-0.5,3.5) ;
457 fOutputList->Add(fNoEvents);
458
fb87d707 459 Int_t binsE[5] = {250, 100, 1000, 100, 10};
460 Double_t xminE[5] = {1.0, -1, 0.0, 0, -0.5};
461 Double_t xmaxE[5] = {3.5, 1, 100.0, 100, 9.5};
462 fEMCAccE = new THnSparseD("fEMCAccE","EMC acceptance & E;#eta;#phi;Energy;Centrality;trugCondition;",5,binsE,xminE,xmaxE);
bfc7c23b 463 fOutputList->Add(fEMCAccE);
464
465 fTrkpt = new TH1F("fTrkpt","track pt",100,0,50);
466 fOutputList->Add(fTrkpt);
467
468 fTrackPtBefTrkCuts = new TH1F("fTrackPtBefTrkCuts","track pt before track cuts",100,0,50);
469 fOutputList->Add(fTrackPtBefTrkCuts);
470
471 fTrackPtAftTrkCuts = new TH1F("fTrackPtAftTrkCuts","track pt after track cuts",100,0,50);
472 fOutputList->Add(fTrackPtAftTrkCuts);
473
474 fTPCnsigma = new TH2F("fTPCnsigma", "TPC - n sigma",100,0,50,200,-10,10);
475 fOutputList->Add(fTPCnsigma);
476
477 fTrkEovPBef = new TH2F("fTrkEovPBef","track E/p before HFE pid",100,0,50,100,0,2);
478 fOutputList->Add(fTrkEovPBef);
479
480 fTrkEovPAft = new TH2F("fTrkEovPAft","track E/p after HFE pid",100,0,50,100,0,2);
481 fOutputList->Add(fTrkEovPAft);
482
483 fdEdxBef = new TH2F("fdEdxBef","track dEdx vs p before HFE pid",100,0,50,150,0,150);
484 fOutputList->Add(fdEdxBef);
485
486 fdEdxAft = new TH2F("fdEdxAft","track dEdx vs p after HFE pid",100,0,50,150,0,150);
487 fOutputList->Add(fdEdxAft);
488
489 fIncpT = new TH2F("fIncpT","HFE pid electro vs. centrality",100,0,100,100,0,50);
490 fOutputList->Add(fIncpT);
491
fb87d707 492 fIncpTM20 = new TH2F("fIncpTM20","HFE pid electro vs. centrality with M20",100,0,100,100,0,50);
493 fOutputList->Add(fIncpTM20);
f09766dd 494
495 Int_t nBinspho[3] = { 100, 100, 500};
496 Double_t minpho[3] = { 0., 0., 0.};
497 Double_t maxpho[3] = {100., 50., 0.5};
498
499 fInvmassLS = new THnSparseD("fInvmassLS", "Inv mass of LS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2);", 3, nBinspho,minpho, maxpho);
bfc7c23b 500 fOutputList->Add(fInvmassLS);
501
f09766dd 502 fInvmassULS = new THnSparseD("fInvmassULS", "Inv mass of ULS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2);", 3, nBinspho,minpho, maxpho);
bfc7c23b 503 fOutputList->Add(fInvmassULS);
504
505 fOpeningAngleLS = new TH1F("fOpeningAngleLS","Opening angle for LS pairs",100,0,1);
506 fOutputList->Add(fOpeningAngleLS);
507
508 fOpeningAngleULS = new TH1F("fOpeningAngleULS","Opening angle for ULS pairs",100,0,1);
509 fOutputList->Add(fOpeningAngleULS);
510
511 fPhotoElecPt = new TH1F("fPhotoElecPt", "photonic electron pt",100,0,50);
512 fOutputList->Add(fPhotoElecPt);
513
f09766dd 514 fPhoElecPt = new TH2F("fPhoElecPt", "Pho-inclusive electron pt",100,0,100,100,0,50);
bfc7c23b 515 fOutputList->Add(fPhoElecPt);
516
fb87d707 517 fPhoElecPtM20 = new TH2F("fPhoElecPtM20", "Pho-inclusive electron pt with M20",100,0,100,100,0,50);
518 fOutputList->Add(fPhoElecPtM20);
519
f09766dd 520 fSameElecPt = new TH2F("fSameElecPt", "Same-inclusive electron pt",100,0,100,100,0,50);
521 fOutputList->Add(fSameElecPt);
522
fb87d707 523 fSameElecPtM20 = new TH2F("fSameElecPtM20", "Same-inclusive electron pt with M20",100,0,100,100,0,50);
524 fOutputList->Add(fSameElecPtM20);
525
bfc7c23b 526 fCent = new TH1F("fCent","Centrality",100,0,100) ;
527 fOutputList->Add(fCent);
528
529 // Make common binning
f09766dd 530 const Double_t kMinP = 2.;
bfc7c23b 531 const Double_t kMaxP = 50.;
532 const Double_t kTPCSigMim = 40.;
533 const Double_t kTPCSigMax = 140.;
534
535 // 1st histogram: TPC dEdx with/without EMCAL (p, pT, TPC Signal, phi, eta, Sig, e/p, ,match, cell, M02, M20, Disp, Centrality, select)
fb87d707 536 Int_t nBins[16] = { 480, 200, 60, 20, 600, 300, 100, 40, 200, 200, 200, 100, 3, 3, 3, 10};
537 Double_t min[16] = {kMinP, kTPCSigMim, 1.0, -1.0, -8.0, 0, 0, 0, 0.0, 0.0, 0.0, 0, -1.5, -0.5, -0.5, -0.5};
538 Double_t max[16] = {kMaxP, kTPCSigMax, 4.0, 1.0, 4.0, 3.0, 0.1, 40, 2.0, 2.0, 2.0, 100, 1.5, 2.5, 2.5, 9.5};
539 fEleInfo = new THnSparseD("fEleInfo", "Electron Info; pT [GeV/c]; TPC signal;phi;eta;nSig; E/p;Rmatch;Ncell;M02;M20;Disp;Centrality;charge;opp;same;trigCond;", 16, nBins, min, max);
bfc7c23b 540 fOutputList->Add(fEleInfo);
541
fb87d707 542 //<--- Trigger info
543
544 fClsEBftTrigCut = new TH1F("fClsEBftTrigCut","cluster E before trigger selection",1000,0,100);
545 fOutputList->Add(fClsEBftTrigCut);
546
547 fClsEAftTrigCut = new TH1F("fClsEAftTrigCut","cluster E if cls has 0 trigcut channel",1000,0,100);
548 fOutputList->Add(fClsEAftTrigCut);
549
550 fClsEAftTrigCut1 = new TH1F("fClsEAftTrigCut1","cluster E if cls with trig channel",1000,0,100);
551 fOutputList->Add(fClsEAftTrigCut1);
552
553 fClsEAftTrigCut2 = new TH1F("fClsEAftTrigCut2","cluster E if cls with trigcut channel",1000,0,100);
554 fOutputList->Add(fClsEAftTrigCut2);
555
556 fClsEAftTrigCut3 = new TH1F("fClsEAftTrigCut3","cluster E if cls with trigcut channel + nCell>Ecorrect",1000,0,100);
557 fOutputList->Add(fClsEAftTrigCut3);
558
559 fClsEAftTrigCut4 = new TH1F("fClsEAftTrigCut4","cluster E if cls with trigcut channel + nCell>Ecorrect + cls time cut",1000,0,100);
560 fOutputList->Add(fClsEAftTrigCut4);
561
562 fClsETime = new TH2F("fClsETime", "Cls time vs E; E; time;",1000,0,100,1000,-0.0000002,0.0000009);
563 fOutputList->Add(fClsETime);
564
565 fClsETime1 = new TH2F("fClsETime1", "Cls time vs E if cls contains trigger channel; E; time;",1000,0,100,1000,-0.0000002,0.0000009);
566 fOutputList->Add(fClsETime1);
567
568 fTrigTimes = new TH1F("fTrigTimes", "Trigger time; time; N;",25,0,25);
569 fOutputList->Add(fTrigTimes);
570
bfc7c23b 571 PostData(1,fOutputList);
572}
573
574//________________________________________________________________________
575void AliAnalysisTaskHFECal::Terminate(Option_t *)
576{
577 // Info("Terminate");
578 AliAnalysisTaskSE::Terminate();
579}
580
581//________________________________________________________________________
582Bool_t AliAnalysisTaskHFECal::ProcessCutStep(Int_t cutStep, AliVParticle *track)
583{
584 // Check single track cuts for a given cut step
585 const Int_t kMCOffset = AliHFEcuts::kNcutStepsMCTrack;
586 if(!fCFM->CheckParticleCuts(cutStep + kMCOffset, track)) return kFALSE;
587 return kTRUE;
588}
589//_________________________________________
f09766dd 590//void AliAnalysisTaskHFECal::SelectPhotonicElectron(Int_t itrack, Double_t cent, AliESDtrack *track, Bool_t &fFlagPhotonicElec)
591void AliAnalysisTaskHFECal::SelectPhotonicElectron(Int_t itrack, Double_t cent, AliESDtrack *track, Bool_t &fFlagPhotonicElec, Bool_t &fFlagConvinatElec)
bfc7c23b 592{
593 //Identify non-heavy flavour electrons using Invariant mass method
594
fb87d707 595 fTrackCuts->SetAcceptKinkDaughters(kFALSE);
f09766dd 596 //fTrackCuts->SetRequireTPCRefit(kTRUE);
597 //fTrackCuts->SetEtaRange(-0.7,0.7);
598 //fTrackCuts->SetRequireSigmaToVertex(kTRUE);
bfc7c23b 599 fTrackCuts->SetMaxChi2PerClusterTPC(3.5);
f09766dd 600 fTrackCuts->SetMinNClustersTPC(70);
bfc7c23b 601
602 const AliESDVertex *pVtx = fESD->GetPrimaryVertex();
603
604 Bool_t flagPhotonicElec = kFALSE;
f09766dd 605 Bool_t flagConvinatElec = kFALSE;
bfc7c23b 606
5e65985c 607 //for(Int_t jTracks = itrack+1; jTracks<fESD->GetNumberOfTracks(); jTracks++){
608 for(Int_t jTracks = 0; jTracks<fESD->GetNumberOfTracks(); jTracks++){
bfc7c23b 609 AliESDtrack* trackAsso = fESD->GetTrack(jTracks);
610 if (!trackAsso) {
611 printf("ERROR: Could not receive track %d\n", jTracks);
612 continue;
613 }
5e65985c 614 if(itrack==jTracks)continue;
615
f09766dd 616 Double_t dEdxAsso = -999., ptPrim=-999., ptAsso=-999., openingAngle = -999.;
bfc7c23b 617 Double_t mass=999., width = -999;
618 Bool_t fFlagLS=kFALSE, fFlagULS=kFALSE;
619
f09766dd 620 ptPrim = track->Pt();
621
bfc7c23b 622 dEdxAsso = trackAsso->GetTPCsignal();
623 ptAsso = trackAsso->Pt();
624 Int_t chargeAsso = trackAsso->Charge();
625 Int_t charge = track->Charge();
626
5e65985c 627 //if(ptAsso <0.3) continue;
628 if(ptAsso <0.5) continue;
bfc7c23b 629 if(!fTrackCuts->AcceptTrack(trackAsso)) continue;
f09766dd 630 if(dEdxAsso <65 || dEdxAsso>100) continue; //11a pass1
bfc7c23b 631
632 Int_t fPDGe1 = 11; Int_t fPDGe2 = 11;
633 if(charge>0) fPDGe1 = -11;
634 if(chargeAsso>0) fPDGe2 = -11;
635
636 if(charge == chargeAsso) fFlagLS = kTRUE;
637 if(charge != chargeAsso) fFlagULS = kTRUE;
638
639 AliKFParticle ge1(*track, fPDGe1);
640 AliKFParticle ge2(*trackAsso, fPDGe2);
641 AliKFParticle recg(ge1, ge2);
642
643 if(recg.GetNDF()<1) continue;
644 Double_t chi2recg = recg.GetChi2()/recg.GetNDF();
645 if(TMath::Sqrt(TMath::Abs(chi2recg))>3.) continue;
646
647 AliKFVertex primV(*pVtx);
648 primV += recg;
649 recg.SetProductionVertex(primV);
650
651 recg.SetMassConstraint(0,0.0001);
652
653 openingAngle = ge1.GetAngle(ge2);
654 if(fFlagLS) fOpeningAngleLS->Fill(openingAngle);
655 if(fFlagULS) fOpeningAngleULS->Fill(openingAngle);
656
657 if(openingAngle > fOpeningAngleCut) continue;
658
659 recg.GetMass(mass,width);
660
f09766dd 661 double phoinfo[3];
662 phoinfo[0] = cent;
663 phoinfo[1] = ptPrim;
664 phoinfo[2] = mass;
665
666 if(fFlagLS) fInvmassLS->Fill(phoinfo);
667 if(fFlagULS) fInvmassULS->Fill(phoinfo);
bfc7c23b 668
669 if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec){
670 flagPhotonicElec = kTRUE;
671 }
f09766dd 672 if(mass<fInvmassCut && fFlagLS && !flagConvinatElec){
673 flagConvinatElec = kTRUE;
674 }
bfc7c23b 675
676 }
677 fFlagPhotonicElec = flagPhotonicElec;
f09766dd 678 fFlagConvinatElec = flagConvinatElec;
bfc7c23b 679
680}
681
fb87d707 682
683//_________________________________________
684void AliAnalysisTaskHFECal::FindTriggerClusters()
685{
686 // constants
687 const int nModuleCols = 2;
688 const int nModuleRows = 5;
689 const int nColsFeeModule = 48;
690 const int nRowsFeeModule = 24;
691 const int nColsFaltroModule = 24;
692 const int nRowsFaltroModule = 12;
693 //const int faltroWidthMax = 20;
694
695 // part 1, trigger extraction -------------------------------------
696 Int_t globCol, globRow;
697 //Int_t ntimes=0, nTrigChannel=0, nTrigChannelCut=0, trigInCut=0;
698 Int_t ntimes=0, nTrigChannel=0, nTrigChannelCut=0;
699
700 //Int_t trigtimes[faltroWidthMax];
701 Double_t cellTime[nColsFeeModule*nModuleCols][nRowsFeeModule*nModuleRows];
702 Double_t cellEnergy[nColsFeeModule*nModuleCols][nRowsFeeModule*nModuleRows];
703 //Double_t fTrigCutLow = 6;
704 //Double_t fTrigCutHigh = 10;
705 Double_t fTimeCutLow = 469e-09;
706 Double_t fTimeCutHigh = 715e-09;
707
708 AliESDCaloTrigger * fCaloTrigger = fESD->GetCaloTrigger( "EMCAL" );
709
710 // erase trigger maps
711 for(Int_t i = 0; i < nColsFaltroModule*nModuleCols; i++ )
712 {
713 for(Int_t j = 0; j < nRowsFaltroModule*nModuleRows; j++ )
714 {
715 ftriggersCut[i][j] = 0;
716 ftriggers[i][j] = 0;
717 ftriggersTime[i][j] = 0;
718 }
719 }
720
721 Int_t iglobCol=0, iglobRow=0;
722 // go through triggers
723 if( fCaloTrigger->GetEntries() > 0 )
724 {
725 // needs reset
726 fCaloTrigger->Reset();
727 while( fCaloTrigger->Next() )
728 {
729 fCaloTrigger->GetPosition( globCol, globRow );
730 fCaloTrigger->GetNL0Times( ntimes );
731 /*
732 // no L0s
733 if( ntimes < 1 ) continue;
734 // get precise timings
735 fCaloTrigger->GetL0Times( trigtimes );
736 trigInCut = 0;
737 for(Int_t i = 0; i < ntimes; i++ )
738 {
739 // save the first trigger time in channel
740 if( i == 0 || triggersTime[globCol][globRow] > trigtimes[i] )
741 triggersTime[globCol][globRow] = trigtimes[i];
742 //printf("trigger times: %d\n",trigtimes[i]);
743 // check if in cut
744 if(trigtimes[i] > fTrigCutLow && trigtimes[i] < fTrigCutHigh )
745 trigInCut = 1;
746
747 fTrigTimes->Fill(trigtimes[i]);
748 }
749 */
750
751 //L1 analysis from AliAnalysisTaskEMCALTriggerQA
752 Int_t bit = 0;
753 fCaloTrigger->GetTriggerBits(bit);
754
755 Int_t ts = 0;
756 fCaloTrigger->GetL1TimeSum(ts);
757 if (ts > 0)ftriggers[globCol][globRow] = 1;
758 // number of triggered channels in event
759 nTrigChannel++;
760 // ... inside cut
761 if(ts>0 && (bit >> 6 & 0x1))
762 {
763 iglobCol = globCol;
764 iglobRow = globRow;
765 nTrigChannelCut++;
766 ftriggersCut[globCol][globRow] = 1;
767 }
768
769 } // calo trigger entries
770 } // has calo trigger entries
771
772 // part 2 go through the clusters here -----------------------------------
773 Int_t nCluster=0, nCell=0, iCell=0, gCell=0;
774 Short_t cellAddr, nSACell, mclabel;
775 //Int_t nSACell, iSACell, mclabel;
776 Int_t iSACell;
777 Double_t cellAmp=0, cellTimeT=0, clusterTime=0, efrac;
778 Int_t nSupMod, nModule, nIphi, nIeta, iphi, ieta, gphi, geta, feta, fphi;
779
780 TRefArray *fCaloClusters = new TRefArray();
781 fESD->GetEMCALClusters( fCaloClusters );
782 nCluster = fCaloClusters->GetEntries();
783
784
785 // save all cells times if there are clusters
786 if( nCluster > 0 ){
787 // erase time map
788 for(Int_t i = 0; i < nColsFeeModule*nModuleCols; i++ ){
789 for(Int_t j = 0; j < nRowsFeeModule*nModuleRows; j++ ){
790 cellTime[i][j] = 0.;
791 cellEnergy[i][j] = 0.;
792 }
793 }
794
795 // get cells
796 AliESDCaloCells *fCaloCells = fESD->GetEMCALCells();
797 //AliVCaloCells fCaloCells = fESD->GetEMCALCells();
798 nSACell = fCaloCells->GetNumberOfCells();
799 for(iSACell = 0; iSACell < nSACell; iSACell++ ){
800 // get the cell info *fCal
801 fCaloCells->GetCell( iSACell, cellAddr, cellAmp, cellTimeT , mclabel, efrac);
802
803 // get cell position
804 fGeom->GetCellIndex( cellAddr, nSupMod, nModule, nIphi, nIeta );
805 fGeom->GetCellPhiEtaIndexInSModule( nSupMod,nModule, nIphi, nIeta, iphi, ieta);
806
807 // convert co global phi eta
808 gphi = iphi + nRowsFeeModule*(nSupMod/2);
809 geta = ieta + nColsFeeModule*(nSupMod%2);
810
811 // save cell time and energy
812 cellTime[geta][gphi] = cellTimeT;
813 cellEnergy[geta][gphi] = cellAmp;
814
815 }
816 }
817
818 Int_t nClusterTrig, nClusterTrigCut;
819 UShort_t *cellAddrs;
820 Double_t clsE=-999, clsEta=-999, clsPhi=-999;
821 Float_t clsPos[3] = {0.,0.,0.};
822
823 for(Int_t icl=0; icl<fESD->GetNumberOfCaloClusters(); icl++)
824 {
825 AliESDCaloCluster *cluster = fESD->GetCaloCluster(icl);
826 if(!cluster || !cluster->IsEMCAL()) continue;
827
828 // get cluster cells
829 nCell = cluster->GetNCells();
830
831 // get cluster energy
832 clsE = cluster->E();
833
834 // get cluster position
835 cluster->GetPosition(clsPos);
836 TVector3 clsPosVec(clsPos[0],clsPos[1],clsPos[2]);
837 clsEta = clsPosVec.Eta();
838 clsPhi = clsPosVec.Phi();
839
840 // get the cell addresses
841 cellAddrs = cluster->GetCellsAbsId();
842
843 // check if the cluster contains cell, that was marked as triggered
844 nClusterTrig = 0;
845 nClusterTrigCut = 0;
846
847 // loop the cells to check, if cluser in acceptance
848 // any cluster with a cell outside acceptance is not considered
849 for( iCell = 0; iCell < nCell; iCell++ )
850 {
851 // get cell position
852 fGeom->GetCellIndex( cellAddrs[iCell], nSupMod, nModule, nIphi, nIeta );
853 fGeom->GetCellPhiEtaIndexInSModule( nSupMod,nModule, nIphi, nIeta, iphi, ieta);
854
855 // convert co global phi eta
856 gphi = iphi + nRowsFeeModule*(nSupMod/2);
857 geta = ieta + nColsFeeModule*(nSupMod%2);
858
859 if( cellTime[geta][gphi] > 0. ){
860 clusterTime += cellTime[geta][gphi];
861 gCell++;
862 }
863
864 // get corresponding FALTRO
865 fphi = gphi / 2;
866 feta = geta / 2;
867
868 // try to match with a triggered
869 if( ftriggers[feta][fphi]==1)
870 { nClusterTrig++;
871 }
872 if( ftriggersCut[feta][fphi]==1)
873 { nClusterTrigCut++;
874 }
875
876 } // cells
877
878
879 if( gCell > 0 )
880 clusterTime = clusterTime / (Double_t)gCell;
881 // fix the reconstriction code time 100ns jumps
882 if( fESD->GetBunchCrossNumber() % 4 < 2 )
883 clusterTime -= 0.0000001;
884
885 fClsETime->Fill(clsE,clusterTime);
886 fClsEBftTrigCut->Fill(clsE);
887
888 if(nClusterTrig>0){
889 fClsETime1->Fill(clsE,clusterTime);
890 }
891
892 if(nClusterTrig>0){
893 cluster->SetChi2(1);
894 fClsEAftTrigCut1->Fill(clsE);
895 }
896
897 if(nClusterTrigCut>0){
898 cluster->SetChi2(2);
899 fClsEAftTrigCut2->Fill(clsE);
900 }
901
902 if(nClusterTrigCut>0 && ( nCell > (1 + clsE / 3)))
903 {
904 cluster->SetChi2(3);
905 fClsEAftTrigCut3->Fill(clsE);
906 }
907
908 if(nClusterTrigCut>0 && (nCell > (1 + clsE / 3) )&&( clusterTime > fTimeCutLow && clusterTime < fTimeCutHigh ))
909 {
910 // cluster->SetChi2(4);
911 fClsEAftTrigCut4->Fill(clsE);
912 }
913 if(nClusterTrigCut<1)
914 {
915 cluster->SetChi2(0);
916
917 fClsEAftTrigCut->Fill(clsE);
918 }
919
920 } // clusters
921}
922
923
924
925
926
927
928