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