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