git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
added AddTask_GammaConvV1.C for Legotrain
[u/mrichter/AliRoot.git] / PWGHF / hfe / AliAnalysisTaskHFECal.cxx
CommitLineData
bfc7c23b 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16// Class for heavy-flavour electron with EMCal triggered events
17// Author: Shingo Sakai
18
19
20#include "TChain.h"
21#include "TTree.h"
22#include "TH2F.h"
23#include "TMath.h"
24#include "TCanvas.h"
25#include "THnSparse.h"
26#include "TLorentzVector.h"
27#include "TString.h"
28#include "TFile.h"
29
30#include "AliAnalysisTask.h"
31#include "AliAnalysisManager.h"
32
33#include "AliESDEvent.h"
34#include "AliESDHandler.h"
35#include "AliAODEvent.h"
36#include "AliAODHandler.h"
37
e7d87aef 38#include "AliMCEventHandler.h"
39#include "AliMCEvent.h"
40#include "AliMCParticle.h"
41
bfc7c23b 42#include "AliAnalysisTaskHFECal.h"
43#include "TGeoGlobalMagField.h"
44#include "AliLog.h"
45#include "AliAnalysisTaskSE.h"
46#include "TRefArray.h"
47#include "TVector.h"
48#include "AliESDInputHandler.h"
49#include "AliESDpid.h"
50#include "AliESDtrackCuts.h"
51#include "AliPhysicsSelection.h"
52#include "AliESDCaloCluster.h"
53#include "AliAODCaloCluster.h"
54#include "AliEMCALRecoUtils.h"
55#include "AliEMCALGeometry.h"
56#include "AliGeomManager.h"
57#include "stdio.h"
58#include "TGeoManager.h"
59#include "iostream"
60#include "fstream"
61
62#include "AliEMCALTrack.h"
63#include "AliMagF.h"
64
65#include "AliKFParticle.h"
66#include "AliKFVertex.h"
67
68#include "AliPID.h"
69#include "AliPIDResponse.h"
70#include "AliHFEcontainer.h"
71#include "AliHFEcuts.h"
72#include "AliHFEpid.h"
73#include "AliHFEpidBase.h"
74#include "AliHFEpidQAmanager.h"
75#include "AliHFEtools.h"
76#include "AliCFContainer.h"
77#include "AliCFManager.h"
78
e7d87aef 79#include "AliStack.h"
80
bfc7c23b 81#include "AliCentrality.h"
82
83ClassImp(AliAnalysisTaskHFECal)
84//________________________________________________________________________
85AliAnalysisTaskHFECal::AliAnalysisTaskHFECal(const char *name)
86 : AliAnalysisTaskSE(name)
87 ,fESD(0)
e7d87aef 88 ,fMC(0)
fb87d707 89 ,fGeom(0)
bfc7c23b 90 ,fOutputList(0)
85985bb0 91 ,fqahist(1)
bfc7c23b 92 ,fTrackCuts(0)
93 ,fCuts(0)
94 ,fIdentifiedAsOutInz(kFALSE)
95 ,fPassTheEventCut(kFALSE)
96 ,fRejectKinkMother(kFALSE)
e7d87aef 97 ,fmcData(kFALSE)
bfc7c23b 98 ,fVz(0.0)
99 ,fCFM(0)
100 ,fPID(0)
101 ,fPIDqa(0)
102 ,fOpeningAngleCut(0.1)
103 ,fInvmassCut(0.01)
104 ,fNoEvents(0)
105 ,fEMCAccE(0)
106 ,fTrkpt(0)
107 ,fTrkEovPBef(0)
108 ,fTrkEovPAft(0)
109 ,fdEdxBef(0)
110 ,fdEdxAft(0)
111 ,fIncpT(0)
fb87d707 112 ,fIncpTM20(0)
bfc7c23b 113 ,fInvmassLS(0)
114 ,fInvmassULS(0)
115 ,fOpeningAngleLS(0)
116 ,fOpeningAngleULS(0)
117 ,fPhotoElecPt(0)
118 ,fPhoElecPt(0)
fb87d707 119 ,fPhoElecPtM20(0)
f09766dd 120 ,fSameElecPt(0)
fb87d707 121 ,fSameElecPtM20(0)
bfc7c23b 122 ,fTrackPtBefTrkCuts(0)
123 ,fTrackPtAftTrkCuts(0)
124 ,fTPCnsigma(0)
125 ,fCent(0)
126 ,fEleInfo(0)
feffe705 127 /*
fb87d707 128 ,fClsEBftTrigCut(0)
129 ,fClsEAftTrigCut(0)
130 ,fClsEAftTrigCut1(0)
131 ,fClsEAftTrigCut2(0)
132 ,fClsEAftTrigCut3(0)
133 ,fClsEAftTrigCut4(0)
134 ,fClsETime(0)
135 ,fClsETime1(0)
136 ,fTrigTimes(0)
42c75692 137 ,fCellCheck(0)
feffe705 138 */
e7d87aef 139 ,fInputHFEMC(0)
feffe705 140 ,fInputAlle(0)
e7d87aef 141 ,fIncpTMChfe(0)
3db00c72 142 ,fIncpTMChfeAll(0)
e7d87aef 143 ,fIncpTMCM20hfe(0)
3db00c72 144 ,fIncpTMCM20hfeAll(0)
e7d87aef 145 ,fIncpTMCpho(0)
146 ,fIncpTMCM20pho(0)
147 ,fPhoElecPtMC(0)
148 ,fPhoElecPtMCM20(0)
149 ,fSameElecPtMC(0)
150 ,fSameElecPtMCM20(0)
a6df418c 151 ,CheckNclust(0)
152 ,CheckNits(0)
d113d7cd 153 ,Hpi0pTcheck(0)
bfc7c23b 154{
155 //Named constructor
156
157 fPID = new AliHFEpid("hfePid");
158 fTrackCuts = new AliESDtrackCuts();
159
160 // Define input and output slots here
161 // Input slot #0 works with a TChain
162 DefineInput(0, TChain::Class());
163 // Output slot #0 id reserved by the base class for AOD
164 // Output slot #1 writes into a TH1 container
165 // DefineOutput(1, TH1I::Class());
166 DefineOutput(1, TList::Class());
167 // DefineOutput(3, TTree::Class());
168}
169
170//________________________________________________________________________
171AliAnalysisTaskHFECal::AliAnalysisTaskHFECal()
172 : AliAnalysisTaskSE("DefaultAnalysis_AliAnalysisTaskHFECal")
173 ,fESD(0)
e7d87aef 174 ,fMC(0)
fb87d707 175 ,fGeom(0)
bfc7c23b 176 ,fOutputList(0)
85985bb0 177 ,fqahist(1)
bfc7c23b 178 ,fTrackCuts(0)
179 ,fCuts(0)
180 ,fIdentifiedAsOutInz(kFALSE)
181 ,fPassTheEventCut(kFALSE)
182 ,fRejectKinkMother(kFALSE)
e7d87aef 183 ,fmcData(kFALSE)
bfc7c23b 184 ,fVz(0.0)
185 ,fCFM(0)
186 ,fPID(0)
187 ,fPIDqa(0)
188 ,fOpeningAngleCut(0.1)
189 ,fInvmassCut(0.01)
190 ,fNoEvents(0)
191 ,fEMCAccE(0)
192 ,fTrkpt(0)
193 ,fTrkEovPBef(0)
194 ,fTrkEovPAft(0)
195 ,fdEdxBef(0)
196 ,fdEdxAft(0)
197 ,fIncpT(0)
fb87d707 198 ,fIncpTM20(0)
bfc7c23b 199 ,fInvmassLS(0)
200 ,fInvmassULS(0)
201 ,fOpeningAngleLS(0)
202 ,fOpeningAngleULS(0)
203 ,fPhotoElecPt(0)
204 ,fPhoElecPt(0)
fb87d707 205 ,fPhoElecPtM20(0)
f09766dd 206 ,fSameElecPt(0)
fb87d707 207 ,fSameElecPtM20(0)
bfc7c23b 208 ,fTrackPtBefTrkCuts(0)
209 ,fTrackPtAftTrkCuts(0)
210 ,fTPCnsigma(0)
211 ,fCent(0)
212 ,fEleInfo(0)
feffe705 213 /*
fb87d707 214 ,fClsEBftTrigCut(0)
215 ,fClsEAftTrigCut(0)
216 ,fClsEAftTrigCut1(0)
217 ,fClsEAftTrigCut2(0)
218 ,fClsEAftTrigCut3(0)
219 ,fClsEAftTrigCut4(0)
220 ,fClsETime(0)
221 ,fClsETime1(0)
222 ,fTrigTimes(0)
42c75692 223 ,fCellCheck(0)
feffe705 224 */
e7d87aef 225 ,fInputHFEMC(0)
feffe705 226 ,fInputAlle(0)
e7d87aef 227 ,fIncpTMChfe(0)
3db00c72 228 ,fIncpTMChfeAll(0)
e7d87aef 229 ,fIncpTMCM20hfe(0)
3db00c72 230 ,fIncpTMCM20hfeAll(0)
e7d87aef 231 ,fIncpTMCpho(0)
232 ,fIncpTMCM20pho(0)
233 ,fPhoElecPtMC(0)
234 ,fPhoElecPtMCM20(0)
235 ,fSameElecPtMC(0)
236 ,fSameElecPtMCM20(0)
a6df418c 237 ,CheckNclust(0)
238 ,CheckNits(0)
d113d7cd 239 ,Hpi0pTcheck(0)
bfc7c23b 240{
241 //Default constructor
242 fPID = new AliHFEpid("hfePid");
243
244 fTrackCuts = new AliESDtrackCuts();
245
246 // Constructor
247 // Define input and output slots here
248 // Input slot #0 works with a TChain
249 DefineInput(0, TChain::Class());
250 // Output slot #0 id reserved by the base class for AOD
251 // Output slot #1 writes into a TH1 container
252 // DefineOutput(1, TH1I::Class());
253 DefineOutput(1, TList::Class());
254 //DefineOutput(3, TTree::Class());
255}
256//_________________________________________
257
258AliAnalysisTaskHFECal::~AliAnalysisTaskHFECal()
259{
260 //Destructor
261
262 delete fOutputList;
fb87d707 263 delete fGeom;
bfc7c23b 264 delete fPID;
a6df418c 265 delete fCuts;
bfc7c23b 266 delete fCFM;
267 delete fPIDqa;
268 delete fTrackCuts;
269}
270//_________________________________________
271
272void AliAnalysisTaskHFECal::UserExec(Option_t*)
273{
274 //Main loop
275 //Called for each event
276
277 // create pointer to event
278 fESD = dynamic_cast<AliESDEvent*>(InputEvent());
279 if (!fESD) {
280 printf("ERROR: fESD not available\n");
281 return;
282 }
283
284 if(!fCuts){
285 AliError("HFE cuts not available");
286 return;
287 }
288
289 if(!fPID->IsInitialized()){
290 // Initialize PID with the given run number
291 AliWarning("PID not initialised, get from Run no");
292 fPID->InitializePID(fESD->GetRunNumber());
293 }
294
e7d87aef 295 if(fmcData)fMC = MCEvent();
296 AliStack* stack = NULL;
297 if(fmcData && fMC)stack = fMC->Stack();
298
299 Float_t cent = -1.;
300 AliCentrality *centrality = fESD->GetCentrality();
301 cent = centrality->GetCentralityPercentile("V0M");
302
303 //---- fill MC track info
304 if(fmcData && fMC)
305 {
306 Int_t nParticles = stack->GetNtrack();
307 for (Int_t iParticle = 0; iParticle < nParticles; iParticle++) {
308 TParticle* particle = stack->Particle(iParticle);
309 int fPDG = particle->GetPdgCode();
310 double mcZvertex = fMC->GetPrimaryVertex()->GetZ();
311 double pTMC = particle->Pt();
feffe705 312 double proR = particle->R();
9b3495ff 313 double etaMC = particle->Eta();
314 if(fabs(etaMC)>0.6)continue;
e7d87aef 315 Bool_t mcInDtoE= kFALSE;
316 Bool_t mcInBtoE= kFALSE;
317
d113d7cd 318 Bool_t MChijing = fMC->IsFromBGEvent(iParticle);
319 if(!MChijing)printf("not MC hijing");
320 int iHijing = 1;
321 if(!MChijing)iHijing = 0;
322 if(fPDG==111)Hpi0pTcheck->Fill(pTMC,iHijing);
323
e7d87aef 324 if(particle->GetFirstMother()>-1 && fabs(fPDG)==11)
325 {
326 int parentPID = stack->Particle(particle->GetFirstMother())->GetPdgCode();
327 if((fabs(parentPID)==411 || fabs(parentPID)==413 || fabs(parentPID)==421 || fabs(parentPID)==423 || fabs(parentPID)==431)&& fabs(fPDG)==11)mcInDtoE = kTRUE;
328 if((fabs(parentPID)==511 || fabs(parentPID)==513 || fabs(parentPID)==521 || fabs(parentPID)==523 || fabs(parentPID)==531)&& fabs(fPDG)==11)mcInBtoE = kTRUE;
329 if((mcInBtoE || mcInDtoE) && fabs(mcZvertex)<10.0)fInputHFEMC->Fill(cent,pTMC);
d113d7cd 330
331
e7d87aef 332 }
333
feffe705 334 if(proR<7.0 && fabs(fPDG)==11)fInputAlle->Fill(cent,pTMC);
335
e7d87aef 336 }
337 }
338
bfc7c23b 339 fNoEvents->Fill(0);
340
341 Int_t fNOtrks = fESD->GetNumberOfTracks();
342 const AliESDVertex *pVtx = fESD->GetPrimaryVertex();
343
344 Double_t pVtxZ = -999;
345 pVtxZ = pVtx->GetZ();
346
347 if(TMath::Abs(pVtxZ)>10) return;
348 fNoEvents->Fill(1);
349
350 if(fNOtrks<2) return;
351 fNoEvents->Fill(2);
352
353 AliPIDResponse *pidResponse = fInputHandler->GetPIDResponse();
354 if(!pidResponse){
355 AliDebug(1, "Using default PID Response");
356 pidResponse = AliHFEtools::GetDefaultPID(kFALSE, fInputEvent->IsA() == AliAODEvent::Class());
357 }
358
359 fPID->SetPIDResponse(pidResponse);
360
361 fCFM->SetRecEventInfo(fESD);
362
e7d87aef 363 //Float_t cent = -1.;
364 //AliCentrality *centrality = fESD->GetCentrality();
365 //cent = centrality->GetCentralityPercentile("V0M");
bfc7c23b 366 fCent->Fill(cent);
367
b12bc641 368 //if(cent>90.) return;
bfc7c23b 369
370 // Calorimeter info.
371
42c75692 372 FindTriggerClusters();
373
bfc7c23b 374 // make EMCAL array
375 for(Int_t iCluster=0; iCluster<fESD->GetNumberOfCaloClusters(); iCluster++)
376 {
377 AliESDCaloCluster *clust = fESD->GetCaloCluster(iCluster);
378 if(clust->IsEMCAL())
379 {
380 double clustE = clust->E();
381 float emcx[3]; // cluster pos
382 clust->GetPosition(emcx);
383 TVector3 clustpos(emcx[0],emcx[1],emcx[2]);
384 double emcphi = clustpos.Phi();
385 double emceta = clustpos.Eta();
fb87d707 386 double calInfo[5];
387 calInfo[0] = emcphi; calInfo[1] = emceta; calInfo[2] = clustE; calInfo[3] = cent; calInfo[4] = clust->Chi2();
85985bb0 388 //if(clustE>1.5)fEMCAccE->Fill(calInfo);
44be9e1d 389 //if(fqahist==1 && clustE>1.5)fEMCAccE->Fill(calInfo);
bfc7c23b 390 }
391 }
392
393 // Track loop
394 for (Int_t iTracks = 0; iTracks < fESD->GetNumberOfTracks(); iTracks++) {
395 AliESDtrack* track = fESD->GetTrack(iTracks);
396 if (!track) {
397 printf("ERROR: Could not receive track %d\n", iTracks);
398 continue;
399 }
e7d87aef 400
401 Bool_t mcPho = kFALSE;
402 Bool_t mcDtoE= kFALSE;
403 Bool_t mcBtoE= kFALSE;
404 double mcele = -1.;
44be9e1d 405 double mcpT = 0.0;
b12bc641 406 double mcMompT = 0.0;
d113d7cd 407
408 int iHijing = 1;
409
e7d87aef 410 if(fmcData && fMC && stack)
411 {
feffe705 412 Int_t label = TMath::Abs(track->GetLabel());
d113d7cd 413
414 Bool_t MChijing = fMC->IsFromBGEvent(label);
415 if(!MChijing)iHijing = 0;
416
e7d87aef 417 TParticle* particle = stack->Particle(label);
418 int mcpid = particle->GetPdgCode();
44be9e1d 419 mcpT = particle->Pt();
9b3495ff 420 //printf("MCpid = %d",mcpid);
e7d87aef 421 if(particle->GetFirstMother()>-1)
422 {
423 int parentPID = stack->Particle(particle->GetFirstMother())->GetPdgCode();
b12bc641 424 mcMompT = stack->Particle(particle->GetFirstMother())->Pt();
e7d87aef 425 if((parentPID==22 || parentPID==111 || parentPID==221)&& fabs(mcpid)==11)mcPho = kTRUE;
426 if((fabs(parentPID)==411 || fabs(parentPID)==413 || fabs(parentPID)==421 || fabs(parentPID)==423 || fabs(parentPID)==431)&& fabs(mcpid)==11)mcDtoE = kTRUE;
427 if((fabs(parentPID)==511 || fabs(parentPID)==513 || fabs(parentPID)==521 || fabs(parentPID)==523 || fabs(parentPID)==531)&& fabs(mcpid)==11)mcBtoE = kTRUE;
428 }
9b3495ff 429 if(fabs(mcpid)==11 && mcDtoE)mcele= 1.;
430 if(fabs(mcpid)==11 && mcBtoE)mcele= 2.;
431 if(fabs(mcpid)==11 && mcPho)mcele= 3.;
e7d87aef 432 }
433
9b3495ff 434 if(TMath::Abs(track->Eta())>0.6) continue;
f09766dd 435 if(TMath::Abs(track->Pt()<2.0)) continue;
bfc7c23b 436
437 fTrackPtBefTrkCuts->Fill(track->Pt());
438 // RecKine: ITSTPC cuts
439 if(!ProcessCutStep(AliHFEcuts::kStepRecKineITSTPC, track)) continue;
440
441 //RecKink
442 if(fRejectKinkMother) { // Quick and dirty fix to reject both kink mothers and daughters
443 if(track->GetKinkIndex(0) != 0) continue;
444 }
445
446 // RecPrim
447 if(!ProcessCutStep(AliHFEcuts::kStepRecPrim, track)) continue;
448
449 // HFEcuts: ITS layers cuts
450 if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsITS, track)) continue;
451
f09766dd 452 // HFE cuts: TPC PID cleanup
453 if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTPC, track)) continue;
454
a6df418c 455 int nTPCcl = track->GetTPCNcls();
6e92c428 456 int nTPCclF = track->GetTPCNclsF();
a6df418c 457 int nITS = track->GetNcls(0);
bfc7c23b 458
459 fTrackPtAftTrkCuts->Fill(track->Pt());
460
461 Double_t mom = -999., eop=-999., pt = -999., dEdx=-999., fTPCnSigma=-10, phi=-999., eta=-999.;
feffe705 462 pt = track->Pt();
463 if(pt<2.0)continue;
bfc7c23b 464
465 // Track extrapolation
466
bfc7c23b 467 Int_t charge = track->Charge();
468 fTrkpt->Fill(pt);
469 mom = track->P();
470 phi = track->Phi();
471 eta = track->Eta();
472 dEdx = track->GetTPCsignal();
473 fTPCnSigma = fPID->GetPIDResponse() ? fPID->GetPIDResponse()->NumberOfSigmasTPC(track, AliPID::kElectron) : 1000;
474
475 double ncells = -1.0;
476 double m20 = -1.0;
477 double m02 = -1.0;
478 double disp = -1.0;
479 double rmatch = -1.0;
480 double nmatch = -1.0;
f09766dd 481 double oppstatus = 0.0;
bfc7c23b 482
f09766dd 483 Bool_t fFlagPhotonicElec = kFALSE;
484 Bool_t fFlagConvinatElec = kFALSE;
bfc7c23b 485
486 Int_t clsId = track->GetEMCALcluster();
487 if (clsId>0){
488 AliESDCaloCluster *clust = fESD->GetCaloCluster(clsId);
489 if(clust && clust->IsEMCAL()){
490
491 double clustE = clust->E();
492 eop = clustE/fabs(mom);
493 //double clustT = clust->GetTOF();
494 ncells = clust->GetNCells();
495 m02 = clust->GetM02();
496 m20 = clust->GetM20();
497 disp = clust->GetDispersion();
498 double delphi = clust->GetTrackDx();
499 double deleta = clust->GetTrackDz();
500 rmatch = sqrt(pow(delphi,2)+pow(deleta,2));
501 nmatch = clust->GetNTracksMatched();
502
c2801a73 503 if(fTPCnSigma>-1.5 && fTPCnSigma<3.0)
504 {
505 SelectPhotonicElectron(iTracks,cent,track,fFlagPhotonicElec,fFlagConvinatElec,fTPCnSigma,m20,eop,mcele);
506 }
507 if(fFlagPhotonicElec)oppstatus = 1.0;
508 if(fFlagConvinatElec)oppstatus = 2.0;
509 if(fFlagPhotonicElec && fFlagConvinatElec)oppstatus = 3.0;
510
feffe705 511 double valdedx[16];
a6df418c 512 valdedx[0] = pt; valdedx[1] = nITS; valdedx[2] = phi; valdedx[3] = eta; valdedx[4] = fTPCnSigma;
6e92c428 513 valdedx[5] = eop; valdedx[6] = rmatch; valdedx[7] = ncells, valdedx[8] = nTPCclF; valdedx[9] = m20; valdedx[10] = mcpT;
a6df418c 514 valdedx[11] = cent; valdedx[12] = charge; valdedx[13] = oppstatus; valdedx[14] = nTPCcl;
feffe705 515 valdedx[15] = mcele;
85985bb0 516 if(fqahist==1)fEleInfo->Fill(valdedx);
bfc7c23b 517
518
519 }
520 }
a6df418c 521
522 if(nITS<2.5)continue;
523 if(nTPCcl<100)continue;
524
525 CheckNclust->Fill(nTPCcl);
526 CheckNits->Fill(nITS);
527
42c75692 528 fdEdxBef->Fill(mom,fTPCnSigma);
bfc7c23b 529 fTPCnsigma->Fill(mom,fTPCnSigma);
f09766dd 530 if(fTPCnSigma >= -1.0 && fTPCnSigma <= 3)fTrkEovPBef->Fill(pt,eop);
bfc7c23b 531
bfc7c23b 532 Int_t pidpassed = 1;
533
534
535 //--- track accepted
536 AliHFEpidObject hfetrack;
537 hfetrack.SetAnalysisType(AliHFEpidObject::kESDanalysis);
538 hfetrack.SetRecTrack(track);
3db00c72 539 double binct = 10.5;
540 if((0.0< cent) && (cent<5.0)) binct = 0.5;
541 if((5.0< cent) && (cent<10.0)) binct = 1.5;
542 if((10.0< cent) && (cent<20.0)) binct = 2.5;
543 if((20.0< cent) && (cent<30.0)) binct = 3.5;
544 if((30.0< cent) && (cent<40.0)) binct = 4.5;
545 if((40.0< cent) && (cent<50.0)) binct = 5.5;
546 if((50.0< cent) && (cent<60.0)) binct = 6.5;
547 if((60.0< cent) && (cent<70.0)) binct = 7.5;
548 if((70.0< cent) && (cent<80.0)) binct = 8.5;
549 if((80.0< cent) && (cent<90.0)) binct = 9.5;
550 if((90.0< cent) && (cent<100.0)) binct = 10.5;
551
552 hfetrack.SetCentrality((int)binct); //added
bfc7c23b 553 hfetrack.SetPbPb();
554 if(!fPID->IsSelected(&hfetrack, NULL, "", fPIDqa)) pidpassed = 0;
555
556 if(pidpassed==0) continue;
557
558 fTrkEovPAft->Fill(pt,eop);
42c75692 559 fdEdxAft->Fill(mom,fTPCnSigma);
bfc7c23b 560
e7d87aef 561 // Fill real data
fb87d707 562 fIncpT->Fill(cent,pt);
bfc7c23b 563 if(fFlagPhotonicElec) fPhoElecPt->Fill(cent,pt);
f09766dd 564 if(fFlagConvinatElec) fSameElecPt->Fill(cent,pt);
fb87d707 565
566 if(m20>0.0 && m20<0.3)
567 {
568 fIncpTM20->Fill(cent,pt);
569 if(fFlagPhotonicElec) fPhoElecPtM20->Fill(cent,pt);
570 if(fFlagConvinatElec) fSameElecPtM20->Fill(cent,pt);
571 }
feffe705 572
e7d87aef 573 // MC
2f37d8dc 574 if(mcele>0)
e7d87aef 575 {
3db00c72 576
577 fIncpTMChfeAll->Fill(cent,pt);
578 if(m20>0.0 && m20<0.3)fIncpTMCM20hfeAll->Fill(cent,pt);
579
e7d87aef 580 if(mcBtoE || mcDtoE)
581 {
582 fIncpTMChfe->Fill(cent,pt);
583 if(m20>0.0 && m20<0.3)fIncpTMCM20hfe->Fill(cent,pt);
584 }
585 if(mcPho)
586 {
d113d7cd 587 double phoval[4];
44be9e1d 588 phoval[0] = cent;
589 phoval[1] = pt;
590 phoval[2] = fTPCnSigma;
d113d7cd 591 phoval[3] = iHijing;
44be9e1d 592
593 fIncpTMCpho->Fill(phoval);
594 if(fFlagPhotonicElec) fPhoElecPtMC->Fill(phoval);
595 if(fFlagConvinatElec) fSameElecPtMC->Fill(phoval);
e7d87aef 596
597 if(m20>0.0 && m20<0.3)
598 {
44be9e1d 599 fIncpTMCM20pho->Fill(phoval);
600 if(fFlagPhotonicElec) fPhoElecPtMCM20->Fill(phoval);
601 if(fFlagConvinatElec) fSameElecPtMCM20->Fill(phoval);
e7d87aef 602 }
603 }
604 }
bfc7c23b 605 }
606 PostData(1, fOutputList);
607}
608//_________________________________________
609void AliAnalysisTaskHFECal::UserCreateOutputObjects()
610{
42c75692 611 //--- Check MC
612
e7d87aef 613 //Bool_t mcData = kFALSE;
42c75692 614 if(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler())
615 {
e7d87aef 616 fmcData = kTRUE;
42c75692 617 printf("+++++ MC Data available");
618 }
e7d87aef 619 if(fmcData)
42c75692 620 {
621 printf("++++++++= MC analysis \n");
622 }
623 else
624 {
625 printf("++++++++ real data analysis \n");
626 }
627
85985bb0 628 printf("+++++++ QA hist %d",fqahist);
629
fb87d707 630 //---- Geometry
631 fGeom = AliEMCALGeometry::GetInstance("EMCAL_COMPLETEV1");
632
bfc7c23b 633 //--------Initialize PID
42c75692 634 //fPID->SetHasMCData(kFALSE);
e7d87aef 635 fPID->SetHasMCData(fmcData);
bfc7c23b 636 if(!fPID->GetNumberOfPIDdetectors())
637 {
638 fPID->AddDetector("TPC", 0);
639 fPID->AddDetector("EMCAL", 1);
640 }
641
3db00c72 642 Double_t params[4];
78ff1095 643 const char *cutmodel;
3db00c72 644 cutmodel = "pol0";
645 params[0] = -1.0; //sigma min
9b3495ff 646 double maxnSig = 3.0;
647 if(fmcData)
648 {
649 params[0] = -5.0; //sigma min
650 maxnSig = 5.0;
651 }
652
653 for(Int_t a=0;a<11;a++)fPID->ConfigureTPCcentralityCut(a,cutmodel,params,maxnSig);
3db00c72 654
bfc7c23b 655 fPID->SortDetectors();
656 fPIDqa = new AliHFEpidQAmanager();
657 fPIDqa->Initialize(fPID);
a6df418c 658
659 //------- fcut --------------
660 fCuts = new AliHFEcuts();
661 fCuts->CreateStandardCuts();
662 //fCuts->SetMinNClustersTPC(100);
663 fCuts->SetMinNClustersTPC(90);
664 fCuts->SetMinRatioTPCclusters(0.6);
665 fCuts->SetTPCmodes(AliHFEextraCuts::kFound, AliHFEextraCuts::kFoundOverFindable);
666 //fCuts->SetMinNClustersITS(3);
667 fCuts->SetMinNClustersITS(2);
668 fCuts->SetCutITSpixel(AliHFEextraCuts::kAny);
669 fCuts->SetCheckITSLayerStatus(kFALSE);
670 fCuts->SetVertexRange(10.);
671 fCuts->SetTOFPIDStep(kFALSE);
672 fCuts->SetPtRange(2, 50);
673 fCuts->SetMaxImpactParam(3.,3.);
674
bfc7c23b 675 //--------Initialize correction Framework and Cuts
676 fCFM = new AliCFManager;
677 const Int_t kNcutSteps = AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kNcutStepsRecTrack + AliHFEcuts::kNcutStepsDETrack;
678 fCFM->SetNStepParticle(kNcutSteps);
679 for(Int_t istep = 0; istep < kNcutSteps; istep++)
680 fCFM->SetParticleCutsList(istep, NULL);
681
682 if(!fCuts){
683 AliWarning("Cuts not available. Default cuts will be used");
684 fCuts = new AliHFEcuts;
685 fCuts->CreateStandardCuts();
686 }
687 fCuts->Initialize(fCFM);
688
689 //---------Output Tlist
690 fOutputList = new TList();
691 fOutputList->SetOwner();
692 fOutputList->Add(fPIDqa->MakeList("PIDQA"));
693
694 fNoEvents = new TH1F("fNoEvents","",4,-0.5,3.5) ;
695 fOutputList->Add(fNoEvents);
696
02cd82a0 697 Int_t binsE[5] = {250, 100, 1000, 200, 10};
fb87d707 698 Double_t xminE[5] = {1.0, -1, 0.0, 0, -0.5};
699 Double_t xmaxE[5] = {3.5, 1, 100.0, 100, 9.5};
700 fEMCAccE = new THnSparseD("fEMCAccE","EMC acceptance & E;#eta;#phi;Energy;Centrality;trugCondition;",5,binsE,xminE,xmaxE);
85985bb0 701 if(fqahist==1)fOutputList->Add(fEMCAccE);
bfc7c23b 702
703 fTrkpt = new TH1F("fTrkpt","track pt",100,0,50);
704 fOutputList->Add(fTrkpt);
705
706 fTrackPtBefTrkCuts = new TH1F("fTrackPtBefTrkCuts","track pt before track cuts",100,0,50);
707 fOutputList->Add(fTrackPtBefTrkCuts);
708
709 fTrackPtAftTrkCuts = new TH1F("fTrackPtAftTrkCuts","track pt after track cuts",100,0,50);
710 fOutputList->Add(fTrackPtAftTrkCuts);
711
712 fTPCnsigma = new TH2F("fTPCnsigma", "TPC - n sigma",100,0,50,200,-10,10);
713 fOutputList->Add(fTPCnsigma);
714
715 fTrkEovPBef = new TH2F("fTrkEovPBef","track E/p before HFE pid",100,0,50,100,0,2);
716 fOutputList->Add(fTrkEovPBef);
717
718 fTrkEovPAft = new TH2F("fTrkEovPAft","track E/p after HFE pid",100,0,50,100,0,2);
719 fOutputList->Add(fTrkEovPAft);
720
42c75692 721 fdEdxBef = new TH2F("fdEdxBef","track dEdx vs p before HFE pid",100,0,50,200,-10,10);
bfc7c23b 722 fOutputList->Add(fdEdxBef);
723
42c75692 724 fdEdxAft = new TH2F("fdEdxAft","track dEdx vs p after HFE pid",100,0,50,200,-10,10);
bfc7c23b 725 fOutputList->Add(fdEdxAft);
726
02cd82a0 727 fIncpT = new TH2F("fIncpT","HFE pid electro vs. centrality",200,0,100,100,0,50);
bfc7c23b 728 fOutputList->Add(fIncpT);
729
02cd82a0 730 fIncpTM20 = new TH2F("fIncpTM20","HFE pid electro vs. centrality with M20",200,0,100,100,0,50);
fb87d707 731 fOutputList->Add(fIncpTM20);
4e01b68c 732
02cd82a0 733 Int_t nBinspho[9] = { 200, 100, 500, 12, 50, 4, 200, 8, 100};
b12bc641 734 Double_t minpho[9] = { 0., 0., 0., -2.5, 0, -0.5, 0,-1.5, 0};
735 Double_t maxpho[9] = {100., 50., 0.5, 3.5, 1, 3.5, 2, 6.5, 50};
f09766dd 736
b12bc641 737 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);
bfc7c23b 738 fOutputList->Add(fInvmassLS);
739
b12bc641 740 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);
bfc7c23b 741 fOutputList->Add(fInvmassULS);
742
743 fOpeningAngleLS = new TH1F("fOpeningAngleLS","Opening angle for LS pairs",100,0,1);
744 fOutputList->Add(fOpeningAngleLS);
745
746 fOpeningAngleULS = new TH1F("fOpeningAngleULS","Opening angle for ULS pairs",100,0,1);
747 fOutputList->Add(fOpeningAngleULS);
748
749 fPhotoElecPt = new TH1F("fPhotoElecPt", "photonic electron pt",100,0,50);
750 fOutputList->Add(fPhotoElecPt);
751
02cd82a0 752 fPhoElecPt = new TH2F("fPhoElecPt", "Pho-inclusive electron pt",200,0,100,100,0,50);
bfc7c23b 753 fOutputList->Add(fPhoElecPt);
754
02cd82a0 755 fPhoElecPtM20 = new TH2F("fPhoElecPtM20", "Pho-inclusive electron pt with M20",200,0,100,100,0,50);
fb87d707 756 fOutputList->Add(fPhoElecPtM20);
757
02cd82a0 758 fSameElecPt = new TH2F("fSameElecPt", "Same-inclusive electron pt",200,0,100,100,0,50);
f09766dd 759 fOutputList->Add(fSameElecPt);
760
02cd82a0 761 fSameElecPtM20 = new TH2F("fSameElecPtM20", "Same-inclusive electron pt with M20",200,0,100,100,0,50);
fb87d707 762 fOutputList->Add(fSameElecPtM20);
763
02cd82a0 764 fCent = new TH1F("fCent","Centrality",200,0,100) ;
bfc7c23b 765 fOutputList->Add(fCent);
766
767 // Make common binning
44be9e1d 768 const Double_t kMinP = 0.;
bfc7c23b 769 const Double_t kMaxP = 50.;
a6df418c 770 //const Double_t kTPCSigMim = 40.;
771 //const Double_t kTPCSigMax = 140.;
bfc7c23b 772
773 // 1st histogram: TPC dEdx with/without EMCAL (p, pT, TPC Signal, phi, eta, Sig, e/p, ,match, cell, M02, M20, Disp, Centrality, select)
a6df418c 774 Int_t nBins[16] = { 250, 10, 60, 20, 100, 300, 50, 40, 200, 200, 250, 200, 3, 5, 100, 8};
775 Double_t min[16] = {kMinP, -0.5, 1.0, -1.0, -6.0, 0, 0, 0, 0.0, 0.0, 0.0, 0, -1.5, -0.5, 80, -1.5};
6e92c428 776 Double_t max[16] = {kMaxP, 9.5, 4.0, 1.0, 4.0, 3.0, 0.1, 40, 200, 2.0, 50.0, 100, 1.5, 4.5, 180, 6.5};
777 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 778 if(fqahist==1)fOutputList->Add(fEleInfo);
bfc7c23b 779
fb87d707 780 //<--- Trigger info
feffe705 781 /*
fb87d707 782 fClsEBftTrigCut = new TH1F("fClsEBftTrigCut","cluster E before trigger selection",1000,0,100);
783 fOutputList->Add(fClsEBftTrigCut);
784
785 fClsEAftTrigCut = new TH1F("fClsEAftTrigCut","cluster E if cls has 0 trigcut channel",1000,0,100);
786 fOutputList->Add(fClsEAftTrigCut);
787
788 fClsEAftTrigCut1 = new TH1F("fClsEAftTrigCut1","cluster E if cls with trig channel",1000,0,100);
789 fOutputList->Add(fClsEAftTrigCut1);
790
791 fClsEAftTrigCut2 = new TH1F("fClsEAftTrigCut2","cluster E if cls with trigcut channel",1000,0,100);
792 fOutputList->Add(fClsEAftTrigCut2);
793
794 fClsEAftTrigCut3 = new TH1F("fClsEAftTrigCut3","cluster E if cls with trigcut channel + nCell>Ecorrect",1000,0,100);
795 fOutputList->Add(fClsEAftTrigCut3);
796
797 fClsEAftTrigCut4 = new TH1F("fClsEAftTrigCut4","cluster E if cls with trigcut channel + nCell>Ecorrect + cls time cut",1000,0,100);
798 fOutputList->Add(fClsEAftTrigCut4);
799
800 fClsETime = new TH2F("fClsETime", "Cls time vs E; E; time;",1000,0,100,1000,-0.0000002,0.0000009);
801 fOutputList->Add(fClsETime);
802
803 fClsETime1 = new TH2F("fClsETime1", "Cls time vs E if cls contains trigger channel; E; time;",1000,0,100,1000,-0.0000002,0.0000009);
804 fOutputList->Add(fClsETime1);
805
806 fTrigTimes = new TH1F("fTrigTimes", "Trigger time; time; N;",25,0,25);
807 fOutputList->Add(fTrigTimes);
808
42c75692 809 fCellCheck = new TH2F("fCellCheck", "Cell vs E; E GeV; Cell ID",10,6,26,12000,0,12000);
810 fOutputList->Add(fCellCheck);
feffe705 811 */
e7d87aef 812 //<---------- MC
813
02cd82a0 814 fInputHFEMC = new TH2F("fInputHFEMC","Input MC HFE pid electro vs. centrality",200,0,100,100,0,50);
e7d87aef 815 fOutputList->Add(fInputHFEMC);
816
02cd82a0 817 fInputAlle = new TH2F("fInputAlle","Input MC electro vs. centrality",200,0,100,100,0,50);
feffe705 818 fOutputList->Add(fInputAlle);
819
02cd82a0 820 fIncpTMChfe = new TH2F("fIncpTMChfe","MC HFE pid electro vs. centrality",200,0,100,100,0,50);
e7d87aef 821 fOutputList->Add(fIncpTMChfe);
822
02cd82a0 823 fIncpTMChfeAll = new TH2F("fIncpTMChfeAll","MC Alle pid electro vs. centrality",200,0,100,100,0,50);
3db00c72 824 fOutputList->Add(fIncpTMChfeAll);
825
02cd82a0 826 fIncpTMCM20hfe = new TH2F("fIncpTMCM20hfe","MC HFE pid electro vs. centrality with M20",200,0,100,100,0,50);
e7d87aef 827 fOutputList->Add(fIncpTMCM20hfe);
828
02cd82a0 829 fIncpTMCM20hfeAll = new TH2F("fIncpTMCM20hfeAll","MC Alle pid electro vs. centrality with M20",200,0,100,100,0,50);
3db00c72 830 fOutputList->Add(fIncpTMCM20hfeAll);
831
44be9e1d 832
d113d7cd 833 Int_t nBinspho2[4] = { 200, 100, 7, 3};
834 Double_t minpho2[4] = { 0., 0., -2.5, -0.5};
835 Double_t maxpho2[4] = {100., 50., 4.5, 2.5};
44be9e1d 836
d113d7cd 837 fIncpTMCpho = new THnSparseD("fIncpTMCpho","MC Pho pid electro vs. centrality",4,nBinspho2,minpho2,maxpho2);
e7d87aef 838 fOutputList->Add(fIncpTMCpho);
839
d113d7cd 840 fIncpTMCM20pho = new THnSparseD("fIncpTMCM20pho","MC Pho pid electro vs. centrality with M20",4,nBinspho2,minpho2,maxpho2);
e7d87aef 841 fOutputList->Add(fIncpTMCM20pho);
842
d113d7cd 843 fPhoElecPtMC = new THnSparseD("fPhoElecPtMC", "MC Pho-inclusive electron pt",4,nBinspho2,minpho2,maxpho2);
e7d87aef 844 fOutputList->Add(fPhoElecPtMC);
845
d113d7cd 846 fPhoElecPtMCM20 = new THnSparseD("fPhoElecPtMCM20", "MC Pho-inclusive electron pt with M20",4,nBinspho2,minpho2,maxpho2);
feffe705 847 fOutputList->Add(fPhoElecPtMCM20);
e7d87aef 848
d113d7cd 849 fSameElecPtMC = new THnSparseD("fSameElecPtMC", "MC Same-inclusive electron pt",4,nBinspho2,minpho2,maxpho2);
e7d87aef 850 fOutputList->Add(fSameElecPtMC);
851
d113d7cd 852 fSameElecPtMCM20 = new THnSparseD("fSameElecPtMCM20", "MC Same-inclusive electron pt with M20",4,nBinspho2,minpho2,maxpho2);
e7d87aef 853 fOutputList->Add(fSameElecPtMCM20);
854
a6df418c 855 CheckNclust = new TH1D("CheckNclust","cluster check",200,0,200);
856 fOutputList->Add(CheckNclust);
857
858 CheckNits = new TH1D("CheckNits","ITS cluster check",8,-0.5,7.5);
859 fOutputList->Add(CheckNits);
e7d87aef 860
d113d7cd 861 Hpi0pTcheck = new TH2D("Hpi0pTcheck","Pi0 pT from Hijing",100,0,50,3,-0.5,2.5);
862 fOutputList->Add(Hpi0pTcheck);
863
864
bfc7c23b 865 PostData(1,fOutputList);
866}
867
868//________________________________________________________________________
869void AliAnalysisTaskHFECal::Terminate(Option_t *)
870{
871 // Info("Terminate");
872 AliAnalysisTaskSE::Terminate();
873}
874
875//________________________________________________________________________
876Bool_t AliAnalysisTaskHFECal::ProcessCutStep(Int_t cutStep, AliVParticle *track)
877{
878 // Check single track cuts for a given cut step
879 const Int_t kMCOffset = AliHFEcuts::kNcutStepsMCTrack;
880 if(!fCFM->CheckParticleCuts(cutStep + kMCOffset, track)) return kFALSE;
881 return kTRUE;
882}
883//_________________________________________
f09766dd 884//void AliAnalysisTaskHFECal::SelectPhotonicElectron(Int_t itrack, Double_t cent, AliESDtrack *track, Bool_t &fFlagPhotonicElec)
c2801a73 885//void AliAnalysisTaskHFECal::SelectPhotonicElectron(Int_t itrack, Double_t cent, AliESDtrack *track, Bool_t &fFlagPhotonicElec, Bool_t &fFlagConvinatElec, Double_t nSig)
886void 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)
bfc7c23b 887{
888 //Identify non-heavy flavour electrons using Invariant mass method
889
fb87d707 890 fTrackCuts->SetAcceptKinkDaughters(kFALSE);
feffe705 891 fTrackCuts->SetRequireTPCRefit(kTRUE);
85985bb0 892 fTrackCuts->SetRequireITSRefit(kTRUE);
feffe705 893 fTrackCuts->SetEtaRange(-0.9,0.9);
f09766dd 894 //fTrackCuts->SetRequireSigmaToVertex(kTRUE);
bfc7c23b 895 fTrackCuts->SetMaxChi2PerClusterTPC(3.5);
feffe705 896 fTrackCuts->SetMinNClustersTPC(90);
bfc7c23b 897
898 const AliESDVertex *pVtx = fESD->GetPrimaryVertex();
899
900 Bool_t flagPhotonicElec = kFALSE;
f09766dd 901 Bool_t flagConvinatElec = kFALSE;
bfc7c23b 902
5e65985c 903 //for(Int_t jTracks = itrack+1; jTracks<fESD->GetNumberOfTracks(); jTracks++){
904 for(Int_t jTracks = 0; jTracks<fESD->GetNumberOfTracks(); jTracks++){
bfc7c23b 905 AliESDtrack* trackAsso = fESD->GetTrack(jTracks);
906 if (!trackAsso) {
907 printf("ERROR: Could not receive track %d\n", jTracks);
908 continue;
909 }
5e65985c 910 if(itrack==jTracks)continue;
911
f09766dd 912 Double_t dEdxAsso = -999., ptPrim=-999., ptAsso=-999., openingAngle = -999.;
bfc7c23b 913 Double_t mass=999., width = -999;
914 Bool_t fFlagLS=kFALSE, fFlagULS=kFALSE;
915
f09766dd 916 ptPrim = track->Pt();
917
bfc7c23b 918 dEdxAsso = trackAsso->GetTPCsignal();
919 ptAsso = trackAsso->Pt();
920 Int_t chargeAsso = trackAsso->Charge();
921 Int_t charge = track->Charge();
922
b12bc641 923
3db00c72 924 if(ptAsso <0.5) continue;
bfc7c23b 925 if(!fTrackCuts->AcceptTrack(trackAsso)) continue;
f09766dd 926 if(dEdxAsso <65 || dEdxAsso>100) continue; //11a pass1
bfc7c23b 927
928 Int_t fPDGe1 = 11; Int_t fPDGe2 = 11;
929 if(charge>0) fPDGe1 = -11;
930 if(chargeAsso>0) fPDGe2 = -11;
b12bc641 931
932 //printf("chargeAsso = %d\n",chargeAsso);
933 //printf("charge = %d\n",charge);
bfc7c23b 934 if(charge == chargeAsso) fFlagLS = kTRUE;
935 if(charge != chargeAsso) fFlagULS = kTRUE;
936
b12bc641 937 //printf("fFlagLS = %d\n",fFlagLS);
938 //printf("fFlagULS = %d\n",fFlagULS);
939 //printf("\n");
940
bfc7c23b 941 AliKFParticle ge1(*track, fPDGe1);
942 AliKFParticle ge2(*trackAsso, fPDGe2);
943 AliKFParticle recg(ge1, ge2);
944
945 if(recg.GetNDF()<1) continue;
946 Double_t chi2recg = recg.GetChi2()/recg.GetNDF();
947 if(TMath::Sqrt(TMath::Abs(chi2recg))>3.) continue;
948
949 AliKFVertex primV(*pVtx);
950 primV += recg;
951 recg.SetProductionVertex(primV);
952
d113d7cd 953 recg.SetMassConstraint(0,0.0001);
bfc7c23b 954
955 openingAngle = ge1.GetAngle(ge2);
956 if(fFlagLS) fOpeningAngleLS->Fill(openingAngle);
957 if(fFlagULS) fOpeningAngleULS->Fill(openingAngle);
958
bfc7c23b 959
960 recg.GetMass(mass,width);
961
c2801a73 962 double ishower = 0;
963 if(shower>0.0 && shower<0.3)ishower = 1;
964
b12bc641 965 double phoinfo[9];
f09766dd 966 phoinfo[0] = cent;
967 phoinfo[1] = ptPrim;
968 phoinfo[2] = mass;
3db00c72 969 phoinfo[3] = nSig;
4e01b68c 970 phoinfo[4] = openingAngle;
c2801a73 971 phoinfo[5] = ishower;
972 phoinfo[6] = ep;
973 phoinfo[7] = mce;
b12bc641 974 phoinfo[8] = ptAsso;
f09766dd 975
976 if(fFlagLS) fInvmassLS->Fill(phoinfo);
977 if(fFlagULS) fInvmassULS->Fill(phoinfo);
ffa14dda 978
b12bc641 979 //printf("fInvmassCut %f\n",fInvmassCut);
980 //printf("openingAngle %f\n",fOpeningAngleCut);
981
ffa14dda 982 if(openingAngle > fOpeningAngleCut) continue;
983
bfc7c23b 984 if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec){
985 flagPhotonicElec = kTRUE;
986 }
f09766dd 987 if(mass<fInvmassCut && fFlagLS && !flagConvinatElec){
988 flagConvinatElec = kTRUE;
989 }
bfc7c23b 990
991 }
992 fFlagPhotonicElec = flagPhotonicElec;
f09766dd 993 fFlagConvinatElec = flagConvinatElec;
bfc7c23b 994
995}
996
fb87d707 997
998//_________________________________________
999void AliAnalysisTaskHFECal::FindTriggerClusters()
1000{
1001 // constants
1002 const int nModuleCols = 2;
1003 const int nModuleRows = 5;
1004 const int nColsFeeModule = 48;
1005 const int nRowsFeeModule = 24;
1006 const int nColsFaltroModule = 24;
1007 const int nRowsFaltroModule = 12;
1008 //const int faltroWidthMax = 20;
1009
1010 // part 1, trigger extraction -------------------------------------
1011 Int_t globCol, globRow;
1012 //Int_t ntimes=0, nTrigChannel=0, nTrigChannelCut=0, trigInCut=0;
1013 Int_t ntimes=0, nTrigChannel=0, nTrigChannelCut=0;
1014
1015 //Int_t trigtimes[faltroWidthMax];
1016 Double_t cellTime[nColsFeeModule*nModuleCols][nRowsFeeModule*nModuleRows];
1017 Double_t cellEnergy[nColsFeeModule*nModuleCols][nRowsFeeModule*nModuleRows];
1018 //Double_t fTrigCutLow = 6;
1019 //Double_t fTrigCutHigh = 10;
1020 Double_t fTimeCutLow = 469e-09;
1021 Double_t fTimeCutHigh = 715e-09;
1022
1023 AliESDCaloTrigger * fCaloTrigger = fESD->GetCaloTrigger( "EMCAL" );
1024
1025 // erase trigger maps
1026 for(Int_t i = 0; i < nColsFaltroModule*nModuleCols; i++ )
1027 {
1028 for(Int_t j = 0; j < nRowsFaltroModule*nModuleRows; j++ )
1029 {
1030 ftriggersCut[i][j] = 0;
1031 ftriggers[i][j] = 0;
1032 ftriggersTime[i][j] = 0;
1033 }
1034 }
1035
1036 Int_t iglobCol=0, iglobRow=0;
1037 // go through triggers
1038 if( fCaloTrigger->GetEntries() > 0 )
1039 {
1040 // needs reset
1041 fCaloTrigger->Reset();
1042 while( fCaloTrigger->Next() )
1043 {
1044 fCaloTrigger->GetPosition( globCol, globRow );
1045 fCaloTrigger->GetNL0Times( ntimes );
1046 /*
1047 // no L0s
1048 if( ntimes < 1 ) continue;
1049 // get precise timings
1050 fCaloTrigger->GetL0Times( trigtimes );
1051 trigInCut = 0;
1052 for(Int_t i = 0; i < ntimes; i++ )
1053 {
1054 // save the first trigger time in channel
1055 if( i == 0 || triggersTime[globCol][globRow] > trigtimes[i] )
1056 triggersTime[globCol][globRow] = trigtimes[i];
1057 //printf("trigger times: %d\n",trigtimes[i]);
1058 // check if in cut
1059 if(trigtimes[i] > fTrigCutLow && trigtimes[i] < fTrigCutHigh )
1060 trigInCut = 1;
1061
1062 fTrigTimes->Fill(trigtimes[i]);
1063 }
1064 */
1065
1066 //L1 analysis from AliAnalysisTaskEMCALTriggerQA
1067 Int_t bit = 0;
1068 fCaloTrigger->GetTriggerBits(bit);
1069
1070 Int_t ts = 0;
1071 fCaloTrigger->GetL1TimeSum(ts);
1072 if (ts > 0)ftriggers[globCol][globRow] = 1;
1073 // number of triggered channels in event
1074 nTrigChannel++;
1075 // ... inside cut
1076 if(ts>0 && (bit >> 6 & 0x1))
1077 {
1078 iglobCol = globCol;
1079 iglobRow = globRow;
1080 nTrigChannelCut++;
1081 ftriggersCut[globCol][globRow] = 1;
1082 }
1083
1084 } // calo trigger entries
1085 } // has calo trigger entries
1086
1087 // part 2 go through the clusters here -----------------------------------
1088 Int_t nCluster=0, nCell=0, iCell=0, gCell=0;
1089 Short_t cellAddr, nSACell, mclabel;
1090 //Int_t nSACell, iSACell, mclabel;
1091 Int_t iSACell;
42c75692 1092 Double_t cellAmp=0, cellTimeT=0, clusterTime=0, efrac=0;
fb87d707 1093 Int_t nSupMod, nModule, nIphi, nIeta, iphi, ieta, gphi, geta, feta, fphi;
1094
1095 TRefArray *fCaloClusters = new TRefArray();
1096 fESD->GetEMCALClusters( fCaloClusters );
1097 nCluster = fCaloClusters->GetEntries();
1098
1099
1100 // save all cells times if there are clusters
1101 if( nCluster > 0 ){
1102 // erase time map
1103 for(Int_t i = 0; i < nColsFeeModule*nModuleCols; i++ ){
1104 for(Int_t j = 0; j < nRowsFeeModule*nModuleRows; j++ ){
1105 cellTime[i][j] = 0.;
1106 cellEnergy[i][j] = 0.;
1107 }
1108 }
1109
1110 // get cells
1111 AliESDCaloCells *fCaloCells = fESD->GetEMCALCells();
1112 //AliVCaloCells fCaloCells = fESD->GetEMCALCells();
1113 nSACell = fCaloCells->GetNumberOfCells();
1114 for(iSACell = 0; iSACell < nSACell; iSACell++ ){
1115 // get the cell info *fCal
1116 fCaloCells->GetCell( iSACell, cellAddr, cellAmp, cellTimeT , mclabel, efrac);
1117
1118 // get cell position
1119 fGeom->GetCellIndex( cellAddr, nSupMod, nModule, nIphi, nIeta );
1120 fGeom->GetCellPhiEtaIndexInSModule( nSupMod,nModule, nIphi, nIeta, iphi, ieta);
1121
1122 // convert co global phi eta
1123 gphi = iphi + nRowsFeeModule*(nSupMod/2);
1124 geta = ieta + nColsFeeModule*(nSupMod%2);
1125
1126 // save cell time and energy
1127 cellTime[geta][gphi] = cellTimeT;
1128 cellEnergy[geta][gphi] = cellAmp;
1129
1130 }
1131 }
1132
1133 Int_t nClusterTrig, nClusterTrigCut;
1134 UShort_t *cellAddrs;
1135 Double_t clsE=-999, clsEta=-999, clsPhi=-999;
1136 Float_t clsPos[3] = {0.,0.,0.};
1137
1138 for(Int_t icl=0; icl<fESD->GetNumberOfCaloClusters(); icl++)
1139 {
1140 AliESDCaloCluster *cluster = fESD->GetCaloCluster(icl);
1141 if(!cluster || !cluster->IsEMCAL()) continue;
1142
1143 // get cluster cells
1144 nCell = cluster->GetNCells();
1145
1146 // get cluster energy
1147 clsE = cluster->E();
1148
1149 // get cluster position
1150 cluster->GetPosition(clsPos);
1151 TVector3 clsPosVec(clsPos[0],clsPos[1],clsPos[2]);
1152 clsEta = clsPosVec.Eta();
1153 clsPhi = clsPosVec.Phi();
1154
1155 // get the cell addresses
1156 cellAddrs = cluster->GetCellsAbsId();
1157
1158 // check if the cluster contains cell, that was marked as triggered
1159 nClusterTrig = 0;
1160 nClusterTrigCut = 0;
1161
1162 // loop the cells to check, if cluser in acceptance
1163 // any cluster with a cell outside acceptance is not considered
1164 for( iCell = 0; iCell < nCell; iCell++ )
1165 {
42c75692 1166 // check hot cell
feffe705 1167 //if(clsE>6.0)fCellCheck->Fill(clsE,cellAddrs[iCell]);
42c75692 1168
fb87d707 1169 // get cell position
1170 fGeom->GetCellIndex( cellAddrs[iCell], nSupMod, nModule, nIphi, nIeta );
1171 fGeom->GetCellPhiEtaIndexInSModule( nSupMod,nModule, nIphi, nIeta, iphi, ieta);
1172
1173 // convert co global phi eta
1174 gphi = iphi + nRowsFeeModule*(nSupMod/2);
1175 geta = ieta + nColsFeeModule*(nSupMod%2);
1176
1177 if( cellTime[geta][gphi] > 0. ){
1178 clusterTime += cellTime[geta][gphi];
1179 gCell++;
1180 }
1181
1182 // get corresponding FALTRO
1183 fphi = gphi / 2;
1184 feta = geta / 2;
1185
1186 // try to match with a triggered
1187 if( ftriggers[feta][fphi]==1)
1188 { nClusterTrig++;
1189 }
1190 if( ftriggersCut[feta][fphi]==1)
1191 { nClusterTrigCut++;
1192 }
1193
1194 } // cells
1195
1196
1197 if( gCell > 0 )
1198 clusterTime = clusterTime / (Double_t)gCell;
1199 // fix the reconstriction code time 100ns jumps
1200 if( fESD->GetBunchCrossNumber() % 4 < 2 )
1201 clusterTime -= 0.0000001;
1202
feffe705 1203 //fClsETime->Fill(clsE,clusterTime);
1204 //fClsEBftTrigCut->Fill(clsE);
fb87d707 1205
1206 if(nClusterTrig>0){
feffe705 1207 //fClsETime1->Fill(clsE,clusterTime);
fb87d707 1208 }
1209
1210 if(nClusterTrig>0){
1211 cluster->SetChi2(1);
feffe705 1212 //fClsEAftTrigCut1->Fill(clsE);
fb87d707 1213 }
1214
1215 if(nClusterTrigCut>0){
1216 cluster->SetChi2(2);
feffe705 1217 //fClsEAftTrigCut2->Fill(clsE);
fb87d707 1218 }
1219
1220 if(nClusterTrigCut>0 && ( nCell > (1 + clsE / 3)))
1221 {
1222 cluster->SetChi2(3);
feffe705 1223 //fClsEAftTrigCut3->Fill(clsE);
fb87d707 1224 }
1225
1226 if(nClusterTrigCut>0 && (nCell > (1 + clsE / 3) )&&( clusterTime > fTimeCutLow && clusterTime < fTimeCutHigh ))
1227 {
1228 // cluster->SetChi2(4);
feffe705 1229 //fClsEAftTrigCut4->Fill(clsE);
fb87d707 1230 }
1231 if(nClusterTrigCut<1)
1232 {
1233 cluster->SetChi2(0);
1234
feffe705 1235 //fClsEAftTrigCut->Fill(clsE);
fb87d707 1236 }
1237
1238 } // clusters
1239}
1240
1241
1242
1243
1244
1245
1246