cahnges from mconnors
[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);
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.;
390 if(fmcData && fMC && stack)
391 {
feffe705 392 Int_t label = TMath::Abs(track->GetLabel());
e7d87aef 393 TParticle* particle = stack->Particle(label);
394 int mcpid = particle->GetPdgCode();
9b3495ff 395 //printf("MCpid = %d",mcpid);
e7d87aef 396 if(particle->GetFirstMother()>-1)
397 {
398 int parentPID = stack->Particle(particle->GetFirstMother())->GetPdgCode();
399 if((parentPID==22 || parentPID==111 || parentPID==221)&& fabs(mcpid)==11)mcPho = kTRUE;
400 if((fabs(parentPID)==411 || fabs(parentPID)==413 || fabs(parentPID)==421 || fabs(parentPID)==423 || fabs(parentPID)==431)&& fabs(mcpid)==11)mcDtoE = kTRUE;
401 if((fabs(parentPID)==511 || fabs(parentPID)==513 || fabs(parentPID)==521 || fabs(parentPID)==523 || fabs(parentPID)==531)&& fabs(mcpid)==11)mcBtoE = kTRUE;
402 }
9b3495ff 403 if(fabs(mcpid)==11 && mcDtoE)mcele= 1.;
404 if(fabs(mcpid)==11 && mcBtoE)mcele= 2.;
405 if(fabs(mcpid)==11 && mcPho)mcele= 3.;
e7d87aef 406 }
407
9b3495ff 408 if(TMath::Abs(track->Eta())>0.6) continue;
f09766dd 409 if(TMath::Abs(track->Pt()<2.0)) continue;
bfc7c23b 410
411 fTrackPtBefTrkCuts->Fill(track->Pt());
412 // RecKine: ITSTPC cuts
413 if(!ProcessCutStep(AliHFEcuts::kStepRecKineITSTPC, track)) continue;
414
415 //RecKink
416 if(fRejectKinkMother) { // Quick and dirty fix to reject both kink mothers and daughters
417 if(track->GetKinkIndex(0) != 0) continue;
418 }
419
420 // RecPrim
421 if(!ProcessCutStep(AliHFEcuts::kStepRecPrim, track)) continue;
422
423 // HFEcuts: ITS layers cuts
424 if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsITS, track)) continue;
425
f09766dd 426 // HFE cuts: TPC PID cleanup
427 if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTPC, track)) continue;
428
bfc7c23b 429
430 fTrackPtAftTrkCuts->Fill(track->Pt());
431
432 Double_t mom = -999., eop=-999., pt = -999., dEdx=-999., fTPCnSigma=-10, phi=-999., eta=-999.;
feffe705 433 pt = track->Pt();
434 if(pt<2.0)continue;
bfc7c23b 435
436 // Track extrapolation
437
bfc7c23b 438 Int_t charge = track->Charge();
439 fTrkpt->Fill(pt);
440 mom = track->P();
441 phi = track->Phi();
442 eta = track->Eta();
443 dEdx = track->GetTPCsignal();
444 fTPCnSigma = fPID->GetPIDResponse() ? fPID->GetPIDResponse()->NumberOfSigmasTPC(track, AliPID::kElectron) : 1000;
445
446 double ncells = -1.0;
447 double m20 = -1.0;
448 double m02 = -1.0;
449 double disp = -1.0;
450 double rmatch = -1.0;
451 double nmatch = -1.0;
f09766dd 452 double oppstatus = 0.0;
bfc7c23b 453
f09766dd 454 Bool_t fFlagPhotonicElec = kFALSE;
455 Bool_t fFlagConvinatElec = kFALSE;
3db00c72 456 if(fTPCnSigma>-2.0 && fTPCnSigma<3.0)
457 {
458 SelectPhotonicElectron(iTracks,cent,track,fFlagPhotonicElec,fFlagConvinatElec,fTPCnSigma);
459 }
f09766dd 460 if(fFlagPhotonicElec)oppstatus = 1.0;
feffe705 461 if(fFlagConvinatElec)oppstatus = 2.0;
bfc7c23b 462
463 Int_t clsId = track->GetEMCALcluster();
464 if (clsId>0){
465 AliESDCaloCluster *clust = fESD->GetCaloCluster(clsId);
466 if(clust && clust->IsEMCAL()){
467
468 double clustE = clust->E();
469 eop = clustE/fabs(mom);
470 //double clustT = clust->GetTOF();
471 ncells = clust->GetNCells();
472 m02 = clust->GetM02();
473 m20 = clust->GetM20();
474 disp = clust->GetDispersion();
475 double delphi = clust->GetTrackDx();
476 double deleta = clust->GetTrackDz();
477 rmatch = sqrt(pow(delphi,2)+pow(deleta,2));
478 nmatch = clust->GetNTracksMatched();
479
feffe705 480 double valdedx[16];
fb87d707 481 valdedx[0] = pt; valdedx[1] = dEdx; valdedx[2] = phi; valdedx[3] = eta; valdedx[4] = fTPCnSigma;
482 valdedx[5] = eop; valdedx[6] = rmatch; valdedx[7] = ncells, valdedx[8] = m02; valdedx[9] = m20; valdedx[10] = disp;
feffe705 483 valdedx[11] = cent; valdedx[12] = charge; valdedx[13] = oppstatus; valdedx[14] = clust->Chi2();
484 valdedx[15] = mcele;
85985bb0 485 if(fqahist==1)fEleInfo->Fill(valdedx);
bfc7c23b 486
487
488 }
489 }
490
42c75692 491 fdEdxBef->Fill(mom,fTPCnSigma);
bfc7c23b 492 fTPCnsigma->Fill(mom,fTPCnSigma);
f09766dd 493 if(fTPCnSigma >= -1.0 && fTPCnSigma <= 3)fTrkEovPBef->Fill(pt,eop);
bfc7c23b 494
bfc7c23b 495 Int_t pidpassed = 1;
496
497
498 //--- track accepted
499 AliHFEpidObject hfetrack;
500 hfetrack.SetAnalysisType(AliHFEpidObject::kESDanalysis);
501 hfetrack.SetRecTrack(track);
3db00c72 502 double binct = 10.5;
503 if((0.0< cent) && (cent<5.0)) binct = 0.5;
504 if((5.0< cent) && (cent<10.0)) binct = 1.5;
505 if((10.0< cent) && (cent<20.0)) binct = 2.5;
506 if((20.0< cent) && (cent<30.0)) binct = 3.5;
507 if((30.0< cent) && (cent<40.0)) binct = 4.5;
508 if((40.0< cent) && (cent<50.0)) binct = 5.5;
509 if((50.0< cent) && (cent<60.0)) binct = 6.5;
510 if((60.0< cent) && (cent<70.0)) binct = 7.5;
511 if((70.0< cent) && (cent<80.0)) binct = 8.5;
512 if((80.0< cent) && (cent<90.0)) binct = 9.5;
513 if((90.0< cent) && (cent<100.0)) binct = 10.5;
514
515 hfetrack.SetCentrality((int)binct); //added
bfc7c23b 516 hfetrack.SetPbPb();
517 if(!fPID->IsSelected(&hfetrack, NULL, "", fPIDqa)) pidpassed = 0;
518
519 if(pidpassed==0) continue;
520
521 fTrkEovPAft->Fill(pt,eop);
42c75692 522 fdEdxAft->Fill(mom,fTPCnSigma);
bfc7c23b 523
e7d87aef 524 // Fill real data
fb87d707 525 fIncpT->Fill(cent,pt);
bfc7c23b 526 if(fFlagPhotonicElec) fPhoElecPt->Fill(cent,pt);
f09766dd 527 if(fFlagConvinatElec) fSameElecPt->Fill(cent,pt);
fb87d707 528
529 if(m20>0.0 && m20<0.3)
530 {
531 fIncpTM20->Fill(cent,pt);
532 if(fFlagPhotonicElec) fPhoElecPtM20->Fill(cent,pt);
533 if(fFlagConvinatElec) fSameElecPtM20->Fill(cent,pt);
534 }
feffe705 535
e7d87aef 536 // MC
2f37d8dc 537 if(mcele>0)
e7d87aef 538 {
3db00c72 539
540 fIncpTMChfeAll->Fill(cent,pt);
541 if(m20>0.0 && m20<0.3)fIncpTMCM20hfeAll->Fill(cent,pt);
542
e7d87aef 543 if(mcBtoE || mcDtoE)
544 {
545 fIncpTMChfe->Fill(cent,pt);
546 if(m20>0.0 && m20<0.3)fIncpTMCM20hfe->Fill(cent,pt);
547 }
548 if(mcPho)
549 {
550 fIncpTMCpho->Fill(cent,pt);
551 if(fFlagPhotonicElec) fPhoElecPtMC->Fill(cent,pt);
552 if(fFlagConvinatElec) fSameElecPtMC->Fill(cent,pt);
553
554 if(m20>0.0 && m20<0.3)
555 {
556 fIncpTMCM20pho->Fill(cent,pt);
557 if(fFlagPhotonicElec) fPhoElecPtMCM20->Fill(cent,pt);
558 if(fFlagConvinatElec) fSameElecPtMCM20->Fill(cent,pt);
559 }
560 }
561 }
bfc7c23b 562 }
563 PostData(1, fOutputList);
564}
565//_________________________________________
566void AliAnalysisTaskHFECal::UserCreateOutputObjects()
567{
42c75692 568 //--- Check MC
569
e7d87aef 570 //Bool_t mcData = kFALSE;
42c75692 571 if(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler())
572 {
e7d87aef 573 fmcData = kTRUE;
42c75692 574 printf("+++++ MC Data available");
575 }
e7d87aef 576 if(fmcData)
42c75692 577 {
578 printf("++++++++= MC analysis \n");
579 }
580 else
581 {
582 printf("++++++++ real data analysis \n");
583 }
584
85985bb0 585 printf("+++++++ QA hist %d",fqahist);
586
fb87d707 587 //---- Geometry
588 fGeom = AliEMCALGeometry::GetInstance("EMCAL_COMPLETEV1");
589
bfc7c23b 590 //--------Initialize PID
42c75692 591 //fPID->SetHasMCData(kFALSE);
e7d87aef 592 fPID->SetHasMCData(fmcData);
bfc7c23b 593 if(!fPID->GetNumberOfPIDdetectors())
594 {
595 fPID->AddDetector("TPC", 0);
596 fPID->AddDetector("EMCAL", 1);
597 }
598
3db00c72 599 Double_t params[4];
78ff1095 600 const char *cutmodel;
3db00c72 601 cutmodel = "pol0";
602 params[0] = -1.0; //sigma min
9b3495ff 603 double maxnSig = 3.0;
604 if(fmcData)
605 {
606 params[0] = -5.0; //sigma min
607 maxnSig = 5.0;
608 }
609
610 for(Int_t a=0;a<11;a++)fPID->ConfigureTPCcentralityCut(a,cutmodel,params,maxnSig);
3db00c72 611
bfc7c23b 612 fPID->SortDetectors();
613 fPIDqa = new AliHFEpidQAmanager();
614 fPIDqa->Initialize(fPID);
615
616 //--------Initialize correction Framework and Cuts
617 fCFM = new AliCFManager;
618 const Int_t kNcutSteps = AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kNcutStepsRecTrack + AliHFEcuts::kNcutStepsDETrack;
619 fCFM->SetNStepParticle(kNcutSteps);
620 for(Int_t istep = 0; istep < kNcutSteps; istep++)
621 fCFM->SetParticleCutsList(istep, NULL);
622
623 if(!fCuts){
624 AliWarning("Cuts not available. Default cuts will be used");
625 fCuts = new AliHFEcuts;
626 fCuts->CreateStandardCuts();
627 }
628 fCuts->Initialize(fCFM);
629
630 //---------Output Tlist
631 fOutputList = new TList();
632 fOutputList->SetOwner();
633 fOutputList->Add(fPIDqa->MakeList("PIDQA"));
634
635 fNoEvents = new TH1F("fNoEvents","",4,-0.5,3.5) ;
636 fOutputList->Add(fNoEvents);
637
fb87d707 638 Int_t binsE[5] = {250, 100, 1000, 100, 10};
639 Double_t xminE[5] = {1.0, -1, 0.0, 0, -0.5};
640 Double_t xmaxE[5] = {3.5, 1, 100.0, 100, 9.5};
641 fEMCAccE = new THnSparseD("fEMCAccE","EMC acceptance & E;#eta;#phi;Energy;Centrality;trugCondition;",5,binsE,xminE,xmaxE);
85985bb0 642 if(fqahist==1)fOutputList->Add(fEMCAccE);
bfc7c23b 643
644 fTrkpt = new TH1F("fTrkpt","track pt",100,0,50);
645 fOutputList->Add(fTrkpt);
646
647 fTrackPtBefTrkCuts = new TH1F("fTrackPtBefTrkCuts","track pt before track cuts",100,0,50);
648 fOutputList->Add(fTrackPtBefTrkCuts);
649
650 fTrackPtAftTrkCuts = new TH1F("fTrackPtAftTrkCuts","track pt after track cuts",100,0,50);
651 fOutputList->Add(fTrackPtAftTrkCuts);
652
653 fTPCnsigma = new TH2F("fTPCnsigma", "TPC - n sigma",100,0,50,200,-10,10);
654 fOutputList->Add(fTPCnsigma);
655
656 fTrkEovPBef = new TH2F("fTrkEovPBef","track E/p before HFE pid",100,0,50,100,0,2);
657 fOutputList->Add(fTrkEovPBef);
658
659 fTrkEovPAft = new TH2F("fTrkEovPAft","track E/p after HFE pid",100,0,50,100,0,2);
660 fOutputList->Add(fTrkEovPAft);
661
42c75692 662 fdEdxBef = new TH2F("fdEdxBef","track dEdx vs p before HFE pid",100,0,50,200,-10,10);
bfc7c23b 663 fOutputList->Add(fdEdxBef);
664
42c75692 665 fdEdxAft = new TH2F("fdEdxAft","track dEdx vs p after HFE pid",100,0,50,200,-10,10);
bfc7c23b 666 fOutputList->Add(fdEdxAft);
667
668 fIncpT = new TH2F("fIncpT","HFE pid electro vs. centrality",100,0,100,100,0,50);
669 fOutputList->Add(fIncpT);
670
fb87d707 671 fIncpTM20 = new TH2F("fIncpTM20","HFE pid electro vs. centrality with M20",100,0,100,100,0,50);
672 fOutputList->Add(fIncpTM20);
4e01b68c 673
674 Int_t nBinspho[5] = { 100, 100, 500, 12, 50};
675 Double_t minpho[5] = { 0., 0., 0., -2.5, 0};
676 Double_t maxpho[5] = {100., 50., 0.5, 3.5, 1};
f09766dd 677
4e01b68c 678 fInvmassLS = new THnSparseD("fInvmassLS", "Inv mass of LS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2); nSigma; angle;", 5, nBinspho,minpho, maxpho);
bfc7c23b 679 fOutputList->Add(fInvmassLS);
680
4e01b68c 681 fInvmassULS = new THnSparseD("fInvmassULS", "Inv mass of ULS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2); nSigma; angle; ", 5, nBinspho,minpho, maxpho);
bfc7c23b 682 fOutputList->Add(fInvmassULS);
683
684 fOpeningAngleLS = new TH1F("fOpeningAngleLS","Opening angle for LS pairs",100,0,1);
685 fOutputList->Add(fOpeningAngleLS);
686
687 fOpeningAngleULS = new TH1F("fOpeningAngleULS","Opening angle for ULS pairs",100,0,1);
688 fOutputList->Add(fOpeningAngleULS);
689
690 fPhotoElecPt = new TH1F("fPhotoElecPt", "photonic electron pt",100,0,50);
691 fOutputList->Add(fPhotoElecPt);
692
f09766dd 693 fPhoElecPt = new TH2F("fPhoElecPt", "Pho-inclusive electron pt",100,0,100,100,0,50);
bfc7c23b 694 fOutputList->Add(fPhoElecPt);
695
fb87d707 696 fPhoElecPtM20 = new TH2F("fPhoElecPtM20", "Pho-inclusive electron pt with M20",100,0,100,100,0,50);
697 fOutputList->Add(fPhoElecPtM20);
698
f09766dd 699 fSameElecPt = new TH2F("fSameElecPt", "Same-inclusive electron pt",100,0,100,100,0,50);
700 fOutputList->Add(fSameElecPt);
701
fb87d707 702 fSameElecPtM20 = new TH2F("fSameElecPtM20", "Same-inclusive electron pt with M20",100,0,100,100,0,50);
703 fOutputList->Add(fSameElecPtM20);
704
bfc7c23b 705 fCent = new TH1F("fCent","Centrality",100,0,100) ;
706 fOutputList->Add(fCent);
707
708 // Make common binning
f09766dd 709 const Double_t kMinP = 2.;
bfc7c23b 710 const Double_t kMaxP = 50.;
711 const Double_t kTPCSigMim = 40.;
712 const Double_t kTPCSigMax = 140.;
713
714 // 1st histogram: TPC dEdx with/without EMCAL (p, pT, TPC Signal, phi, eta, Sig, e/p, ,match, cell, M02, M20, Disp, Centrality, select)
9b3495ff 715 Int_t nBins[16] = { 480, 200, 60, 20, 600, 300, 100, 40, 200, 200, 200, 100, 3, 5, 10, 8};
feffe705 716 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, -1.5};
9b3495ff 717 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, 4.5, 9.5, 6.5};
feffe705 718 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;MCele", 16, nBins, min, max);
85985bb0 719 if(fqahist==1)fOutputList->Add(fEleInfo);
bfc7c23b 720
fb87d707 721 //<--- Trigger info
feffe705 722 /*
fb87d707 723 fClsEBftTrigCut = new TH1F("fClsEBftTrigCut","cluster E before trigger selection",1000,0,100);
724 fOutputList->Add(fClsEBftTrigCut);
725
726 fClsEAftTrigCut = new TH1F("fClsEAftTrigCut","cluster E if cls has 0 trigcut channel",1000,0,100);
727 fOutputList->Add(fClsEAftTrigCut);
728
729 fClsEAftTrigCut1 = new TH1F("fClsEAftTrigCut1","cluster E if cls with trig channel",1000,0,100);
730 fOutputList->Add(fClsEAftTrigCut1);
731
732 fClsEAftTrigCut2 = new TH1F("fClsEAftTrigCut2","cluster E if cls with trigcut channel",1000,0,100);
733 fOutputList->Add(fClsEAftTrigCut2);
734
735 fClsEAftTrigCut3 = new TH1F("fClsEAftTrigCut3","cluster E if cls with trigcut channel + nCell>Ecorrect",1000,0,100);
736 fOutputList->Add(fClsEAftTrigCut3);
737
738 fClsEAftTrigCut4 = new TH1F("fClsEAftTrigCut4","cluster E if cls with trigcut channel + nCell>Ecorrect + cls time cut",1000,0,100);
739 fOutputList->Add(fClsEAftTrigCut4);
740
741 fClsETime = new TH2F("fClsETime", "Cls time vs E; E; time;",1000,0,100,1000,-0.0000002,0.0000009);
742 fOutputList->Add(fClsETime);
743
744 fClsETime1 = new TH2F("fClsETime1", "Cls time vs E if cls contains trigger channel; E; time;",1000,0,100,1000,-0.0000002,0.0000009);
745 fOutputList->Add(fClsETime1);
746
747 fTrigTimes = new TH1F("fTrigTimes", "Trigger time; time; N;",25,0,25);
748 fOutputList->Add(fTrigTimes);
749
42c75692 750 fCellCheck = new TH2F("fCellCheck", "Cell vs E; E GeV; Cell ID",10,6,26,12000,0,12000);
751 fOutputList->Add(fCellCheck);
feffe705 752 */
e7d87aef 753 //<---------- MC
754
755 fInputHFEMC = new TH2F("fInputHFEMC","Input MC HFE pid electro vs. centrality",100,0,100,100,0,50);
756 fOutputList->Add(fInputHFEMC);
757
feffe705 758 fInputAlle = new TH2F("fInputAlle","Input MC electro vs. centrality",100,0,100,100,0,50);
759 fOutputList->Add(fInputAlle);
760
e7d87aef 761 fIncpTMChfe = new TH2F("fIncpTMChfe","MC HFE pid electro vs. centrality",100,0,100,100,0,50);
762 fOutputList->Add(fIncpTMChfe);
763
3db00c72 764 fIncpTMChfeAll = new TH2F("fIncpTMChfeAll","MC Alle pid electro vs. centrality",100,0,100,100,0,50);
765 fOutputList->Add(fIncpTMChfeAll);
766
e7d87aef 767 fIncpTMCM20hfe = new TH2F("fIncpTMCM20hfe","MC HFE pid electro vs. centrality with M20",100,0,100,100,0,50);
768 fOutputList->Add(fIncpTMCM20hfe);
769
3db00c72 770 fIncpTMCM20hfeAll = new TH2F("fIncpTMCM20hfeAll","MC Alle pid electro vs. centrality with M20",100,0,100,100,0,50);
771 fOutputList->Add(fIncpTMCM20hfeAll);
772
e7d87aef 773 fIncpTMCpho = new TH2F("fIncpTMCpho","MC Pho pid electro vs. centrality",100,0,100,100,0,50);
774 fOutputList->Add(fIncpTMCpho);
775
776 fIncpTMCM20pho = new TH2F("fIncpTMCM20pho","MC Pho pid electro vs. centrality with M20",100,0,100,100,0,50);
777 fOutputList->Add(fIncpTMCM20pho);
778
779 fPhoElecPtMC = new TH2F("fPhoElecPtMC", "MC Pho-inclusive electron pt",100,0,100,100,0,50);
780 fOutputList->Add(fPhoElecPtMC);
781
782 fPhoElecPtMCM20 = new TH2F("fPhoElecPtMCM20", "MC Pho-inclusive electron pt with M20",100,0,100,100,0,50);
feffe705 783 fOutputList->Add(fPhoElecPtMCM20);
e7d87aef 784
785 fSameElecPtMC = new TH2F("fSameElecPtMC", "MC Same-inclusive electron pt",100,0,100,100,0,50);
786 fOutputList->Add(fSameElecPtMC);
787
788 fSameElecPtMCM20 = new TH2F("fSameElecPtMCM20", "MC Same-inclusive electron pt with M20",100,0,100,100,0,50);
789 fOutputList->Add(fSameElecPtMCM20);
790
791
bfc7c23b 792 PostData(1,fOutputList);
793}
794
795//________________________________________________________________________
796void AliAnalysisTaskHFECal::Terminate(Option_t *)
797{
798 // Info("Terminate");
799 AliAnalysisTaskSE::Terminate();
800}
801
802//________________________________________________________________________
803Bool_t AliAnalysisTaskHFECal::ProcessCutStep(Int_t cutStep, AliVParticle *track)
804{
805 // Check single track cuts for a given cut step
806 const Int_t kMCOffset = AliHFEcuts::kNcutStepsMCTrack;
807 if(!fCFM->CheckParticleCuts(cutStep + kMCOffset, track)) return kFALSE;
808 return kTRUE;
809}
810//_________________________________________
f09766dd 811//void AliAnalysisTaskHFECal::SelectPhotonicElectron(Int_t itrack, Double_t cent, AliESDtrack *track, Bool_t &fFlagPhotonicElec)
3db00c72 812void AliAnalysisTaskHFECal::SelectPhotonicElectron(Int_t itrack, Double_t cent, AliESDtrack *track, Bool_t &fFlagPhotonicElec, Bool_t &fFlagConvinatElec, Double_t nSig)
bfc7c23b 813{
814 //Identify non-heavy flavour electrons using Invariant mass method
815
fb87d707 816 fTrackCuts->SetAcceptKinkDaughters(kFALSE);
feffe705 817 fTrackCuts->SetRequireTPCRefit(kTRUE);
85985bb0 818 fTrackCuts->SetRequireITSRefit(kTRUE);
feffe705 819 fTrackCuts->SetEtaRange(-0.9,0.9);
f09766dd 820 //fTrackCuts->SetRequireSigmaToVertex(kTRUE);
bfc7c23b 821 fTrackCuts->SetMaxChi2PerClusterTPC(3.5);
feffe705 822 fTrackCuts->SetMinNClustersTPC(90);
bfc7c23b 823
824 const AliESDVertex *pVtx = fESD->GetPrimaryVertex();
825
826 Bool_t flagPhotonicElec = kFALSE;
f09766dd 827 Bool_t flagConvinatElec = kFALSE;
bfc7c23b 828
5e65985c 829 //for(Int_t jTracks = itrack+1; jTracks<fESD->GetNumberOfTracks(); jTracks++){
830 for(Int_t jTracks = 0; jTracks<fESD->GetNumberOfTracks(); jTracks++){
bfc7c23b 831 AliESDtrack* trackAsso = fESD->GetTrack(jTracks);
832 if (!trackAsso) {
833 printf("ERROR: Could not receive track %d\n", jTracks);
834 continue;
835 }
5e65985c 836 if(itrack==jTracks)continue;
837
f09766dd 838 Double_t dEdxAsso = -999., ptPrim=-999., ptAsso=-999., openingAngle = -999.;
bfc7c23b 839 Double_t mass=999., width = -999;
840 Bool_t fFlagLS=kFALSE, fFlagULS=kFALSE;
841
f09766dd 842 ptPrim = track->Pt();
843
bfc7c23b 844 dEdxAsso = trackAsso->GetTPCsignal();
845 ptAsso = trackAsso->Pt();
846 Int_t chargeAsso = trackAsso->Charge();
847 Int_t charge = track->Charge();
848
5e65985c 849 //if(ptAsso <0.3) continue;
3db00c72 850 if(ptAsso <0.5) continue;
bfc7c23b 851 if(!fTrackCuts->AcceptTrack(trackAsso)) continue;
f09766dd 852 if(dEdxAsso <65 || dEdxAsso>100) continue; //11a pass1
bfc7c23b 853
854 Int_t fPDGe1 = 11; Int_t fPDGe2 = 11;
855 if(charge>0) fPDGe1 = -11;
856 if(chargeAsso>0) fPDGe2 = -11;
857
858 if(charge == chargeAsso) fFlagLS = kTRUE;
859 if(charge != chargeAsso) fFlagULS = kTRUE;
860
861 AliKFParticle ge1(*track, fPDGe1);
862 AliKFParticle ge2(*trackAsso, fPDGe2);
863 AliKFParticle recg(ge1, ge2);
864
865 if(recg.GetNDF()<1) continue;
866 Double_t chi2recg = recg.GetChi2()/recg.GetNDF();
867 if(TMath::Sqrt(TMath::Abs(chi2recg))>3.) continue;
868
869 AliKFVertex primV(*pVtx);
870 primV += recg;
871 recg.SetProductionVertex(primV);
872
873 recg.SetMassConstraint(0,0.0001);
874
875 openingAngle = ge1.GetAngle(ge2);
876 if(fFlagLS) fOpeningAngleLS->Fill(openingAngle);
877 if(fFlagULS) fOpeningAngleULS->Fill(openingAngle);
878
bfc7c23b 879
880 recg.GetMass(mass,width);
881
4e01b68c 882 double phoinfo[5];
f09766dd 883 phoinfo[0] = cent;
884 phoinfo[1] = ptPrim;
885 phoinfo[2] = mass;
3db00c72 886 phoinfo[3] = nSig;
4e01b68c 887 phoinfo[4] = openingAngle;
f09766dd 888
889 if(fFlagLS) fInvmassLS->Fill(phoinfo);
890 if(fFlagULS) fInvmassULS->Fill(phoinfo);
ffa14dda 891
892 if(openingAngle > fOpeningAngleCut) continue;
893
bfc7c23b 894 if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec){
895 flagPhotonicElec = kTRUE;
896 }
f09766dd 897 if(mass<fInvmassCut && fFlagLS && !flagConvinatElec){
898 flagConvinatElec = kTRUE;
899 }
bfc7c23b 900
901 }
902 fFlagPhotonicElec = flagPhotonicElec;
f09766dd 903 fFlagConvinatElec = flagConvinatElec;
bfc7c23b 904
905}
906
fb87d707 907
908//_________________________________________
909void AliAnalysisTaskHFECal::FindTriggerClusters()
910{
911 // constants
912 const int nModuleCols = 2;
913 const int nModuleRows = 5;
914 const int nColsFeeModule = 48;
915 const int nRowsFeeModule = 24;
916 const int nColsFaltroModule = 24;
917 const int nRowsFaltroModule = 12;
918 //const int faltroWidthMax = 20;
919
920 // part 1, trigger extraction -------------------------------------
921 Int_t globCol, globRow;
922 //Int_t ntimes=0, nTrigChannel=0, nTrigChannelCut=0, trigInCut=0;
923 Int_t ntimes=0, nTrigChannel=0, nTrigChannelCut=0;
924
925 //Int_t trigtimes[faltroWidthMax];
926 Double_t cellTime[nColsFeeModule*nModuleCols][nRowsFeeModule*nModuleRows];
927 Double_t cellEnergy[nColsFeeModule*nModuleCols][nRowsFeeModule*nModuleRows];
928 //Double_t fTrigCutLow = 6;
929 //Double_t fTrigCutHigh = 10;
930 Double_t fTimeCutLow = 469e-09;
931 Double_t fTimeCutHigh = 715e-09;
932
933 AliESDCaloTrigger * fCaloTrigger = fESD->GetCaloTrigger( "EMCAL" );
934
935 // erase trigger maps
936 for(Int_t i = 0; i < nColsFaltroModule*nModuleCols; i++ )
937 {
938 for(Int_t j = 0; j < nRowsFaltroModule*nModuleRows; j++ )
939 {
940 ftriggersCut[i][j] = 0;
941 ftriggers[i][j] = 0;
942 ftriggersTime[i][j] = 0;
943 }
944 }
945
946 Int_t iglobCol=0, iglobRow=0;
947 // go through triggers
948 if( fCaloTrigger->GetEntries() > 0 )
949 {
950 // needs reset
951 fCaloTrigger->Reset();
952 while( fCaloTrigger->Next() )
953 {
954 fCaloTrigger->GetPosition( globCol, globRow );
955 fCaloTrigger->GetNL0Times( ntimes );
956 /*
957 // no L0s
958 if( ntimes < 1 ) continue;
959 // get precise timings
960 fCaloTrigger->GetL0Times( trigtimes );
961 trigInCut = 0;
962 for(Int_t i = 0; i < ntimes; i++ )
963 {
964 // save the first trigger time in channel
965 if( i == 0 || triggersTime[globCol][globRow] > trigtimes[i] )
966 triggersTime[globCol][globRow] = trigtimes[i];
967 //printf("trigger times: %d\n",trigtimes[i]);
968 // check if in cut
969 if(trigtimes[i] > fTrigCutLow && trigtimes[i] < fTrigCutHigh )
970 trigInCut = 1;
971
972 fTrigTimes->Fill(trigtimes[i]);
973 }
974 */
975
976 //L1 analysis from AliAnalysisTaskEMCALTriggerQA
977 Int_t bit = 0;
978 fCaloTrigger->GetTriggerBits(bit);
979
980 Int_t ts = 0;
981 fCaloTrigger->GetL1TimeSum(ts);
982 if (ts > 0)ftriggers[globCol][globRow] = 1;
983 // number of triggered channels in event
984 nTrigChannel++;
985 // ... inside cut
986 if(ts>0 && (bit >> 6 & 0x1))
987 {
988 iglobCol = globCol;
989 iglobRow = globRow;
990 nTrigChannelCut++;
991 ftriggersCut[globCol][globRow] = 1;
992 }
993
994 } // calo trigger entries
995 } // has calo trigger entries
996
997 // part 2 go through the clusters here -----------------------------------
998 Int_t nCluster=0, nCell=0, iCell=0, gCell=0;
999 Short_t cellAddr, nSACell, mclabel;
1000 //Int_t nSACell, iSACell, mclabel;
1001 Int_t iSACell;
42c75692 1002 Double_t cellAmp=0, cellTimeT=0, clusterTime=0, efrac=0;
fb87d707 1003 Int_t nSupMod, nModule, nIphi, nIeta, iphi, ieta, gphi, geta, feta, fphi;
1004
1005 TRefArray *fCaloClusters = new TRefArray();
1006 fESD->GetEMCALClusters( fCaloClusters );
1007 nCluster = fCaloClusters->GetEntries();
1008
1009
1010 // save all cells times if there are clusters
1011 if( nCluster > 0 ){
1012 // erase time map
1013 for(Int_t i = 0; i < nColsFeeModule*nModuleCols; i++ ){
1014 for(Int_t j = 0; j < nRowsFeeModule*nModuleRows; j++ ){
1015 cellTime[i][j] = 0.;
1016 cellEnergy[i][j] = 0.;
1017 }
1018 }
1019
1020 // get cells
1021 AliESDCaloCells *fCaloCells = fESD->GetEMCALCells();
1022 //AliVCaloCells fCaloCells = fESD->GetEMCALCells();
1023 nSACell = fCaloCells->GetNumberOfCells();
1024 for(iSACell = 0; iSACell < nSACell; iSACell++ ){
1025 // get the cell info *fCal
1026 fCaloCells->GetCell( iSACell, cellAddr, cellAmp, cellTimeT , mclabel, efrac);
1027
1028 // get cell position
1029 fGeom->GetCellIndex( cellAddr, nSupMod, nModule, nIphi, nIeta );
1030 fGeom->GetCellPhiEtaIndexInSModule( nSupMod,nModule, nIphi, nIeta, iphi, ieta);
1031
1032 // convert co global phi eta
1033 gphi = iphi + nRowsFeeModule*(nSupMod/2);
1034 geta = ieta + nColsFeeModule*(nSupMod%2);
1035
1036 // save cell time and energy
1037 cellTime[geta][gphi] = cellTimeT;
1038 cellEnergy[geta][gphi] = cellAmp;
1039
1040 }
1041 }
1042
1043 Int_t nClusterTrig, nClusterTrigCut;
1044 UShort_t *cellAddrs;
1045 Double_t clsE=-999, clsEta=-999, clsPhi=-999;
1046 Float_t clsPos[3] = {0.,0.,0.};
1047
1048 for(Int_t icl=0; icl<fESD->GetNumberOfCaloClusters(); icl++)
1049 {
1050 AliESDCaloCluster *cluster = fESD->GetCaloCluster(icl);
1051 if(!cluster || !cluster->IsEMCAL()) continue;
1052
1053 // get cluster cells
1054 nCell = cluster->GetNCells();
1055
1056 // get cluster energy
1057 clsE = cluster->E();
1058
1059 // get cluster position
1060 cluster->GetPosition(clsPos);
1061 TVector3 clsPosVec(clsPos[0],clsPos[1],clsPos[2]);
1062 clsEta = clsPosVec.Eta();
1063 clsPhi = clsPosVec.Phi();
1064
1065 // get the cell addresses
1066 cellAddrs = cluster->GetCellsAbsId();
1067
1068 // check if the cluster contains cell, that was marked as triggered
1069 nClusterTrig = 0;
1070 nClusterTrigCut = 0;
1071
1072 // loop the cells to check, if cluser in acceptance
1073 // any cluster with a cell outside acceptance is not considered
1074 for( iCell = 0; iCell < nCell; iCell++ )
1075 {
42c75692 1076 // check hot cell
feffe705 1077 //if(clsE>6.0)fCellCheck->Fill(clsE,cellAddrs[iCell]);
42c75692 1078
fb87d707 1079 // get cell position
1080 fGeom->GetCellIndex( cellAddrs[iCell], nSupMod, nModule, nIphi, nIeta );
1081 fGeom->GetCellPhiEtaIndexInSModule( nSupMod,nModule, nIphi, nIeta, iphi, ieta);
1082
1083 // convert co global phi eta
1084 gphi = iphi + nRowsFeeModule*(nSupMod/2);
1085 geta = ieta + nColsFeeModule*(nSupMod%2);
1086
1087 if( cellTime[geta][gphi] > 0. ){
1088 clusterTime += cellTime[geta][gphi];
1089 gCell++;
1090 }
1091
1092 // get corresponding FALTRO
1093 fphi = gphi / 2;
1094 feta = geta / 2;
1095
1096 // try to match with a triggered
1097 if( ftriggers[feta][fphi]==1)
1098 { nClusterTrig++;
1099 }
1100 if( ftriggersCut[feta][fphi]==1)
1101 { nClusterTrigCut++;
1102 }
1103
1104 } // cells
1105
1106
1107 if( gCell > 0 )
1108 clusterTime = clusterTime / (Double_t)gCell;
1109 // fix the reconstriction code time 100ns jumps
1110 if( fESD->GetBunchCrossNumber() % 4 < 2 )
1111 clusterTime -= 0.0000001;
1112
feffe705 1113 //fClsETime->Fill(clsE,clusterTime);
1114 //fClsEBftTrigCut->Fill(clsE);
fb87d707 1115
1116 if(nClusterTrig>0){
feffe705 1117 //fClsETime1->Fill(clsE,clusterTime);
fb87d707 1118 }
1119
1120 if(nClusterTrig>0){
1121 cluster->SetChi2(1);
feffe705 1122 //fClsEAftTrigCut1->Fill(clsE);
fb87d707 1123 }
1124
1125 if(nClusterTrigCut>0){
1126 cluster->SetChi2(2);
feffe705 1127 //fClsEAftTrigCut2->Fill(clsE);
fb87d707 1128 }
1129
1130 if(nClusterTrigCut>0 && ( nCell > (1 + clsE / 3)))
1131 {
1132 cluster->SetChi2(3);
feffe705 1133 //fClsEAftTrigCut3->Fill(clsE);
fb87d707 1134 }
1135
1136 if(nClusterTrigCut>0 && (nCell > (1 + clsE / 3) )&&( clusterTime > fTimeCutLow && clusterTime < fTimeCutHigh ))
1137 {
1138 // cluster->SetChi2(4);
feffe705 1139 //fClsEAftTrigCut4->Fill(clsE);
fb87d707 1140 }
1141 if(nClusterTrigCut<1)
1142 {
1143 cluster->SetChi2(0);
1144
feffe705 1145 //fClsEAftTrigCut->Fill(clsE);
fb87d707 1146 }
1147
1148 } // clusters
1149}
1150
1151
1152
1153
1154
1155
1156