on Proof dont delete the output objects in the destructor
[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)
5cc5103f 91 ,fqahist(0)
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)
8806ce6c 115 ,fInvmassLSmc(0)
116 ,fInvmassULSmc(0)
93c189c5 117 ,fInvmassLSmc0(0)
118 ,fInvmassLSmc1(0)
119 ,fInvmassLSmc2(0)
120 ,fInvmassLSmc3(0)
121 ,fInvmassULSmc0(0)
122 ,fInvmassULSmc1(0)
123 ,fInvmassULSmc2(0)
124 ,fInvmassULSmc3(0)
bfc7c23b 125 ,fOpeningAngleLS(0)
126 ,fOpeningAngleULS(0)
127 ,fPhotoElecPt(0)
128 ,fPhoElecPt(0)
fb87d707 129 ,fPhoElecPtM20(0)
f09766dd 130 ,fSameElecPt(0)
fb87d707 131 ,fSameElecPtM20(0)
bfc7c23b 132 ,fTrackPtBefTrkCuts(0)
133 ,fTrackPtAftTrkCuts(0)
134 ,fTPCnsigma(0)
135 ,fCent(0)
136 ,fEleInfo(0)
feffe705 137 /*
fb87d707 138 ,fClsEBftTrigCut(0)
139 ,fClsEAftTrigCut(0)
140 ,fClsEAftTrigCut1(0)
141 ,fClsEAftTrigCut2(0)
142 ,fClsEAftTrigCut3(0)
143 ,fClsEAftTrigCut4(0)
144 ,fClsETime(0)
145 ,fClsETime1(0)
146 ,fTrigTimes(0)
42c75692 147 ,fCellCheck(0)
feffe705 148 */
e7d87aef 149 ,fInputHFEMC(0)
feffe705 150 ,fInputAlle(0)
e7d87aef 151 ,fIncpTMChfe(0)
3db00c72 152 ,fIncpTMChfeAll(0)
e7d87aef 153 ,fIncpTMCM20hfe(0)
3db00c72 154 ,fIncpTMCM20hfeAll(0)
60544aea 155 ,fIncpTMCM20hfeCheck(0)
e7d87aef 156 ,fIncpTMCpho(0)
157 ,fIncpTMCM20pho(0)
158 ,fPhoElecPtMC(0)
159 ,fPhoElecPtMCM20(0)
160 ,fSameElecPtMC(0)
161 ,fSameElecPtMCM20(0)
e4b0faf2 162 ,fIncpTMCM20pho_pi0e(0)
163 ,fPhoElecPtMCM20_pi0e(0)
164 ,fSameElecPtMCM20_pi0e(0)
93c189c5 165 ,fIncpTMCM20pho_eta(0)
166 ,fPhoElecPtMCM20_eta(0)
167 ,fSameElecPtMCM20_eta(0)
a6df418c 168 ,CheckNclust(0)
169 ,CheckNits(0)
d113d7cd 170 ,Hpi0pTcheck(0)
19325033 171 ,HETApTcheck(0)
e4b0faf2 172 ,HphopTcheck(0)
93c189c5 173 ,fpTCheck(0)
50919258 174 ,fMomDtoE(0)
bfc7c23b 175{
176 //Named constructor
177
178 fPID = new AliHFEpid("hfePid");
179 fTrackCuts = new AliESDtrackCuts();
180
181 // Define input and output slots here
182 // Input slot #0 works with a TChain
183 DefineInput(0, TChain::Class());
184 // Output slot #0 id reserved by the base class for AOD
185 // Output slot #1 writes into a TH1 container
186 // DefineOutput(1, TH1I::Class());
187 DefineOutput(1, TList::Class());
188 // DefineOutput(3, TTree::Class());
189}
190
191//________________________________________________________________________
192AliAnalysisTaskHFECal::AliAnalysisTaskHFECal()
193 : AliAnalysisTaskSE("DefaultAnalysis_AliAnalysisTaskHFECal")
194 ,fESD(0)
e7d87aef 195 ,fMC(0)
fb87d707 196 ,fGeom(0)
bfc7c23b 197 ,fOutputList(0)
5cc5103f 198 ,fqahist(0)
bfc7c23b 199 ,fTrackCuts(0)
200 ,fCuts(0)
201 ,fIdentifiedAsOutInz(kFALSE)
202 ,fPassTheEventCut(kFALSE)
203 ,fRejectKinkMother(kFALSE)
e7d87aef 204 ,fmcData(kFALSE)
bfc7c23b 205 ,fVz(0.0)
206 ,fCFM(0)
207 ,fPID(0)
208 ,fPIDqa(0)
209 ,fOpeningAngleCut(0.1)
210 ,fInvmassCut(0.01)
211 ,fNoEvents(0)
212 ,fEMCAccE(0)
213 ,fTrkpt(0)
214 ,fTrkEovPBef(0)
215 ,fTrkEovPAft(0)
216 ,fdEdxBef(0)
217 ,fdEdxAft(0)
218 ,fIncpT(0)
fb87d707 219 ,fIncpTM20(0)
bfc7c23b 220 ,fInvmassLS(0)
221 ,fInvmassULS(0)
8806ce6c 222 ,fInvmassLSmc(0)
223 ,fInvmassULSmc(0)
93c189c5 224 ,fInvmassLSmc0(0)
225 ,fInvmassLSmc1(0)
226 ,fInvmassLSmc2(0)
227 ,fInvmassLSmc3(0)
228 ,fInvmassULSmc0(0)
229 ,fInvmassULSmc1(0)
230 ,fInvmassULSmc2(0)
231 ,fInvmassULSmc3(0)
bfc7c23b 232 ,fOpeningAngleLS(0)
233 ,fOpeningAngleULS(0)
234 ,fPhotoElecPt(0)
235 ,fPhoElecPt(0)
fb87d707 236 ,fPhoElecPtM20(0)
f09766dd 237 ,fSameElecPt(0)
fb87d707 238 ,fSameElecPtM20(0)
bfc7c23b 239 ,fTrackPtBefTrkCuts(0)
240 ,fTrackPtAftTrkCuts(0)
241 ,fTPCnsigma(0)
242 ,fCent(0)
243 ,fEleInfo(0)
feffe705 244 /*
fb87d707 245 ,fClsEBftTrigCut(0)
246 ,fClsEAftTrigCut(0)
247 ,fClsEAftTrigCut1(0)
248 ,fClsEAftTrigCut2(0)
249 ,fClsEAftTrigCut3(0)
250 ,fClsEAftTrigCut4(0)
251 ,fClsETime(0)
252 ,fClsETime1(0)
253 ,fTrigTimes(0)
42c75692 254 ,fCellCheck(0)
feffe705 255 */
e7d87aef 256 ,fInputHFEMC(0)
feffe705 257 ,fInputAlle(0)
e7d87aef 258 ,fIncpTMChfe(0)
3db00c72 259 ,fIncpTMChfeAll(0)
e7d87aef 260 ,fIncpTMCM20hfe(0)
3db00c72 261 ,fIncpTMCM20hfeAll(0)
60544aea 262 ,fIncpTMCM20hfeCheck(0)
e7d87aef 263 ,fIncpTMCpho(0)
264 ,fIncpTMCM20pho(0)
265 ,fPhoElecPtMC(0)
266 ,fPhoElecPtMCM20(0)
267 ,fSameElecPtMC(0)
268 ,fSameElecPtMCM20(0)
e4b0faf2 269 ,fIncpTMCM20pho_pi0e(0)
270 ,fPhoElecPtMCM20_pi0e(0)
271 ,fSameElecPtMCM20_pi0e(0)
93c189c5 272 ,fIncpTMCM20pho_eta(0)
273 ,fPhoElecPtMCM20_eta(0)
274 ,fSameElecPtMCM20_eta(0)
a6df418c 275 ,CheckNclust(0)
276 ,CheckNits(0)
d113d7cd 277 ,Hpi0pTcheck(0)
19325033 278 ,HETApTcheck(0)
e4b0faf2 279 ,HphopTcheck(0)
93c189c5 280 ,fpTCheck(0)
50919258 281 ,fMomDtoE(0)
bfc7c23b 282{
283 //Default constructor
284 fPID = new AliHFEpid("hfePid");
285
286 fTrackCuts = new AliESDtrackCuts();
287
288 // Constructor
289 // Define input and output slots here
290 // Input slot #0 works with a TChain
291 DefineInput(0, TChain::Class());
292 // Output slot #0 id reserved by the base class for AOD
293 // Output slot #1 writes into a TH1 container
294 // DefineOutput(1, TH1I::Class());
295 DefineOutput(1, TList::Class());
296 //DefineOutput(3, TTree::Class());
297}
298//_________________________________________
299
300AliAnalysisTaskHFECal::~AliAnalysisTaskHFECal()
301{
302 //Destructor
303
304 delete fOutputList;
fb87d707 305 delete fGeom;
bfc7c23b 306 delete fPID;
a6df418c 307 delete fCuts;
bfc7c23b 308 delete fCFM;
309 delete fPIDqa;
310 delete fTrackCuts;
311}
312//_________________________________________
313
314void AliAnalysisTaskHFECal::UserExec(Option_t*)
315{
316 //Main loop
317 //Called for each event
318
319 // create pointer to event
320 fESD = dynamic_cast<AliESDEvent*>(InputEvent());
321 if (!fESD) {
322 printf("ERROR: fESD not available\n");
323 return;
324 }
325
326 if(!fCuts){
327 AliError("HFE cuts not available");
328 return;
329 }
330
331 if(!fPID->IsInitialized()){
332 // Initialize PID with the given run number
333 AliWarning("PID not initialised, get from Run no");
334 fPID->InitializePID(fESD->GetRunNumber());
335 }
336
e7d87aef 337 if(fmcData)fMC = MCEvent();
338 AliStack* stack = NULL;
339 if(fmcData && fMC)stack = fMC->Stack();
340
341 Float_t cent = -1.;
342 AliCentrality *centrality = fESD->GetCentrality();
343 cent = centrality->GetCentralityPercentile("V0M");
344
345 //---- fill MC track info
346 if(fmcData && fMC)
347 {
348 Int_t nParticles = stack->GetNtrack();
349 for (Int_t iParticle = 0; iParticle < nParticles; iParticle++) {
350 TParticle* particle = stack->Particle(iParticle);
351 int fPDG = particle->GetPdgCode();
352 double mcZvertex = fMC->GetPrimaryVertex()->GetZ();
353 double pTMC = particle->Pt();
feffe705 354 double proR = particle->R();
9b3495ff 355 double etaMC = particle->Eta();
356 if(fabs(etaMC)>0.6)continue;
e7d87aef 357 Bool_t mcInDtoE= kFALSE;
358 Bool_t mcInBtoE= kFALSE;
359
d113d7cd 360 Bool_t MChijing = fMC->IsFromBGEvent(iParticle);
fd426741 361 //if(!MChijing)printf("not MC hijing");
d113d7cd 362 int iHijing = 1;
363 if(!MChijing)iHijing = 0;
364 if(fPDG==111)Hpi0pTcheck->Fill(pTMC,iHijing);
19325033 365 if(fPDG==221)HETApTcheck->Fill(pTMC,iHijing);
d113d7cd 366
e4b0faf2 367 if(particle->GetFirstMother()>-1 && fPDG==22)
368 {
369 int parentPID = stack->Particle(particle->GetFirstMother())->GetPdgCode();
370 if(parentPID==111 || parentPID==221)HphopTcheck->Fill(pTMC,iHijing); // pi0->g & eta->g
371 }
372
e7d87aef 373 if(particle->GetFirstMother()>-1 && fabs(fPDG)==11)
374 {
375 int parentPID = stack->Particle(particle->GetFirstMother())->GetPdgCode();
376 if((fabs(parentPID)==411 || fabs(parentPID)==413 || fabs(parentPID)==421 || fabs(parentPID)==423 || fabs(parentPID)==431)&& fabs(fPDG)==11)mcInDtoE = kTRUE;
377 if((fabs(parentPID)==511 || fabs(parentPID)==513 || fabs(parentPID)==521 || fabs(parentPID)==523 || fabs(parentPID)==531)&& fabs(fPDG)==11)mcInBtoE = kTRUE;
378 if((mcInBtoE || mcInDtoE) && fabs(mcZvertex)<10.0)fInputHFEMC->Fill(cent,pTMC);
379 }
380
feffe705 381 if(proR<7.0 && fabs(fPDG)==11)fInputAlle->Fill(cent,pTMC);
382
e7d87aef 383 }
384 }
385
bfc7c23b 386 fNoEvents->Fill(0);
387
388 Int_t fNOtrks = fESD->GetNumberOfTracks();
389 const AliESDVertex *pVtx = fESD->GetPrimaryVertex();
390
391 Double_t pVtxZ = -999;
392 pVtxZ = pVtx->GetZ();
393
394 if(TMath::Abs(pVtxZ)>10) return;
395 fNoEvents->Fill(1);
396
397 if(fNOtrks<2) return;
398 fNoEvents->Fill(2);
399
400 AliPIDResponse *pidResponse = fInputHandler->GetPIDResponse();
401 if(!pidResponse){
402 AliDebug(1, "Using default PID Response");
403 pidResponse = AliHFEtools::GetDefaultPID(kFALSE, fInputEvent->IsA() == AliAODEvent::Class());
404 }
405
406 fPID->SetPIDResponse(pidResponse);
407
408 fCFM->SetRecEventInfo(fESD);
409
e7d87aef 410 //Float_t cent = -1.;
411 //AliCentrality *centrality = fESD->GetCentrality();
412 //cent = centrality->GetCentralityPercentile("V0M");
bfc7c23b 413 fCent->Fill(cent);
414
b12bc641 415 //if(cent>90.) return;
bfc7c23b 416
417 // Calorimeter info.
418
42c75692 419 FindTriggerClusters();
420
bfc7c23b 421 // make EMCAL array
422 for(Int_t iCluster=0; iCluster<fESD->GetNumberOfCaloClusters(); iCluster++)
423 {
424 AliESDCaloCluster *clust = fESD->GetCaloCluster(iCluster);
425 if(clust->IsEMCAL())
426 {
427 double clustE = clust->E();
428 float emcx[3]; // cluster pos
429 clust->GetPosition(emcx);
430 TVector3 clustpos(emcx[0],emcx[1],emcx[2]);
431 double emcphi = clustpos.Phi();
432 double emceta = clustpos.Eta();
fb87d707 433 double calInfo[5];
434 calInfo[0] = emcphi; calInfo[1] = emceta; calInfo[2] = clustE; calInfo[3] = cent; calInfo[4] = clust->Chi2();
85985bb0 435 //if(clustE>1.5)fEMCAccE->Fill(calInfo);
44be9e1d 436 //if(fqahist==1 && clustE>1.5)fEMCAccE->Fill(calInfo);
bfc7c23b 437 }
438 }
439
440 // Track loop
441 for (Int_t iTracks = 0; iTracks < fESD->GetNumberOfTracks(); iTracks++) {
442 AliESDtrack* track = fESD->GetTrack(iTracks);
443 if (!track) {
444 printf("ERROR: Could not receive track %d\n", iTracks);
445 continue;
446 }
e7d87aef 447
448 Bool_t mcPho = kFALSE;
449 Bool_t mcDtoE= kFALSE;
450 Bool_t mcBtoE= kFALSE;
93c189c5 451 Bool_t mcOrgPi0 = kFALSE;
452 Bool_t mcOrgEta = kFALSE;
e7d87aef 453 double mcele = -1.;
44be9e1d 454 double mcpT = 0.0;
b12bc641 455 double mcMompT = 0.0;
e4b0faf2 456 //double mcGrandMompT = 0.0;
457 double mcWeight = -10.0;
d113d7cd 458
459 int iHijing = 1;
460
e7d87aef 461 if(fmcData && fMC && stack)
462 {
feffe705 463 Int_t label = TMath::Abs(track->GetLabel());
d113d7cd 464
465 Bool_t MChijing = fMC->IsFromBGEvent(label);
466 if(!MChijing)iHijing = 0;
467
e7d87aef 468 TParticle* particle = stack->Particle(label);
469 int mcpid = particle->GetPdgCode();
44be9e1d 470 mcpT = particle->Pt();
9b3495ff 471 //printf("MCpid = %d",mcpid);
e7d87aef 472 if(particle->GetFirstMother()>-1)
473 {
e4b0faf2 474 int parentlabel = particle->GetFirstMother();
e7d87aef 475 int parentPID = stack->Particle(particle->GetFirstMother())->GetPdgCode();
b12bc641 476 mcMompT = stack->Particle(particle->GetFirstMother())->Pt();
e7d87aef 477 if((parentPID==22 || parentPID==111 || parentPID==221)&& fabs(mcpid)==11)mcPho = kTRUE;
478 if((fabs(parentPID)==411 || fabs(parentPID)==413 || fabs(parentPID)==421 || fabs(parentPID)==423 || fabs(parentPID)==431)&& fabs(mcpid)==11)mcDtoE = kTRUE;
479 if((fabs(parentPID)==511 || fabs(parentPID)==513 || fabs(parentPID)==521 || fabs(parentPID)==523 || fabs(parentPID)==531)&& fabs(mcpid)==11)mcBtoE = kTRUE;
e4b0faf2 480
50919258 481 // make D->e pT correlation
482 if(mcDtoE)fMomDtoE->Fill(mcpT,mcMompT);
483
484 // pi->e (Dalitz)
485 if(parentPID==111 && fabs(mcpid)==11 && mcMompT>0.0)
e4b0faf2 486 {
93c189c5 487 mcOrgPi0 = kTRUE;
8806ce6c 488 mcWeight = GetMCweight(mcMompT);
e4b0faf2 489 }
93c189c5 490 // eta->e (Dalitz)
491 if(parentPID==221 && fabs(mcpid)==11 && mcMompT>0.0)
492 {
493 mcOrgEta = kTRUE;
494 mcWeight = GetMCweightEta(mcMompT);
495 }
496
e4b0faf2 497
498 // access grand parent pi0->g->e
499 TParticle* particle_parent = stack->Particle(parentlabel); // get parent pointer
500 if(particle_parent->GetFirstMother()>-1 && parentPID==22 && fabs(mcpid)==11) // get grand parent g->e
501 {
502 int grand_parentPID = stack->Particle(particle_parent->GetFirstMother())->GetPdgCode();
93c189c5 503 double pTtmp = stack->Particle(particle_parent->GetFirstMother())->Pt();
504 if(grand_parentPID==111 && pTtmp>0.0)
e4b0faf2 505 {
93c189c5 506 mcOrgPi0 = kTRUE;
11704995 507 mcWeight = GetMCweight(pTtmp);
93c189c5 508 }
509 if(grand_parentPID==221 && pTtmp>0.0)
510 {
511 mcOrgEta = kTRUE;
512 mcWeight = GetMCweightEta(pTtmp);
e4b0faf2 513 }
514 }
e7d87aef 515 }
9b3495ff 516 if(fabs(mcpid)==11 && mcDtoE)mcele= 1.;
517 if(fabs(mcpid)==11 && mcBtoE)mcele= 2.;
518 if(fabs(mcpid)==11 && mcPho)mcele= 3.;
e7d87aef 519 }
520
50919258 521 //printf("weight = %f\n",mcWeight);
522
9b3495ff 523 if(TMath::Abs(track->Eta())>0.6) continue;
f09766dd 524 if(TMath::Abs(track->Pt()<2.0)) continue;
bfc7c23b 525
526 fTrackPtBefTrkCuts->Fill(track->Pt());
527 // RecKine: ITSTPC cuts
528 if(!ProcessCutStep(AliHFEcuts::kStepRecKineITSTPC, track)) continue;
529
530 //RecKink
531 if(fRejectKinkMother) { // Quick and dirty fix to reject both kink mothers and daughters
532 if(track->GetKinkIndex(0) != 0) continue;
533 }
534
535 // RecPrim
536 if(!ProcessCutStep(AliHFEcuts::kStepRecPrim, track)) continue;
537
538 // HFEcuts: ITS layers cuts
539 if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsITS, track)) continue;
540
f09766dd 541 // HFE cuts: TPC PID cleanup
542 if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTPC, track)) continue;
543
a6df418c 544 int nTPCcl = track->GetTPCNcls();
6e92c428 545 int nTPCclF = track->GetTPCNclsF();
a6df418c 546 int nITS = track->GetNcls(0);
bfc7c23b 547
548 fTrackPtAftTrkCuts->Fill(track->Pt());
549
550 Double_t mom = -999., eop=-999., pt = -999., dEdx=-999., fTPCnSigma=-10, phi=-999., eta=-999.;
feffe705 551 pt = track->Pt();
552 if(pt<2.0)continue;
bfc7c23b 553
554 // Track extrapolation
555
bfc7c23b 556 Int_t charge = track->Charge();
557 fTrkpt->Fill(pt);
558 mom = track->P();
559 phi = track->Phi();
560 eta = track->Eta();
561 dEdx = track->GetTPCsignal();
562 fTPCnSigma = fPID->GetPIDResponse() ? fPID->GetPIDResponse()->NumberOfSigmasTPC(track, AliPID::kElectron) : 1000;
563
564 double ncells = -1.0;
565 double m20 = -1.0;
566 double m02 = -1.0;
567 double disp = -1.0;
568 double rmatch = -1.0;
569 double nmatch = -1.0;
f09766dd 570 double oppstatus = 0.0;
bfc7c23b 571
f09766dd 572 Bool_t fFlagPhotonicElec = kFALSE;
573 Bool_t fFlagConvinatElec = kFALSE;
bfc7c23b 574
575 Int_t clsId = track->GetEMCALcluster();
576 if (clsId>0){
577 AliESDCaloCluster *clust = fESD->GetCaloCluster(clsId);
578 if(clust && clust->IsEMCAL()){
579
580 double clustE = clust->E();
581 eop = clustE/fabs(mom);
582 //double clustT = clust->GetTOF();
583 ncells = clust->GetNCells();
584 m02 = clust->GetM02();
585 m20 = clust->GetM20();
586 disp = clust->GetDispersion();
587 double delphi = clust->GetTrackDx();
588 double deleta = clust->GetTrackDz();
589 rmatch = sqrt(pow(delphi,2)+pow(deleta,2));
590 nmatch = clust->GetNTracksMatched();
591
c2801a73 592 if(fTPCnSigma>-1.5 && fTPCnSigma<3.0)
593 {
8806ce6c 594 SelectPhotonicElectron(iTracks,cent,track,stack,fFlagPhotonicElec,fFlagConvinatElec,fTPCnSigma,m20,eop,mcele,mcWeight,iHijing);
c2801a73 595 }
596 if(fFlagPhotonicElec)oppstatus = 1.0;
597 if(fFlagConvinatElec)oppstatus = 2.0;
598 if(fFlagPhotonicElec && fFlagConvinatElec)oppstatus = 3.0;
599
feffe705 600 double valdedx[16];
a6df418c 601 valdedx[0] = pt; valdedx[1] = nITS; valdedx[2] = phi; valdedx[3] = eta; valdedx[4] = fTPCnSigma;
6e92c428 602 valdedx[5] = eop; valdedx[6] = rmatch; valdedx[7] = ncells, valdedx[8] = nTPCclF; valdedx[9] = m20; valdedx[10] = mcpT;
a6df418c 603 valdedx[11] = cent; valdedx[12] = charge; valdedx[13] = oppstatus; valdedx[14] = nTPCcl;
feffe705 604 valdedx[15] = mcele;
85985bb0 605 if(fqahist==1)fEleInfo->Fill(valdedx);
bfc7c23b 606
607
608 }
609 }
a6df418c 610
611 if(nITS<2.5)continue;
612 if(nTPCcl<100)continue;
613
614 CheckNclust->Fill(nTPCcl);
615 CheckNits->Fill(nITS);
616
42c75692 617 fdEdxBef->Fill(mom,fTPCnSigma);
bfc7c23b 618 fTPCnsigma->Fill(mom,fTPCnSigma);
f09766dd 619 if(fTPCnSigma >= -1.0 && fTPCnSigma <= 3)fTrkEovPBef->Fill(pt,eop);
bfc7c23b 620
bfc7c23b 621 Int_t pidpassed = 1;
622
623
624 //--- track accepted
625 AliHFEpidObject hfetrack;
626 hfetrack.SetAnalysisType(AliHFEpidObject::kESDanalysis);
627 hfetrack.SetRecTrack(track);
3db00c72 628 double binct = 10.5;
629 if((0.0< cent) && (cent<5.0)) binct = 0.5;
630 if((5.0< cent) && (cent<10.0)) binct = 1.5;
631 if((10.0< cent) && (cent<20.0)) binct = 2.5;
632 if((20.0< cent) && (cent<30.0)) binct = 3.5;
633 if((30.0< cent) && (cent<40.0)) binct = 4.5;
634 if((40.0< cent) && (cent<50.0)) binct = 5.5;
635 if((50.0< cent) && (cent<60.0)) binct = 6.5;
636 if((60.0< cent) && (cent<70.0)) binct = 7.5;
637 if((70.0< cent) && (cent<80.0)) binct = 8.5;
638 if((80.0< cent) && (cent<90.0)) binct = 9.5;
639 if((90.0< cent) && (cent<100.0)) binct = 10.5;
640
641 hfetrack.SetCentrality((int)binct); //added
bfc7c23b 642 hfetrack.SetPbPb();
643 if(!fPID->IsSelected(&hfetrack, NULL, "", fPIDqa)) pidpassed = 0;
644
645 if(pidpassed==0) continue;
646
647 fTrkEovPAft->Fill(pt,eop);
42c75692 648 fdEdxAft->Fill(mom,fTPCnSigma);
bfc7c23b 649
e7d87aef 650 // Fill real data
fb87d707 651 fIncpT->Fill(cent,pt);
bfc7c23b 652 if(fFlagPhotonicElec) fPhoElecPt->Fill(cent,pt);
f09766dd 653 if(fFlagConvinatElec) fSameElecPt->Fill(cent,pt);
fb87d707 654
655 if(m20>0.0 && m20<0.3)
656 {
657 fIncpTM20->Fill(cent,pt);
658 if(fFlagPhotonicElec) fPhoElecPtM20->Fill(cent,pt);
659 if(fFlagConvinatElec) fSameElecPtM20->Fill(cent,pt);
660 }
feffe705 661
e7d87aef 662 // MC
e4b0faf2 663 if(mcele>0) // select MC electrons
e7d87aef 664 {
3db00c72 665
666 fIncpTMChfeAll->Fill(cent,pt);
667 if(m20>0.0 && m20<0.3)fIncpTMCM20hfeAll->Fill(cent,pt);
668
e4b0faf2 669 if(mcBtoE || mcDtoE) // select B->e & D->e
e7d87aef 670 {
671 fIncpTMChfe->Fill(cent,pt);
672 if(m20>0.0 && m20<0.3)fIncpTMCM20hfe->Fill(cent,pt);
60544aea 673 if(m20>0.0 && m20<0.3)fIncpTMCM20hfeCheck->Fill(cent,mcpT);
e7d87aef 674 }
e4b0faf2 675
676 if(mcPho) // select photonic electrons
e7d87aef 677 {
987053ce 678 double phoval[5];
44be9e1d 679 phoval[0] = cent;
680 phoval[1] = pt;
681 phoval[2] = fTPCnSigma;
d113d7cd 682 phoval[3] = iHijing;
987053ce 683 phoval[4] = mcMompT;
44be9e1d 684
685 fIncpTMCpho->Fill(phoval);
686 if(fFlagPhotonicElec) fPhoElecPtMC->Fill(phoval);
687 if(fFlagConvinatElec) fSameElecPtMC->Fill(phoval);
e7d87aef 688
689 if(m20>0.0 && m20<0.3)
690 {
44be9e1d 691 fIncpTMCM20pho->Fill(phoval);
692 if(fFlagPhotonicElec) fPhoElecPtMCM20->Fill(phoval);
693 if(fFlagConvinatElec) fSameElecPtMCM20->Fill(phoval);
e4b0faf2 694 // pi0->g->e
695 if(mcWeight>-1)
696 {
987053ce 697 if(iHijing==1)mcWeight = 1.0;
93c189c5 698 if(mcOrgPi0)
699 {
700 fIncpTMCM20pho_pi0e->Fill(phoval,mcWeight);
701 if(fFlagPhotonicElec) fPhoElecPtMCM20_pi0e->Fill(phoval,mcWeight);
702 if(fFlagConvinatElec) fSameElecPtMCM20_pi0e->Fill(phoval,mcWeight);
703 }
704 if(mcOrgEta)
705 {
706 fIncpTMCM20pho_eta->Fill(phoval,mcWeight);
707 if(fFlagPhotonicElec) fPhoElecPtMCM20_eta->Fill(phoval,mcWeight);
708 if(fFlagConvinatElec) fSameElecPtMCM20_eta->Fill(phoval,mcWeight);
709 }
710 // --- eta
e4b0faf2 711 }
e7d87aef 712 }
713 }
714 }
bfc7c23b 715 }
716 PostData(1, fOutputList);
717}
718//_________________________________________
719void AliAnalysisTaskHFECal::UserCreateOutputObjects()
720{
42c75692 721 //--- Check MC
722
e7d87aef 723 //Bool_t mcData = kFALSE;
42c75692 724 if(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler())
725 {
e7d87aef 726 fmcData = kTRUE;
42c75692 727 printf("+++++ MC Data available");
728 }
e7d87aef 729 if(fmcData)
42c75692 730 {
731 printf("++++++++= MC analysis \n");
732 }
733 else
734 {
735 printf("++++++++ real data analysis \n");
736 }
737
85985bb0 738 printf("+++++++ QA hist %d",fqahist);
739
fb87d707 740 //---- Geometry
741 fGeom = AliEMCALGeometry::GetInstance("EMCAL_COMPLETEV1");
742
bfc7c23b 743 //--------Initialize PID
42c75692 744 //fPID->SetHasMCData(kFALSE);
e7d87aef 745 fPID->SetHasMCData(fmcData);
bfc7c23b 746 if(!fPID->GetNumberOfPIDdetectors())
747 {
748 fPID->AddDetector("TPC", 0);
749 fPID->AddDetector("EMCAL", 1);
750 }
751
3db00c72 752 Double_t params[4];
78ff1095 753 const char *cutmodel;
3db00c72 754 cutmodel = "pol0";
755 params[0] = -1.0; //sigma min
9b3495ff 756 double maxnSig = 3.0;
757 if(fmcData)
758 {
759 params[0] = -5.0; //sigma min
760 maxnSig = 5.0;
761 }
762
763 for(Int_t a=0;a<11;a++)fPID->ConfigureTPCcentralityCut(a,cutmodel,params,maxnSig);
3db00c72 764
bfc7c23b 765 fPID->SortDetectors();
766 fPIDqa = new AliHFEpidQAmanager();
767 fPIDqa->Initialize(fPID);
a6df418c 768
769 //------- fcut --------------
770 fCuts = new AliHFEcuts();
771 fCuts->CreateStandardCuts();
772 //fCuts->SetMinNClustersTPC(100);
773 fCuts->SetMinNClustersTPC(90);
774 fCuts->SetMinRatioTPCclusters(0.6);
775 fCuts->SetTPCmodes(AliHFEextraCuts::kFound, AliHFEextraCuts::kFoundOverFindable);
776 //fCuts->SetMinNClustersITS(3);
777 fCuts->SetMinNClustersITS(2);
778 fCuts->SetCutITSpixel(AliHFEextraCuts::kAny);
779 fCuts->SetCheckITSLayerStatus(kFALSE);
780 fCuts->SetVertexRange(10.);
781 fCuts->SetTOFPIDStep(kFALSE);
782 fCuts->SetPtRange(2, 50);
783 fCuts->SetMaxImpactParam(3.,3.);
784
bfc7c23b 785 //--------Initialize correction Framework and Cuts
786 fCFM = new AliCFManager;
787 const Int_t kNcutSteps = AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kNcutStepsRecTrack + AliHFEcuts::kNcutStepsDETrack;
788 fCFM->SetNStepParticle(kNcutSteps);
789 for(Int_t istep = 0; istep < kNcutSteps; istep++)
790 fCFM->SetParticleCutsList(istep, NULL);
791
792 if(!fCuts){
793 AliWarning("Cuts not available. Default cuts will be used");
794 fCuts = new AliHFEcuts;
795 fCuts->CreateStandardCuts();
796 }
797 fCuts->Initialize(fCFM);
798
799 //---------Output Tlist
800 fOutputList = new TList();
801 fOutputList->SetOwner();
802 fOutputList->Add(fPIDqa->MakeList("PIDQA"));
803
804 fNoEvents = new TH1F("fNoEvents","",4,-0.5,3.5) ;
805 fOutputList->Add(fNoEvents);
806
02cd82a0 807 Int_t binsE[5] = {250, 100, 1000, 200, 10};
fb87d707 808 Double_t xminE[5] = {1.0, -1, 0.0, 0, -0.5};
809 Double_t xmaxE[5] = {3.5, 1, 100.0, 100, 9.5};
810 fEMCAccE = new THnSparseD("fEMCAccE","EMC acceptance & E;#eta;#phi;Energy;Centrality;trugCondition;",5,binsE,xminE,xmaxE);
85985bb0 811 if(fqahist==1)fOutputList->Add(fEMCAccE);
bfc7c23b 812
813 fTrkpt = new TH1F("fTrkpt","track pt",100,0,50);
814 fOutputList->Add(fTrkpt);
815
816 fTrackPtBefTrkCuts = new TH1F("fTrackPtBefTrkCuts","track pt before track cuts",100,0,50);
817 fOutputList->Add(fTrackPtBefTrkCuts);
818
819 fTrackPtAftTrkCuts = new TH1F("fTrackPtAftTrkCuts","track pt after track cuts",100,0,50);
820 fOutputList->Add(fTrackPtAftTrkCuts);
821
822 fTPCnsigma = new TH2F("fTPCnsigma", "TPC - n sigma",100,0,50,200,-10,10);
823 fOutputList->Add(fTPCnsigma);
824
825 fTrkEovPBef = new TH2F("fTrkEovPBef","track E/p before HFE pid",100,0,50,100,0,2);
826 fOutputList->Add(fTrkEovPBef);
827
828 fTrkEovPAft = new TH2F("fTrkEovPAft","track E/p after HFE pid",100,0,50,100,0,2);
829 fOutputList->Add(fTrkEovPAft);
830
42c75692 831 fdEdxBef = new TH2F("fdEdxBef","track dEdx vs p before HFE pid",100,0,50,200,-10,10);
bfc7c23b 832 fOutputList->Add(fdEdxBef);
833
42c75692 834 fdEdxAft = new TH2F("fdEdxAft","track dEdx vs p after HFE pid",100,0,50,200,-10,10);
bfc7c23b 835 fOutputList->Add(fdEdxAft);
836
02cd82a0 837 fIncpT = new TH2F("fIncpT","HFE pid electro vs. centrality",200,0,100,100,0,50);
bfc7c23b 838 fOutputList->Add(fIncpT);
839
02cd82a0 840 fIncpTM20 = new TH2F("fIncpTM20","HFE pid electro vs. centrality with M20",200,0,100,100,0,50);
fb87d707 841 fOutputList->Add(fIncpTM20);
4e01b68c 842
02cd82a0 843 Int_t nBinspho[9] = { 200, 100, 500, 12, 50, 4, 200, 8, 100};
b12bc641 844 Double_t minpho[9] = { 0., 0., 0., -2.5, 0, -0.5, 0,-1.5, 0};
845 Double_t maxpho[9] = {100., 50., 0.5, 3.5, 1, 3.5, 2, 6.5, 50};
f09766dd 846
b12bc641 847 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);
5cc5103f 848 if(fqahist==1)fOutputList->Add(fInvmassLS);
bfc7c23b 849
b12bc641 850 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);
5cc5103f 851 if(fqahist==1)fOutputList->Add(fInvmassULS);
bfc7c23b 852
8806ce6c 853 fInvmassLSmc = new THnSparseD("fInvmassLSmc", "Inv mass of LS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2); nSigma; angle; m20cut; eop; Mcele;", 9, nBinspho,minpho, maxpho);
5cc5103f 854 if(fqahist==1)fOutputList->Add(fInvmassLSmc);
8806ce6c 855
856 fInvmassULSmc = new THnSparseD("fInvmassULSmc", "Inv mass of ULS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2); nSigma; angle; m20cut; eop; MCele", 9, nBinspho,minpho, maxpho);
5cc5103f 857 if(fqahist==1)fOutputList->Add(fInvmassULSmc);
8806ce6c 858
93c189c5 859 fInvmassLSmc0 = new TH2D("fInvmassLSmc0", "Inv mass of LS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2)",20,0,20,500,0,0.5 );
860 fInvmassLSmc0->Sumw2();
861 fOutputList->Add(fInvmassLSmc0);
862
863 fInvmassLSmc1 = new TH2D("fInvmassLSmc1", "Inv mass of LS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2)",20,0,20,500,0,0.5 );
864 fInvmassLSmc1->Sumw2();
865 fOutputList->Add(fInvmassLSmc1);
866
867 fInvmassLSmc2 = new TH2D("fInvmassLSmc2", "Inv mass of LS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2)",20,0,20,500,0,0.5 );
868 fInvmassLSmc2->Sumw2();
869 fOutputList->Add(fInvmassLSmc2);
870
871 fInvmassLSmc3 = new TH2D("fInvmassLSmc3", "Inv mass of LS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2)",20,0,20,500,0,0.5 );
872 fInvmassLSmc3->Sumw2();
873 fOutputList->Add(fInvmassLSmc3);
874
875 fInvmassULSmc0 = new TH2D("fInvmassULSmc0", "Inv mass of ULS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2)",20,0,20,500,0,0.5 );
876 fInvmassULSmc0->Sumw2();
877 fOutputList->Add(fInvmassULSmc0);
878
879 fInvmassULSmc1 = new TH2D("fInvmassULSmc1", "Inv mass of ULS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2)",20,0,20,500,0,0.5 );
880 fInvmassULSmc1->Sumw2();
881 fOutputList->Add(fInvmassULSmc1);
882
883 fInvmassULSmc2 = new TH2D("fInvmassULSmc2", "Inv mass of ULS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2)",20,0,20,500,0,0.5 );
884 fInvmassULSmc2->Sumw2();
885 fOutputList->Add(fInvmassULSmc2);
886
887 fInvmassULSmc3 = new TH2D("fInvmassULSmc3", "Inv mass of ULS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2)",20,0,20,500,0,0.5 );
888 fInvmassULSmc3->Sumw2();
93c189c5 889 fOutputList->Add(fInvmassULSmc3);
890
bfc7c23b 891 fOpeningAngleLS = new TH1F("fOpeningAngleLS","Opening angle for LS pairs",100,0,1);
892 fOutputList->Add(fOpeningAngleLS);
893
894 fOpeningAngleULS = new TH1F("fOpeningAngleULS","Opening angle for ULS pairs",100,0,1);
895 fOutputList->Add(fOpeningAngleULS);
896
897 fPhotoElecPt = new TH1F("fPhotoElecPt", "photonic electron pt",100,0,50);
898 fOutputList->Add(fPhotoElecPt);
899
02cd82a0 900 fPhoElecPt = new TH2F("fPhoElecPt", "Pho-inclusive electron pt",200,0,100,100,0,50);
bfc7c23b 901 fOutputList->Add(fPhoElecPt);
902
02cd82a0 903 fPhoElecPtM20 = new TH2F("fPhoElecPtM20", "Pho-inclusive electron pt with M20",200,0,100,100,0,50);
fb87d707 904 fOutputList->Add(fPhoElecPtM20);
905
02cd82a0 906 fSameElecPt = new TH2F("fSameElecPt", "Same-inclusive electron pt",200,0,100,100,0,50);
f09766dd 907 fOutputList->Add(fSameElecPt);
908
02cd82a0 909 fSameElecPtM20 = new TH2F("fSameElecPtM20", "Same-inclusive electron pt with M20",200,0,100,100,0,50);
fb87d707 910 fOutputList->Add(fSameElecPtM20);
911
02cd82a0 912 fCent = new TH1F("fCent","Centrality",200,0,100) ;
bfc7c23b 913 fOutputList->Add(fCent);
914
915 // Make common binning
44be9e1d 916 const Double_t kMinP = 0.;
bfc7c23b 917 const Double_t kMaxP = 50.;
a6df418c 918 //const Double_t kTPCSigMim = 40.;
919 //const Double_t kTPCSigMax = 140.;
bfc7c23b 920
921 // 1st histogram: TPC dEdx with/without EMCAL (p, pT, TPC Signal, phi, eta, Sig, e/p, ,match, cell, M02, M20, Disp, Centrality, select)
a6df418c 922 Int_t nBins[16] = { 250, 10, 60, 20, 100, 300, 50, 40, 200, 200, 250, 200, 3, 5, 100, 8};
923 Double_t min[16] = {kMinP, -0.5, 1.0, -1.0, -6.0, 0, 0, 0, 0.0, 0.0, 0.0, 0, -1.5, -0.5, 80, -1.5};
6e92c428 924 Double_t max[16] = {kMaxP, 9.5, 4.0, 1.0, 4.0, 3.0, 0.1, 40, 200, 2.0, 50.0, 100, 1.5, 4.5, 180, 6.5};
925 fEleInfo = new THnSparseD("fEleInfo", "Electron Info; pT [GeV/c]; TPC signal;phi;eta;nSig; E/p;Rmatch;Ncell;clsF;M20;mcpT;Centrality;charge;opp;same;trigCond;MCele", 16, nBins, min, max);
85985bb0 926 if(fqahist==1)fOutputList->Add(fEleInfo);
bfc7c23b 927
fb87d707 928 //<--- Trigger info
feffe705 929 /*
fb87d707 930 fClsEBftTrigCut = new TH1F("fClsEBftTrigCut","cluster E before trigger selection",1000,0,100);
931 fOutputList->Add(fClsEBftTrigCut);
932
933 fClsEAftTrigCut = new TH1F("fClsEAftTrigCut","cluster E if cls has 0 trigcut channel",1000,0,100);
934 fOutputList->Add(fClsEAftTrigCut);
935
936 fClsEAftTrigCut1 = new TH1F("fClsEAftTrigCut1","cluster E if cls with trig channel",1000,0,100);
937 fOutputList->Add(fClsEAftTrigCut1);
938
939 fClsEAftTrigCut2 = new TH1F("fClsEAftTrigCut2","cluster E if cls with trigcut channel",1000,0,100);
940 fOutputList->Add(fClsEAftTrigCut2);
941
942 fClsEAftTrigCut3 = new TH1F("fClsEAftTrigCut3","cluster E if cls with trigcut channel + nCell>Ecorrect",1000,0,100);
943 fOutputList->Add(fClsEAftTrigCut3);
944
945 fClsEAftTrigCut4 = new TH1F("fClsEAftTrigCut4","cluster E if cls with trigcut channel + nCell>Ecorrect + cls time cut",1000,0,100);
946 fOutputList->Add(fClsEAftTrigCut4);
947
948 fClsETime = new TH2F("fClsETime", "Cls time vs E; E; time;",1000,0,100,1000,-0.0000002,0.0000009);
949 fOutputList->Add(fClsETime);
950
951 fClsETime1 = new TH2F("fClsETime1", "Cls time vs E if cls contains trigger channel; E; time;",1000,0,100,1000,-0.0000002,0.0000009);
952 fOutputList->Add(fClsETime1);
953
954 fTrigTimes = new TH1F("fTrigTimes", "Trigger time; time; N;",25,0,25);
955 fOutputList->Add(fTrigTimes);
956
42c75692 957 fCellCheck = new TH2F("fCellCheck", "Cell vs E; E GeV; Cell ID",10,6,26,12000,0,12000);
958 fOutputList->Add(fCellCheck);
feffe705 959 */
e7d87aef 960 //<---------- MC
961
02cd82a0 962 fInputHFEMC = new TH2F("fInputHFEMC","Input MC HFE pid electro vs. centrality",200,0,100,100,0,50);
e7d87aef 963 fOutputList->Add(fInputHFEMC);
964
02cd82a0 965 fInputAlle = new TH2F("fInputAlle","Input MC electro vs. centrality",200,0,100,100,0,50);
feffe705 966 fOutputList->Add(fInputAlle);
967
02cd82a0 968 fIncpTMChfe = new TH2F("fIncpTMChfe","MC HFE pid electro vs. centrality",200,0,100,100,0,50);
e7d87aef 969 fOutputList->Add(fIncpTMChfe);
970
02cd82a0 971 fIncpTMChfeAll = new TH2F("fIncpTMChfeAll","MC Alle pid electro vs. centrality",200,0,100,100,0,50);
3db00c72 972 fOutputList->Add(fIncpTMChfeAll);
973
02cd82a0 974 fIncpTMCM20hfe = new TH2F("fIncpTMCM20hfe","MC HFE pid electro vs. centrality with M20",200,0,100,100,0,50);
e7d87aef 975 fOutputList->Add(fIncpTMCM20hfe);
976
02cd82a0 977 fIncpTMCM20hfeAll = new TH2F("fIncpTMCM20hfeAll","MC Alle pid electro vs. centrality with M20",200,0,100,100,0,50);
3db00c72 978 fOutputList->Add(fIncpTMCM20hfeAll);
979
60544aea 980 fIncpTMCM20hfeCheck = new TH2F("fIncpTMCM20hfeCheck","MC HFE pid electro vs. centrality with M20 Check",200,0,100,100,0,50);
981 fOutputList->Add(fIncpTMCM20hfeCheck);
44be9e1d 982
987053ce 983 Int_t nBinspho2[5] = { 200, 100, 7, 3, 100};
984 Double_t minpho2[5] = { 0., 0., -2.5, -0.5, 0.};
985 Double_t maxpho2[5] = {100., 50., 4.5, 2.5, 50.};
44be9e1d 986
987053ce 987 fIncpTMCpho = new THnSparseD("fIncpTMCpho","MC Pho pid electro vs. centrality",5,nBinspho2,minpho2,maxpho2);
e7d87aef 988 fOutputList->Add(fIncpTMCpho);
989
987053ce 990 fIncpTMCM20pho = new THnSparseD("fIncpTMCM20pho","MC Pho pid electro vs. centrality with M20",5,nBinspho2,minpho2,maxpho2);
e7d87aef 991 fOutputList->Add(fIncpTMCM20pho);
992
987053ce 993 fPhoElecPtMC = new THnSparseD("fPhoElecPtMC", "MC Pho-inclusive electron pt",5,nBinspho2,minpho2,maxpho2);
e7d87aef 994 fOutputList->Add(fPhoElecPtMC);
995
987053ce 996 fPhoElecPtMCM20 = new THnSparseD("fPhoElecPtMCM20", "MC Pho-inclusive electron pt with M20",5,nBinspho2,minpho2,maxpho2);
feffe705 997 fOutputList->Add(fPhoElecPtMCM20);
e7d87aef 998
987053ce 999 fSameElecPtMC = new THnSparseD("fSameElecPtMC", "MC Same-inclusive electron pt",5,nBinspho2,minpho2,maxpho2);
e7d87aef 1000 fOutputList->Add(fSameElecPtMC);
1001
987053ce 1002 fSameElecPtMCM20 = new THnSparseD("fSameElecPtMCM20", "MC Same-inclusive electron pt with M20",5,nBinspho2,minpho2,maxpho2);
e7d87aef 1003 fOutputList->Add(fSameElecPtMCM20);
1004
987053ce 1005 fIncpTMCM20pho_pi0e = new THnSparseD("fIncpTMCM20pho_pi0e","MC Pho pi0->e pid electro vs. centrality with M20",5,nBinspho2,minpho2,maxpho2);
1006 fIncpTMCM20pho_pi0e->Sumw2();
e4b0faf2 1007 fOutputList->Add(fIncpTMCM20pho_pi0e);
1008
987053ce 1009 fPhoElecPtMCM20_pi0e = new THnSparseD("fPhoElecPtMCM20_pi0e", "MC Pho-inclusive electron pt with M20 pi0->e",5,nBinspho2,minpho2,maxpho2);
1010 fPhoElecPtMCM20_pi0e->Sumw2();
e4b0faf2 1011 fOutputList->Add(fPhoElecPtMCM20_pi0e);
1012
987053ce 1013 fSameElecPtMCM20_pi0e = new THnSparseD("fSameElecPtMCM20_pi0e", "MC Same-inclusive electron pt pi0->e",5,nBinspho2,minpho2,maxpho2);
1014 fSameElecPtMCM20_pi0e->Sumw2();
e4b0faf2 1015 fOutputList->Add(fSameElecPtMCM20_pi0e);
93c189c5 1016 //
1017 fIncpTMCM20pho_eta = new THnSparseD("fIncpTMCM20pho_eta","MC Pho pi0->e pid electro vs. centrality with M20",5,nBinspho2,minpho2,maxpho2);
1018 fIncpTMCM20pho_eta->Sumw2();
1019 fOutputList->Add(fIncpTMCM20pho_eta);
1020
1021 fPhoElecPtMCM20_eta = new THnSparseD("fPhoElecPtMCM20_eta", "MC Pho-inclusive electron pt with M20 pi0->e",5,nBinspho2,minpho2,maxpho2);
1022 fPhoElecPtMCM20_eta->Sumw2();
1023 fOutputList->Add(fPhoElecPtMCM20_eta);
1024
1025 fSameElecPtMCM20_eta = new THnSparseD("fSameElecPtMCM20_eta", "MC Same-inclusive electron pt pi0->e",5,nBinspho2,minpho2,maxpho2);
1026 fSameElecPtMCM20_eta->Sumw2();
1027 fOutputList->Add(fSameElecPtMCM20_eta);
1028
e4b0faf2 1029
a6df418c 1030 CheckNclust = new TH1D("CheckNclust","cluster check",200,0,200);
1031 fOutputList->Add(CheckNclust);
1032
1033 CheckNits = new TH1D("CheckNits","ITS cluster check",8,-0.5,7.5);
1034 fOutputList->Add(CheckNits);
e7d87aef 1035
d113d7cd 1036 Hpi0pTcheck = new TH2D("Hpi0pTcheck","Pi0 pT from Hijing",100,0,50,3,-0.5,2.5);
1037 fOutputList->Add(Hpi0pTcheck);
1038
19325033 1039 HETApTcheck = new TH2D("HETApTcheck","Eta pT from Hijing",100,0,50,3,-0.5,2.5);
1040 fOutputList->Add(HETApTcheck);
1041
e4b0faf2 1042 HphopTcheck = new TH2D("HphopTcheck","Pho pT from Hijing",100,0,50,3,-0.5,2.5);
1043 fOutputList->Add(HphopTcheck);
1044
93c189c5 1045 fpTCheck = new TH1D("fpTCheck","pT check",500,0,50);
1046 fOutputList->Add(fpTCheck);
1047
50919258 1048 fMomDtoE = new TH2D("fMomDtoE","D->E pT correlations;e p_{T} GeV/c;D p_{T} GeV/c",400,0,40,400,0,40);
1049 fOutputList->Add(fMomDtoE);
d113d7cd 1050
bfc7c23b 1051 PostData(1,fOutputList);
1052}
1053
1054//________________________________________________________________________
1055void AliAnalysisTaskHFECal::Terminate(Option_t *)
1056{
1057 // Info("Terminate");
1058 AliAnalysisTaskSE::Terminate();
1059}
1060
1061//________________________________________________________________________
1062Bool_t AliAnalysisTaskHFECal::ProcessCutStep(Int_t cutStep, AliVParticle *track)
1063{
1064 // Check single track cuts for a given cut step
1065 const Int_t kMCOffset = AliHFEcuts::kNcutStepsMCTrack;
1066 if(!fCFM->CheckParticleCuts(cutStep + kMCOffset, track)) return kFALSE;
1067 return kTRUE;
1068}
1069//_________________________________________
f09766dd 1070//void AliAnalysisTaskHFECal::SelectPhotonicElectron(Int_t itrack, Double_t cent, AliESDtrack *track, Bool_t &fFlagPhotonicElec)
c2801a73 1071//void AliAnalysisTaskHFECal::SelectPhotonicElectron(Int_t itrack, Double_t cent, AliESDtrack *track, Bool_t &fFlagPhotonicElec, Bool_t &fFlagConvinatElec, Double_t nSig)
8806ce6c 1072void AliAnalysisTaskHFECal::SelectPhotonicElectron(Int_t itrack, Double_t cent, AliESDtrack *track, AliStack *stack ,Bool_t &fFlagPhotonicElec, Bool_t &fFlagConvinatElec, Double_t nSig, Double_t shower, Double_t ep, Double_t mce, Double_t w, Int_t ibgevent)
bfc7c23b 1073{
1074 //Identify non-heavy flavour electrons using Invariant mass method
1075
fb87d707 1076 fTrackCuts->SetAcceptKinkDaughters(kFALSE);
feffe705 1077 fTrackCuts->SetRequireTPCRefit(kTRUE);
85985bb0 1078 fTrackCuts->SetRequireITSRefit(kTRUE);
feffe705 1079 fTrackCuts->SetEtaRange(-0.9,0.9);
f09766dd 1080 //fTrackCuts->SetRequireSigmaToVertex(kTRUE);
bfc7c23b 1081 fTrackCuts->SetMaxChi2PerClusterTPC(3.5);
feffe705 1082 fTrackCuts->SetMinNClustersTPC(90);
bfc7c23b 1083
1084 const AliESDVertex *pVtx = fESD->GetPrimaryVertex();
1085
93c189c5 1086 double ptEle = track->Pt(); //add
1087 if(ibgevent==0 && w > 0.0)
1088 {
1089 fpTCheck->Fill(ptEle,w);
1090 }
1091
bfc7c23b 1092 Bool_t flagPhotonicElec = kFALSE;
f09766dd 1093 Bool_t flagConvinatElec = kFALSE;
bfc7c23b 1094
8806ce6c 1095 int p1 = 0;
1096 if(mce==3)
1097 {
1098 Int_t label = TMath::Abs(track->GetLabel());
1099 TParticle* particle = stack->Particle(label);
1100 p1 = particle->GetFirstMother();
1101 }
1102
5e65985c 1103 //for(Int_t jTracks = itrack+1; jTracks<fESD->GetNumberOfTracks(); jTracks++){
1104 for(Int_t jTracks = 0; jTracks<fESD->GetNumberOfTracks(); jTracks++){
bfc7c23b 1105 AliESDtrack* trackAsso = fESD->GetTrack(jTracks);
1106 if (!trackAsso) {
1107 printf("ERROR: Could not receive track %d\n", jTracks);
1108 continue;
1109 }
5e65985c 1110 if(itrack==jTracks)continue;
d86a815c 1111 int jbgevent = 0;
5e65985c 1112
8806ce6c 1113 int p2 = 0;
1114 if(mce==3)
1115 {
1116 Int_t label2 = TMath::Abs(trackAsso->GetLabel());
1117 TParticle* particle2 = stack->Particle(label2);
d86a815c 1118 Bool_t MChijing_ass = fMC->IsFromBGEvent(label2);
1119 if(MChijing_ass)jbgevent =1;
8806ce6c 1120 if(particle2->GetFirstMother()>-1)
1121 p2 = particle2->GetFirstMother();
1122 }
1123
f09766dd 1124 Double_t dEdxAsso = -999., ptPrim=-999., ptAsso=-999., openingAngle = -999.;
bfc7c23b 1125 Double_t mass=999., width = -999;
1126 Bool_t fFlagLS=kFALSE, fFlagULS=kFALSE;
1127
93c189c5 1128 //ptPrim = track->Pt();
1129 ptPrim = ptEle;
f09766dd 1130
bfc7c23b 1131 dEdxAsso = trackAsso->GetTPCsignal();
1132 ptAsso = trackAsso->Pt();
1133 Int_t chargeAsso = trackAsso->Charge();
1134 Int_t charge = track->Charge();
1135
b12bc641 1136
3db00c72 1137 if(ptAsso <0.5) continue;
bfc7c23b 1138 if(!fTrackCuts->AcceptTrack(trackAsso)) continue;
f09766dd 1139 if(dEdxAsso <65 || dEdxAsso>100) continue; //11a pass1
bfc7c23b 1140
1141 Int_t fPDGe1 = 11; Int_t fPDGe2 = 11;
1142 if(charge>0) fPDGe1 = -11;
1143 if(chargeAsso>0) fPDGe2 = -11;
b12bc641 1144
1145 //printf("chargeAsso = %d\n",chargeAsso);
1146 //printf("charge = %d\n",charge);
bfc7c23b 1147 if(charge == chargeAsso) fFlagLS = kTRUE;
1148 if(charge != chargeAsso) fFlagULS = kTRUE;
1149
b12bc641 1150 //printf("fFlagLS = %d\n",fFlagLS);
1151 //printf("fFlagULS = %d\n",fFlagULS);
1152 //printf("\n");
1153
bfc7c23b 1154 AliKFParticle ge1(*track, fPDGe1);
1155 AliKFParticle ge2(*trackAsso, fPDGe2);
1156 AliKFParticle recg(ge1, ge2);
1157
1158 if(recg.GetNDF()<1) continue;
1159 Double_t chi2recg = recg.GetChi2()/recg.GetNDF();
1160 if(TMath::Sqrt(TMath::Abs(chi2recg))>3.) continue;
1161
d86a815c 1162 // for v5-03-70-AN comment
bfc7c23b 1163 AliKFVertex primV(*pVtx);
1164 primV += recg;
1165 recg.SetProductionVertex(primV);
1166
d113d7cd 1167 recg.SetMassConstraint(0,0.0001);
d86a815c 1168
bfc7c23b 1169 openingAngle = ge1.GetAngle(ge2);
1170 if(fFlagLS) fOpeningAngleLS->Fill(openingAngle);
1171 if(fFlagULS) fOpeningAngleULS->Fill(openingAngle);
1172
bfc7c23b 1173
1174 recg.GetMass(mass,width);
1175
c2801a73 1176 double ishower = 0;
1177 if(shower>0.0 && shower<0.3)ishower = 1;
1178
b12bc641 1179 double phoinfo[9];
f09766dd 1180 phoinfo[0] = cent;
1181 phoinfo[1] = ptPrim;
1182 phoinfo[2] = mass;
3db00c72 1183 phoinfo[3] = nSig;
4e01b68c 1184 phoinfo[4] = openingAngle;
c2801a73 1185 phoinfo[5] = ishower;
1186 phoinfo[6] = ep;
1187 phoinfo[7] = mce;
b12bc641 1188 phoinfo[8] = ptAsso;
f09766dd 1189
1190 if(fFlagLS) fInvmassLS->Fill(phoinfo);
1191 if(fFlagULS) fInvmassULS->Fill(phoinfo);
d86a815c 1192 if(fFlagLS && ibgevent==0 && jbgevent==0) fInvmassLSmc->Fill(phoinfo,w);
93c189c5 1193 if(fFlagULS && ibgevent==0 && jbgevent==0)
1194 {
1195 fInvmassULSmc->Fill(phoinfo,w);
1196 }
b12bc641 1197 //printf("fInvmassCut %f\n",fInvmassCut);
1198 //printf("openingAngle %f\n",fOpeningAngleCut);
1199
ffa14dda 1200 if(openingAngle > fOpeningAngleCut) continue;
bfc7c23b 1201
d86a815c 1202 // for real data
93c189c5 1203 //printf("mce =%f\n",mce);
d86a815c 1204 if(mce<-0.5) // mce==-1. is real
1205 {
93c189c5 1206 //printf("Real data\n");
d86a815c 1207 if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec){
1208 //if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec && (p1==p2)){ <--- only MC train (55,56) v5-03-68-AN & 69 for check
1209 flagPhotonicElec = kTRUE;
1210 }
1211 if(mass<fInvmassCut && fFlagLS && !flagConvinatElec){
1212 flagConvinatElec = kTRUE;
1213 }
1214 }
1215 // for MC data
1216 else
1217 {
93c189c5 1218 //printf("MC data\n");
1219
1220 if(w>0.0)
1221 {
1222 if(fFlagLS && ibgevent==0 && jbgevent==0) fInvmassLSmc0->Fill(ptPrim,mass,w);
1223 if(fFlagULS && ibgevent==0 && jbgevent==0) fInvmassULSmc0->Fill(ptPrim,mass,w);
1224 if(fFlagLS && ibgevent==0 && jbgevent==0) fInvmassLSmc1->Fill(ptPrim,mass);
1225 if(fFlagULS && ibgevent==0 && jbgevent==0) fInvmassULSmc1->Fill(ptPrim,mass);
1226 if(fFlagLS && ibgevent==0 && jbgevent==0 && (p1==p2)) fInvmassLSmc2->Fill(ptPrim,mass,w);
1227 if(fFlagULS && ibgevent==0 && jbgevent==0 && (p1==p2)) fInvmassULSmc2->Fill(ptPrim,mass,w);
1228 if(fFlagLS && ibgevent==0 && jbgevent==0 && (p1==p2)) fInvmassLSmc3->Fill(ptPrim,mass);
1229 if(fFlagULS && ibgevent==0 && jbgevent==0 && (p1==p2)) fInvmassULSmc3->Fill(ptPrim,mass);
1230 }
d86a815c 1231 if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec && (ibgevent==jbgevent)){
1232 //if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec && (p1==p2)){ <--- only MC train (55,56) v5-03-68-AN & 69 for check
1233 flagPhotonicElec = kTRUE;
1234 }
1235 if(mass<fInvmassCut && fFlagLS && !flagConvinatElec && (ibgevent==jbgevent)){
1236 flagConvinatElec = kTRUE;
1237 }
1238 }
1239
bfc7c23b 1240 }
1241 fFlagPhotonicElec = flagPhotonicElec;
f09766dd 1242 fFlagConvinatElec = flagConvinatElec;
bfc7c23b 1243
1244}
8806ce6c 1245//-------------------------------------------
1246//void AliAnalysisTaskHFECal::GetMCweight(double mcPi0pT)
1247double AliAnalysisTaskHFECal::GetMCweight(double mcPi0pT)
1248{
1249 double weight = 1.0;
1250
1251 if(mcPi0pT>0.0 && mcPi0pT<5.0)
1252 {
1253 weight = 0.323*mcPi0pT/(TMath::Exp(-1.6+0.767*mcPi0pT+0.0285*mcPi0pT*mcPi0pT));
1254 }
1255 else
1256 {
1257 weight = 115.0/(0.718*mcPi0pT*TMath::Power(mcPi0pT,3.65));
1258 }
1259 return weight;
1260}
fb87d707 1261
93c189c5 1262double AliAnalysisTaskHFECal::GetMCweightEta(double mcEtapT)
1263{
1264 double weight = 1.0;
1265
1266 weight = 223.3/TMath::Power((TMath::Exp(-0.17*mcEtapT-0.0322*mcEtapT*mcEtapT)+mcEtapT/1.69),5.65);
1267 return weight;
1268}
1269
1270
fb87d707 1271//_________________________________________
1272void AliAnalysisTaskHFECal::FindTriggerClusters()
1273{
1274 // constants
1275 const int nModuleCols = 2;
1276 const int nModuleRows = 5;
1277 const int nColsFeeModule = 48;
1278 const int nRowsFeeModule = 24;
1279 const int nColsFaltroModule = 24;
1280 const int nRowsFaltroModule = 12;
1281 //const int faltroWidthMax = 20;
1282
1283 // part 1, trigger extraction -------------------------------------
1284 Int_t globCol, globRow;
1285 //Int_t ntimes=0, nTrigChannel=0, nTrigChannelCut=0, trigInCut=0;
1286 Int_t ntimes=0, nTrigChannel=0, nTrigChannelCut=0;
1287
1288 //Int_t trigtimes[faltroWidthMax];
1289 Double_t cellTime[nColsFeeModule*nModuleCols][nRowsFeeModule*nModuleRows];
1290 Double_t cellEnergy[nColsFeeModule*nModuleCols][nRowsFeeModule*nModuleRows];
1291 //Double_t fTrigCutLow = 6;
1292 //Double_t fTrigCutHigh = 10;
1293 Double_t fTimeCutLow = 469e-09;
1294 Double_t fTimeCutHigh = 715e-09;
1295
1296 AliESDCaloTrigger * fCaloTrigger = fESD->GetCaloTrigger( "EMCAL" );
1297
1298 // erase trigger maps
1299 for(Int_t i = 0; i < nColsFaltroModule*nModuleCols; i++ )
1300 {
1301 for(Int_t j = 0; j < nRowsFaltroModule*nModuleRows; j++ )
1302 {
1303 ftriggersCut[i][j] = 0;
1304 ftriggers[i][j] = 0;
1305 ftriggersTime[i][j] = 0;
1306 }
1307 }
1308
1309 Int_t iglobCol=0, iglobRow=0;
1310 // go through triggers
1311 if( fCaloTrigger->GetEntries() > 0 )
1312 {
1313 // needs reset
1314 fCaloTrigger->Reset();
1315 while( fCaloTrigger->Next() )
1316 {
1317 fCaloTrigger->GetPosition( globCol, globRow );
1318 fCaloTrigger->GetNL0Times( ntimes );
1319 /*
1320 // no L0s
1321 if( ntimes < 1 ) continue;
1322 // get precise timings
1323 fCaloTrigger->GetL0Times( trigtimes );
1324 trigInCut = 0;
1325 for(Int_t i = 0; i < ntimes; i++ )
1326 {
1327 // save the first trigger time in channel
1328 if( i == 0 || triggersTime[globCol][globRow] > trigtimes[i] )
1329 triggersTime[globCol][globRow] = trigtimes[i];
1330 //printf("trigger times: %d\n",trigtimes[i]);
1331 // check if in cut
1332 if(trigtimes[i] > fTrigCutLow && trigtimes[i] < fTrigCutHigh )
1333 trigInCut = 1;
1334
1335 fTrigTimes->Fill(trigtimes[i]);
1336 }
1337 */
1338
1339 //L1 analysis from AliAnalysisTaskEMCALTriggerQA
1340 Int_t bit = 0;
1341 fCaloTrigger->GetTriggerBits(bit);
1342
1343 Int_t ts = 0;
1344 fCaloTrigger->GetL1TimeSum(ts);
1345 if (ts > 0)ftriggers[globCol][globRow] = 1;
1346 // number of triggered channels in event
1347 nTrigChannel++;
1348 // ... inside cut
1349 if(ts>0 && (bit >> 6 & 0x1))
1350 {
1351 iglobCol = globCol;
1352 iglobRow = globRow;
1353 nTrigChannelCut++;
1354 ftriggersCut[globCol][globRow] = 1;
1355 }
1356
1357 } // calo trigger entries
1358 } // has calo trigger entries
1359
1360 // part 2 go through the clusters here -----------------------------------
1361 Int_t nCluster=0, nCell=0, iCell=0, gCell=0;
1362 Short_t cellAddr, nSACell, mclabel;
1363 //Int_t nSACell, iSACell, mclabel;
1364 Int_t iSACell;
42c75692 1365 Double_t cellAmp=0, cellTimeT=0, clusterTime=0, efrac=0;
fb87d707 1366 Int_t nSupMod, nModule, nIphi, nIeta, iphi, ieta, gphi, geta, feta, fphi;
1367
1368 TRefArray *fCaloClusters = new TRefArray();
1369 fESD->GetEMCALClusters( fCaloClusters );
1370 nCluster = fCaloClusters->GetEntries();
1371
1372
1373 // save all cells times if there are clusters
1374 if( nCluster > 0 ){
1375 // erase time map
1376 for(Int_t i = 0; i < nColsFeeModule*nModuleCols; i++ ){
1377 for(Int_t j = 0; j < nRowsFeeModule*nModuleRows; j++ ){
1378 cellTime[i][j] = 0.;
1379 cellEnergy[i][j] = 0.;
1380 }
1381 }
1382
1383 // get cells
1384 AliESDCaloCells *fCaloCells = fESD->GetEMCALCells();
1385 //AliVCaloCells fCaloCells = fESD->GetEMCALCells();
1386 nSACell = fCaloCells->GetNumberOfCells();
1387 for(iSACell = 0; iSACell < nSACell; iSACell++ ){
1388 // get the cell info *fCal
1389 fCaloCells->GetCell( iSACell, cellAddr, cellAmp, cellTimeT , mclabel, efrac);
1390
1391 // get cell position
1392 fGeom->GetCellIndex( cellAddr, nSupMod, nModule, nIphi, nIeta );
1393 fGeom->GetCellPhiEtaIndexInSModule( nSupMod,nModule, nIphi, nIeta, iphi, ieta);
1394
1395 // convert co global phi eta
1396 gphi = iphi + nRowsFeeModule*(nSupMod/2);
1397 geta = ieta + nColsFeeModule*(nSupMod%2);
1398
1399 // save cell time and energy
1400 cellTime[geta][gphi] = cellTimeT;
1401 cellEnergy[geta][gphi] = cellAmp;
1402
1403 }
1404 }
1405
1406 Int_t nClusterTrig, nClusterTrigCut;
1407 UShort_t *cellAddrs;
1408 Double_t clsE=-999, clsEta=-999, clsPhi=-999;
1409 Float_t clsPos[3] = {0.,0.,0.};
1410
1411 for(Int_t icl=0; icl<fESD->GetNumberOfCaloClusters(); icl++)
1412 {
1413 AliESDCaloCluster *cluster = fESD->GetCaloCluster(icl);
1414 if(!cluster || !cluster->IsEMCAL()) continue;
1415
1416 // get cluster cells
1417 nCell = cluster->GetNCells();
1418
1419 // get cluster energy
1420 clsE = cluster->E();
1421
1422 // get cluster position
1423 cluster->GetPosition(clsPos);
1424 TVector3 clsPosVec(clsPos[0],clsPos[1],clsPos[2]);
1425 clsEta = clsPosVec.Eta();
1426 clsPhi = clsPosVec.Phi();
1427
1428 // get the cell addresses
1429 cellAddrs = cluster->GetCellsAbsId();
1430
1431 // check if the cluster contains cell, that was marked as triggered
1432 nClusterTrig = 0;
1433 nClusterTrigCut = 0;
1434
1435 // loop the cells to check, if cluser in acceptance
1436 // any cluster with a cell outside acceptance is not considered
1437 for( iCell = 0; iCell < nCell; iCell++ )
1438 {
42c75692 1439 // check hot cell
feffe705 1440 //if(clsE>6.0)fCellCheck->Fill(clsE,cellAddrs[iCell]);
42c75692 1441
fb87d707 1442 // get cell position
1443 fGeom->GetCellIndex( cellAddrs[iCell], nSupMod, nModule, nIphi, nIeta );
1444 fGeom->GetCellPhiEtaIndexInSModule( nSupMod,nModule, nIphi, nIeta, iphi, ieta);
1445
1446 // convert co global phi eta
1447 gphi = iphi + nRowsFeeModule*(nSupMod/2);
1448 geta = ieta + nColsFeeModule*(nSupMod%2);
1449
1450 if( cellTime[geta][gphi] > 0. ){
1451 clusterTime += cellTime[geta][gphi];
1452 gCell++;
1453 }
1454
1455 // get corresponding FALTRO
1456 fphi = gphi / 2;
1457 feta = geta / 2;
1458
1459 // try to match with a triggered
1460 if( ftriggers[feta][fphi]==1)
1461 { nClusterTrig++;
1462 }
1463 if( ftriggersCut[feta][fphi]==1)
1464 { nClusterTrigCut++;
1465 }
1466
1467 } // cells
1468
1469
1470 if( gCell > 0 )
1471 clusterTime = clusterTime / (Double_t)gCell;
1472 // fix the reconstriction code time 100ns jumps
1473 if( fESD->GetBunchCrossNumber() % 4 < 2 )
1474 clusterTime -= 0.0000001;
1475
feffe705 1476 //fClsETime->Fill(clsE,clusterTime);
1477 //fClsEBftTrigCut->Fill(clsE);
fb87d707 1478
1479 if(nClusterTrig>0){
feffe705 1480 //fClsETime1->Fill(clsE,clusterTime);
fb87d707 1481 }
1482
1483 if(nClusterTrig>0){
1484 cluster->SetChi2(1);
feffe705 1485 //fClsEAftTrigCut1->Fill(clsE);
fb87d707 1486 }
1487
1488 if(nClusterTrigCut>0){
1489 cluster->SetChi2(2);
feffe705 1490 //fClsEAftTrigCut2->Fill(clsE);
fb87d707 1491 }
1492
1493 if(nClusterTrigCut>0 && ( nCell > (1 + clsE / 3)))
1494 {
1495 cluster->SetChi2(3);
feffe705 1496 //fClsEAftTrigCut3->Fill(clsE);
fb87d707 1497 }
1498
1499 if(nClusterTrigCut>0 && (nCell > (1 + clsE / 3) )&&( clusterTime > fTimeCutLow && clusterTime < fTimeCutHigh ))
1500 {
1501 // cluster->SetChi2(4);
feffe705 1502 //fClsEAftTrigCut4->Fill(clsE);
fb87d707 1503 }
1504 if(nClusterTrigCut<1)
1505 {
1506 cluster->SetChi2(0);
1507
feffe705 1508 //fClsEAftTrigCut->Fill(clsE);
fb87d707 1509 }
1510
1511 } // clusters
1512}
1513
1514
1515
1516
1517
1518
1519