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