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