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