]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGHF/hfe/AliAnalysisTaskHFECal.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
#101318: Patch for various problems in AliROOT
[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"
8efacc22 23#include "TF1.h"
bfc7c23b 24#include "TMath.h"
25#include "TCanvas.h"
26#include "THnSparse.h"
27#include "TLorentzVector.h"
28#include "TString.h"
29#include "TFile.h"
09022877 30#include "TGraphErrors.h"
bfc7c23b 31
32#include "AliAnalysisTask.h"
33#include "AliAnalysisManager.h"
34
35#include "AliESDEvent.h"
36#include "AliESDHandler.h"
37#include "AliAODEvent.h"
38#include "AliAODHandler.h"
39
e7d87aef 40#include "AliMCEventHandler.h"
41#include "AliMCEvent.h"
42#include "AliMCParticle.h"
43
bfc7c23b 44#include "AliAnalysisTaskHFECal.h"
45#include "TGeoGlobalMagField.h"
46#include "AliLog.h"
47#include "AliAnalysisTaskSE.h"
48#include "TRefArray.h"
49#include "TVector.h"
50#include "AliESDInputHandler.h"
51#include "AliESDpid.h"
52#include "AliESDtrackCuts.h"
53#include "AliPhysicsSelection.h"
54#include "AliESDCaloCluster.h"
55#include "AliAODCaloCluster.h"
56#include "AliEMCALRecoUtils.h"
57#include "AliEMCALGeometry.h"
58#include "AliGeomManager.h"
59#include "stdio.h"
60#include "TGeoManager.h"
61#include "iostream"
62#include "fstream"
63
64#include "AliEMCALTrack.h"
65#include "AliMagF.h"
66
67#include "AliKFParticle.h"
68#include "AliKFVertex.h"
69
70#include "AliPID.h"
71#include "AliPIDResponse.h"
72#include "AliHFEcontainer.h"
73#include "AliHFEcuts.h"
74#include "AliHFEpid.h"
75#include "AliHFEpidBase.h"
76#include "AliHFEpidQAmanager.h"
77#include "AliHFEtools.h"
78#include "AliCFContainer.h"
79#include "AliCFManager.h"
80
e7d87aef 81#include "AliStack.h"
82
bfc7c23b 83#include "AliCentrality.h"
84
c2aad3ae 85using namespace std;
86
bfc7c23b 87ClassImp(AliAnalysisTaskHFECal)
88//________________________________________________________________________
89AliAnalysisTaskHFECal::AliAnalysisTaskHFECal(const char *name)
90 : AliAnalysisTaskSE(name)
91 ,fESD(0)
e7d87aef 92 ,fMC(0)
46305ed6 93 ,stack(0)
fb87d707 94 ,fGeom(0)
bfc7c23b 95 ,fOutputList(0)
5e0e45b3 96 ,fqahist(1)
bfc7c23b 97 ,fTrackCuts(0)
98 ,fCuts(0)
99 ,fIdentifiedAsOutInz(kFALSE)
100 ,fPassTheEventCut(kFALSE)
101 ,fRejectKinkMother(kFALSE)
e7d87aef 102 ,fmcData(kFALSE)
bfc7c23b 103 ,fVz(0.0)
104 ,fCFM(0)
105 ,fPID(0)
106 ,fPIDqa(0)
107 ,fOpeningAngleCut(0.1)
f9447b3b 108 ,fInvmassCut(0.05)
109 //,fInvmassCut(0.1) // no mass
bfc7c23b 110 ,fNoEvents(0)
111 ,fEMCAccE(0)
f4e0d2d5 112 ,hEMCAccE(0)
bfc7c23b 113 ,fTrkpt(0)
114 ,fTrkEovPBef(0)
115 ,fTrkEovPAft(0)
116 ,fdEdxBef(0)
117 ,fdEdxAft(0)
118 ,fIncpT(0)
fb87d707 119 ,fIncpTM20(0)
bfc7c23b 120 ,fInvmassLS(0)
121 ,fInvmassULS(0)
8806ce6c 122 ,fInvmassLSmc(0)
123 ,fInvmassULSmc(0)
93c189c5 124 ,fInvmassLSmc0(0)
125 ,fInvmassLSmc1(0)
126 ,fInvmassLSmc2(0)
127 ,fInvmassLSmc3(0)
128 ,fInvmassULSmc0(0)
129 ,fInvmassULSmc1(0)
130 ,fInvmassULSmc2(0)
131 ,fInvmassULSmc3(0)
bfc7c23b 132 ,fOpeningAngleLS(0)
133 ,fOpeningAngleULS(0)
134 ,fPhotoElecPt(0)
135 ,fPhoElecPt(0)
fb87d707 136 ,fPhoElecPtM20(0)
f09766dd 137 ,fSameElecPt(0)
fb87d707 138 ,fSameElecPtM20(0)
bfc7c23b 139 ,fTrackPtBefTrkCuts(0)
140 ,fTrackPtAftTrkCuts(0)
141 ,fTPCnsigma(0)
142 ,fCent(0)
143 ,fEleInfo(0)
e1f0fb74 144 ,fElenSigma(0)
feffe705 145 /*
fb87d707 146 ,fClsEBftTrigCut(0)
147 ,fClsEAftTrigCut(0)
148 ,fClsEAftTrigCut1(0)
149 ,fClsEAftTrigCut2(0)
150 ,fClsEAftTrigCut3(0)
151 ,fClsEAftTrigCut4(0)
152 ,fClsETime(0)
153 ,fClsETime1(0)
154 ,fTrigTimes(0)
42c75692 155 ,fCellCheck(0)
feffe705 156 */
e7d87aef 157 ,fInputHFEMC(0)
feffe705 158 ,fInputAlle(0)
e7d87aef 159 ,fIncpTMChfe(0)
3db00c72 160 ,fIncpTMChfeAll(0)
e7d87aef 161 ,fIncpTMCM20hfe(0)
3db00c72 162 ,fIncpTMCM20hfeAll(0)
60544aea 163 ,fIncpTMCM20hfeCheck(0)
bd6ee6dd 164 ,fInputHFEMC_weight(0)
165 ,fIncpTMCM20hfeCheck_weight(0)
e7d87aef 166 ,fIncpTMCpho(0)
167 ,fIncpTMCM20pho(0)
168 ,fPhoElecPtMC(0)
169 ,fPhoElecPtMCM20(0)
170 ,fSameElecPtMC(0)
171 ,fSameElecPtMCM20(0)
e4b0faf2 172 ,fIncpTMCM20pho_pi0e(0)
173 ,fPhoElecPtMCM20_pi0e(0)
174 ,fSameElecPtMCM20_pi0e(0)
93c189c5 175 ,fIncpTMCM20pho_eta(0)
176 ,fPhoElecPtMCM20_eta(0)
177 ,fSameElecPtMCM20_eta(0)
697ecf6b 178 ,fIncpTMCpho_pi0e_TPC(0)
179 ,fPhoElecPtMC_pi0e_TPC(0)
180 ,fSameElecPtMC_pi0e_TPC(0)
a6df418c 181 ,CheckNclust(0)
182 ,CheckNits(0)
d113d7cd 183 ,Hpi0pTcheck(0)
19325033 184 ,HETApTcheck(0)
bd6ee6dd 185 ,HphopTcheck(0)
5e0e45b3 186 ,HDpTcheck(0)
187 ,HBpTcheck(0)
93c189c5 188 ,fpTCheck(0)
50919258 189 ,fMomDtoE(0)
70448e3b 190 ,fLabelCheck(0)
191 ,fgeoFake(0)
d19af75d 192 ,fFakeTrk0(0)
193 ,fFakeTrk1(0)
59d998de 194 ,ftimingEle(0)
09022877 195 //,fnSigEtaCorr(NULL)
bfc7c23b 196{
197 //Named constructor
198
199 fPID = new AliHFEpid("hfePid");
200 fTrackCuts = new AliESDtrackCuts();
201
09022877 202 for(int i=0; i<7; i++)fnSigEtaCorr[i] = 0;
203
bfc7c23b 204 // Define input and output slots here
205 // Input slot #0 works with a TChain
206 DefineInput(0, TChain::Class());
207 // Output slot #0 id reserved by the base class for AOD
208 // Output slot #1 writes into a TH1 container
209 // DefineOutput(1, TH1I::Class());
210 DefineOutput(1, TList::Class());
211 // DefineOutput(3, TTree::Class());
212}
213
214//________________________________________________________________________
215AliAnalysisTaskHFECal::AliAnalysisTaskHFECal()
216 : AliAnalysisTaskSE("DefaultAnalysis_AliAnalysisTaskHFECal")
217 ,fESD(0)
e7d87aef 218 ,fMC(0)
46305ed6 219 ,stack(0)
fb87d707 220 ,fGeom(0)
bfc7c23b 221 ,fOutputList(0)
5e0e45b3 222 ,fqahist(1)
bfc7c23b 223 ,fTrackCuts(0)
224 ,fCuts(0)
225 ,fIdentifiedAsOutInz(kFALSE)
226 ,fPassTheEventCut(kFALSE)
227 ,fRejectKinkMother(kFALSE)
e7d87aef 228 ,fmcData(kFALSE)
bfc7c23b 229 ,fVz(0.0)
230 ,fCFM(0)
231 ,fPID(0)
232 ,fPIDqa(0)
233 ,fOpeningAngleCut(0.1)
f9447b3b 234 ,fInvmassCut(0.05)
235 //,fInvmassCut(0.1)
bfc7c23b 236 ,fNoEvents(0)
237 ,fEMCAccE(0)
f4e0d2d5 238 ,hEMCAccE(0)
bfc7c23b 239 ,fTrkpt(0)
240 ,fTrkEovPBef(0)
241 ,fTrkEovPAft(0)
242 ,fdEdxBef(0)
243 ,fdEdxAft(0)
244 ,fIncpT(0)
fb87d707 245 ,fIncpTM20(0)
bfc7c23b 246 ,fInvmassLS(0)
247 ,fInvmassULS(0)
8806ce6c 248 ,fInvmassLSmc(0)
249 ,fInvmassULSmc(0)
93c189c5 250 ,fInvmassLSmc0(0)
251 ,fInvmassLSmc1(0)
252 ,fInvmassLSmc2(0)
253 ,fInvmassLSmc3(0)
254 ,fInvmassULSmc0(0)
255 ,fInvmassULSmc1(0)
256 ,fInvmassULSmc2(0)
257 ,fInvmassULSmc3(0)
bfc7c23b 258 ,fOpeningAngleLS(0)
259 ,fOpeningAngleULS(0)
260 ,fPhotoElecPt(0)
261 ,fPhoElecPt(0)
fb87d707 262 ,fPhoElecPtM20(0)
f09766dd 263 ,fSameElecPt(0)
fb87d707 264 ,fSameElecPtM20(0)
bfc7c23b 265 ,fTrackPtBefTrkCuts(0)
266 ,fTrackPtAftTrkCuts(0)
267 ,fTPCnsigma(0)
268 ,fCent(0)
269 ,fEleInfo(0)
e1f0fb74 270 ,fElenSigma(0)
feffe705 271 /*
fb87d707 272 ,fClsEBftTrigCut(0)
273 ,fClsEAftTrigCut(0)
274 ,fClsEAftTrigCut1(0)
275 ,fClsEAftTrigCut2(0)
276 ,fClsEAftTrigCut3(0)
277 ,fClsEAftTrigCut4(0)
278 ,fClsETime(0)
279 ,fClsETime1(0)
280 ,fTrigTimes(0)
42c75692 281 ,fCellCheck(0)
feffe705 282 */
e7d87aef 283 ,fInputHFEMC(0)
feffe705 284 ,fInputAlle(0)
e7d87aef 285 ,fIncpTMChfe(0)
3db00c72 286 ,fIncpTMChfeAll(0)
e7d87aef 287 ,fIncpTMCM20hfe(0)
3db00c72 288 ,fIncpTMCM20hfeAll(0)
60544aea 289 ,fIncpTMCM20hfeCheck(0)
bd6ee6dd 290 ,fInputHFEMC_weight(0)
291 ,fIncpTMCM20hfeCheck_weight(0)
e7d87aef 292 ,fIncpTMCpho(0)
293 ,fIncpTMCM20pho(0)
294 ,fPhoElecPtMC(0)
295 ,fPhoElecPtMCM20(0)
296 ,fSameElecPtMC(0)
297 ,fSameElecPtMCM20(0)
e4b0faf2 298 ,fIncpTMCM20pho_pi0e(0)
299 ,fPhoElecPtMCM20_pi0e(0)
300 ,fSameElecPtMCM20_pi0e(0)
93c189c5 301 ,fIncpTMCM20pho_eta(0)
302 ,fPhoElecPtMCM20_eta(0)
303 ,fSameElecPtMCM20_eta(0)
697ecf6b 304 ,fIncpTMCpho_pi0e_TPC(0)
305 ,fPhoElecPtMC_pi0e_TPC(0)
306 ,fSameElecPtMC_pi0e_TPC(0)
a6df418c 307 ,CheckNclust(0)
308 ,CheckNits(0)
d113d7cd 309 ,Hpi0pTcheck(0)
19325033 310 ,HETApTcheck(0)
e4b0faf2 311 ,HphopTcheck(0)
5e0e45b3 312 ,HDpTcheck(0)
313 ,HBpTcheck(0)
93c189c5 314 ,fpTCheck(0)
50919258 315 ,fMomDtoE(0)
70448e3b 316 ,fLabelCheck(0)
317 ,fgeoFake(0)
d19af75d 318 ,fFakeTrk0(0)
319 ,fFakeTrk1(0)
59d998de 320 ,ftimingEle(0)
09022877 321 //,fnSigEtaCorr(NULL)
bfc7c23b 322{
323 //Default constructor
324 fPID = new AliHFEpid("hfePid");
325
326 fTrackCuts = new AliESDtrackCuts();
327
09022877 328 for(int i=0; i<7; i++)fnSigEtaCorr[i] = 0;
329
bfc7c23b 330 // Constructor
331 // Define input and output slots here
332 // Input slot #0 works with a TChain
333 DefineInput(0, TChain::Class());
334 // Output slot #0 id reserved by the base class for AOD
335 // Output slot #1 writes into a TH1 container
336 // DefineOutput(1, TH1I::Class());
337 DefineOutput(1, TList::Class());
338 //DefineOutput(3, TTree::Class());
339}
340//_________________________________________
341
342AliAnalysisTaskHFECal::~AliAnalysisTaskHFECal()
343{
344 //Destructor
345
346 delete fOutputList;
fb87d707 347 delete fGeom;
bfc7c23b 348 delete fPID;
a6df418c 349 delete fCuts;
bfc7c23b 350 delete fCFM;
351 delete fPIDqa;
352 delete fTrackCuts;
353}
354//_________________________________________
355
356void AliAnalysisTaskHFECal::UserExec(Option_t*)
357{
358 //Main loop
359 //Called for each event
360
361 // create pointer to event
362 fESD = dynamic_cast<AliESDEvent*>(InputEvent());
363 if (!fESD) {
364 printf("ERROR: fESD not available\n");
365 return;
366 }
367
368 if(!fCuts){
369 AliError("HFE cuts not available");
370 return;
371 }
372
373 if(!fPID->IsInitialized()){
374 // Initialize PID with the given run number
375 AliWarning("PID not initialised, get from Run no");
376 fPID->InitializePID(fESD->GetRunNumber());
377 }
378
e7d87aef 379 if(fmcData)fMC = MCEvent();
46305ed6 380 //AliStack* stack = NULL;
e7d87aef 381 if(fmcData && fMC)stack = fMC->Stack();
382
383 Float_t cent = -1.;
384 AliCentrality *centrality = fESD->GetCentrality();
385 cent = centrality->GetCentralityPercentile("V0M");
386
387 //---- fill MC track info
388 if(fmcData && fMC)
389 {
390 Int_t nParticles = stack->GetNtrack();
391 for (Int_t iParticle = 0; iParticle < nParticles; iParticle++) {
392 TParticle* particle = stack->Particle(iParticle);
393 int fPDG = particle->GetPdgCode();
394 double mcZvertex = fMC->GetPrimaryVertex()->GetZ();
395 double pTMC = particle->Pt();
feffe705 396 double proR = particle->R();
9b3495ff 397 double etaMC = particle->Eta();
398 if(fabs(etaMC)>0.6)continue;
e7d87aef 399 Bool_t mcInDtoE= kFALSE;
400 Bool_t mcInBtoE= kFALSE;
401
d113d7cd 402 Bool_t MChijing = fMC->IsFromBGEvent(iParticle);
fd426741 403 //if(!MChijing)printf("not MC hijing");
d113d7cd 404 int iHijing = 1;
405 if(!MChijing)iHijing = 0;
68d718e7 406 double mcphoinfo[5];
407 mcphoinfo[0] = cent;
408 mcphoinfo[1] = pTMC;
409 mcphoinfo[2] = iHijing;
410 //if(fPDG==111)Hpi0pTcheck->Fill(pTMC,iHijing);
411 //if(fPDG==221)HETApTcheck->Fill(pTMC,iHijing);
412 if(fPDG==111)Hpi0pTcheck->Fill(mcphoinfo);
413 if(fPDG==221)HETApTcheck->Fill(mcphoinfo);
bd6ee6dd 414 if(fabs(fPDG)==411 || fabs(fPDG)==413 || fabs(fPDG)==421 || fabs(fPDG)==423 || fabs(fPDG)==431)HDpTcheck->Fill(pTMC,iHijing);
415 if(fabs(fPDG)==511 || fabs(fPDG)==513 || fabs(fPDG)==521 || fabs(fPDG)==523 || fabs(fPDG)==531)HBpTcheck->Fill(pTMC,iHijing);
d113d7cd 416
e4b0faf2 417 if(particle->GetFirstMother()>-1 && fPDG==22)
418 {
419 int parentPID = stack->Particle(particle->GetFirstMother())->GetPdgCode();
420 if(parentPID==111 || parentPID==221)HphopTcheck->Fill(pTMC,iHijing); // pi0->g & eta->g
421 }
422
e7d87aef 423 if(particle->GetFirstMother()>-1 && fabs(fPDG)==11)
424 {
425 int parentPID = stack->Particle(particle->GetFirstMother())->GetPdgCode();
bd6ee6dd 426 double pTMCparent = stack->Particle(particle->GetFirstMother())->Pt();
e7d87aef 427 if((fabs(parentPID)==411 || fabs(parentPID)==413 || fabs(parentPID)==421 || fabs(parentPID)==423 || fabs(parentPID)==431)&& fabs(fPDG)==11)mcInDtoE = kTRUE;
428 if((fabs(parentPID)==511 || fabs(parentPID)==513 || fabs(parentPID)==521 || fabs(parentPID)==523 || fabs(parentPID)==531)&& fabs(fPDG)==11)mcInBtoE = kTRUE;
bd6ee6dd 429 if((mcInBtoE || mcInDtoE) && fabs(mcZvertex)<10.0)
430 {
431 fInputHFEMC->Fill(cent,pTMC);
432 double mcinfo[5];
433 mcinfo[0] = cent;
434 mcinfo[1] = pTMC;
435 mcinfo[2] = 0.0;
436 mcinfo[3] = iHijing;
437 mcinfo[4] = pTMCparent;
438 fInputHFEMC_weight->Fill(mcinfo);
439 }
e7d87aef 440 }
441
5e0e45b3 442
feffe705 443 if(proR<7.0 && fabs(fPDG)==11)fInputAlle->Fill(cent,pTMC);
444
e7d87aef 445 }
446 }
447
bfc7c23b 448 fNoEvents->Fill(0);
449
450 Int_t fNOtrks = fESD->GetNumberOfTracks();
451 const AliESDVertex *pVtx = fESD->GetPrimaryVertex();
452
453 Double_t pVtxZ = -999;
454 pVtxZ = pVtx->GetZ();
455
456 if(TMath::Abs(pVtxZ)>10) return;
457 fNoEvents->Fill(1);
458
f4e0d2d5 459 if(fNOtrks<1) return;
bfc7c23b 460 fNoEvents->Fill(2);
461
462 AliPIDResponse *pidResponse = fInputHandler->GetPIDResponse();
463 if(!pidResponse){
464 AliDebug(1, "Using default PID Response");
465 pidResponse = AliHFEtools::GetDefaultPID(kFALSE, fInputEvent->IsA() == AliAODEvent::Class());
466 }
467
468 fPID->SetPIDResponse(pidResponse);
469
470 fCFM->SetRecEventInfo(fESD);
471
e7d87aef 472 //Float_t cent = -1.;
473 //AliCentrality *centrality = fESD->GetCentrality();
474 //cent = centrality->GetCentralityPercentile("V0M");
bfc7c23b 475 fCent->Fill(cent);
476
b12bc641 477 //if(cent>90.) return;
bfc7c23b 478
479 // Calorimeter info.
480
42c75692 481 FindTriggerClusters();
482
bfc7c23b 483 // make EMCAL array
484 for(Int_t iCluster=0; iCluster<fESD->GetNumberOfCaloClusters(); iCluster++)
485 {
486 AliESDCaloCluster *clust = fESD->GetCaloCluster(iCluster);
68d718e7 487 if(clust && clust->IsEMCAL())
bfc7c23b 488 {
489 double clustE = clust->E();
490 float emcx[3]; // cluster pos
491 clust->GetPosition(emcx);
492 TVector3 clustpos(emcx[0],emcx[1],emcx[2]);
493 double emcphi = clustpos.Phi();
494 double emceta = clustpos.Eta();
fb87d707 495 double calInfo[5];
496 calInfo[0] = emcphi; calInfo[1] = emceta; calInfo[2] = clustE; calInfo[3] = cent; calInfo[4] = clust->Chi2();
f4e0d2d5 497 //fEMCAccE->Fill(calInfo);
498 //if(clustE>3.0)fEMCAccE->Fill(calInfo);
44be9e1d 499 //if(fqahist==1 && clustE>1.5)fEMCAccE->Fill(calInfo);
f4e0d2d5 500 hEMCAccE->Fill(cent,clustE);
bfc7c23b 501 }
502 }
503
504 // Track loop
505 for (Int_t iTracks = 0; iTracks < fESD->GetNumberOfTracks(); iTracks++) {
506 AliESDtrack* track = fESD->GetTrack(iTracks);
507 if (!track) {
508 printf("ERROR: Could not receive track %d\n", iTracks);
509 continue;
510 }
e7d87aef 511
46305ed6 512 int parentlabel = 99999;
513 int parentPID = 99999;
514 int grand_parentlabel = 99999;
515 int grand_parentPID = 99999;
e7d87aef 516 Bool_t mcPho = kFALSE;
517 Bool_t mcDtoE= kFALSE;
518 Bool_t mcBtoE= kFALSE;
93c189c5 519 Bool_t mcOrgPi0 = kFALSE;
520 Bool_t mcOrgEta = kFALSE;
e7d87aef 521 double mcele = -1.;
44be9e1d 522 double mcpT = 0.0;
b12bc641 523 double mcMompT = 0.0;
e4b0faf2 524 //double mcGrandMompT = 0.0;
525 double mcWeight = -10.0;
d113d7cd 526
527 int iHijing = 1;
70448e3b 528 int mcLabel = -1;
d113d7cd 529
e7d87aef 530 if(fmcData && fMC && stack)
531 {
feffe705 532 Int_t label = TMath::Abs(track->GetLabel());
70448e3b 533 mcLabel = track->GetLabel();
bd6ee6dd 534
8efacc22 535 if(mcLabel>-1)
bd6ee6dd 536 {
537
538 Bool_t MChijing = fMC->IsFromBGEvent(label);
539 if(!MChijing)iHijing = 0;
540
541 TParticle* particle = stack->Particle(label);
542 int mcpid = particle->GetPdgCode();
543 mcpT = particle->Pt();
544 //printf("MCpid = %d",mcpid);
545 if(particle->GetFirstMother()>-1)
546 {
547 //int parentlabel = particle->GetFirstMother();
548 //int parentPID = stack->Particle(particle->GetFirstMother())->GetPdgCode();
549 //mcMompT = stack->Particle(particle->GetFirstMother())->Pt();
550 FindMother(particle, parentlabel, parentPID);
551 mcMompT = stack->Particle(parentlabel)->Pt();
552 if((parentPID==22 || parentPID==111 || parentPID==221)&& fabs(mcpid)==11)mcPho = kTRUE;
553 if((fabs(parentPID)==411 || fabs(parentPID)==413 || fabs(parentPID)==421 || fabs(parentPID)==423 || fabs(parentPID)==431)&& fabs(mcpid)==11)mcDtoE = kTRUE;
554 if((fabs(parentPID)==511 || fabs(parentPID)==513 || fabs(parentPID)==521 || fabs(parentPID)==523 || fabs(parentPID)==531)&& fabs(mcpid)==11)mcBtoE = kTRUE;
555
556 // make D->e pT correlation
557 if(mcDtoE)fMomDtoE->Fill(mcpT,mcMompT);
558
559 //cout << "check PID = " << parentPID << endl;
560 //cout << "check pho = " << mcPho << endl;
561 //cout << "check D or B = " << mcDtoE << endl;
562 // pi->e (Dalitz)
563 if(parentPID==111 && fabs(mcpid)==11 && mcMompT>0.0)
564 {
565 //cout << "find pi0->e " << endl;
566 mcOrgPi0 = kTRUE;
567 mcWeight = GetMCweight(mcMompT);
568 }
569 // eta->e (Dalitz)
570 if(parentPID==221 && fabs(mcpid)==11 && mcMompT>0.0)
571 {
572 //cout << "find Eta->e " << endl;
573 mcOrgEta = kTRUE;
574 mcWeight = GetMCweightEta(mcMompT);
575 }
576
577 // access grand parent
578 TParticle* particle_parent = stack->Particle(parentlabel); // get parent pointer
579 //if(particle_parent->GetFirstMother()>-1 && parentPID==22 && fabs(mcpid)==11) // get grand parent g->e
580 if(particle_parent->GetFirstMother()>-1 && (parentPID==22 || parentPID==111) && fabs(mcpid)==11) // get grand parent g->e
581 {
582 //int grand_parentPID = stack->Particle(particle_parent->GetFirstMother())->GetPdgCode();
583 //double pTtmp = stack->Particle(particle_parent->GetFirstMother())->Pt();
584 FindMother(particle_parent, grand_parentlabel, grand_parentPID);
585 double mcGrandpT = stack->Particle(grand_parentlabel)->Pt();
586 if(grand_parentPID==111 && mcGrandpT>0.0)
587 {
588 // check eta->pi0 decay !
589 int grand2_parentlabel = 99999; int grand2_parentPID = 99999;
590 TParticle* particle_grand = stack->Particle(grand_parentlabel); // get parent pointer
591 FindMother(particle_grand, grand2_parentlabel, grand2_parentPID);
592 if(grand2_parentPID==221)
593 {
594 //cout << "find Eta->e " << endl;
595 double mcGrandpT2 = stack->Particle(grand2_parentlabel)->Pt();
596 mcOrgEta = kTRUE;
597 mcWeight = GetMCweight(mcGrandpT2);
598 mcMompT = mcGrandpT2;
599 }
600 else
601 {
602 //cout << "find pi0->e " << endl;
603 mcOrgPi0 = kTRUE;
604 mcWeight = GetMCweight(mcGrandpT);
605 mcMompT = mcGrandpT;
606 }
607 }
608
609 if(grand_parentPID==221 && mcGrandpT>0.0)
610 {
611 //cout << "find Eta->e " << endl;
612 mcOrgEta = kTRUE;
613 mcOrgPi0 = kFALSE;
614 mcWeight = GetMCweightEta(mcGrandpT);
615 mcMompT = mcGrandpT;
616 }
617 }
618 }
619
cc1605fe 620 //cout << "===================="<<endl;
621 //cout << "mcDtoE : " << mcDtoE << endl;
622 //cout << "mcBtoE : " << mcBtoE << endl;
623 //cout << "mcPho : " << mcPho << endl;
624
625 if(fabs(mcpid)==11)mcele= 0.;
626 //cout << "check e: " << mcele << endl;
bd6ee6dd 627 if(fabs(mcpid)==11 && mcDtoE)mcele= 1.;
cc1605fe 628 //cout << "check D->e: " << mcele << endl;
bd6ee6dd 629 if(fabs(mcpid)==11 && mcBtoE)mcele= 2.;
cc1605fe 630 //cout << "check B->e: " << mcele << endl;
bd6ee6dd 631 if(fabs(mcpid)==11 && mcPho)mcele= 3.;
cc1605fe 632 //cout << "check Pho->e: " << mcele << endl;
bd6ee6dd 633
09022877 634 //cout << "check PID " << endl;
fd2126aa 635 if(fabs(mcpid)!=11)
636 {
09022877 637 //cout << "!= 11" << endl;
638 //cout << mcpid << endl;
fd2126aa 639 }
640 if(mcele==-1)
641 {
09022877 642 //cout << "mcele==-1" << endl;
643 //cout << mcele << endl;
644 //cout << mcpid << endl;
fd2126aa 645 }
646
bd6ee6dd 647 } // end of mcLabel>-1
648
649 } // end of MC info.
e7d87aef 650
2b4460a6 651 //cout << "Pi0 = " << mcOrgPi0 << " ; Eta = " << mcOrgEta << endl;
50919258 652 //printf("weight = %f\n",mcWeight);
653
9b3495ff 654 if(TMath::Abs(track->Eta())>0.6) continue;
f09766dd 655 if(TMath::Abs(track->Pt()<2.0)) continue;
bfc7c23b 656
657 fTrackPtBefTrkCuts->Fill(track->Pt());
658 // RecKine: ITSTPC cuts
659 if(!ProcessCutStep(AliHFEcuts::kStepRecKineITSTPC, track)) continue;
660
661 //RecKink
662 if(fRejectKinkMother) { // Quick and dirty fix to reject both kink mothers and daughters
663 if(track->GetKinkIndex(0) != 0) continue;
664 }
665
666 // RecPrim
667 if(!ProcessCutStep(AliHFEcuts::kStepRecPrim, track)) continue;
668
669 // HFEcuts: ITS layers cuts
670 if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsITS, track)) continue;
671
f09766dd 672 // HFE cuts: TPC PID cleanup
673 if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTPC, track)) continue;
674
a6df418c 675 int nTPCcl = track->GetTPCNcls();
96167a04 676 //int nTPCclF = track->GetTPCNclsF(); // warnings
a6df418c 677 int nITS = track->GetNcls(0);
bfc7c23b 678
679 fTrackPtAftTrkCuts->Fill(track->Pt());
680
681 Double_t mom = -999., eop=-999., pt = -999., dEdx=-999., fTPCnSigma=-10, phi=-999., eta=-999.;
feffe705 682 pt = track->Pt();
683 if(pt<2.0)continue;
bfc7c23b 684
685 // Track extrapolation
686
bfc7c23b 687 Int_t charge = track->Charge();
688 fTrkpt->Fill(pt);
689 mom = track->P();
690 phi = track->Phi();
691 eta = track->Eta();
692 dEdx = track->GetTPCsignal();
693 fTPCnSigma = fPID->GetPIDResponse() ? fPID->GetPIDResponse()->NumberOfSigmasTPC(track, AliPID::kElectron) : 1000;
8af58b0d 694
09022877 695 //cout << "nSigma correctoon-----" << endl;
696 //cout << "org = " << fTPCnSigma << endl;
697 if(!fmcData) // nsigma eta correction
698 {
699 double nSigexpCorr = NsigmaCorrection(eta,cent);
700 fTPCnSigma -= nSigexpCorr;
701 }
702
703 //cout << "correction = " << fTPCnSigma << endl;
704
bfc7c23b 705 double ncells = -1.0;
706 double m20 = -1.0;
707 double m02 = -1.0;
708 double disp = -1.0;
709 double rmatch = -1.0;
710 double nmatch = -1.0;
f09766dd 711 double oppstatus = 0.0;
59d998de 712 double emctof = 0.0;
bfc7c23b 713
f09766dd 714 Bool_t fFlagPhotonicElec = kFALSE;
715 Bool_t fFlagConvinatElec = kFALSE;
bfc7c23b 716
717 Int_t clsId = track->GetEMCALcluster();
718 if (clsId>0){
719 AliESDCaloCluster *clust = fESD->GetCaloCluster(clsId);
720 if(clust && clust->IsEMCAL()){
721
722 double clustE = clust->E();
723 eop = clustE/fabs(mom);
09022877 724 //cout << "eop org = "<< eop << endl;
8efacc22 725 if(mcLabel>-1.0)
726 {
727 double mceopcorr = MCEopMeanCorrection(pt,cent);
728 eop += mceopcorr;
729 }
09022877 730 //cout << "eop corr = " << eop << endl;
8efacc22 731
bfc7c23b 732 //double clustT = clust->GetTOF();
733 ncells = clust->GetNCells();
734 m02 = clust->GetM02();
735 m20 = clust->GetM20();
736 disp = clust->GetDispersion();
737 double delphi = clust->GetTrackDx();
738 double deleta = clust->GetTrackDz();
739 rmatch = sqrt(pow(delphi,2)+pow(deleta,2));
740 nmatch = clust->GetNTracksMatched();
59d998de 741 emctof = clust->GetTOF();
742 //cout << "emctof = " << emctof << endl;
bfc7c23b 743
c2801a73 744 if(fTPCnSigma>-1.5 && fTPCnSigma<3.0)
745 {
2b4460a6 746 SelectPhotonicElectron(iTracks,cent,track,fFlagPhotonicElec,fFlagConvinatElec,fTPCnSigma,m20,eop,mcele,mcWeight,iHijing,mcOrgPi0,mcOrgEta);
c2801a73 747 }
748 if(fFlagPhotonicElec)oppstatus = 1.0;
749 if(fFlagConvinatElec)oppstatus = 2.0;
750 if(fFlagPhotonicElec && fFlagConvinatElec)oppstatus = 3.0;
751
feffe705 752 double valdedx[16];
a6df418c 753 valdedx[0] = pt; valdedx[1] = nITS; valdedx[2] = phi; valdedx[3] = eta; valdedx[4] = fTPCnSigma;
615ffe11 754 //valdedx[5] = eop; valdedx[6] = rmatch; valdedx[7] = ncells, valdedx[8] = nTPCclF; valdedx[9] = m20; valdedx[10] = mcpT;
755 valdedx[5] = eop; valdedx[6] = rmatch; valdedx[7] = ncells, valdedx[8] = nmatch; valdedx[9] = m20; valdedx[10] = mcpT;
d19af75d 756 valdedx[11] = cent; valdedx[12] = dEdx; valdedx[13] = oppstatus; valdedx[14] = nTPCcl;
feffe705 757 valdedx[15] = mcele;
85985bb0 758 if(fqahist==1)fEleInfo->Fill(valdedx);
bfc7c23b 759
760
761 }
762 }
e1f0fb74 763
764 //Get Cal info PID response
765 double eop2;
766 double ss[4];
767 Double_t nSigmaEop = fPID->GetPIDResponse()->NumberOfSigmasEMCAL(track,AliPID::kElectron,eop2,ss);
768 if(fTPCnSigma>-1.5 && fTPCnSigma<3.0 && nITS>2.5 && nTPCcl>100)
769 {
770 double valEop[3];
771 valEop[0] = cent;
772 valEop[1] = pt;
773 valEop[2] = nSigmaEop;
774 fElenSigma->Fill(valEop);
775 }
776
777 // ============ PID
778
a6df418c 779 if(nITS<2.5)continue;
780 if(nTPCcl<100)continue;
781
782 CheckNclust->Fill(nTPCcl);
783 CheckNits->Fill(nITS);
784
42c75692 785 fdEdxBef->Fill(mom,fTPCnSigma);
bfc7c23b 786 fTPCnsigma->Fill(mom,fTPCnSigma);
f09766dd 787 if(fTPCnSigma >= -1.0 && fTPCnSigma <= 3)fTrkEovPBef->Fill(pt,eop);
bfc7c23b 788
bfc7c23b 789 Int_t pidpassed = 1;
697ecf6b 790
791 // check reco eff. with TPC
792
793 double phoval[5];
794 phoval[0] = cent;
795 phoval[1] = pt;
796 phoval[2] = fTPCnSigma;
797 phoval[3] = iHijing;
798 phoval[4] = mcMompT;
799
d19af75d 800 if((fTPCnSigma >= -1.0 && fTPCnSigma <= 3) && mcele>-1 && mcPho && mcOrgPi0)
697ecf6b 801 {
802 if(iHijing==1)mcWeight = 1.0;
803 fIncpTMCpho_pi0e_TPC->Fill(phoval,mcWeight);
804 if(fFlagPhotonicElec) fPhoElecPtMC_pi0e_TPC->Fill(phoval,mcWeight);
805 if(fFlagConvinatElec) fSameElecPtMC_pi0e_TPC->Fill(phoval,mcWeight);
806 }
bfc7c23b 807
808 //--- track accepted
809 AliHFEpidObject hfetrack;
810 hfetrack.SetAnalysisType(AliHFEpidObject::kESDanalysis);
811 hfetrack.SetRecTrack(track);
3db00c72 812 double binct = 10.5;
813 if((0.0< cent) && (cent<5.0)) binct = 0.5;
814 if((5.0< cent) && (cent<10.0)) binct = 1.5;
815 if((10.0< cent) && (cent<20.0)) binct = 2.5;
816 if((20.0< cent) && (cent<30.0)) binct = 3.5;
817 if((30.0< cent) && (cent<40.0)) binct = 4.5;
818 if((40.0< cent) && (cent<50.0)) binct = 5.5;
819 if((50.0< cent) && (cent<60.0)) binct = 6.5;
820 if((60.0< cent) && (cent<70.0)) binct = 7.5;
821 if((70.0< cent) && (cent<80.0)) binct = 8.5;
822 if((80.0< cent) && (cent<90.0)) binct = 9.5;
823 if((90.0< cent) && (cent<100.0)) binct = 10.5;
824
825 hfetrack.SetCentrality((int)binct); //added
bfc7c23b 826 hfetrack.SetPbPb();
827 if(!fPID->IsSelected(&hfetrack, NULL, "", fPIDqa)) pidpassed = 0;
828
829 if(pidpassed==0) continue;
830
831 fTrkEovPAft->Fill(pt,eop);
42c75692 832 fdEdxAft->Fill(mom,fTPCnSigma);
bfc7c23b 833
e7d87aef 834 // Fill real data
fb87d707 835 fIncpT->Fill(cent,pt);
bfc7c23b 836 if(fFlagPhotonicElec) fPhoElecPt->Fill(cent,pt);
f09766dd 837 if(fFlagConvinatElec) fSameElecPt->Fill(cent,pt);
fb87d707 838
839 if(m20>0.0 && m20<0.3)
840 {
59d998de 841 fIncpTM20->Fill(cent,pt);
842 ftimingEle->Fill(pt,emctof);
fb87d707 843 if(fFlagPhotonicElec) fPhoElecPtM20->Fill(cent,pt);
844 if(fFlagConvinatElec) fSameElecPtM20->Fill(cent,pt);
845 }
feffe705 846
e7d87aef 847 // MC
70448e3b 848 // check label for electron candidiates
849
850 int idlabel = 1;
851 if(mcLabel==0)idlabel = 0;
852 fLabelCheck->Fill(pt,idlabel);
853 if(mcLabel==0)fgeoFake->Fill(phi,eta);
854
d19af75d 855 if(mcLabel<0 && m20>0.0 && m20<0.3 && fTPCnSigma>-1 && fTPCnSigma<3)
856 {
857 fFakeTrk0->Fill(cent,pt);
858 }
859
860 if(mcele>-1) // select MC electrons
e7d87aef 861 {
3db00c72 862
863 fIncpTMChfeAll->Fill(cent,pt);
864 if(m20>0.0 && m20<0.3)fIncpTMCM20hfeAll->Fill(cent,pt);
d19af75d 865 if(m20>0.0 && m20<0.3 && fTPCnSigma>-1 && fTPCnSigma<3)fFakeTrk1->Fill(cent,pt);
3db00c72 866
e4b0faf2 867 if(mcBtoE || mcDtoE) // select B->e & D->e
e7d87aef 868 {
869 fIncpTMChfe->Fill(cent,pt);
bd6ee6dd 870 //if(m20>0.0 && m20<0.3)fIncpTMCM20hfe->Fill(cent,pt);
871 //if(m20>0.0 && m20<0.3)fIncpTMCM20hfeCheck->Fill(cent,mcpT);
872 if(m20>0.0 && m20<0.3)
873 {
09022877 874 //cout << "MC label = " << mcLabel << endl;
bd6ee6dd 875 fIncpTMCM20hfe->Fill(cent,pt);
876 fIncpTMCM20hfeCheck->Fill(cent,mcpT);
877 fIncpTMCM20hfeCheck_weight->Fill(phoval);
878 }
e7d87aef 879 }
e4b0faf2 880
881 if(mcPho) // select photonic electrons
e7d87aef 882 {
44be9e1d 883
884 fIncpTMCpho->Fill(phoval);
885 if(fFlagPhotonicElec) fPhoElecPtMC->Fill(phoval);
886 if(fFlagConvinatElec) fSameElecPtMC->Fill(phoval);
e7d87aef 887
888 if(m20>0.0 && m20<0.3)
889 {
44be9e1d 890 fIncpTMCM20pho->Fill(phoval);
891 if(fFlagPhotonicElec) fPhoElecPtMCM20->Fill(phoval);
892 if(fFlagConvinatElec) fSameElecPtMCM20->Fill(phoval);
e4b0faf2 893 // pi0->g->e
894 if(mcWeight>-1)
895 {
987053ce 896 if(iHijing==1)mcWeight = 1.0;
93c189c5 897 if(mcOrgPi0)
898 {
697ecf6b 899 fIncpTMCM20pho_pi0e->Fill(phoval,mcWeight);
900 if(fFlagPhotonicElec) fPhoElecPtMCM20_pi0e->Fill(phoval,mcWeight);
901 if(fFlagConvinatElec) fSameElecPtMCM20_pi0e->Fill(phoval,mcWeight);
902 //fIncpTMCM20pho_pi0e->Fill(phoval); // v5-04-02-AN & v5-04-06-AN
903 //if(fFlagPhotonicElec) fPhoElecPtMCM20_pi0e->Fill(phoval);
904 //if(fFlagConvinatElec) fSameElecPtMCM20_pi0e->Fill(phoval);
93c189c5 905 }
906 if(mcOrgEta)
907 {
697ecf6b 908 fIncpTMCM20pho_eta->Fill(phoval,mcWeight);
909 if(fFlagPhotonicElec) fPhoElecPtMCM20_eta->Fill(phoval,mcWeight);
910 if(fFlagConvinatElec) fSameElecPtMCM20_eta->Fill(phoval,mcWeight);
911 //fIncpTMCM20pho_eta->Fill(phoval);
912 //if(fFlagPhotonicElec) fPhoElecPtMCM20_eta->Fill(phoval);
913 //if(fFlagConvinatElec) fSameElecPtMCM20_eta->Fill(phoval);
93c189c5 914 }
915 // --- eta
e4b0faf2 916 }
e7d87aef 917 }
918 }
919 }
bfc7c23b 920 }
921 PostData(1, fOutputList);
922}
923//_________________________________________
924void AliAnalysisTaskHFECal::UserCreateOutputObjects()
925{
42c75692 926 //--- Check MC
927
e7d87aef 928 //Bool_t mcData = kFALSE;
42c75692 929 if(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler())
930 {
e7d87aef 931 fmcData = kTRUE;
42c75692 932 printf("+++++ MC Data available");
933 }
e7d87aef 934 if(fmcData)
42c75692 935 {
936 printf("++++++++= MC analysis \n");
937 }
938 else
939 {
940 printf("++++++++ real data analysis \n");
941 }
942
85985bb0 943 printf("+++++++ QA hist %d",fqahist);
944
fb87d707 945 //---- Geometry
946 fGeom = AliEMCALGeometry::GetInstance("EMCAL_COMPLETEV1");
947
bfc7c23b 948 //--------Initialize PID
42c75692 949 //fPID->SetHasMCData(kFALSE);
e7d87aef 950 fPID->SetHasMCData(fmcData);
bfc7c23b 951 if(!fPID->GetNumberOfPIDdetectors())
952 {
953 fPID->AddDetector("TPC", 0);
954 fPID->AddDetector("EMCAL", 1);
955 }
956
3db00c72 957 Double_t params[4];
78ff1095 958 const char *cutmodel;
3db00c72 959 cutmodel = "pol0";
960 params[0] = -1.0; //sigma min
9b3495ff 961 double maxnSig = 3.0;
962 if(fmcData)
963 {
964 params[0] = -5.0; //sigma min
965 maxnSig = 5.0;
966 }
967
968 for(Int_t a=0;a<11;a++)fPID->ConfigureTPCcentralityCut(a,cutmodel,params,maxnSig);
3db00c72 969
bfc7c23b 970 fPID->SortDetectors();
971 fPIDqa = new AliHFEpidQAmanager();
972 fPIDqa->Initialize(fPID);
a6df418c 973
974 //------- fcut --------------
975 fCuts = new AliHFEcuts();
976 fCuts->CreateStandardCuts();
977 //fCuts->SetMinNClustersTPC(100);
978 fCuts->SetMinNClustersTPC(90);
979 fCuts->SetMinRatioTPCclusters(0.6);
980 fCuts->SetTPCmodes(AliHFEextraCuts::kFound, AliHFEextraCuts::kFoundOverFindable);
981 //fCuts->SetMinNClustersITS(3);
982 fCuts->SetMinNClustersITS(2);
983 fCuts->SetCutITSpixel(AliHFEextraCuts::kAny);
984 fCuts->SetCheckITSLayerStatus(kFALSE);
985 fCuts->SetVertexRange(10.);
986 fCuts->SetTOFPIDStep(kFALSE);
987 fCuts->SetPtRange(2, 50);
988 fCuts->SetMaxImpactParam(3.,3.);
989
bfc7c23b 990 //--------Initialize correction Framework and Cuts
991 fCFM = new AliCFManager;
992 const Int_t kNcutSteps = AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kNcutStepsRecTrack + AliHFEcuts::kNcutStepsDETrack;
993 fCFM->SetNStepParticle(kNcutSteps);
994 for(Int_t istep = 0; istep < kNcutSteps; istep++)
995 fCFM->SetParticleCutsList(istep, NULL);
996
997 if(!fCuts){
998 AliWarning("Cuts not available. Default cuts will be used");
999 fCuts = new AliHFEcuts;
1000 fCuts->CreateStandardCuts();
1001 }
1002 fCuts->Initialize(fCFM);
1003
1004 //---------Output Tlist
1005 fOutputList = new TList();
1006 fOutputList->SetOwner();
1007 fOutputList->Add(fPIDqa->MakeList("PIDQA"));
1008
1009 fNoEvents = new TH1F("fNoEvents","",4,-0.5,3.5) ;
1010 fOutputList->Add(fNoEvents);
1011
02cd82a0 1012 Int_t binsE[5] = {250, 100, 1000, 200, 10};
fb87d707 1013 Double_t xminE[5] = {1.0, -1, 0.0, 0, -0.5};
1014 Double_t xmaxE[5] = {3.5, 1, 100.0, 100, 9.5};
1015 fEMCAccE = new THnSparseD("fEMCAccE","EMC acceptance & E;#eta;#phi;Energy;Centrality;trugCondition;",5,binsE,xminE,xmaxE);
85985bb0 1016 if(fqahist==1)fOutputList->Add(fEMCAccE);
bfc7c23b 1017
f4e0d2d5 1018 hEMCAccE = new TH2F("hEMCAccE","Cluster Energy",200,0,100,100,0,20);
1019 fOutputList->Add(hEMCAccE);
1020
bfc7c23b 1021 fTrkpt = new TH1F("fTrkpt","track pt",100,0,50);
1022 fOutputList->Add(fTrkpt);
1023
1024 fTrackPtBefTrkCuts = new TH1F("fTrackPtBefTrkCuts","track pt before track cuts",100,0,50);
1025 fOutputList->Add(fTrackPtBefTrkCuts);
1026
1027 fTrackPtAftTrkCuts = new TH1F("fTrackPtAftTrkCuts","track pt after track cuts",100,0,50);
1028 fOutputList->Add(fTrackPtAftTrkCuts);
1029
1030 fTPCnsigma = new TH2F("fTPCnsigma", "TPC - n sigma",100,0,50,200,-10,10);
1031 fOutputList->Add(fTPCnsigma);
1032
1033 fTrkEovPBef = new TH2F("fTrkEovPBef","track E/p before HFE pid",100,0,50,100,0,2);
1034 fOutputList->Add(fTrkEovPBef);
1035
1036 fTrkEovPAft = new TH2F("fTrkEovPAft","track E/p after HFE pid",100,0,50,100,0,2);
1037 fOutputList->Add(fTrkEovPAft);
1038
42c75692 1039 fdEdxBef = new TH2F("fdEdxBef","track dEdx vs p before HFE pid",100,0,50,200,-10,10);
bfc7c23b 1040 fOutputList->Add(fdEdxBef);
1041
42c75692 1042 fdEdxAft = new TH2F("fdEdxAft","track dEdx vs p after HFE pid",100,0,50,200,-10,10);
bfc7c23b 1043 fOutputList->Add(fdEdxAft);
1044
02cd82a0 1045 fIncpT = new TH2F("fIncpT","HFE pid electro vs. centrality",200,0,100,100,0,50);
bfc7c23b 1046 fOutputList->Add(fIncpT);
1047
02cd82a0 1048 fIncpTM20 = new TH2F("fIncpTM20","HFE pid electro vs. centrality with M20",200,0,100,100,0,50);
fb87d707 1049 fOutputList->Add(fIncpTM20);
4e01b68c 1050
02cd82a0 1051 Int_t nBinspho[9] = { 200, 100, 500, 12, 50, 4, 200, 8, 100};
b12bc641 1052 Double_t minpho[9] = { 0., 0., 0., -2.5, 0, -0.5, 0,-1.5, 0};
1053 Double_t maxpho[9] = {100., 50., 0.5, 3.5, 1, 3.5, 2, 6.5, 50};
f09766dd 1054
b12bc641 1055 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 1056 if(fqahist==1)fOutputList->Add(fInvmassLS);
bfc7c23b 1057
b12bc641 1058 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 1059 if(fqahist==1)fOutputList->Add(fInvmassULS);
bfc7c23b 1060
8806ce6c 1061 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 1062 if(fqahist==1)fOutputList->Add(fInvmassLSmc);
8806ce6c 1063
1064 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 1065 if(fqahist==1)fOutputList->Add(fInvmassULSmc);
8806ce6c 1066
93c189c5 1067 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 );
1068 fInvmassLSmc0->Sumw2();
1069 fOutputList->Add(fInvmassLSmc0);
1070
1071 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 );
1072 fInvmassLSmc1->Sumw2();
1073 fOutputList->Add(fInvmassLSmc1);
1074
1075 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 );
1076 fInvmassLSmc2->Sumw2();
1077 fOutputList->Add(fInvmassLSmc2);
1078
1079 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 );
1080 fInvmassLSmc3->Sumw2();
1081 fOutputList->Add(fInvmassLSmc3);
1082
1083 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 );
1084 fInvmassULSmc0->Sumw2();
1085 fOutputList->Add(fInvmassULSmc0);
1086
1087 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 );
1088 fInvmassULSmc1->Sumw2();
1089 fOutputList->Add(fInvmassULSmc1);
1090
1091 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 );
1092 fInvmassULSmc2->Sumw2();
1093 fOutputList->Add(fInvmassULSmc2);
1094
1095 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 );
1096 fInvmassULSmc3->Sumw2();
93c189c5 1097 fOutputList->Add(fInvmassULSmc3);
1098
bfc7c23b 1099 fOpeningAngleLS = new TH1F("fOpeningAngleLS","Opening angle for LS pairs",100,0,1);
1100 fOutputList->Add(fOpeningAngleLS);
1101
1102 fOpeningAngleULS = new TH1F("fOpeningAngleULS","Opening angle for ULS pairs",100,0,1);
1103 fOutputList->Add(fOpeningAngleULS);
1104
1105 fPhotoElecPt = new TH1F("fPhotoElecPt", "photonic electron pt",100,0,50);
1106 fOutputList->Add(fPhotoElecPt);
1107
02cd82a0 1108 fPhoElecPt = new TH2F("fPhoElecPt", "Pho-inclusive electron pt",200,0,100,100,0,50);
bfc7c23b 1109 fOutputList->Add(fPhoElecPt);
1110
02cd82a0 1111 fPhoElecPtM20 = new TH2F("fPhoElecPtM20", "Pho-inclusive electron pt with M20",200,0,100,100,0,50);
fb87d707 1112 fOutputList->Add(fPhoElecPtM20);
1113
02cd82a0 1114 fSameElecPt = new TH2F("fSameElecPt", "Same-inclusive electron pt",200,0,100,100,0,50);
f09766dd 1115 fOutputList->Add(fSameElecPt);
1116
02cd82a0 1117 fSameElecPtM20 = new TH2F("fSameElecPtM20", "Same-inclusive electron pt with M20",200,0,100,100,0,50);
fb87d707 1118 fOutputList->Add(fSameElecPtM20);
1119
02cd82a0 1120 fCent = new TH1F("fCent","Centrality",200,0,100) ;
bfc7c23b 1121 fOutputList->Add(fCent);
1122
1123 // Make common binning
44be9e1d 1124 const Double_t kMinP = 0.;
d19af75d 1125 const Double_t kMaxP = 20.;
bfc7c23b 1126
1127 // 1st histogram: TPC dEdx with/without EMCAL (p, pT, TPC Signal, phi, eta, Sig, e/p, ,match, cell, M02, M20, Disp, Centrality, select)
d19af75d 1128 Int_t nBins[16] = { 100, 7, 60, 20, 90, 250, 25, 40, 10, 200, 100, 100, 500, 5, 100, 8};
1129 Double_t min[16] = {kMinP, -0.5, 1.0, -1.0, -5.0, 0, 0, 0, 0.0, 0.0, 0.0, 0, 0, -0.5, 80, -1.5};
1130 Double_t max[16] = {kMaxP, 6.5, 4.0, 1.0, 4.0, 2.5, 0.05, 40, 10, 2.0, 20.0, 100, 100, 4.5, 180, 6.5};
6e92c428 1131 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 1132 if(fqahist==1)fOutputList->Add(fEleInfo);
bfc7c23b 1133
e1f0fb74 1134 // Make common binning
1135 Int_t nBinsEop[3] = { 10, 50, 100};
1136 Double_t minEop[3] = { 0, 0, -5};
1137 Double_t maxEop[3] = {100, 50, 5};
1138 fElenSigma= new THnSparseD("fElenSigma", "Electron nSigma; cent; pT [GeV/c]; nSigma", 3, nBinsEop, minEop, maxEop);
1139 fOutputList->Add(fElenSigma);
1140
1141
fb87d707 1142 //<--- Trigger info
feffe705 1143 /*
fb87d707 1144 fClsEBftTrigCut = new TH1F("fClsEBftTrigCut","cluster E before trigger selection",1000,0,100);
1145 fOutputList->Add(fClsEBftTrigCut);
1146
1147 fClsEAftTrigCut = new TH1F("fClsEAftTrigCut","cluster E if cls has 0 trigcut channel",1000,0,100);
1148 fOutputList->Add(fClsEAftTrigCut);
1149
1150 fClsEAftTrigCut1 = new TH1F("fClsEAftTrigCut1","cluster E if cls with trig channel",1000,0,100);
1151 fOutputList->Add(fClsEAftTrigCut1);
1152
1153 fClsEAftTrigCut2 = new TH1F("fClsEAftTrigCut2","cluster E if cls with trigcut channel",1000,0,100);
1154 fOutputList->Add(fClsEAftTrigCut2);
1155
1156 fClsEAftTrigCut3 = new TH1F("fClsEAftTrigCut3","cluster E if cls with trigcut channel + nCell>Ecorrect",1000,0,100);
1157 fOutputList->Add(fClsEAftTrigCut3);
1158
1159 fClsEAftTrigCut4 = new TH1F("fClsEAftTrigCut4","cluster E if cls with trigcut channel + nCell>Ecorrect + cls time cut",1000,0,100);
1160 fOutputList->Add(fClsEAftTrigCut4);
1161
1162 fClsETime = new TH2F("fClsETime", "Cls time vs E; E; time;",1000,0,100,1000,-0.0000002,0.0000009);
1163 fOutputList->Add(fClsETime);
1164
1165 fClsETime1 = new TH2F("fClsETime1", "Cls time vs E if cls contains trigger channel; E; time;",1000,0,100,1000,-0.0000002,0.0000009);
1166 fOutputList->Add(fClsETime1);
1167
1168 fTrigTimes = new TH1F("fTrigTimes", "Trigger time; time; N;",25,0,25);
1169 fOutputList->Add(fTrigTimes);
1170
42c75692 1171 fCellCheck = new TH2F("fCellCheck", "Cell vs E; E GeV; Cell ID",10,6,26,12000,0,12000);
1172 fOutputList->Add(fCellCheck);
feffe705 1173 */
e7d87aef 1174 //<---------- MC
1175
02cd82a0 1176 fInputHFEMC = new TH2F("fInputHFEMC","Input MC HFE pid electro vs. centrality",200,0,100,100,0,50);
e7d87aef 1177 fOutputList->Add(fInputHFEMC);
1178
02cd82a0 1179 fInputAlle = new TH2F("fInputAlle","Input MC electro vs. centrality",200,0,100,100,0,50);
feffe705 1180 fOutputList->Add(fInputAlle);
1181
02cd82a0 1182 fIncpTMChfe = new TH2F("fIncpTMChfe","MC HFE pid electro vs. centrality",200,0,100,100,0,50);
e7d87aef 1183 fOutputList->Add(fIncpTMChfe);
1184
02cd82a0 1185 fIncpTMChfeAll = new TH2F("fIncpTMChfeAll","MC Alle pid electro vs. centrality",200,0,100,100,0,50);
3db00c72 1186 fOutputList->Add(fIncpTMChfeAll);
1187
02cd82a0 1188 fIncpTMCM20hfe = new TH2F("fIncpTMCM20hfe","MC HFE pid electro vs. centrality with M20",200,0,100,100,0,50);
e7d87aef 1189 fOutputList->Add(fIncpTMCM20hfe);
1190
02cd82a0 1191 fIncpTMCM20hfeAll = new TH2F("fIncpTMCM20hfeAll","MC Alle pid electro vs. centrality with M20",200,0,100,100,0,50);
3db00c72 1192 fOutputList->Add(fIncpTMCM20hfeAll);
1193
60544aea 1194 fIncpTMCM20hfeCheck = new TH2F("fIncpTMCM20hfeCheck","MC HFE pid electro vs. centrality with M20 Check",200,0,100,100,0,50);
1195 fOutputList->Add(fIncpTMCM20hfeCheck);
44be9e1d 1196
2b4460a6 1197 Int_t nBinspho2[5] = { 200, 100, 7, 3, 700};
987053ce 1198 Double_t minpho2[5] = { 0., 0., -2.5, -0.5, 0.};
2b4460a6 1199 Double_t maxpho2[5] = {100., 50., 4.5, 2.5, 70.};
bd6ee6dd 1200
1201 fInputHFEMC_weight = new THnSparseD("fInputHFEMC_weight", "MC HFE electron pt",5,nBinspho2,minpho2,maxpho2);
1202 fOutputList->Add(fInputHFEMC_weight);
1203
1204 fIncpTMCM20hfeCheck_weight = new THnSparseD("fIncpTMCM20hfeCheck_weight", "HFE electron pt with M20",5,nBinspho2,minpho2,maxpho2);
1205 fOutputList->Add(fIncpTMCM20hfeCheck_weight);
1206
987053ce 1207 fIncpTMCpho = new THnSparseD("fIncpTMCpho","MC Pho pid electro vs. centrality",5,nBinspho2,minpho2,maxpho2);
e7d87aef 1208 fOutputList->Add(fIncpTMCpho);
1209
987053ce 1210 fIncpTMCM20pho = new THnSparseD("fIncpTMCM20pho","MC Pho pid electro vs. centrality with M20",5,nBinspho2,minpho2,maxpho2);
e7d87aef 1211 fOutputList->Add(fIncpTMCM20pho);
1212
987053ce 1213 fPhoElecPtMC = new THnSparseD("fPhoElecPtMC", "MC Pho-inclusive electron pt",5,nBinspho2,minpho2,maxpho2);
e7d87aef 1214 fOutputList->Add(fPhoElecPtMC);
1215
987053ce 1216 fPhoElecPtMCM20 = new THnSparseD("fPhoElecPtMCM20", "MC Pho-inclusive electron pt with M20",5,nBinspho2,minpho2,maxpho2);
feffe705 1217 fOutputList->Add(fPhoElecPtMCM20);
e7d87aef 1218
987053ce 1219 fSameElecPtMC = new THnSparseD("fSameElecPtMC", "MC Same-inclusive electron pt",5,nBinspho2,minpho2,maxpho2);
e7d87aef 1220 fOutputList->Add(fSameElecPtMC);
1221
987053ce 1222 fSameElecPtMCM20 = new THnSparseD("fSameElecPtMCM20", "MC Same-inclusive electron pt with M20",5,nBinspho2,minpho2,maxpho2);
e7d87aef 1223 fOutputList->Add(fSameElecPtMCM20);
1224
987053ce 1225 fIncpTMCM20pho_pi0e = new THnSparseD("fIncpTMCM20pho_pi0e","MC Pho pi0->e pid electro vs. centrality with M20",5,nBinspho2,minpho2,maxpho2);
1226 fIncpTMCM20pho_pi0e->Sumw2();
e4b0faf2 1227 fOutputList->Add(fIncpTMCM20pho_pi0e);
1228
987053ce 1229 fPhoElecPtMCM20_pi0e = new THnSparseD("fPhoElecPtMCM20_pi0e", "MC Pho-inclusive electron pt with M20 pi0->e",5,nBinspho2,minpho2,maxpho2);
1230 fPhoElecPtMCM20_pi0e->Sumw2();
e4b0faf2 1231 fOutputList->Add(fPhoElecPtMCM20_pi0e);
1232
987053ce 1233 fSameElecPtMCM20_pi0e = new THnSparseD("fSameElecPtMCM20_pi0e", "MC Same-inclusive electron pt pi0->e",5,nBinspho2,minpho2,maxpho2);
1234 fSameElecPtMCM20_pi0e->Sumw2();
e4b0faf2 1235 fOutputList->Add(fSameElecPtMCM20_pi0e);
697ecf6b 1236 //
93c189c5 1237 fIncpTMCM20pho_eta = new THnSparseD("fIncpTMCM20pho_eta","MC Pho pi0->e pid electro vs. centrality with M20",5,nBinspho2,minpho2,maxpho2);
1238 fIncpTMCM20pho_eta->Sumw2();
1239 fOutputList->Add(fIncpTMCM20pho_eta);
1240
1241 fPhoElecPtMCM20_eta = new THnSparseD("fPhoElecPtMCM20_eta", "MC Pho-inclusive electron pt with M20 pi0->e",5,nBinspho2,minpho2,maxpho2);
1242 fPhoElecPtMCM20_eta->Sumw2();
1243 fOutputList->Add(fPhoElecPtMCM20_eta);
1244
1245 fSameElecPtMCM20_eta = new THnSparseD("fSameElecPtMCM20_eta", "MC Same-inclusive electron pt pi0->e",5,nBinspho2,minpho2,maxpho2);
1246 fSameElecPtMCM20_eta->Sumw2();
1247 fOutputList->Add(fSameElecPtMCM20_eta);
697ecf6b 1248 // ------------
1249 fIncpTMCpho_pi0e_TPC = new THnSparseD("fIncpTMCpho_pi0e_TPC","MC Pho pi0->e pid electro vs. centrality with M20",5,nBinspho2,minpho2,maxpho2);
1250 fIncpTMCpho_pi0e_TPC->Sumw2();
1251 fOutputList->Add(fIncpTMCpho_pi0e_TPC);
1252
1253 fPhoElecPtMC_pi0e_TPC = new THnSparseD("fPhoElecPtMC_pi0e_TPC", "MC Pho-inclusive electron pt with pi0->e",5,nBinspho2,minpho2,maxpho2);
1254 fPhoElecPtMC_pi0e_TPC->Sumw2();
1255 fOutputList->Add(fPhoElecPtMC_pi0e_TPC);
1256
1257 fSameElecPtMC_pi0e_TPC = new THnSparseD("fSameElecPtMC_pi0e_TPC", "MC Same-inclusive electron pt pi0->e",5,nBinspho2,minpho2,maxpho2);
1258 fSameElecPtMC_pi0e_TPC->Sumw2();
1259 fOutputList->Add(fSameElecPtMC_pi0e_TPC);
1260 //-------------
e4b0faf2 1261
a6df418c 1262 CheckNclust = new TH1D("CheckNclust","cluster check",200,0,200);
1263 fOutputList->Add(CheckNclust);
1264
1265 CheckNits = new TH1D("CheckNits","ITS cluster check",8,-0.5,7.5);
1266 fOutputList->Add(CheckNits);
68d718e7 1267 /*
d113d7cd 1268 Hpi0pTcheck = new TH2D("Hpi0pTcheck","Pi0 pT from Hijing",100,0,50,3,-0.5,2.5);
1269 fOutputList->Add(Hpi0pTcheck);
1270
19325033 1271 HETApTcheck = new TH2D("HETApTcheck","Eta pT from Hijing",100,0,50,3,-0.5,2.5);
1272 fOutputList->Add(HETApTcheck);
68d718e7 1273 */
1274
1275 Int_t nBinspho3[3] = { 200, 100, 3};
1276 Double_t minpho3[3] = { 0., 0., -0.5};
1277 Double_t maxpho3[3] = {100., 50., 2.5};
19325033 1278
68d718e7 1279 Hpi0pTcheck = new THnSparseD("Hpi0pTcheck","Pi0 pT from Hijing",3,nBinspho3,minpho3,maxpho3);
1280 fOutputList->Add(Hpi0pTcheck);
1281
1282 HETApTcheck = new THnSparseD("HETApTcheck","Eta pT from Hijing",3,nBinspho3,minpho3,maxpho3);
1283 fOutputList->Add(HETApTcheck);
1284 //--
e4b0faf2 1285 HphopTcheck = new TH2D("HphopTcheck","Pho pT from Hijing",100,0,50,3,-0.5,2.5);
1286 fOutputList->Add(HphopTcheck);
5e0e45b3 1287 //
1288 HDpTcheck = new TH2D("HDpTcheck","D pT from Hijing",100,0,50,3,-0.5,2.5);
1289 fOutputList->Add(HDpTcheck);
1290
1291 HBpTcheck = new TH2D("HBpTcheck","B pT from Hijing",100,0,50,3,-0.5,2.5);
1292 fOutputList->Add(HBpTcheck);
1293 //
e4b0faf2 1294
93c189c5 1295 fpTCheck = new TH1D("fpTCheck","pT check",500,0,50);
1296 fOutputList->Add(fpTCheck);
1297
50919258 1298 fMomDtoE = new TH2D("fMomDtoE","D->E pT correlations;e p_{T} GeV/c;D p_{T} GeV/c",400,0,40,400,0,40);
1299 fOutputList->Add(fMomDtoE);
d113d7cd 1300
70448e3b 1301 fLabelCheck = new TH2D("fLabelCheck","MC label",50,0,50,5,-1.5,3.5);
1302 fOutputList->Add(fLabelCheck);
1303
1304 fgeoFake = new TH2D("fgeoFake","Label==0 eta and phi",628,0,6.28,200,-1,1);
1305 fOutputList->Add(fgeoFake);
1306
d19af75d 1307 fFakeTrk0 = new TH2D("fFakeTrk0","fake trakcs",10,0,100,20,0,20);
1308 fOutputList->Add(fFakeTrk0);
1309
1310 fFakeTrk1 = new TH2D("fFakeTrk1","true all e a.f. eID",10,0,100,20,0,20);
1311 fOutputList->Add(fFakeTrk1);
1312
59d998de 1313 ftimingEle = new TH2D("ftimingEle","electron TOF",100,0,20,100,1e-7,1e-6);
1314 fOutputList->Add(ftimingEle);
1315
09022877 1316 // eta correction
1317 // note: parameters 01/31new.TPCnSigmaEtaDep
1318 // 70-90 delta_eta = 0.2
1319
1320 double etaval[12] = {-0.55,-0.45,-0.35,-0.25,-0.15,-0.05,0.05,0.15,0.25,0.35,0.45,0.55};
1321 double corr0[12]= {-0.569177,-0.528844,-0.391979,-0.165494,0.0283495,0.156171,0.266353,0.13103,-0.0250842,-0.274089,-0.45481,-0.536291}; // 0-10 (done)
1322 double corr1[12]= {-0.404742,-0.278953,-0.218069,0.00139927,0.191412,0.354403,0.524594,0.341778,0.244199,-0.112146,-0.160692,-0.352832}; // 10-20 (done)
1323 double corr2[12] = {-0.306007,-0.16821,-0.0248635,0.202233,0.447051,0.497197,0.712251,0.433482,0.337907,0.168426,-0.0693229,-0.0728351}; // 20-30 (done)
1324 double corr3[12] = {-0.13884,-0.0503553,0.104403,0.389773,0.50697,0.539048,0.751642,0.655636,0.518563,0.308156,0.0361159,-0.0491439}; // 30-40 (done)
1325 double corr4[12] = {-0.0319431,0.0808711,0.208774,0.443217,0.557762,0.61453,0.889519,0.808282,0.620394,0.267092,0.15241,-0.0458664}; // 40-50 (done)
1326 double corr5[12] = {-0.130625,0.0189124,0.190344,0.467431,0.546353,0.672251,0.731541,0.802101,0.437108,0.294081,0.193682,0.159074}; // 50-70(done)
1327 double corr6[12] = {0.0600197,0.0600197,0.358366,0.358366,0.973734,0.973734,0.759812,0.759812,0.667861,0.667861,0.415635,0.415635}; // 70-90(done)
1328
1329 fnSigEtaCorr[0] = new TGraphErrors(12,etaval,corr0); // 0-10
1330 fnSigEtaCorr[1] = new TGraphErrors(12,etaval,corr1); // 10-20
1331 fnSigEtaCorr[2] = new TGraphErrors(12,etaval,corr2); // 20-30
1332 fnSigEtaCorr[3] = new TGraphErrors(12,etaval,corr3); // 30-40
1333 fnSigEtaCorr[4] = new TGraphErrors(12,etaval,corr4); // 40-50
1334 fnSigEtaCorr[5] = new TGraphErrors(12,etaval,corr5); // 50-70
1335 fnSigEtaCorr[6] = new TGraphErrors(12,etaval,corr6); // 70-90
1336
1337
bfc7c23b 1338 PostData(1,fOutputList);
1339}
1340
1341//________________________________________________________________________
1342void AliAnalysisTaskHFECal::Terminate(Option_t *)
1343{
1344 // Info("Terminate");
1345 AliAnalysisTaskSE::Terminate();
1346}
1347
1348//________________________________________________________________________
1349Bool_t AliAnalysisTaskHFECal::ProcessCutStep(Int_t cutStep, AliVParticle *track)
1350{
1351 // Check single track cuts for a given cut step
1352 const Int_t kMCOffset = AliHFEcuts::kNcutStepsMCTrack;
1353 if(!fCFM->CheckParticleCuts(cutStep + kMCOffset, track)) return kFALSE;
1354 return kTRUE;
1355}
1356//_________________________________________
f09766dd 1357//void AliAnalysisTaskHFECal::SelectPhotonicElectron(Int_t itrack, Double_t cent, AliESDtrack *track, Bool_t &fFlagPhotonicElec)
c2801a73 1358//void AliAnalysisTaskHFECal::SelectPhotonicElectron(Int_t itrack, Double_t cent, AliESDtrack *track, Bool_t &fFlagPhotonicElec, Bool_t &fFlagConvinatElec, Double_t nSig)
2b4460a6 1359void 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 1360{
1361 //Identify non-heavy flavour electrons using Invariant mass method
1362
fb87d707 1363 fTrackCuts->SetAcceptKinkDaughters(kFALSE);
feffe705 1364 fTrackCuts->SetRequireTPCRefit(kTRUE);
85985bb0 1365 fTrackCuts->SetRequireITSRefit(kTRUE);
feffe705 1366 fTrackCuts->SetEtaRange(-0.9,0.9);
f09766dd 1367 //fTrackCuts->SetRequireSigmaToVertex(kTRUE);
bfc7c23b 1368 fTrackCuts->SetMaxChi2PerClusterTPC(3.5);
feffe705 1369 fTrackCuts->SetMinNClustersTPC(90);
bfc7c23b 1370
1371 const AliESDVertex *pVtx = fESD->GetPrimaryVertex();
1372
93c189c5 1373 double ptEle = track->Pt(); //add
1374 if(ibgevent==0 && w > 0.0)
1375 {
1376 fpTCheck->Fill(ptEle,w);
1377 }
1378
bfc7c23b 1379 Bool_t flagPhotonicElec = kFALSE;
f09766dd 1380 Bool_t flagConvinatElec = kFALSE;
bfc7c23b 1381
8806ce6c 1382 int p1 = 0;
1383 if(mce==3)
1384 {
1385 Int_t label = TMath::Abs(track->GetLabel());
1386 TParticle* particle = stack->Particle(label);
1387 p1 = particle->GetFirstMother();
1388 }
1389
5e65985c 1390 //for(Int_t jTracks = itrack+1; jTracks<fESD->GetNumberOfTracks(); jTracks++){
1391 for(Int_t jTracks = 0; jTracks<fESD->GetNumberOfTracks(); jTracks++){
bfc7c23b 1392 AliESDtrack* trackAsso = fESD->GetTrack(jTracks);
1393 if (!trackAsso) {
1394 printf("ERROR: Could not receive track %d\n", jTracks);
1395 continue;
1396 }
5e65985c 1397 if(itrack==jTracks)continue;
d86a815c 1398 int jbgevent = 0;
5e65985c 1399
8806ce6c 1400 int p2 = 0;
1401 if(mce==3)
1402 {
1403 Int_t label2 = TMath::Abs(trackAsso->GetLabel());
1404 TParticle* particle2 = stack->Particle(label2);
d86a815c 1405 Bool_t MChijing_ass = fMC->IsFromBGEvent(label2);
1406 if(MChijing_ass)jbgevent =1;
8806ce6c 1407 if(particle2->GetFirstMother()>-1)
1408 p2 = particle2->GetFirstMother();
1409 }
1410
f09766dd 1411 Double_t dEdxAsso = -999., ptPrim=-999., ptAsso=-999., openingAngle = -999.;
bfc7c23b 1412 Double_t mass=999., width = -999;
1413 Bool_t fFlagLS=kFALSE, fFlagULS=kFALSE;
1414
93c189c5 1415 //ptPrim = track->Pt();
1416 ptPrim = ptEle;
f09766dd 1417
bfc7c23b 1418 dEdxAsso = trackAsso->GetTPCsignal();
1419 ptAsso = trackAsso->Pt();
1420 Int_t chargeAsso = trackAsso->Charge();
1421 Int_t charge = track->Charge();
1422
b12bc641 1423
3db00c72 1424 if(ptAsso <0.5) continue;
bfc7c23b 1425 if(!fTrackCuts->AcceptTrack(trackAsso)) continue;
f09766dd 1426 if(dEdxAsso <65 || dEdxAsso>100) continue; //11a pass1
bfc7c23b 1427
1428 Int_t fPDGe1 = 11; Int_t fPDGe2 = 11;
1429 if(charge>0) fPDGe1 = -11;
1430 if(chargeAsso>0) fPDGe2 = -11;
b12bc641 1431
1432 //printf("chargeAsso = %d\n",chargeAsso);
1433 //printf("charge = %d\n",charge);
bfc7c23b 1434 if(charge == chargeAsso) fFlagLS = kTRUE;
1435 if(charge != chargeAsso) fFlagULS = kTRUE;
1436
b12bc641 1437 //printf("fFlagLS = %d\n",fFlagLS);
1438 //printf("fFlagULS = %d\n",fFlagULS);
7edfba87 1439 printf("\n");
b12bc641 1440
bfc7c23b 1441 AliKFParticle ge1(*track, fPDGe1);
1442 AliKFParticle ge2(*trackAsso, fPDGe2);
1443 AliKFParticle recg(ge1, ge2);
1444
d19af75d 1445 // vertex
bfc7c23b 1446 AliKFVertex primV(*pVtx);
1447 primV += recg;
1448 recg.SetProductionVertex(primV);
1449
d19af75d 1450 // check chi2
1451 if(recg.GetNDF()<1) continue;
f9447b3b 1452
1453 // mass const.
1454 recg.SetMassConstraint(0,0.0001);
1455 // v5-04-50-AN no constrain
1456 recg.GetMass(mass,width);
1457
c3e32eae 1458 Double_t chi2recg = recg.GetChi2()/recg.GetNDF();
1459 if(TMath::Sqrt(TMath::Abs(chi2recg))>30.) continue;
1460
f9447b3b 1461 // angle
bfc7c23b 1462 openingAngle = ge1.GetAngle(ge2);
1463 if(fFlagLS) fOpeningAngleLS->Fill(openingAngle);
1464 if(fFlagULS) fOpeningAngleULS->Fill(openingAngle);
1465
c2801a73 1466 double ishower = 0;
1467 if(shower>0.0 && shower<0.3)ishower = 1;
1468
b12bc641 1469 double phoinfo[9];
f09766dd 1470 phoinfo[0] = cent;
1471 phoinfo[1] = ptPrim;
1472 phoinfo[2] = mass;
3db00c72 1473 phoinfo[3] = nSig;
4e01b68c 1474 phoinfo[4] = openingAngle;
c2801a73 1475 phoinfo[5] = ishower;
1476 phoinfo[6] = ep;
1477 phoinfo[7] = mce;
b12bc641 1478 phoinfo[8] = ptAsso;
f09766dd 1479
1480 if(fFlagLS) fInvmassLS->Fill(phoinfo);
1481 if(fFlagULS) fInvmassULS->Fill(phoinfo);
d86a815c 1482 if(fFlagLS && ibgevent==0 && jbgevent==0) fInvmassLSmc->Fill(phoinfo,w);
93c189c5 1483 if(fFlagULS && ibgevent==0 && jbgevent==0)
1484 {
1485 fInvmassULSmc->Fill(phoinfo,w);
1486 }
b12bc641 1487 //printf("fInvmassCut %f\n",fInvmassCut);
1488 //printf("openingAngle %f\n",fOpeningAngleCut);
1489
f9447b3b 1490 // angle cut
ffa14dda 1491 if(openingAngle > fOpeningAngleCut) continue;
bfc7c23b 1492
d86a815c 1493 // for real data
93c189c5 1494 //printf("mce =%f\n",mce);
d86a815c 1495 if(mce<-0.5) // mce==-1. is real
1496 {
93c189c5 1497 //printf("Real data\n");
d86a815c 1498 if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec){
1499 //if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec && (p1==p2)){ <--- only MC train (55,56) v5-03-68-AN & 69 for check
1500 flagPhotonicElec = kTRUE;
1501 }
1502 if(mass<fInvmassCut && fFlagLS && !flagConvinatElec){
1503 flagConvinatElec = kTRUE;
1504 }
1505 }
1506 // for MC data
1507 else
1508 {
93c189c5 1509 //printf("MC data\n");
1510
1511 if(w>0.0)
1512 {
2b4460a6 1513 //cout << "tagpi0 = " << tagpi0 << " ; tageta = " << tageta << endl;
1514 if(fFlagLS && ibgevent==0 && jbgevent==0 && tagpi0) fInvmassLSmc0->Fill(ptPrim,mass,w);
1515 if(fFlagULS && ibgevent==0 && jbgevent==0 && tagpi0) fInvmassULSmc0->Fill(ptPrim,mass,w);
1516 if(fFlagLS && ibgevent==0 && jbgevent==0 && tageta) fInvmassLSmc1->Fill(ptPrim,mass,w);
1517 if(fFlagULS && ibgevent==0 && jbgevent==0 && tageta) fInvmassULSmc1->Fill(ptPrim,mass,w);
1518 if(fFlagLS && ibgevent==0 && jbgevent==0 && (p1==p2) && tagpi0) fInvmassLSmc2->Fill(ptPrim,mass,w);
1519 if(fFlagULS && ibgevent==0 && jbgevent==0 && (p1==p2) && tagpi0) fInvmassULSmc2->Fill(ptPrim,mass,w);
1520 if(fFlagLS && ibgevent==0 && jbgevent==0 && (p1==p2) && tageta) fInvmassLSmc3->Fill(ptPrim,mass,w);
1521 if(fFlagULS && ibgevent==0 && jbgevent==0 && (p1==p2) && tageta) fInvmassULSmc3->Fill(ptPrim,mass,w);
93c189c5 1522 }
d86a815c 1523 if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec && (ibgevent==jbgevent)){
1524 //if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec && (p1==p2)){ <--- only MC train (55,56) v5-03-68-AN & 69 for check
1525 flagPhotonicElec = kTRUE;
1526 }
1527 if(mass<fInvmassCut && fFlagLS && !flagConvinatElec && (ibgevent==jbgevent)){
1528 flagConvinatElec = kTRUE;
1529 }
1530 }
1531
bfc7c23b 1532 }
1533 fFlagPhotonicElec = flagPhotonicElec;
f09766dd 1534 fFlagConvinatElec = flagConvinatElec;
bfc7c23b 1535
1536}
8806ce6c 1537//-------------------------------------------
46305ed6 1538
1539void AliAnalysisTaskHFECal::FindMother(TParticle* part, int &label, int &pid)
1540{
1541 //int label = 99999;
1542 //int pid = 99999;
1543
1544 if(part->GetFirstMother()>-1)
1545 {
1546 label = part->GetFirstMother();
1547 pid = stack->Particle(label)->GetPdgCode();
1548 }
70448e3b 1549 //cout << "Find Mother : label = " << label << " ; pid" << pid << endl;
46305ed6 1550}
1551
8806ce6c 1552double AliAnalysisTaskHFECal::GetMCweight(double mcPi0pT)
1553{
1554 double weight = 1.0;
1555
1556 if(mcPi0pT>0.0 && mcPi0pT<5.0)
1557 {
1558 weight = 0.323*mcPi0pT/(TMath::Exp(-1.6+0.767*mcPi0pT+0.0285*mcPi0pT*mcPi0pT));
1559 }
1560 else
1561 {
1562 weight = 115.0/(0.718*mcPi0pT*TMath::Power(mcPi0pT,3.65));
1563 }
1564 return weight;
1565}
fb87d707 1566
93c189c5 1567double AliAnalysisTaskHFECal::GetMCweightEta(double mcEtapT)
1568{
1569 double weight = 1.0;
1570
1571 weight = 223.3/TMath::Power((TMath::Exp(-0.17*mcEtapT-0.0322*mcEtapT*mcEtapT)+mcEtapT/1.69),5.65);
1572 return weight;
1573}
1574
1575
fb87d707 1576//_________________________________________
1577void AliAnalysisTaskHFECal::FindTriggerClusters()
1578{
bd6ee6dd 1579 //cout << "finding trigger patch" << endl;
fb87d707 1580 // constants
1581 const int nModuleCols = 2;
1582 const int nModuleRows = 5;
1583 const int nColsFeeModule = 48;
1584 const int nRowsFeeModule = 24;
1585 const int nColsFaltroModule = 24;
1586 const int nRowsFaltroModule = 12;
1587 //const int faltroWidthMax = 20;
1588
1589 // part 1, trigger extraction -------------------------------------
1590 Int_t globCol, globRow;
1591 //Int_t ntimes=0, nTrigChannel=0, nTrigChannelCut=0, trigInCut=0;
1592 Int_t ntimes=0, nTrigChannel=0, nTrigChannelCut=0;
1593
1594 //Int_t trigtimes[faltroWidthMax];
1595 Double_t cellTime[nColsFeeModule*nModuleCols][nRowsFeeModule*nModuleRows];
1596 Double_t cellEnergy[nColsFeeModule*nModuleCols][nRowsFeeModule*nModuleRows];
1597 //Double_t fTrigCutLow = 6;
1598 //Double_t fTrigCutHigh = 10;
1599 Double_t fTimeCutLow = 469e-09;
1600 Double_t fTimeCutHigh = 715e-09;
1601
1602 AliESDCaloTrigger * fCaloTrigger = fESD->GetCaloTrigger( "EMCAL" );
1603
1604 // erase trigger maps
1605 for(Int_t i = 0; i < nColsFaltroModule*nModuleCols; i++ )
1606 {
1607 for(Int_t j = 0; j < nRowsFaltroModule*nModuleRows; j++ )
1608 {
1609 ftriggersCut[i][j] = 0;
1610 ftriggers[i][j] = 0;
1611 ftriggersTime[i][j] = 0;
1612 }
1613 }
1614
1615 Int_t iglobCol=0, iglobRow=0;
1616 // go through triggers
1617 if( fCaloTrigger->GetEntries() > 0 )
1618 {
1619 // needs reset
1620 fCaloTrigger->Reset();
1621 while( fCaloTrigger->Next() )
1622 {
1623 fCaloTrigger->GetPosition( globCol, globRow );
1624 fCaloTrigger->GetNL0Times( ntimes );
1625 /*
1626 // no L0s
1627 if( ntimes < 1 ) continue;
1628 // get precise timings
1629 fCaloTrigger->GetL0Times( trigtimes );
1630 trigInCut = 0;
1631 for(Int_t i = 0; i < ntimes; i++ )
1632 {
1633 // save the first trigger time in channel
1634 if( i == 0 || triggersTime[globCol][globRow] > trigtimes[i] )
1635 triggersTime[globCol][globRow] = trigtimes[i];
1636 //printf("trigger times: %d\n",trigtimes[i]);
1637 // check if in cut
1638 if(trigtimes[i] > fTrigCutLow && trigtimes[i] < fTrigCutHigh )
1639 trigInCut = 1;
1640
1641 fTrigTimes->Fill(trigtimes[i]);
1642 }
1643 */
1644
1645 //L1 analysis from AliAnalysisTaskEMCALTriggerQA
1646 Int_t bit = 0;
1647 fCaloTrigger->GetTriggerBits(bit);
bd6ee6dd 1648 //cout << "bit = " << bit << endl;
fb87d707 1649
1650 Int_t ts = 0;
1651 fCaloTrigger->GetL1TimeSum(ts);
bd6ee6dd 1652 //cout << "ts = " << ts << endl;
fb87d707 1653 if (ts > 0)ftriggers[globCol][globRow] = 1;
1654 // number of triggered channels in event
1655 nTrigChannel++;
1656 // ... inside cut
1657 if(ts>0 && (bit >> 6 & 0x1))
1658 {
1659 iglobCol = globCol;
1660 iglobRow = globRow;
1661 nTrigChannelCut++;
bd6ee6dd 1662 //cout << "ts cut = " << ts << endl;
1663 //cout << "globCol = " << globCol << endl;
1664 //cout << "globRow = " << globRow << endl;
fb87d707 1665 ftriggersCut[globCol][globRow] = 1;
1666 }
1667
1668 } // calo trigger entries
1669 } // has calo trigger entries
1670
1671 // part 2 go through the clusters here -----------------------------------
bd6ee6dd 1672 //cout << " part 2 go through the clusters here ----------------------------------- " << endl;
fb87d707 1673 Int_t nCluster=0, nCell=0, iCell=0, gCell=0;
60d77596 1674 Short_t cellAddr, nSACell;
1675 Int_t mclabel;
fb87d707 1676 //Int_t nSACell, iSACell, mclabel;
1677 Int_t iSACell;
42c75692 1678 Double_t cellAmp=0, cellTimeT=0, clusterTime=0, efrac=0;
fb87d707 1679 Int_t nSupMod, nModule, nIphi, nIeta, iphi, ieta, gphi, geta, feta, fphi;
1680
1681 TRefArray *fCaloClusters = new TRefArray();
1682 fESD->GetEMCALClusters( fCaloClusters );
1683 nCluster = fCaloClusters->GetEntries();
1684
1685
1686 // save all cells times if there are clusters
1687 if( nCluster > 0 ){
1688 // erase time map
1689 for(Int_t i = 0; i < nColsFeeModule*nModuleCols; i++ ){
1690 for(Int_t j = 0; j < nRowsFeeModule*nModuleRows; j++ ){
1691 cellTime[i][j] = 0.;
1692 cellEnergy[i][j] = 0.;
1693 }
1694 }
1695
1696 // get cells
1697 AliESDCaloCells *fCaloCells = fESD->GetEMCALCells();
1698 //AliVCaloCells fCaloCells = fESD->GetEMCALCells();
1699 nSACell = fCaloCells->GetNumberOfCells();
1700 for(iSACell = 0; iSACell < nSACell; iSACell++ ){
1701 // get the cell info *fCal
1702 fCaloCells->GetCell( iSACell, cellAddr, cellAmp, cellTimeT , mclabel, efrac);
1703
1704 // get cell position
1705 fGeom->GetCellIndex( cellAddr, nSupMod, nModule, nIphi, nIeta );
1706 fGeom->GetCellPhiEtaIndexInSModule( nSupMod,nModule, nIphi, nIeta, iphi, ieta);
1707
1708 // convert co global phi eta
1709 gphi = iphi + nRowsFeeModule*(nSupMod/2);
1710 geta = ieta + nColsFeeModule*(nSupMod%2);
1711
1712 // save cell time and energy
1713 cellTime[geta][gphi] = cellTimeT;
1714 cellEnergy[geta][gphi] = cellAmp;
1715
1716 }
1717 }
1718
1719 Int_t nClusterTrig, nClusterTrigCut;
1720 UShort_t *cellAddrs;
1721 Double_t clsE=-999, clsEta=-999, clsPhi=-999;
1722 Float_t clsPos[3] = {0.,0.,0.};
1723
1724 for(Int_t icl=0; icl<fESD->GetNumberOfCaloClusters(); icl++)
1725 {
1726 AliESDCaloCluster *cluster = fESD->GetCaloCluster(icl);
1727 if(!cluster || !cluster->IsEMCAL()) continue;
1728
1729 // get cluster cells
1730 nCell = cluster->GetNCells();
1731
1732 // get cluster energy
1733 clsE = cluster->E();
1734
1735 // get cluster position
1736 cluster->GetPosition(clsPos);
1737 TVector3 clsPosVec(clsPos[0],clsPos[1],clsPos[2]);
1738 clsEta = clsPosVec.Eta();
1739 clsPhi = clsPosVec.Phi();
1740
1741 // get the cell addresses
1742 cellAddrs = cluster->GetCellsAbsId();
1743
1744 // check if the cluster contains cell, that was marked as triggered
1745 nClusterTrig = 0;
1746 nClusterTrigCut = 0;
1747
1748 // loop the cells to check, if cluser in acceptance
1749 // any cluster with a cell outside acceptance is not considered
1750 for( iCell = 0; iCell < nCell; iCell++ )
1751 {
42c75692 1752 // check hot cell
feffe705 1753 //if(clsE>6.0)fCellCheck->Fill(clsE,cellAddrs[iCell]);
42c75692 1754
fb87d707 1755 // get cell position
1756 fGeom->GetCellIndex( cellAddrs[iCell], nSupMod, nModule, nIphi, nIeta );
1757 fGeom->GetCellPhiEtaIndexInSModule( nSupMod,nModule, nIphi, nIeta, iphi, ieta);
1758
1759 // convert co global phi eta
1760 gphi = iphi + nRowsFeeModule*(nSupMod/2);
1761 geta = ieta + nColsFeeModule*(nSupMod%2);
1762
1763 if( cellTime[geta][gphi] > 0. ){
1764 clusterTime += cellTime[geta][gphi];
1765 gCell++;
1766 }
1767
1768 // get corresponding FALTRO
1769 fphi = gphi / 2;
1770 feta = geta / 2;
1771
bd6ee6dd 1772 //cout << "fphi = " << fphi << endl;
1773 //cout << "feta = " << feta << endl;
1774
fb87d707 1775 // try to match with a triggered
1776 if( ftriggers[feta][fphi]==1)
1777 { nClusterTrig++;
1778 }
1779 if( ftriggersCut[feta][fphi]==1)
1780 { nClusterTrigCut++;
1781 }
1782
bd6ee6dd 1783 //cout << "nClusterTrigCut : " << nClusterTrigCut << endl;
1784
fb87d707 1785 } // cells
1786
1787
1788 if( gCell > 0 )
1789 clusterTime = clusterTime / (Double_t)gCell;
1790 // fix the reconstriction code time 100ns jumps
1791 if( fESD->GetBunchCrossNumber() % 4 < 2 )
1792 clusterTime -= 0.0000001;
1793
feffe705 1794 //fClsETime->Fill(clsE,clusterTime);
1795 //fClsEBftTrigCut->Fill(clsE);
fb87d707 1796
1797 if(nClusterTrig>0){
feffe705 1798 //fClsETime1->Fill(clsE,clusterTime);
fb87d707 1799 }
1800
1801 if(nClusterTrig>0){
1802 cluster->SetChi2(1);
feffe705 1803 //fClsEAftTrigCut1->Fill(clsE);
fb87d707 1804 }
1805
1806 if(nClusterTrigCut>0){
1807 cluster->SetChi2(2);
feffe705 1808 //fClsEAftTrigCut2->Fill(clsE);
fb87d707 1809 }
1810
1811 if(nClusterTrigCut>0 && ( nCell > (1 + clsE / 3)))
1812 {
1813 cluster->SetChi2(3);
feffe705 1814 //fClsEAftTrigCut3->Fill(clsE);
fb87d707 1815 }
1816
1817 if(nClusterTrigCut>0 && (nCell > (1 + clsE / 3) )&&( clusterTime > fTimeCutLow && clusterTime < fTimeCutHigh ))
1818 {
1819 // cluster->SetChi2(4);
feffe705 1820 //fClsEAftTrigCut4->Fill(clsE);
fb87d707 1821 }
1822 if(nClusterTrigCut<1)
1823 {
1824 cluster->SetChi2(0);
1825
feffe705 1826 //fClsEAftTrigCut->Fill(clsE);
fb87d707 1827 }
1828
1829 } // clusters
1830}
1831
09022877 1832// <-------- only MC correction
8efacc22 1833double AliAnalysisTaskHFECal::MCEopMeanCorrection(double pTmc, float central)
8af58b0d 1834{
8efacc22 1835 TF1 *fcorr0 = new TF1("fcorr0","[0]*tanh([1]+[2]*x)");
1836 TF1 *fcorr1 = new TF1("fcorr1","[0]*tanh([1]+[2]*x)");
1837
8af58b0d 1838 double shift = 0.0;
8efacc22 1839
1840 if(central>0 && central<=10)
1841 {
1842 fcorr0->SetParameters(1.045,1.288,3.18e-01); //
1843 fcorr1->SetParameters(9.91e-01,3.466,2.344);
1844 }
09022877 1845 else if(central>10 && central<=20)
8efacc22 1846 {
1847 fcorr0->SetParameters(1.029,8.254e-01,4.07e-01);
1848 fcorr1->SetParameters(0.975,2.276,1.501e-01);
1849 }
1850 else if(central>20 && central<=30)
1851 {
1852 fcorr0->SetParameters(1.01,8.795e-01,3.904e-01);
1853 fcorr1->SetParameters(9.675e-01,1.654,2.583e-01);
1854 }
1855 else if(central>30 && central<=40)
1856 {
1857 fcorr0->SetParameters(1.00,1.466,2.305e-1);
1858 fcorr1->SetParameters(9.637e-01,1.473,2.754e-01);
1859 }
1860 else if(central>40 && central<=50)
1861 {
fd2126aa 1862 fcorr0->SetParameters(1.00,1.422,1.518e-01);
8efacc22 1863 fcorr1->SetParameters(9.59e-01,1.421,2.931e-01);
1864 }
1865
1866 else if(central>50 && central<=70)
1867 {
1868 fcorr0->SetParameters(0.989,2.495,2.167);
e1f0fb74 1869 fcorr1->SetParameters(0.961,1.734,1.438e-01);
8efacc22 1870 }
1871 else if(central>70 && central<=100)
1872 {
1873 fcorr0->SetParameters(0.981,-3.373,3.93327);
fd2126aa 1874 fcorr1->SetParameters(9.574e-01,1.698,1.58e-01);
8efacc22 1875 }
1876
1877
1878 shift = fcorr0->Eval(pTmc)-fcorr1->Eval(pTmc);
1879
8af58b0d 1880 return shift;
1881}
fb87d707 1882
09022877 1883// <-------- only Data correction
1884double AliAnalysisTaskHFECal::NsigmaCorrection(double tmpeta, float central)
1885{
1886 int icent = 0;
1887
1888 if(central>=0 && central<10)
1889 {
1890 icent = 0;
1891 }
1892 else if(central>=10 && central<20)
1893 {
1894 icent = 1;
1895 }
1896 else if(central>=20 && central<30)
1897 {
1898 icent = 2;
1899 }
1900 else if(central>=30 && central<40)
1901 {
1902 icent = 3;
1903 }
1904 else if(central>=40 && central<50)
1905 {
1906 icent = 4;
1907 }
1908 else if(central>=50 && central<70)
1909 {
1910 icent = 5;
1911 }
1912 else
1913 {
1914 icent = 6;
1915 }
1916
1917 double shift = fnSigEtaCorr[icent]->Eval(tmpeta);
1918
1919 //cout << "eta correction"<< endl;
1920 //cout << "cent = "<< central<< endl;
1921 //cout << "icent = "<< icent << endl;
1922 //cout << "shift = "<< shift << endl;
1923
1924 return shift;
fb87d707 1925
09022877 1926}
fb87d707 1927
1928