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