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