git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
-updates on user addtasks
[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
e4b0faf2 510 // access grand parent pi0->g->e
511 TParticle* particle_parent = stack->Particle(parentlabel); // get parent pointer
512 if(particle_parent->GetFirstMother()>-1 && parentPID==22 && fabs(mcpid)==11) // get grand parent g->e
513 {
46305ed6 514 //int grand_parentPID = stack->Particle(particle_parent->GetFirstMother())->GetPdgCode();
515 //double pTtmp = stack->Particle(particle_parent->GetFirstMother())->Pt();
516 FindMother(particle_parent, grand_parentlabel, grand_parentPID);
517 double mcGrandpT = stack->Particle(grand_parentlabel)->Pt();
518 if(grand_parentPID==111 && mcGrandpT>0.0)
e4b0faf2 519 {
46305ed6 520 // check eta->pi0 decay !
521 int grand2_parentlabel = 99999; int grand2_parentPID = 99999;
522 TParticle* particle_grand = stack->Particle(grand_parentlabel); // get parent pointer
523 FindMother(particle_grand, grand2_parentlabel, grand2_parentPID);
524 if(grand2_parentPID==221)
525 {
526 double mcGrandpT2 = stack->Particle(grand2_parentlabel)->Pt();
527 mcOrgEta = kTRUE;
528 mcWeight = GetMCweight(mcGrandpT2);
529 }
530 else
531 {
532 mcOrgPi0 = kTRUE;
533 mcWeight = GetMCweight(mcGrandpT);
534 }
93c189c5 535 }
46305ed6 536 if(grand_parentPID==221 && mcGrandpT>0.0)
93c189c5 537 {
538 mcOrgEta = kTRUE;
46305ed6 539 mcWeight = GetMCweightEta(mcGrandpT);
e4b0faf2 540 }
541 }
e7d87aef 542 }
9b3495ff 543 if(fabs(mcpid)==11 && mcDtoE)mcele= 1.;
544 if(fabs(mcpid)==11 && mcBtoE)mcele= 2.;
545 if(fabs(mcpid)==11 && mcPho)mcele= 3.;
e7d87aef 546 }
547
50919258 548 //printf("weight = %f\n",mcWeight);
549
9b3495ff 550 if(TMath::Abs(track->Eta())>0.6) continue;
f09766dd 551 if(TMath::Abs(track->Pt()<2.0)) continue;
bfc7c23b 552
553 fTrackPtBefTrkCuts->Fill(track->Pt());
554 // RecKine: ITSTPC cuts
555 if(!ProcessCutStep(AliHFEcuts::kStepRecKineITSTPC, track)) continue;
556
557 //RecKink
558 if(fRejectKinkMother) { // Quick and dirty fix to reject both kink mothers and daughters
559 if(track->GetKinkIndex(0) != 0) continue;
560 }
561
562 // RecPrim
563 if(!ProcessCutStep(AliHFEcuts::kStepRecPrim, track)) continue;
564
565 // HFEcuts: ITS layers cuts
566 if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsITS, track)) continue;
567
f09766dd 568 // HFE cuts: TPC PID cleanup
569 if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTPC, track)) continue;
570
a6df418c 571 int nTPCcl = track->GetTPCNcls();
6e92c428 572 int nTPCclF = track->GetTPCNclsF();
a6df418c 573 int nITS = track->GetNcls(0);
bfc7c23b 574
575 fTrackPtAftTrkCuts->Fill(track->Pt());
576
577 Double_t mom = -999., eop=-999., pt = -999., dEdx=-999., fTPCnSigma=-10, phi=-999., eta=-999.;
feffe705 578 pt = track->Pt();
579 if(pt<2.0)continue;
bfc7c23b 580
581 // Track extrapolation
582
bfc7c23b 583 Int_t charge = track->Charge();
584 fTrkpt->Fill(pt);
585 mom = track->P();
586 phi = track->Phi();
587 eta = track->Eta();
588 dEdx = track->GetTPCsignal();
589 fTPCnSigma = fPID->GetPIDResponse() ? fPID->GetPIDResponse()->NumberOfSigmasTPC(track, AliPID::kElectron) : 1000;
590
591 double ncells = -1.0;
592 double m20 = -1.0;
593 double m02 = -1.0;
594 double disp = -1.0;
595 double rmatch = -1.0;
596 double nmatch = -1.0;
f09766dd 597 double oppstatus = 0.0;
bfc7c23b 598
f09766dd 599 Bool_t fFlagPhotonicElec = kFALSE;
600 Bool_t fFlagConvinatElec = kFALSE;
bfc7c23b 601
602 Int_t clsId = track->GetEMCALcluster();
603 if (clsId>0){
604 AliESDCaloCluster *clust = fESD->GetCaloCluster(clsId);
605 if(clust && clust->IsEMCAL()){
606
607 double clustE = clust->E();
608 eop = clustE/fabs(mom);
609 //double clustT = clust->GetTOF();
610 ncells = clust->GetNCells();
611 m02 = clust->GetM02();
612 m20 = clust->GetM20();
613 disp = clust->GetDispersion();
614 double delphi = clust->GetTrackDx();
615 double deleta = clust->GetTrackDz();
616 rmatch = sqrt(pow(delphi,2)+pow(deleta,2));
617 nmatch = clust->GetNTracksMatched();
618
c2801a73 619 if(fTPCnSigma>-1.5 && fTPCnSigma<3.0)
620 {
46305ed6 621 SelectPhotonicElectron(iTracks,cent,track,fFlagPhotonicElec,fFlagConvinatElec,fTPCnSigma,m20,eop,mcele,mcWeight,iHijing);
c2801a73 622 }
623 if(fFlagPhotonicElec)oppstatus = 1.0;
624 if(fFlagConvinatElec)oppstatus = 2.0;
625 if(fFlagPhotonicElec && fFlagConvinatElec)oppstatus = 3.0;
626
feffe705 627 double valdedx[16];
a6df418c 628 valdedx[0] = pt; valdedx[1] = nITS; valdedx[2] = phi; valdedx[3] = eta; valdedx[4] = fTPCnSigma;
6e92c428 629 valdedx[5] = eop; valdedx[6] = rmatch; valdedx[7] = ncells, valdedx[8] = nTPCclF; valdedx[9] = m20; valdedx[10] = mcpT;
a6df418c 630 valdedx[11] = cent; valdedx[12] = charge; valdedx[13] = oppstatus; valdedx[14] = nTPCcl;
feffe705 631 valdedx[15] = mcele;
85985bb0 632 if(fqahist==1)fEleInfo->Fill(valdedx);
bfc7c23b 633
634
635 }
636 }
a6df418c 637
638 if(nITS<2.5)continue;
639 if(nTPCcl<100)continue;
640
641 CheckNclust->Fill(nTPCcl);
642 CheckNits->Fill(nITS);
643
42c75692 644 fdEdxBef->Fill(mom,fTPCnSigma);
bfc7c23b 645 fTPCnsigma->Fill(mom,fTPCnSigma);
f09766dd 646 if(fTPCnSigma >= -1.0 && fTPCnSigma <= 3)fTrkEovPBef->Fill(pt,eop);
bfc7c23b 647
bfc7c23b 648 Int_t pidpassed = 1;
649
650
651 //--- track accepted
652 AliHFEpidObject hfetrack;
653 hfetrack.SetAnalysisType(AliHFEpidObject::kESDanalysis);
654 hfetrack.SetRecTrack(track);
3db00c72 655 double binct = 10.5;
656 if((0.0< cent) && (cent<5.0)) binct = 0.5;
657 if((5.0< cent) && (cent<10.0)) binct = 1.5;
658 if((10.0< cent) && (cent<20.0)) binct = 2.5;
659 if((20.0< cent) && (cent<30.0)) binct = 3.5;
660 if((30.0< cent) && (cent<40.0)) binct = 4.5;
661 if((40.0< cent) && (cent<50.0)) binct = 5.5;
662 if((50.0< cent) && (cent<60.0)) binct = 6.5;
663 if((60.0< cent) && (cent<70.0)) binct = 7.5;
664 if((70.0< cent) && (cent<80.0)) binct = 8.5;
665 if((80.0< cent) && (cent<90.0)) binct = 9.5;
666 if((90.0< cent) && (cent<100.0)) binct = 10.5;
667
668 hfetrack.SetCentrality((int)binct); //added
bfc7c23b 669 hfetrack.SetPbPb();
670 if(!fPID->IsSelected(&hfetrack, NULL, "", fPIDqa)) pidpassed = 0;
671
672 if(pidpassed==0) continue;
673
674 fTrkEovPAft->Fill(pt,eop);
42c75692 675 fdEdxAft->Fill(mom,fTPCnSigma);
bfc7c23b 676
e7d87aef 677 // Fill real data
fb87d707 678 fIncpT->Fill(cent,pt);
bfc7c23b 679 if(fFlagPhotonicElec) fPhoElecPt->Fill(cent,pt);
f09766dd 680 if(fFlagConvinatElec) fSameElecPt->Fill(cent,pt);
fb87d707 681
682 if(m20>0.0 && m20<0.3)
683 {
684 fIncpTM20->Fill(cent,pt);
685 if(fFlagPhotonicElec) fPhoElecPtM20->Fill(cent,pt);
686 if(fFlagConvinatElec) fSameElecPtM20->Fill(cent,pt);
687 }
feffe705 688
e7d87aef 689 // MC
e4b0faf2 690 if(mcele>0) // select MC electrons
e7d87aef 691 {
3db00c72 692
693 fIncpTMChfeAll->Fill(cent,pt);
694 if(m20>0.0 && m20<0.3)fIncpTMCM20hfeAll->Fill(cent,pt);
695
e4b0faf2 696 if(mcBtoE || mcDtoE) // select B->e & D->e
e7d87aef 697 {
698 fIncpTMChfe->Fill(cent,pt);
699 if(m20>0.0 && m20<0.3)fIncpTMCM20hfe->Fill(cent,pt);
60544aea 700 if(m20>0.0 && m20<0.3)fIncpTMCM20hfeCheck->Fill(cent,mcpT);
e7d87aef 701 }
e4b0faf2 702
703 if(mcPho) // select photonic electrons
e7d87aef 704 {
987053ce 705 double phoval[5];
44be9e1d 706 phoval[0] = cent;
707 phoval[1] = pt;
708 phoval[2] = fTPCnSigma;
d113d7cd 709 phoval[3] = iHijing;
987053ce 710 phoval[4] = mcMompT;
44be9e1d 711
712 fIncpTMCpho->Fill(phoval);
713 if(fFlagPhotonicElec) fPhoElecPtMC->Fill(phoval);
714 if(fFlagConvinatElec) fSameElecPtMC->Fill(phoval);
e7d87aef 715
716 if(m20>0.0 && m20<0.3)
717 {
44be9e1d 718 fIncpTMCM20pho->Fill(phoval);
719 if(fFlagPhotonicElec) fPhoElecPtMCM20->Fill(phoval);
720 if(fFlagConvinatElec) fSameElecPtMCM20->Fill(phoval);
e4b0faf2 721 // pi0->g->e
722 if(mcWeight>-1)
723 {
987053ce 724 if(iHijing==1)mcWeight = 1.0;
93c189c5 725 if(mcOrgPi0)
726 {
727 fIncpTMCM20pho_pi0e->Fill(phoval,mcWeight);
728 if(fFlagPhotonicElec) fPhoElecPtMCM20_pi0e->Fill(phoval,mcWeight);
729 if(fFlagConvinatElec) fSameElecPtMCM20_pi0e->Fill(phoval,mcWeight);
730 }
731 if(mcOrgEta)
732 {
733 fIncpTMCM20pho_eta->Fill(phoval,mcWeight);
734 if(fFlagPhotonicElec) fPhoElecPtMCM20_eta->Fill(phoval,mcWeight);
735 if(fFlagConvinatElec) fSameElecPtMCM20_eta->Fill(phoval,mcWeight);
736 }
737 // --- eta
e4b0faf2 738 }
e7d87aef 739 }
740 }
741 }
bfc7c23b 742 }
743 PostData(1, fOutputList);
744}
745//_________________________________________
746void AliAnalysisTaskHFECal::UserCreateOutputObjects()
747{
42c75692 748 //--- Check MC
749
e7d87aef 750 //Bool_t mcData = kFALSE;
42c75692 751 if(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler())
752 {
e7d87aef 753 fmcData = kTRUE;
42c75692 754 printf("+++++ MC Data available");
755 }
e7d87aef 756 if(fmcData)
42c75692 757 {
758 printf("++++++++= MC analysis \n");
759 }
760 else
761 {
762 printf("++++++++ real data analysis \n");
763 }
764
85985bb0 765 printf("+++++++ QA hist %d",fqahist);
766
fb87d707 767 //---- Geometry
768 fGeom = AliEMCALGeometry::GetInstance("EMCAL_COMPLETEV1");
769
bfc7c23b 770 //--------Initialize PID
42c75692 771 //fPID->SetHasMCData(kFALSE);
e7d87aef 772 fPID->SetHasMCData(fmcData);
bfc7c23b 773 if(!fPID->GetNumberOfPIDdetectors())
774 {
775 fPID->AddDetector("TPC", 0);
776 fPID->AddDetector("EMCAL", 1);
777 }
778
3db00c72 779 Double_t params[4];
78ff1095 780 const char *cutmodel;
3db00c72 781 cutmodel = "pol0";
782 params[0] = -1.0; //sigma min
9b3495ff 783 double maxnSig = 3.0;
784 if(fmcData)
785 {
786 params[0] = -5.0; //sigma min
787 maxnSig = 5.0;
788 }
789
790 for(Int_t a=0;a<11;a++)fPID->ConfigureTPCcentralityCut(a,cutmodel,params,maxnSig);
3db00c72 791
bfc7c23b 792 fPID->SortDetectors();
793 fPIDqa = new AliHFEpidQAmanager();
794 fPIDqa->Initialize(fPID);
a6df418c 795
796 //------- fcut --------------
797 fCuts = new AliHFEcuts();
798 fCuts->CreateStandardCuts();
799 //fCuts->SetMinNClustersTPC(100);
800 fCuts->SetMinNClustersTPC(90);
801 fCuts->SetMinRatioTPCclusters(0.6);
802 fCuts->SetTPCmodes(AliHFEextraCuts::kFound, AliHFEextraCuts::kFoundOverFindable);
803 //fCuts->SetMinNClustersITS(3);
804 fCuts->SetMinNClustersITS(2);
805 fCuts->SetCutITSpixel(AliHFEextraCuts::kAny);
806 fCuts->SetCheckITSLayerStatus(kFALSE);
807 fCuts->SetVertexRange(10.);
808 fCuts->SetTOFPIDStep(kFALSE);
809 fCuts->SetPtRange(2, 50);
810 fCuts->SetMaxImpactParam(3.,3.);
811
bfc7c23b 812 //--------Initialize correction Framework and Cuts
813 fCFM = new AliCFManager;
814 const Int_t kNcutSteps = AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kNcutStepsRecTrack + AliHFEcuts::kNcutStepsDETrack;
815 fCFM->SetNStepParticle(kNcutSteps);
816 for(Int_t istep = 0; istep < kNcutSteps; istep++)
817 fCFM->SetParticleCutsList(istep, NULL);
818
819 if(!fCuts){
820 AliWarning("Cuts not available. Default cuts will be used");
821 fCuts = new AliHFEcuts;
822 fCuts->CreateStandardCuts();
823 }
824 fCuts->Initialize(fCFM);
825
826 //---------Output Tlist
827 fOutputList = new TList();
828 fOutputList->SetOwner();
829 fOutputList->Add(fPIDqa->MakeList("PIDQA"));
830
831 fNoEvents = new TH1F("fNoEvents","",4,-0.5,3.5) ;
832 fOutputList->Add(fNoEvents);
833
02cd82a0 834 Int_t binsE[5] = {250, 100, 1000, 200, 10};
fb87d707 835 Double_t xminE[5] = {1.0, -1, 0.0, 0, -0.5};
836 Double_t xmaxE[5] = {3.5, 1, 100.0, 100, 9.5};
837 fEMCAccE = new THnSparseD("fEMCAccE","EMC acceptance & E;#eta;#phi;Energy;Centrality;trugCondition;",5,binsE,xminE,xmaxE);
85985bb0 838 if(fqahist==1)fOutputList->Add(fEMCAccE);
bfc7c23b 839
840 fTrkpt = new TH1F("fTrkpt","track pt",100,0,50);
841 fOutputList->Add(fTrkpt);
842
843 fTrackPtBefTrkCuts = new TH1F("fTrackPtBefTrkCuts","track pt before track cuts",100,0,50);
844 fOutputList->Add(fTrackPtBefTrkCuts);
845
846 fTrackPtAftTrkCuts = new TH1F("fTrackPtAftTrkCuts","track pt after track cuts",100,0,50);
847 fOutputList->Add(fTrackPtAftTrkCuts);
848
849 fTPCnsigma = new TH2F("fTPCnsigma", "TPC - n sigma",100,0,50,200,-10,10);
850 fOutputList->Add(fTPCnsigma);
851
852 fTrkEovPBef = new TH2F("fTrkEovPBef","track E/p before HFE pid",100,0,50,100,0,2);
853 fOutputList->Add(fTrkEovPBef);
854
855 fTrkEovPAft = new TH2F("fTrkEovPAft","track E/p after HFE pid",100,0,50,100,0,2);
856 fOutputList->Add(fTrkEovPAft);
857
42c75692 858 fdEdxBef = new TH2F("fdEdxBef","track dEdx vs p before HFE pid",100,0,50,200,-10,10);
bfc7c23b 859 fOutputList->Add(fdEdxBef);
860
42c75692 861 fdEdxAft = new TH2F("fdEdxAft","track dEdx vs p after HFE pid",100,0,50,200,-10,10);
bfc7c23b 862 fOutputList->Add(fdEdxAft);
863
02cd82a0 864 fIncpT = new TH2F("fIncpT","HFE pid electro vs. centrality",200,0,100,100,0,50);
bfc7c23b 865 fOutputList->Add(fIncpT);
866
02cd82a0 867 fIncpTM20 = new TH2F("fIncpTM20","HFE pid electro vs. centrality with M20",200,0,100,100,0,50);
fb87d707 868 fOutputList->Add(fIncpTM20);
4e01b68c 869
02cd82a0 870 Int_t nBinspho[9] = { 200, 100, 500, 12, 50, 4, 200, 8, 100};
b12bc641 871 Double_t minpho[9] = { 0., 0., 0., -2.5, 0, -0.5, 0,-1.5, 0};
872 Double_t maxpho[9] = {100., 50., 0.5, 3.5, 1, 3.5, 2, 6.5, 50};
f09766dd 873
b12bc641 874 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 875 if(fqahist==1)fOutputList->Add(fInvmassLS);
bfc7c23b 876
b12bc641 877 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 878 if(fqahist==1)fOutputList->Add(fInvmassULS);
bfc7c23b 879
8806ce6c 880 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 881 if(fqahist==1)fOutputList->Add(fInvmassLSmc);
8806ce6c 882
883 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 884 if(fqahist==1)fOutputList->Add(fInvmassULSmc);
8806ce6c 885
93c189c5 886 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 );
887 fInvmassLSmc0->Sumw2();
888 fOutputList->Add(fInvmassLSmc0);
889
890 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 );
891 fInvmassLSmc1->Sumw2();
892 fOutputList->Add(fInvmassLSmc1);
893
894 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 );
895 fInvmassLSmc2->Sumw2();
896 fOutputList->Add(fInvmassLSmc2);
897
898 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 );
899 fInvmassLSmc3->Sumw2();
900 fOutputList->Add(fInvmassLSmc3);
901
902 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 );
903 fInvmassULSmc0->Sumw2();
904 fOutputList->Add(fInvmassULSmc0);
905
906 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 );
907 fInvmassULSmc1->Sumw2();
908 fOutputList->Add(fInvmassULSmc1);
909
910 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 );
911 fInvmassULSmc2->Sumw2();
912 fOutputList->Add(fInvmassULSmc2);
913
914 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 );
915 fInvmassULSmc3->Sumw2();
93c189c5 916 fOutputList->Add(fInvmassULSmc3);
917
bfc7c23b 918 fOpeningAngleLS = new TH1F("fOpeningAngleLS","Opening angle for LS pairs",100,0,1);
919 fOutputList->Add(fOpeningAngleLS);
920
921 fOpeningAngleULS = new TH1F("fOpeningAngleULS","Opening angle for ULS pairs",100,0,1);
922 fOutputList->Add(fOpeningAngleULS);
923
924 fPhotoElecPt = new TH1F("fPhotoElecPt", "photonic electron pt",100,0,50);
925 fOutputList->Add(fPhotoElecPt);
926
02cd82a0 927 fPhoElecPt = new TH2F("fPhoElecPt", "Pho-inclusive electron pt",200,0,100,100,0,50);
bfc7c23b 928 fOutputList->Add(fPhoElecPt);
929
02cd82a0 930 fPhoElecPtM20 = new TH2F("fPhoElecPtM20", "Pho-inclusive electron pt with M20",200,0,100,100,0,50);
fb87d707 931 fOutputList->Add(fPhoElecPtM20);
932
02cd82a0 933 fSameElecPt = new TH2F("fSameElecPt", "Same-inclusive electron pt",200,0,100,100,0,50);
f09766dd 934 fOutputList->Add(fSameElecPt);
935
02cd82a0 936 fSameElecPtM20 = new TH2F("fSameElecPtM20", "Same-inclusive electron pt with M20",200,0,100,100,0,50);
fb87d707 937 fOutputList->Add(fSameElecPtM20);
938
02cd82a0 939 fCent = new TH1F("fCent","Centrality",200,0,100) ;
bfc7c23b 940 fOutputList->Add(fCent);
941
942 // Make common binning
44be9e1d 943 const Double_t kMinP = 0.;
bfc7c23b 944 const Double_t kMaxP = 50.;
a6df418c 945 //const Double_t kTPCSigMim = 40.;
946 //const Double_t kTPCSigMax = 140.;
bfc7c23b 947
948 // 1st histogram: TPC dEdx with/without EMCAL (p, pT, TPC Signal, phi, eta, Sig, e/p, ,match, cell, M02, M20, Disp, Centrality, select)
a6df418c 949 Int_t nBins[16] = { 250, 10, 60, 20, 100, 300, 50, 40, 200, 200, 250, 200, 3, 5, 100, 8};
950 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 951 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};
952 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 953 if(fqahist==1)fOutputList->Add(fEleInfo);
bfc7c23b 954
fb87d707 955 //<--- Trigger info
feffe705 956 /*
fb87d707 957 fClsEBftTrigCut = new TH1F("fClsEBftTrigCut","cluster E before trigger selection",1000,0,100);
958 fOutputList->Add(fClsEBftTrigCut);
959
960 fClsEAftTrigCut = new TH1F("fClsEAftTrigCut","cluster E if cls has 0 trigcut channel",1000,0,100);
961 fOutputList->Add(fClsEAftTrigCut);
962
963 fClsEAftTrigCut1 = new TH1F("fClsEAftTrigCut1","cluster E if cls with trig channel",1000,0,100);
964 fOutputList->Add(fClsEAftTrigCut1);
965
966 fClsEAftTrigCut2 = new TH1F("fClsEAftTrigCut2","cluster E if cls with trigcut channel",1000,0,100);
967 fOutputList->Add(fClsEAftTrigCut2);
968
969 fClsEAftTrigCut3 = new TH1F("fClsEAftTrigCut3","cluster E if cls with trigcut channel + nCell>Ecorrect",1000,0,100);
970 fOutputList->Add(fClsEAftTrigCut3);
971
972 fClsEAftTrigCut4 = new TH1F("fClsEAftTrigCut4","cluster E if cls with trigcut channel + nCell>Ecorrect + cls time cut",1000,0,100);
973 fOutputList->Add(fClsEAftTrigCut4);
974
975 fClsETime = new TH2F("fClsETime", "Cls time vs E; E; time;",1000,0,100,1000,-0.0000002,0.0000009);
976 fOutputList->Add(fClsETime);
977
978 fClsETime1 = new TH2F("fClsETime1", "Cls time vs E if cls contains trigger channel; E; time;",1000,0,100,1000,-0.0000002,0.0000009);
979 fOutputList->Add(fClsETime1);
980
981 fTrigTimes = new TH1F("fTrigTimes", "Trigger time; time; N;",25,0,25);
982 fOutputList->Add(fTrigTimes);
983
42c75692 984 fCellCheck = new TH2F("fCellCheck", "Cell vs E; E GeV; Cell ID",10,6,26,12000,0,12000);
985 fOutputList->Add(fCellCheck);
feffe705 986 */
e7d87aef 987 //<---------- MC
988
02cd82a0 989 fInputHFEMC = new TH2F("fInputHFEMC","Input MC HFE pid electro vs. centrality",200,0,100,100,0,50);
e7d87aef 990 fOutputList->Add(fInputHFEMC);
991
02cd82a0 992 fInputAlle = new TH2F("fInputAlle","Input MC electro vs. centrality",200,0,100,100,0,50);
feffe705 993 fOutputList->Add(fInputAlle);
994
02cd82a0 995 fIncpTMChfe = new TH2F("fIncpTMChfe","MC HFE pid electro vs. centrality",200,0,100,100,0,50);
e7d87aef 996 fOutputList->Add(fIncpTMChfe);
997
02cd82a0 998 fIncpTMChfeAll = new TH2F("fIncpTMChfeAll","MC Alle pid electro vs. centrality",200,0,100,100,0,50);
3db00c72 999 fOutputList->Add(fIncpTMChfeAll);
1000
02cd82a0 1001 fIncpTMCM20hfe = new TH2F("fIncpTMCM20hfe","MC HFE pid electro vs. centrality with M20",200,0,100,100,0,50);
e7d87aef 1002 fOutputList->Add(fIncpTMCM20hfe);
1003
02cd82a0 1004 fIncpTMCM20hfeAll = new TH2F("fIncpTMCM20hfeAll","MC Alle pid electro vs. centrality with M20",200,0,100,100,0,50);
3db00c72 1005 fOutputList->Add(fIncpTMCM20hfeAll);
1006
60544aea 1007 fIncpTMCM20hfeCheck = new TH2F("fIncpTMCM20hfeCheck","MC HFE pid electro vs. centrality with M20 Check",200,0,100,100,0,50);
1008 fOutputList->Add(fIncpTMCM20hfeCheck);
44be9e1d 1009
987053ce 1010 Int_t nBinspho2[5] = { 200, 100, 7, 3, 100};
1011 Double_t minpho2[5] = { 0., 0., -2.5, -0.5, 0.};
1012 Double_t maxpho2[5] = {100., 50., 4.5, 2.5, 50.};
44be9e1d 1013
987053ce 1014 fIncpTMCpho = new THnSparseD("fIncpTMCpho","MC Pho pid electro vs. centrality",5,nBinspho2,minpho2,maxpho2);
e7d87aef 1015 fOutputList->Add(fIncpTMCpho);
1016
987053ce 1017 fIncpTMCM20pho = new THnSparseD("fIncpTMCM20pho","MC Pho pid electro vs. centrality with M20",5,nBinspho2,minpho2,maxpho2);
e7d87aef 1018 fOutputList->Add(fIncpTMCM20pho);
1019
987053ce 1020 fPhoElecPtMC = new THnSparseD("fPhoElecPtMC", "MC Pho-inclusive electron pt",5,nBinspho2,minpho2,maxpho2);
e7d87aef 1021 fOutputList->Add(fPhoElecPtMC);
1022
987053ce 1023 fPhoElecPtMCM20 = new THnSparseD("fPhoElecPtMCM20", "MC Pho-inclusive electron pt with M20",5,nBinspho2,minpho2,maxpho2);
feffe705 1024 fOutputList->Add(fPhoElecPtMCM20);
e7d87aef 1025
987053ce 1026 fSameElecPtMC = new THnSparseD("fSameElecPtMC", "MC Same-inclusive electron pt",5,nBinspho2,minpho2,maxpho2);
e7d87aef 1027 fOutputList->Add(fSameElecPtMC);
1028
987053ce 1029 fSameElecPtMCM20 = new THnSparseD("fSameElecPtMCM20", "MC Same-inclusive electron pt with M20",5,nBinspho2,minpho2,maxpho2);
e7d87aef 1030 fOutputList->Add(fSameElecPtMCM20);
1031
987053ce 1032 fIncpTMCM20pho_pi0e = new THnSparseD("fIncpTMCM20pho_pi0e","MC Pho pi0->e pid electro vs. centrality with M20",5,nBinspho2,minpho2,maxpho2);
1033 fIncpTMCM20pho_pi0e->Sumw2();
e4b0faf2 1034 fOutputList->Add(fIncpTMCM20pho_pi0e);
1035
987053ce 1036 fPhoElecPtMCM20_pi0e = new THnSparseD("fPhoElecPtMCM20_pi0e", "MC Pho-inclusive electron pt with M20 pi0->e",5,nBinspho2,minpho2,maxpho2);
1037 fPhoElecPtMCM20_pi0e->Sumw2();
e4b0faf2 1038 fOutputList->Add(fPhoElecPtMCM20_pi0e);
1039
987053ce 1040 fSameElecPtMCM20_pi0e = new THnSparseD("fSameElecPtMCM20_pi0e", "MC Same-inclusive electron pt pi0->e",5,nBinspho2,minpho2,maxpho2);
1041 fSameElecPtMCM20_pi0e->Sumw2();
e4b0faf2 1042 fOutputList->Add(fSameElecPtMCM20_pi0e);
93c189c5 1043 //
1044 fIncpTMCM20pho_eta = new THnSparseD("fIncpTMCM20pho_eta","MC Pho pi0->e pid electro vs. centrality with M20",5,nBinspho2,minpho2,maxpho2);
1045 fIncpTMCM20pho_eta->Sumw2();
1046 fOutputList->Add(fIncpTMCM20pho_eta);
1047
1048 fPhoElecPtMCM20_eta = new THnSparseD("fPhoElecPtMCM20_eta", "MC Pho-inclusive electron pt with M20 pi0->e",5,nBinspho2,minpho2,maxpho2);
1049 fPhoElecPtMCM20_eta->Sumw2();
1050 fOutputList->Add(fPhoElecPtMCM20_eta);
1051
1052 fSameElecPtMCM20_eta = new THnSparseD("fSameElecPtMCM20_eta", "MC Same-inclusive electron pt pi0->e",5,nBinspho2,minpho2,maxpho2);
1053 fSameElecPtMCM20_eta->Sumw2();
1054 fOutputList->Add(fSameElecPtMCM20_eta);
1055
e4b0faf2 1056
a6df418c 1057 CheckNclust = new TH1D("CheckNclust","cluster check",200,0,200);
1058 fOutputList->Add(CheckNclust);
1059
1060 CheckNits = new TH1D("CheckNits","ITS cluster check",8,-0.5,7.5);
1061 fOutputList->Add(CheckNits);
e7d87aef 1062
d113d7cd 1063 Hpi0pTcheck = new TH2D("Hpi0pTcheck","Pi0 pT from Hijing",100,0,50,3,-0.5,2.5);
1064 fOutputList->Add(Hpi0pTcheck);
1065
19325033 1066 HETApTcheck = new TH2D("HETApTcheck","Eta pT from Hijing",100,0,50,3,-0.5,2.5);
1067 fOutputList->Add(HETApTcheck);
1068
e4b0faf2 1069 HphopTcheck = new TH2D("HphopTcheck","Pho pT from Hijing",100,0,50,3,-0.5,2.5);
1070 fOutputList->Add(HphopTcheck);
1071
93c189c5 1072 fpTCheck = new TH1D("fpTCheck","pT check",500,0,50);
1073 fOutputList->Add(fpTCheck);
1074
50919258 1075 fMomDtoE = new TH2D("fMomDtoE","D->E pT correlations;e p_{T} GeV/c;D p_{T} GeV/c",400,0,40,400,0,40);
1076 fOutputList->Add(fMomDtoE);
d113d7cd 1077
bfc7c23b 1078 PostData(1,fOutputList);
1079}
1080
1081//________________________________________________________________________
1082void AliAnalysisTaskHFECal::Terminate(Option_t *)
1083{
1084 // Info("Terminate");
1085 AliAnalysisTaskSE::Terminate();
1086}
1087
1088//________________________________________________________________________
1089Bool_t AliAnalysisTaskHFECal::ProcessCutStep(Int_t cutStep, AliVParticle *track)
1090{
1091 // Check single track cuts for a given cut step
1092 const Int_t kMCOffset = AliHFEcuts::kNcutStepsMCTrack;
1093 if(!fCFM->CheckParticleCuts(cutStep + kMCOffset, track)) return kFALSE;
1094 return kTRUE;
1095}
1096//_________________________________________
f09766dd 1097//void AliAnalysisTaskHFECal::SelectPhotonicElectron(Int_t itrack, Double_t cent, AliESDtrack *track, Bool_t &fFlagPhotonicElec)
c2801a73 1098//void AliAnalysisTaskHFECal::SelectPhotonicElectron(Int_t itrack, Double_t cent, AliESDtrack *track, Bool_t &fFlagPhotonicElec, Bool_t &fFlagConvinatElec, Double_t nSig)
46305ed6 1099void 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)
bfc7c23b 1100{
1101 //Identify non-heavy flavour electrons using Invariant mass method
1102
fb87d707 1103 fTrackCuts->SetAcceptKinkDaughters(kFALSE);
feffe705 1104 fTrackCuts->SetRequireTPCRefit(kTRUE);
85985bb0 1105 fTrackCuts->SetRequireITSRefit(kTRUE);
feffe705 1106 fTrackCuts->SetEtaRange(-0.9,0.9);
f09766dd 1107 //fTrackCuts->SetRequireSigmaToVertex(kTRUE);
bfc7c23b 1108 fTrackCuts->SetMaxChi2PerClusterTPC(3.5);
feffe705 1109 fTrackCuts->SetMinNClustersTPC(90);
bfc7c23b 1110
1111 const AliESDVertex *pVtx = fESD->GetPrimaryVertex();
1112
93c189c5 1113 double ptEle = track->Pt(); //add
1114 if(ibgevent==0 && w > 0.0)
1115 {
1116 fpTCheck->Fill(ptEle,w);
1117 }
1118
bfc7c23b 1119 Bool_t flagPhotonicElec = kFALSE;
f09766dd 1120 Bool_t flagConvinatElec = kFALSE;
bfc7c23b 1121
8806ce6c 1122 int p1 = 0;
1123 if(mce==3)
1124 {
1125 Int_t label = TMath::Abs(track->GetLabel());
1126 TParticle* particle = stack->Particle(label);
1127 p1 = particle->GetFirstMother();
1128 }
1129
5e65985c 1130 //for(Int_t jTracks = itrack+1; jTracks<fESD->GetNumberOfTracks(); jTracks++){
1131 for(Int_t jTracks = 0; jTracks<fESD->GetNumberOfTracks(); jTracks++){
bfc7c23b 1132 AliESDtrack* trackAsso = fESD->GetTrack(jTracks);
1133 if (!trackAsso) {
1134 printf("ERROR: Could not receive track %d\n", jTracks);
1135 continue;
1136 }
5e65985c 1137 if(itrack==jTracks)continue;
d86a815c 1138 int jbgevent = 0;
5e65985c 1139
8806ce6c 1140 int p2 = 0;
1141 if(mce==3)
1142 {
1143 Int_t label2 = TMath::Abs(trackAsso->GetLabel());
1144 TParticle* particle2 = stack->Particle(label2);
d86a815c 1145 Bool_t MChijing_ass = fMC->IsFromBGEvent(label2);
1146 if(MChijing_ass)jbgevent =1;
8806ce6c 1147 if(particle2->GetFirstMother()>-1)
1148 p2 = particle2->GetFirstMother();
1149 }
1150
f09766dd 1151 Double_t dEdxAsso = -999., ptPrim=-999., ptAsso=-999., openingAngle = -999.;
bfc7c23b 1152 Double_t mass=999., width = -999;
1153 Bool_t fFlagLS=kFALSE, fFlagULS=kFALSE;
1154
93c189c5 1155 //ptPrim = track->Pt();
1156 ptPrim = ptEle;
f09766dd 1157
bfc7c23b 1158 dEdxAsso = trackAsso->GetTPCsignal();
1159 ptAsso = trackAsso->Pt();
1160 Int_t chargeAsso = trackAsso->Charge();
1161 Int_t charge = track->Charge();
1162
b12bc641 1163
3db00c72 1164 if(ptAsso <0.5) continue;
bfc7c23b 1165 if(!fTrackCuts->AcceptTrack(trackAsso)) continue;
f09766dd 1166 if(dEdxAsso <65 || dEdxAsso>100) continue; //11a pass1
bfc7c23b 1167
1168 Int_t fPDGe1 = 11; Int_t fPDGe2 = 11;
1169 if(charge>0) fPDGe1 = -11;
1170 if(chargeAsso>0) fPDGe2 = -11;
b12bc641 1171
1172 //printf("chargeAsso = %d\n",chargeAsso);
1173 //printf("charge = %d\n",charge);
bfc7c23b 1174 if(charge == chargeAsso) fFlagLS = kTRUE;
1175 if(charge != chargeAsso) fFlagULS = kTRUE;
1176
b12bc641 1177 //printf("fFlagLS = %d\n",fFlagLS);
1178 //printf("fFlagULS = %d\n",fFlagULS);
1179 //printf("\n");
1180
bfc7c23b 1181 AliKFParticle ge1(*track, fPDGe1);
1182 AliKFParticle ge2(*trackAsso, fPDGe2);
1183 AliKFParticle recg(ge1, ge2);
1184
1185 if(recg.GetNDF()<1) continue;
1186 Double_t chi2recg = recg.GetChi2()/recg.GetNDF();
1187 if(TMath::Sqrt(TMath::Abs(chi2recg))>3.) continue;
1188
d86a815c 1189 // for v5-03-70-AN comment
bfc7c23b 1190 AliKFVertex primV(*pVtx);
1191 primV += recg;
1192 recg.SetProductionVertex(primV);
1193
d113d7cd 1194 recg.SetMassConstraint(0,0.0001);
d86a815c 1195
bfc7c23b 1196 openingAngle = ge1.GetAngle(ge2);
1197 if(fFlagLS) fOpeningAngleLS->Fill(openingAngle);
1198 if(fFlagULS) fOpeningAngleULS->Fill(openingAngle);
1199
bfc7c23b 1200
1201 recg.GetMass(mass,width);
1202
c2801a73 1203 double ishower = 0;
1204 if(shower>0.0 && shower<0.3)ishower = 1;
1205
b12bc641 1206 double phoinfo[9];
f09766dd 1207 phoinfo[0] = cent;
1208 phoinfo[1] = ptPrim;
1209 phoinfo[2] = mass;
3db00c72 1210 phoinfo[3] = nSig;
4e01b68c 1211 phoinfo[4] = openingAngle;
c2801a73 1212 phoinfo[5] = ishower;
1213 phoinfo[6] = ep;
1214 phoinfo[7] = mce;
b12bc641 1215 phoinfo[8] = ptAsso;
f09766dd 1216
1217 if(fFlagLS) fInvmassLS->Fill(phoinfo);
1218 if(fFlagULS) fInvmassULS->Fill(phoinfo);
d86a815c 1219 if(fFlagLS && ibgevent==0 && jbgevent==0) fInvmassLSmc->Fill(phoinfo,w);
93c189c5 1220 if(fFlagULS && ibgevent==0 && jbgevent==0)
1221 {
1222 fInvmassULSmc->Fill(phoinfo,w);
1223 }
b12bc641 1224 //printf("fInvmassCut %f\n",fInvmassCut);
1225 //printf("openingAngle %f\n",fOpeningAngleCut);
1226
ffa14dda 1227 if(openingAngle > fOpeningAngleCut) continue;
bfc7c23b 1228
d86a815c 1229 // for real data
93c189c5 1230 //printf("mce =%f\n",mce);
d86a815c 1231 if(mce<-0.5) // mce==-1. is real
1232 {
93c189c5 1233 //printf("Real data\n");
d86a815c 1234 if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec){
1235 //if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec && (p1==p2)){ <--- only MC train (55,56) v5-03-68-AN & 69 for check
1236 flagPhotonicElec = kTRUE;
1237 }
1238 if(mass<fInvmassCut && fFlagLS && !flagConvinatElec){
1239 flagConvinatElec = kTRUE;
1240 }
1241 }
1242 // for MC data
1243 else
1244 {
93c189c5 1245 //printf("MC data\n");
1246
1247 if(w>0.0)
1248 {
1249 if(fFlagLS && ibgevent==0 && jbgevent==0) fInvmassLSmc0->Fill(ptPrim,mass,w);
1250 if(fFlagULS && ibgevent==0 && jbgevent==0) fInvmassULSmc0->Fill(ptPrim,mass,w);
1251 if(fFlagLS && ibgevent==0 && jbgevent==0) fInvmassLSmc1->Fill(ptPrim,mass);
1252 if(fFlagULS && ibgevent==0 && jbgevent==0) fInvmassULSmc1->Fill(ptPrim,mass);
1253 if(fFlagLS && ibgevent==0 && jbgevent==0 && (p1==p2)) fInvmassLSmc2->Fill(ptPrim,mass,w);
1254 if(fFlagULS && ibgevent==0 && jbgevent==0 && (p1==p2)) fInvmassULSmc2->Fill(ptPrim,mass,w);
1255 if(fFlagLS && ibgevent==0 && jbgevent==0 && (p1==p2)) fInvmassLSmc3->Fill(ptPrim,mass);
1256 if(fFlagULS && ibgevent==0 && jbgevent==0 && (p1==p2)) fInvmassULSmc3->Fill(ptPrim,mass);
1257 }
d86a815c 1258 if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec && (ibgevent==jbgevent)){
1259 //if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec && (p1==p2)){ <--- only MC train (55,56) v5-03-68-AN & 69 for check
1260 flagPhotonicElec = kTRUE;
1261 }
1262 if(mass<fInvmassCut && fFlagLS && !flagConvinatElec && (ibgevent==jbgevent)){
1263 flagConvinatElec = kTRUE;
1264 }
1265 }
1266
bfc7c23b 1267 }
1268 fFlagPhotonicElec = flagPhotonicElec;
f09766dd 1269 fFlagConvinatElec = flagConvinatElec;
bfc7c23b 1270
1271}
8806ce6c 1272//-------------------------------------------
46305ed6 1273
1274void AliAnalysisTaskHFECal::FindMother(TParticle* part, int &label, int &pid)
1275{
1276 //int label = 99999;
1277 //int pid = 99999;
1278
1279 if(part->GetFirstMother()>-1)
1280 {
1281 label = part->GetFirstMother();
1282 pid = stack->Particle(label)->GetPdgCode();
1283 }
1284 //cout << "Find Mother : label = " << label << " ; pid" << pid << endl;
1285}
1286
8806ce6c 1287double AliAnalysisTaskHFECal::GetMCweight(double mcPi0pT)
1288{
1289 double weight = 1.0;
1290
1291 if(mcPi0pT>0.0 && mcPi0pT<5.0)
1292 {
1293 weight = 0.323*mcPi0pT/(TMath::Exp(-1.6+0.767*mcPi0pT+0.0285*mcPi0pT*mcPi0pT));
1294 }
1295 else
1296 {
1297 weight = 115.0/(0.718*mcPi0pT*TMath::Power(mcPi0pT,3.65));
1298 }
1299 return weight;
1300}
fb87d707 1301
93c189c5 1302double AliAnalysisTaskHFECal::GetMCweightEta(double mcEtapT)
1303{
1304 double weight = 1.0;
1305
1306 weight = 223.3/TMath::Power((TMath::Exp(-0.17*mcEtapT-0.0322*mcEtapT*mcEtapT)+mcEtapT/1.69),5.65);
1307 return weight;
1308}
1309
1310
fb87d707 1311//_________________________________________
1312void AliAnalysisTaskHFECal::FindTriggerClusters()
1313{
1314 // constants
1315 const int nModuleCols = 2;
1316 const int nModuleRows = 5;
1317 const int nColsFeeModule = 48;
1318 const int nRowsFeeModule = 24;
1319 const int nColsFaltroModule = 24;
1320 const int nRowsFaltroModule = 12;
1321 //const int faltroWidthMax = 20;
1322
1323 // part 1, trigger extraction -------------------------------------
1324 Int_t globCol, globRow;
1325 //Int_t ntimes=0, nTrigChannel=0, nTrigChannelCut=0, trigInCut=0;
1326 Int_t ntimes=0, nTrigChannel=0, nTrigChannelCut=0;
1327
1328 //Int_t trigtimes[faltroWidthMax];
1329 Double_t cellTime[nColsFeeModule*nModuleCols][nRowsFeeModule*nModuleRows];
1330 Double_t cellEnergy[nColsFeeModule*nModuleCols][nRowsFeeModule*nModuleRows];
1331 //Double_t fTrigCutLow = 6;
1332 //Double_t fTrigCutHigh = 10;
1333 Double_t fTimeCutLow = 469e-09;
1334 Double_t fTimeCutHigh = 715e-09;
1335
1336 AliESDCaloTrigger * fCaloTrigger = fESD->GetCaloTrigger( "EMCAL" );
1337
1338 // erase trigger maps
1339 for(Int_t i = 0; i < nColsFaltroModule*nModuleCols; i++ )
1340 {
1341 for(Int_t j = 0; j < nRowsFaltroModule*nModuleRows; j++ )
1342 {
1343 ftriggersCut[i][j] = 0;
1344 ftriggers[i][j] = 0;
1345 ftriggersTime[i][j] = 0;
1346 }
1347 }
1348
1349 Int_t iglobCol=0, iglobRow=0;
1350 // go through triggers
1351 if( fCaloTrigger->GetEntries() > 0 )
1352 {
1353 // needs reset
1354 fCaloTrigger->Reset();
1355 while( fCaloTrigger->Next() )
1356 {
1357 fCaloTrigger->GetPosition( globCol, globRow );
1358 fCaloTrigger->GetNL0Times( ntimes );
1359 /*
1360 // no L0s
1361 if( ntimes < 1 ) continue;
1362 // get precise timings
1363 fCaloTrigger->GetL0Times( trigtimes );
1364 trigInCut = 0;
1365 for(Int_t i = 0; i < ntimes; i++ )
1366 {
1367 // save the first trigger time in channel
1368 if( i == 0 || triggersTime[globCol][globRow] > trigtimes[i] )
1369 triggersTime[globCol][globRow] = trigtimes[i];
1370 //printf("trigger times: %d\n",trigtimes[i]);
1371 // check if in cut
1372 if(trigtimes[i] > fTrigCutLow && trigtimes[i] < fTrigCutHigh )
1373 trigInCut = 1;
1374
1375 fTrigTimes->Fill(trigtimes[i]);
1376 }
1377 */
1378
1379 //L1 analysis from AliAnalysisTaskEMCALTriggerQA
1380 Int_t bit = 0;
1381 fCaloTrigger->GetTriggerBits(bit);
1382
1383 Int_t ts = 0;
1384 fCaloTrigger->GetL1TimeSum(ts);
1385 if (ts > 0)ftriggers[globCol][globRow] = 1;
1386 // number of triggered channels in event
1387 nTrigChannel++;
1388 // ... inside cut
1389 if(ts>0 && (bit >> 6 & 0x1))
1390 {
1391 iglobCol = globCol;
1392 iglobRow = globRow;
1393 nTrigChannelCut++;
1394 ftriggersCut[globCol][globRow] = 1;
1395 }
1396
1397 } // calo trigger entries
1398 } // has calo trigger entries
1399
1400 // part 2 go through the clusters here -----------------------------------
1401 Int_t nCluster=0, nCell=0, iCell=0, gCell=0;
1402 Short_t cellAddr, nSACell, mclabel;
1403 //Int_t nSACell, iSACell, mclabel;
1404 Int_t iSACell;
42c75692 1405 Double_t cellAmp=0, cellTimeT=0, clusterTime=0, efrac=0;
fb87d707 1406 Int_t nSupMod, nModule, nIphi, nIeta, iphi, ieta, gphi, geta, feta, fphi;
1407
1408 TRefArray *fCaloClusters = new TRefArray();
1409 fESD->GetEMCALClusters( fCaloClusters );
1410 nCluster = fCaloClusters->GetEntries();
1411
1412
1413 // save all cells times if there are clusters
1414 if( nCluster > 0 ){
1415 // erase time map
1416 for(Int_t i = 0; i < nColsFeeModule*nModuleCols; i++ ){
1417 for(Int_t j = 0; j < nRowsFeeModule*nModuleRows; j++ ){
1418 cellTime[i][j] = 0.;
1419 cellEnergy[i][j] = 0.;
1420 }
1421 }
1422
1423 // get cells
1424 AliESDCaloCells *fCaloCells = fESD->GetEMCALCells();
1425 //AliVCaloCells fCaloCells = fESD->GetEMCALCells();
1426 nSACell = fCaloCells->GetNumberOfCells();
1427 for(iSACell = 0; iSACell < nSACell; iSACell++ ){
1428 // get the cell info *fCal
1429 fCaloCells->GetCell( iSACell, cellAddr, cellAmp, cellTimeT , mclabel, efrac);
1430
1431 // get cell position
1432 fGeom->GetCellIndex( cellAddr, nSupMod, nModule, nIphi, nIeta );
1433 fGeom->GetCellPhiEtaIndexInSModule( nSupMod,nModule, nIphi, nIeta, iphi, ieta);
1434
1435 // convert co global phi eta
1436 gphi = iphi + nRowsFeeModule*(nSupMod/2);
1437 geta = ieta + nColsFeeModule*(nSupMod%2);
1438
1439 // save cell time and energy
1440 cellTime[geta][gphi] = cellTimeT;
1441 cellEnergy[geta][gphi] = cellAmp;
1442
1443 }
1444 }
1445
1446 Int_t nClusterTrig, nClusterTrigCut;
1447 UShort_t *cellAddrs;
1448 Double_t clsE=-999, clsEta=-999, clsPhi=-999;
1449 Float_t clsPos[3] = {0.,0.,0.};
1450
1451 for(Int_t icl=0; icl<fESD->GetNumberOfCaloClusters(); icl++)
1452 {
1453 AliESDCaloCluster *cluster = fESD->GetCaloCluster(icl);
1454 if(!cluster || !cluster->IsEMCAL()) continue;
1455
1456 // get cluster cells
1457 nCell = cluster->GetNCells();
1458
1459 // get cluster energy
1460 clsE = cluster->E();
1461
1462 // get cluster position
1463 cluster->GetPosition(clsPos);
1464 TVector3 clsPosVec(clsPos[0],clsPos[1],clsPos[2]);
1465 clsEta = clsPosVec.Eta();
1466 clsPhi = clsPosVec.Phi();
1467
1468 // get the cell addresses
1469 cellAddrs = cluster->GetCellsAbsId();
1470
1471 // check if the cluster contains cell, that was marked as triggered
1472 nClusterTrig = 0;
1473 nClusterTrigCut = 0;
1474
1475 // loop the cells to check, if cluser in acceptance
1476 // any cluster with a cell outside acceptance is not considered
1477 for( iCell = 0; iCell < nCell; iCell++ )
1478 {
42c75692 1479 // check hot cell
feffe705 1480 //if(clsE>6.0)fCellCheck->Fill(clsE,cellAddrs[iCell]);
42c75692 1481
fb87d707 1482 // get cell position
1483 fGeom->GetCellIndex( cellAddrs[iCell], nSupMod, nModule, nIphi, nIeta );
1484 fGeom->GetCellPhiEtaIndexInSModule( nSupMod,nModule, nIphi, nIeta, iphi, ieta);
1485
1486 // convert co global phi eta
1487 gphi = iphi + nRowsFeeModule*(nSupMod/2);
1488 geta = ieta + nColsFeeModule*(nSupMod%2);
1489
1490 if( cellTime[geta][gphi] > 0. ){
1491 clusterTime += cellTime[geta][gphi];
1492 gCell++;
1493 }
1494
1495 // get corresponding FALTRO
1496 fphi = gphi / 2;
1497 feta = geta / 2;
1498
1499 // try to match with a triggered
1500 if( ftriggers[feta][fphi]==1)
1501 { nClusterTrig++;
1502 }
1503 if( ftriggersCut[feta][fphi]==1)
1504 { nClusterTrigCut++;
1505 }
1506
1507 } // cells
1508
1509
1510 if( gCell > 0 )
1511 clusterTime = clusterTime / (Double_t)gCell;
1512 // fix the reconstriction code time 100ns jumps
1513 if( fESD->GetBunchCrossNumber() % 4 < 2 )
1514 clusterTime -= 0.0000001;
1515
feffe705 1516 //fClsETime->Fill(clsE,clusterTime);
1517 //fClsEBftTrigCut->Fill(clsE);
fb87d707 1518
1519 if(nClusterTrig>0){
feffe705 1520 //fClsETime1->Fill(clsE,clusterTime);
fb87d707 1521 }
1522
1523 if(nClusterTrig>0){
1524 cluster->SetChi2(1);
feffe705 1525 //fClsEAftTrigCut1->Fill(clsE);
fb87d707 1526 }
1527
1528 if(nClusterTrigCut>0){
1529 cluster->SetChi2(2);
feffe705 1530 //fClsEAftTrigCut2->Fill(clsE);
fb87d707 1531 }
1532
1533 if(nClusterTrigCut>0 && ( nCell > (1 + clsE / 3)))
1534 {
1535 cluster->SetChi2(3);
feffe705 1536 //fClsEAftTrigCut3->Fill(clsE);
fb87d707 1537 }
1538
1539 if(nClusterTrigCut>0 && (nCell > (1 + clsE / 3) )&&( clusterTime > fTimeCutLow && clusterTime < fTimeCutHigh ))
1540 {
1541 // cluster->SetChi2(4);
feffe705 1542 //fClsEAftTrigCut4->Fill(clsE);
fb87d707 1543 }
1544 if(nClusterTrigCut<1)
1545 {
1546 cluster->SetChi2(0);
1547
feffe705 1548 //fClsEAftTrigCut->Fill(clsE);
fb87d707 1549 }
1550
1551 } // clusters
1552}
1553
1554
1555
1556
1557
1558
1559