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