1 /*************************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
10 * without fee, provided that the above copyright notice appears in all *
11 * copies and that both the copyright notice and this permission notice *
12 * appear in the supporting documentation. The authors make no claims *
13 * about the suitability of this software for any purpose. It is *
14 * provided "as is" without express or implied warranty. *
16 *************************************************************************************/
18 /*************************************************************************************
20 * Class for the Selection of Non-Heavy-Flavour-Electrons trought *
21 * the invariant mass method. The selection can be done from two *
22 * different algorithms, which can be choosed calling the function *
23 * "SetAlgorithm(TString Algorithm)". *
25 * Authors: R.Bailhache, C.A.Schmidt *
27 *************************************************************************************/
30 #include "THnSparse.h"
32 #include "TLorentzVector.h"
33 #include "TParticle.h"
35 #include "TDatabasePDG.h"
37 #include "AliVEvent.h"
38 #include "AliMCEvent.h"
39 #include "AliESDEvent.h"
40 #include "AliMCParticle.h"
41 #include "AliAODMCParticle.h"
42 #include "AliAODEvent.h"
43 #include "AliAODVertex.h"
44 #include "AliAODTrack.h"
45 #include "AliVTrack.h"
46 #include "AliESDtrack.h"
47 #include "AliESDtrackCuts.h"
48 #include "AliPIDResponse.h"
51 #include "AliKFParticle.h"
52 #include "AliKFVertex.h"
54 #include "AliHFEcuts.h"
55 #include "AliHFEpid.h"
56 #include "AliHFEpidQAmanager.h"
57 #include "AliHFEtools.h"
58 #include "AliHFEmcQA.h"
60 #include "AliHFENonPhotonicElectron.h"
62 ClassImp(AliHFENonPhotonicElectron)
63 //________________________________________________________________________
64 AliHFENonPhotonicElectron::AliHFENonPhotonicElectron(const char *name, const Char_t *title)
68 ,fAODArrayMCInfo (NULL)
70 ,fHFEBackgroundCuts (NULL)
77 ,fChi2OverNDFCut (3.0)
79 // ,fMaxOpeningTheta (0.02)
80 // ,fMaxOpeningPhi (0.1)
81 ,fMaxOpening3D (TMath::Pi())
83 ,fSetMassConstraint (kFALSE)
84 ,fSelectCategory1tracks(kTRUE)
85 ,fSelectCategory2tracks(kFALSE)
91 ,fCounterPoolBackground (0)
101 ,fWeightsSource(NULL)
102 ,fIncElectronRadius(NULL)
103 ,fRecElectronRadius(NULL)
104 // ,fUSignAngle (NULL)
105 // ,fLSignAngle (NULL)
110 fPIDBackground = new AliHFEpid("hfePidBackground");
111 fPIDBackgroundQA = new AliHFEpidQAmanager;
114 //________________________________________________________________________
115 AliHFENonPhotonicElectron::AliHFENonPhotonicElectron()
119 ,fAODArrayMCInfo (NULL)
121 ,fHFEBackgroundCuts (NULL)
122 ,fPIDBackground (0x0)
123 ,fPIDBackgroundQA (0)
127 ,fAlgorithmMA (kTRUE)
128 ,fChi2OverNDFCut (3.0)
130 // ,fMaxOpeningTheta (0.02)
131 // ,fMaxOpeningPhi (0.1)
132 ,fMaxOpening3D (TMath::TwoPi())
134 ,fSetMassConstraint (kFALSE)
135 ,fSelectCategory1tracks(kTRUE)
136 ,fSelectCategory2tracks(kFALSE)
142 ,fCounterPoolBackground (0)
152 ,fWeightsSource(NULL)
153 ,fIncElectronRadius(NULL)
154 ,fRecElectronRadius(NULL)
155 // ,fUSignAngle (NULL)
156 // ,fLSignAngle (NULL)
161 fPIDBackground = new AliHFEpid("hfePidBackground");
162 fPIDBackgroundQA = new AliHFEpidQAmanager;
165 //________________________________________________________________________
166 AliHFENonPhotonicElectron::AliHFENonPhotonicElectron(const AliHFENonPhotonicElectron &ref)
170 ,fAODArrayMCInfo (NULL)
171 ,fLevelBack (ref.fLevelBack)
172 ,fHFEBackgroundCuts (ref.fHFEBackgroundCuts)
173 ,fPIDBackground (ref.fPIDBackground)
174 ,fPIDBackgroundQA (ref.fPIDBackgroundQA)
175 ,fkPIDRespons (ref.fkPIDRespons)
176 ,fPtBinning(ref.fPtBinning)
177 ,fEtaBinning(ref.fEtaBinning)
178 ,fAlgorithmMA (ref.fAlgorithmMA)
179 ,fChi2OverNDFCut (ref.fChi2OverNDFCut)
180 ,fMaxDCA (ref.fMaxDCA)
181 // ,fMaxOpeningTheta (ref.fMaxOpeningTheta)
182 // ,fMaxOpeningPhi (ref.fMaxOpeningPhi)
183 ,fMaxOpening3D (ref.fMaxOpening3D)
184 ,fMaxInvMass (ref.fMaxInvMass)
185 ,fSetMassConstraint (ref.fSetMassConstraint)
186 ,fSelectCategory1tracks(ref.fSelectCategory1tracks)
187 ,fSelectCategory2tracks(ref.fSelectCategory2tracks)
188 ,fITSmeanShift(ref.fITSmeanShift)
189 ,fITSnSigmaHigh(ref.fITSnSigmaHigh)
190 ,fITSnSigmaLow(ref.fITSnSigmaLow)
193 ,fCounterPoolBackground (0)
195 ,fListOutput (ref.fListOutput)
196 ,fAssElectron (ref.fAssElectron)
197 ,fIncElectron (ref.fIncElectron)
200 ,fUSmatches(ref.fUSmatches)
201 ,fLSmatches(ref.fLSmatches)
202 ,fHnsigmaITS(ref.fHnsigmaITS)
203 ,fWeightsSource(ref.fWeightsSource)
204 ,fIncElectronRadius(ref.fIncElectronRadius)
205 ,fRecElectronRadius(ref.fRecElectronRadius)
206 // ,fUSignAngle (ref.fUSignAngle)
207 // ,fLSignAngle (ref.fLSignAngle)
215 //____________________________________________________________
216 AliHFENonPhotonicElectron &AliHFENonPhotonicElectron::operator=(const AliHFENonPhotonicElectron &ref){
218 // Assignment operator
220 if(this == &ref) ref.Copy(*this);
224 //_________________________________________
225 AliHFENonPhotonicElectron::~AliHFENonPhotonicElectron()
230 if(fArraytrack) delete fArraytrack;
231 //if(fHFEBackgroundCuts) delete fHFEBackgroundCuts;
232 if(fPIDBackground) delete fPIDBackground;
233 if(fPIDBackgroundQA) delete fPIDBackgroundQA;
236 //_____________________________________________________________________________________________
237 void AliHFENonPhotonicElectron::Init()
243 //printf("Analysis Mode for AliHFENonPhotonicElectron: %s Analysis\n", fIsAOD ? "AOD" : "ESD");
245 if(!fListOutput) fListOutput = new TList;
246 fListOutput->SetName("HFENonPhotonicElectron");
247 fListOutput->SetOwner();
249 if(!fHFEBackgroundCuts) fHFEBackgroundCuts = new AliHFEcuts();
250 if(fIsAOD) fHFEBackgroundCuts->SetAOD();
251 fHFEBackgroundCuts->Initialize();
252 if(fHFEBackgroundCuts->IsQAOn()) {
253 fListOutput->Add(fHFEBackgroundCuts->GetQAhistograms());
257 if(!fPIDBackground) fPIDBackground = new AliHFEpid("default pid");
258 if(fMCEvent || fAODArrayMCInfo) fPIDBackground->SetHasMCData(kTRUE); // does nothing since the fMCEvent are set afterwards at the moment
259 if(!fPIDBackground->GetNumberOfPIDdetectors())
261 //fPIDBackground->AddDetector("TOF", 0);
262 fPIDBackground->AddDetector("TPC", 0);
264 AliInfo("PID Background QA switched on");
265 fPIDBackgroundQA->Initialize(fPIDBackground);
266 fListOutput->Add(fPIDBackgroundQA->MakeList("HFENP_PID_Background"));
267 fPIDBackground->SortDetectors();
269 const Int_t kBinsPtDefault = 35;
270 Double_t binLimPtDefault[kBinsPtDefault+1] = {0., 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1., 1.1, 1.2, 1.3, 1.4, 1.5, 1.75, 2., 2.25, 2.5, 2.75, 3., 3.5, 4., 4.5, 5., 5.5, 6., 7., 8., 10., 12., 14., 16., 18., 20.};
271 const Int_t kBinsEtaInclusiveDefault = 8;
272 Double_t binLimEtaInclusiveDefault[kBinsEtaInclusiveDefault+1] = {-0.8, -0.6, -0.4, -0.2, 0., 0.2, 0.4, 0.6, 0.8};
273 const Int_t kBinsEtaAssociated = 15;
274 Double_t binLimEtaAssociat[kBinsEtaAssociated+1] = {-1.5,-1.3,-1.1,-0.9,-0.7,-0.5,-0.3,-0.1,0.1,0.3,0.5,0.7,0.9,1.1,1.3,1.5};
276 if(!fPtBinning.GetSize()) fPtBinning.Set(kBinsPtDefault+1, binLimPtDefault);
277 if(!fEtaBinning.GetSize()) fEtaBinning.Set(kBinsEtaInclusiveDefault+1, binLimEtaInclusiveDefault);
279 //Int_t nBinsP = 400;
280 //Double_t minP = 0.0;
281 //Double_t maxP = 20.0;
282 //Double_t binLimP[nBinsP+1];
283 //for(Int_t i=0; i<=nBinsP; i++) binLimP[i]=(Double_t)minP + (maxP-minP)/nBinsP*(Double_t)i ;
287 Double_t maxC = 11.0;
288 Double_t binLimC[nBinsC+1];
289 for(Int_t i=0; i<=nBinsC; i++) binLimC[i]=(Double_t)minC + (maxC-minC)/nBinsC*(Double_t)i ;
291 Int_t nBinsSource = 10;
292 Double_t minSource = 0.;
293 Double_t maxSource = 10.;
294 Double_t binLimSource[nBinsSource+1];
295 for(Int_t i=0; i<=nBinsSource; i++) binLimSource[i]=(Double_t)minSource + (maxSource-minSource)/nBinsSource*(Double_t)i ;
297 Int_t nBinsInvMass = 30;
298 Double_t minInvMass = 0.;
299 Double_t maxInvMass = 0.3;
300 Double_t binLimInvMass[nBinsInvMass+1];
301 for(Int_t i=0; i<=nBinsInvMass; i++) binLimInvMass[i]=(Double_t)minInvMass + (maxInvMass-minInvMass)/nBinsInvMass*(Double_t)i ;
304 Double_t minPhi = 0.0;
305 Double_t maxPhi = TMath::Pi();
306 Double_t binLimPhi[nBinsPhi+1];
307 for(Int_t i=0; i<=nBinsPhi; i++)
309 binLimPhi[i]=(Double_t)minPhi + (maxPhi-minPhi)/nBinsPhi*(Double_t)i ;
310 AliDebug(2,Form("bin phi is %f for %d",binLimPhi[i],i));
313 Int_t nBinsAngle = 72;
314 Double_t minAngle = 0.0;
315 Double_t maxAngle = 0.4;
316 Double_t binLimAngle[nBinsAngle+1];
317 for(Int_t i=0; i<=nBinsAngle; i++)
319 binLimAngle[i]=(Double_t)minAngle + (maxAngle-minAngle)/nBinsAngle*(Double_t)i ;
320 AliDebug(2,Form("bin phi is %f for %d",binLimAngle[i],i));
323 // Constrain histograms
324 const Int_t nDimSingle=4;
325 const Int_t nDimPair=9;
326 Int_t nBinPair[nDimPair] = {nBinsPhi,nBinsC,fPtBinning.GetSize()-1,nBinsInvMass,nBinsSource,nBinsAngle,fPtBinning.GetSize()-1,fEtaBinning.GetSize()-1,kBinsEtaAssociated};
328 // Associated Electron
329 Int_t nBinAssElectron[nDimSingle] = {nBinsC,fPtBinning.GetSize()-1,nBinsSource,kBinsEtaAssociated};
330 fAssElectron = new THnSparseF("fAssElectron","fAssElectron",nDimSingle,nBinAssElectron);
331 fAssElectron->SetBinEdges(0,binLimC);
332 fAssElectron->SetBinEdges(1,fPtBinning.GetArray());
333 fAssElectron->SetBinEdges(2,binLimSource);
334 fAssElectron->SetBinEdges(3,binLimEtaAssociat);
335 fAssElectron->Sumw2();
336 AliDebug(2,"AliHFENonPhotonicElectron: fAssElectron");
338 // Inclusive Electron
339 Int_t nBinIncElectron[nDimSingle] = {nBinsC,fPtBinning.GetSize()-1,nBinsSource,fEtaBinning.GetSize()-1};
340 fIncElectron = new THnSparseF("fIncElectron","fIncElectron",nDimSingle,nBinIncElectron);
341 fIncElectron->SetBinEdges(0,binLimC);
342 fIncElectron->SetBinEdges(1,fPtBinning.GetArray());
343 fIncElectron->SetBinEdges(2,binLimSource);
344 fIncElectron->SetBinEdges(3,fEtaBinning.GetArray());
345 fIncElectron->Sumw2();
346 AliDebug(2,"AliHFENonPhotonicElectron: fIncElectron");
348 // ee invariant mass Unlike Sign
349 fUSign = new THnSparseF("fUSign","fUSign",nDimPair,nBinPair);
350 fUSign->SetBinEdges(0,binLimPhi);
351 fUSign->SetBinEdges(1,binLimC);
352 fUSign->SetBinEdges(2,fPtBinning.GetArray());
353 fUSign->SetBinEdges(3,binLimInvMass);
354 fUSign->SetBinEdges(4,binLimSource);
355 fUSign->SetBinEdges(5,binLimAngle);
356 fUSign->SetBinEdges(6,fPtBinning.GetArray());
357 fUSign->SetBinEdges(7,fEtaBinning.GetArray());
358 fUSign->SetBinEdges(8,binLimEtaAssociat);
360 AliDebug(2,"AliHFENonPhotonicElectron: fUSign");
362 // ee invariant mass Like Sign
363 fLSign = new THnSparseF("fLSign","fLSign",nDimPair,nBinPair);
364 fLSign->SetBinEdges(0,binLimPhi);
365 fLSign->SetBinEdges(1,binLimC);
366 fLSign->SetBinEdges(2,fPtBinning.GetArray());
367 fLSign->SetBinEdges(3,binLimInvMass);
368 fLSign->SetBinEdges(4,binLimSource);
369 fLSign->SetBinEdges(5,binLimAngle);
370 fLSign->SetBinEdges(6,fPtBinning.GetArray());
371 fLSign->SetBinEdges(7,fEtaBinning.GetArray());
372 fLSign->SetBinEdges(8,binLimEtaAssociat);
374 AliDebug(2,"AliHFENonPhotonicElectron: fLSign");
376 // Histograms counting the number of like sign / unlike sign matches per inclusive track
377 const Int_t nBinsMatches = 50;
378 Double_t binLimMatches[nBinsMatches+1];
379 for(int ib = 0; ib <= nBinsMatches; ib++) binLimMatches[ib] = ib;
380 const Int_t nDimMatches = 3; // centrality, pt_inc, number of matches
381 const Int_t nBinsMatchHist[nDimMatches] = {nBinsC, fPtBinning.GetSize()-1, nBinsMatches};
382 fUSmatches = new THnSparseF("fUSmatches", "fUSmatches", nDimMatches, nBinsMatchHist);
383 fUSmatches->SetBinEdges(0,binLimC);
384 fUSmatches->SetBinEdges(1,fPtBinning.GetArray());
385 fUSmatches->SetBinEdges(2,binLimMatches);
387 fLSmatches = new THnSparseF("fLSmatches", "fLSmatches", nDimMatches, nBinsMatchHist);
388 fLSmatches->SetBinEdges(0,binLimC);
389 fLSmatches->SetBinEdges(1,fPtBinning.GetArray());
390 fLSmatches->SetBinEdges(2,binLimMatches);
392 // Histograms for radius studies
393 Int_t nBinsradius = 50;
394 Double_t minradius = 0.0;
395 Double_t maxradius = 100.0;
396 Double_t binLimradius[nBinsradius+1];
397 for(Int_t i=0; i<=nBinsradius; i++) binLimradius[i]=(Double_t)minradius + (maxradius-minradius)/nBinsradius*(Double_t)i ;
398 const Int_t nDimIncElectronRadius = 4; // centrality, pt_inc, radius
399 const Int_t nBinsIncElectronRadius[nDimIncElectronRadius] = {nBinsC, fPtBinning.GetSize()-1, nBinsradius, nBinsSource};
400 fIncElectronRadius = new THnSparseF("fIncElectronRadius", "fIncElectronRadius", nDimIncElectronRadius, nBinsIncElectronRadius);
401 fIncElectronRadius->SetBinEdges(0,binLimC);
402 fIncElectronRadius->SetBinEdges(1,fPtBinning.GetArray());
403 fIncElectronRadius->SetBinEdges(2,binLimradius);
404 fIncElectronRadius->SetBinEdges(3,binLimSource);
406 fRecElectronRadius = new THnSparseF("fRecElectronRadius", "fRecElectronRadius", nDimIncElectronRadius, nBinsIncElectronRadius);
407 fRecElectronRadius->SetBinEdges(0,binLimC);
408 fRecElectronRadius->SetBinEdges(1,fPtBinning.GetArray());
409 fRecElectronRadius->SetBinEdges(2,binLimradius);
410 fRecElectronRadius->SetBinEdges(2,binLimSource);
413 // ee angle Unlike Sign
414 const Int_t nDimUSignAngle=3;
415 Int_t nBinUSignAngle[nDimUSignAngle] = {nBinsAngle,nBinsC,nBinsSource};
416 fUSignAngle = new THnSparseF("fUSignAngle","fUSignAngle",nDimUSignAngle,nBinUSignAngle);
417 fUSignAngle->SetBinEdges(0,binLimAngle);
418 fUSignAngle->SetBinEdges(1,binLimC);
419 fUSignAngle->SetBinEdges(2,binLimSource);
420 fUSignAngle->Sumw2();
421 AliDebug(2,"AliHFENonPhotonicElectron: fUSignAngle");
423 // ee angle Like Sign
424 const Int_t nDimLSignAngle=3;
425 Int_t nBinLSignAngle[nDimLSignAngle] = {nBinsAngle,nBinsC,nBinsSource};
426 fLSignAngle = new THnSparseF("fLSignAngle","fLSignAngle",nDimLSignAngle,nBinLSignAngle);
427 fLSignAngle->SetBinEdges(0,binLimAngle);
428 fLSignAngle->SetBinEdges(1,binLimC);
429 fLSignAngle->SetBinEdges(2,binLimSource);
430 fLSignAngle->Sumw2();
431 AliDebug(2,"AliHFENonPhotonicElectron: fLSignAngle");
434 // control histogram for ITS PID
435 fHnsigmaITS = new TH2F("fHnsigmaITS", "Number of sigmas in the ITS", 30, 0., 0.3, 1200, -10., 10.);
437 // control histogram for weights sources
438 fWeightsSource = new TH2F("fWeightsSource", "Source code for weights", 11, -1.5, 9.5, 29, -1.5, 27.5);
440 fListOutput->Add(fAssElectron);
441 fListOutput->Add(fIncElectron);
442 fListOutput->Add(fUSign);
443 fListOutput->Add(fLSign);
444 fListOutput->Add(fUSmatches);
445 fListOutput->Add(fLSmatches);
446 fListOutput->Add(fHnsigmaITS);
447 fListOutput->Add(fWeightsSource);
448 fListOutput->Add(fIncElectronRadius);
449 fListOutput->Add(fRecElectronRadius);
450 // fListOutput->Add(fUSignAngle);
451 // fListOutput->Add(fLSignAngle);
455 //_____________________________________________________________________________________________
456 void AliHFENonPhotonicElectron::SetWithWeights(Int_t levelBack)
459 // Init the HFE level
461 if(levelBack >= 0) fLevelBack = levelBack;
465 //_____________________________________________________________________________________________
466 void AliHFENonPhotonicElectron::SetMCEvent(AliMCEvent *mcEvent)
476 //_____________________________________________________________________________________________
477 void AliHFENonPhotonicElectron::SetAODArrayMCInfo(TClonesArray *aodArrayMCInfo)
480 // Pass the mcEvent info
483 fAODArrayMCInfo = aodArrayMCInfo;
487 //_____________________________________________________________________________________________
488 void AliHFENonPhotonicElectron::InitRun(const AliVEvent *inputEvent,const AliPIDResponse *pidResponse)
496 AliDebug(1, "Using default PID Response");
497 Bool_t hasmc = kFALSE;
498 if(fMCEvent || fAODArrayMCInfo) hasmc=kTRUE;
499 pidResponse = AliHFEtools::GetDefaultPID(hasmc, inputEvent->IsA() == AliESDEvent::Class());
502 if(!fPIDBackground) return;
503 fPIDBackground->SetPIDResponse(pidResponse);
505 if(!fPIDBackground->IsInitialized())
507 // Initialize PID with the given run number
508 fPIDBackground->InitializePID(inputEvent->GetRunNumber());
513 //_____________________________________________________________________________________________
514 Int_t AliHFENonPhotonicElectron::FillPoolAssociatedTracks(AliVEvent *inputEvent, Int_t binct)
517 // Fill the pool of associated tracks
518 // Return the number of associated tracks
523 fHFEBackgroundCuts->SetRecEvent(inputEvent);
524 Int_t nbtracks = inputEvent->GetNumberOfTracks();
527 fArraytrack->~TArrayI();
528 new(fArraytrack) TArrayI(nbtracks);
530 fArraytrack = new TArrayI(nbtracks);
533 fCounterPoolBackground = 0;
535 Bool_t isSelected(kFALSE);
536 Bool_t isAOD = (dynamic_cast<AliAODEvent *>(inputEvent) != NULL);
537 AliDebug(2, Form("isAOD: %s", isAOD ? "yes" : "no"));
538 for(Int_t k = 0; k < nbtracks; k++) {
539 AliVTrack *track = (AliVTrack *) inputEvent->GetTrack(k);
544 if(fSelectCategory1tracks && FilterCategory1Track(track, isAOD, binct)) isSelected = kTRUE;
545 else if(fSelectCategory2tracks && FilterCategory2Track(track, isAOD)) isSelected = kTRUE;
548 AliDebug(2,Form("fCounterPoolBackground %d, track %d",fCounterPoolBackground,k));
549 fArraytrack->AddAt(k,fCounterPoolBackground);
550 fCounterPoolBackground++;
554 //printf(Form("Associated Pool: Tracks %d, fCounterPoolBackground %d \n", nbtracks, fCounterPoolBackground));
556 return fCounterPoolBackground;
560 //_____________________________________________________________________________________________
561 Int_t AliHFENonPhotonicElectron::CountPoolAssociated(AliVEvent *inputEvent, Int_t binct)
564 // Count the pool of assiocated tracks
568 if(fnumberfound > 0) //!count only events with an inclusive electron
570 Double_t valueAssElectron[4] = { (Double_t)binct, -1, -1}; //Centrality Pt Source
572 Int_t indexmother2 = -1;
573 AliVTrack *track2 = 0x0;
575 for(Int_t ii = 0; ii < fCounterPoolBackground; ii++){
576 iTrack2 = fArraytrack->At(ii);
577 AliDebug(2,Form("track %d",iTrack2));
578 track2 = (AliVTrack *)inputEvent->GetTrack(iTrack2);
581 //printf("ERROR: Could not receive track %d", iTrack2);
586 if(fMCEvent || fAODArrayMCInfo) valueAssElectron[2] = FindMother(TMath::Abs(track2->GetLabel()), indexmother2) ;
588 fkPIDRespons = fPIDBackground->GetPIDResponse();
590 valueAssElectron[1] = track2->Pt() ;
591 valueAssElectron[3] = track2->Eta() ;
593 fAssElectron->Fill( valueAssElectron) ;
595 //printf(Form("Associated Pool: fCounterPoolBackground %d \n", fCounterPoolBackground));
600 //_____________________________________________________________________________________________
601 Int_t AliHFENonPhotonicElectron::LookAtNonHFE(Int_t iTrack1, AliVTrack *track1, AliVEvent *vEvent, Double_t weight, Int_t binct, Double_t deltaphi, Int_t source, Int_t indexmother,Int_t mcQAsource)
607 /***********************************************************************************
609 * iTrack1: index of the tagged electrons in AliVEvent *
610 * track1: tagged electron *
612 * weight: weight in pt if not realistic *
613 * binct: centrality bin *
614 * deltaphi: phi-phi event plane for v2 *
615 * source: MC sources *
616 * indexmother: MC index mother *
619 * return -1 if nothing *
620 * return 2 if opposite charge within the mass range *
621 * return 4 if like charge within the mass range *
622 * return 6 if opposite & like charge within the mass range *
624 ***********************************************************************************/
626 //printf("weight %f and source %d\n",weight,source);
629 AliAODEvent *aodeventu = dynamic_cast<AliAODEvent*>(vEvent);
630 Int_t taggedphotonic = -1;
632 AliDebug(2,Form("fCounterPoolBackground %d in LookAtNonHFE!!!",fCounterPoolBackground));
633 if(!fArraytrack) return taggedphotonic;
634 AliDebug(2,Form("process track %d",iTrack1));
635 AliDebug(1,Form("Inclusive source is %d\n", source));
637 fkPIDRespons = fPIDBackground->GetPIDResponse();
639 //Set Fill-Arrays for THnSparse
640 Double_t valueIncElectron[4] = { (Double_t)binct, track1->Pt(), (Double_t)source, track1->Eta()}; //Centrality Pt Source P
641 Double_t valueSign[9] = { deltaphi, (Double_t)binct, track1->Pt(), -1, (Double_t)source, -1, -1, track1->Eta(), -1}; //DeltaPhi Centrality Pt InvariantMass Source Angle Pt
642 //Double_t valueAngle[3] = { -1, binct, source};
643 Double_t valueradius[4] = { (Double_t)binct, track1->Pt(), 0.,(Double_t)source};
646 Int_t pdg1 = CheckPdg(TMath::Abs(track1->GetLabel()));
647 Double_t radius = Radius(TMath::Abs(track1->GetLabel()));
648 AliKFParticle::SetField(vEvent->GetMagneticField());
649 AliKFVertex primV(*(vEvent->GetPrimaryVertex()));
650 valueradius[2] = radius;
652 AliVTrack *track2(NULL);
654 Int_t indexmother2 = -1;
657 Float_t fCharge2 = 0;
659 // count number of matches with opposite/same sign track in the given mass range
660 Int_t countsMatchLikesign(0),
661 countsMatchUnlikesign(0);
665 Double_t invmass(-1);
667 Float_t fCharge1 = track1->Charge(); //Charge from track1
669 Bool_t kUSignPhotonic = kFALSE;
670 Bool_t kLSignPhotonic = kFALSE;
672 //! FILL Inclusive Electron
673 fWeightsSource->Fill(source,mcQAsource);
674 fIncElectron->Fill(valueIncElectron,weight);
676 fIncElectronRadius->Fill(valueradius,weight);
678 //printf(Form("Inclusive Pool: TrackNr. %d, fnumberfound %d \n", iTrack1, fnumberfound));
680 for(Int_t idex = 0; idex < fCounterPoolBackground; idex++){
681 iTrack2 = fArraytrack->At(idex);
682 AliDebug(2,Form("track %d",iTrack2));
683 track2 = (AliVTrack *)vEvent->GetTrack(iTrack2);
686 //printf("ERROR: Could not receive track %d", iTrack2);
690 fCharge2 = track2->Charge(); //Charge from track2
693 //valueAngle[2] = source;
694 valueSign[4] = source;
695 valueSign[6] = track2->Pt();
696 valueSign[8] = track2->Eta();
698 // track cuts and PID already done
700 // Checking if it is the same Track!
701 if(iTrack2==iTrack1) continue;
702 AliDebug(2,"Different");
705 if(fMCEvent || fAODArrayMCInfo){
706 AliDebug(2, "Checking for source");
707 source2 = FindMother(TMath::Abs(track2->GetLabel()), indexmother2);
708 AliDebug(2, Form("source is %d", source2));
709 pdg2 = CheckPdg(TMath::Abs(track2->GetLabel()));
711 if(source == kElectronfromconversion){
712 AliDebug(2, Form("Electron from conversion (source %d), paired with source %d", source, source2));
713 AliDebug(2, Form("Index of the mothers: incl %d, associated %d", indexmother, indexmother2));
714 AliDebug(2, Form("PDGs: incl %d, associated %d", pdg1, pdg2));
718 if((indexmother2 == indexmother) && (source == source2) && ((pdg1*pdg2)<0.0)){
719 AliDebug(1, "Real pair");
721 case kElectronfromconversion:
722 valueSign[4] = kElectronfromconversionboth;
723 valueradius[3] = kElectronfromconversionboth;
725 case kElectronfrompi0:
726 valueSign[4] = kElectronfrompi0both;
727 valueradius[3] = kElectronfrompi0both;
729 case kElectronfrometa:
730 valueSign[4] = kElectronfrometaboth;
731 valueradius[3] = kElectronfrometaboth;
739 // Use TLorentzVector
740 if(!MakePairDCA(track1, track2, vEvent, (aodeventu != NULL), invmass, angle)) continue;
743 if(!MakePairKF(track1, track2, primV, invmass, angle)) continue;
746 valueSign[3] = invmass;
747 valueSign[5] = angle;
749 //if((fCharge1*fCharge2)>0.0) fLSignAngle->Fill(&valueAngle[0],weight);
750 //else fUSignAngle->Fill(&valueAngle[0],weight);
752 if(angle > fMaxOpening3D) continue; //! Cut on Opening Angle
753 if(invmass > fMaxInvMass) continue; //! Cut on Invariant Mass
755 if((fCharge1*fCharge2)>0.0){
756 if(invmass < 1.0) fLSign->Fill( valueSign, weight);
757 // count like-sign background matched pairs per inclusive based on mass cut
758 if(invmass < 0.14) countsMatchLikesign++;
759 AliDebug(1, "Selected Like sign");
761 if(invmass < 1.0)fUSign->Fill( valueSign, weight);
762 // count unlike-sign matched pairs per inclusive based on mass cut
764 countsMatchUnlikesign++;
765 fRecElectronRadius->Fill(valueradius,weight);
767 AliDebug(1, "Selected Unike sign");
770 if((fCharge1*fCharge2)>0.0) kLSignPhotonic=kTRUE;
771 else kUSignPhotonic=kTRUE;
775 Double_t valCountsLS[3] = {(Double_t)binct, track1->Pt(), (Double_t)countsMatchLikesign},
776 valCountsUS[3] = {(Double_t)binct, track1->Pt(), (Double_t)countsMatchUnlikesign};
777 fUSmatches->Fill(valCountsUS);
778 fLSmatches->Fill(valCountsLS);
780 if( kUSignPhotonic && kLSignPhotonic) taggedphotonic = 6;
781 if(!kUSignPhotonic && kLSignPhotonic) taggedphotonic = 4;
782 if( kUSignPhotonic && !kLSignPhotonic) taggedphotonic = 2;
784 return taggedphotonic;
787 //_________________________________________________________________________
788 Int_t AliHFENonPhotonicElectron::FindMother(Int_t tr, Int_t &indexmother) const {
790 // Find the mother if MC
793 if(!fMCEvent && !fAODArrayMCInfo) return -1;
795 Int_t pdg = CheckPdg(tr);
796 if(TMath::Abs(pdg)!= 11)
802 indexmother = IsMotherGamma(tr);
803 if(indexmother > 0) return kElectronfromconversion;
804 indexmother = IsMotherPi0(tr);
805 if(indexmother > 0) return kElectronfrompi0;
806 indexmother = IsMotherC(tr);
807 if(indexmother > 0) return kElectronfromC;
808 indexmother = IsMotherB(tr);
809 if(indexmother > 0) return kElectronfromB;
810 indexmother = IsMotherEta(tr);
811 if(indexmother > 0) return kElectronfrometa;
813 return kElectronfromother;
816 //________________________________________________________________________________________________
817 Int_t AliHFENonPhotonicElectron::CheckPdg(Int_t tr) const {
820 // Return the pdg of the particle
824 if(tr < 0) return pdgcode;
826 AliMCParticle *mctrackesd = NULL; AliAODMCParticle *mctrackaod = NULL;
828 AliVParticle *mctrack = fMCEvent->GetTrack(tr);
830 if((mctrackesd = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(tr))))) pdgcode = mctrackesd->PdgCode();
831 else if((mctrackaod = dynamic_cast<AliAODMCParticle *>(fMCEvent->GetTrack(TMath::Abs(tr))))) pdgcode = mctrackaod->GetPdgCode();
833 } else if(fAODArrayMCInfo) {
834 if(tr < fAODArrayMCInfo->GetEntriesFast()){
835 mctrackaod = (AliAODMCParticle *) fAODArrayMCInfo->At(tr);
836 if(mctrackaod) return pdgcode = mctrackaod->GetPdgCode();
843 //________________________________________________________________________________________________
844 Double_t AliHFENonPhotonicElectron::Radius(Int_t tr) const {
847 // Return the production vertex radius
850 Double_t radius = 0.;
851 if(tr < 0) return radius;
853 AliMCParticle *mctrackesd = NULL; AliAODMCParticle *mctrackaod = NULL;
855 AliVParticle *mctrack = fMCEvent->GetTrack(tr);
857 if((mctrackesd = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(tr))))) radius = TMath::Sqrt(mctrackesd->Xv()*mctrackesd->Xv()+mctrackesd->Yv()*mctrackesd->Yv());
858 else if((mctrackaod = dynamic_cast<AliAODMCParticle *>(fMCEvent->GetTrack(TMath::Abs(tr))))) radius = TMath::Sqrt(mctrackaod->Xv()*mctrackaod->Xv()+mctrackaod->Yv()*mctrackaod->Yv());
860 } else if(fAODArrayMCInfo) {
861 if(tr < fAODArrayMCInfo->GetEntriesFast()){
862 mctrackaod = (AliAODMCParticle *) fAODArrayMCInfo->At(tr);
863 if(mctrackaod) return radius = TMath::Sqrt(mctrackaod->Xv()*mctrackaod->Xv()+mctrackaod->Yv()*mctrackaod->Yv());
870 //_______________________________________________________________________________________________
871 Int_t AliHFENonPhotonicElectron::GetMotherPDG(Int_t tr, Int_t &motherIndex) const {
873 // Returns the mother PDG of the track (return value) and the index of the
876 if(tr < 0) return -1;
879 AliMCParticle *mctrackesd(NULL); AliAODMCParticle *mctrackaod(NULL);
883 AliDebug(2, "Using MC Event");
884 AliVParticle *mctrack = fMCEvent->GetTrack(tr);
886 if((mctrackesd = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(tr))))){
888 TParticle *particle = mctrackesd->Particle();
892 motherIndex = particle->GetFirstMother();
893 if(motherIndex >= 0){
894 AliMCParticle *mothertrack = NULL;
895 if((mothertrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(motherIndex))))){
896 TParticle * mother = mothertrack->Particle();
897 pdg = mother->GetPdgCode();
901 } else if((mctrackaod = dynamic_cast<AliAODMCParticle *>(fMCEvent->GetTrack(TMath::Abs(tr))))){
904 motherIndex = mctrackaod->GetMother();
905 if(motherIndex >= 0){
906 AliAODMCParticle *mothertrack = dynamic_cast<AliAODMCParticle *>(fMCEvent->GetTrack(motherIndex));
907 if(mothertrack) pdg = mothertrack->GetPdgCode();
911 } else if(fAODArrayMCInfo) {
912 AliDebug(2, "Using AOD list");
913 if(tr < fAODArrayMCInfo->GetEntriesFast()){
914 mctrackaod = (AliAODMCParticle *) fAODArrayMCInfo->At(tr);
918 motherIndex = mctrackaod->GetMother();
919 if(motherIndex >= 0 && motherIndex < fAODArrayMCInfo->GetEntriesFast()){
920 AliAODMCParticle *mothertrack = dynamic_cast<AliAODMCParticle *>(fAODArrayMCInfo->At(TMath::Abs(motherIndex)));
921 if(mothertrack) pdg = mothertrack->GetPdgCode();
929 //_______________________________________________________________________________________________
930 Int_t AliHFENonPhotonicElectron::IsMotherGamma(Int_t tr) const {
933 // Return the lab of gamma mother or -1 if not gamma
936 Int_t imother(-1), pdg(-1);
937 pdg = GetMotherPDG(tr, imother);
941 if(TMath::Abs(pdg) == 22){
942 AliDebug(2, "Gamma Mother selected");
945 if(TMath::Abs(pdg) == 11){
946 AliDebug(2, "Mother is electron - look further in hierarchy");
947 return IsMotherGamma(imother);
949 AliDebug(2, "Nothing selected");
952 AliDebug(2, "Not mother");
956 //________________________________________________________________________________________________
957 Int_t AliHFENonPhotonicElectron::IsMotherPi0(Int_t tr) const {
960 // Return the lab of pi0 mother or -1 if not pi0
963 Int_t imother(-1), pdg(-1);
964 pdg = GetMotherPDG(tr, imother);
968 if(TMath::Abs(pdg) == 111){
969 AliDebug(2, "Pi0 Mother selected");
972 if(TMath::Abs(pdg) == 11){
973 AliDebug(2, "Mother is electron - look further in hierarchy");
974 return IsMotherPi0(imother);
976 AliDebug(2, "Nothing selected");
979 AliDebug(2, "Not mother");
982 //________________________________________________________________________________________________
983 Int_t AliHFENonPhotonicElectron::IsMotherC(Int_t tr) const {
986 // Return the lab of signal mother or -1 if not from C
989 Int_t imother(-1), pdg(-1);
990 pdg = GetMotherPDG(tr, imother);
994 if((TMath::Abs(pdg)==411) || (TMath::Abs(pdg)==421) || (TMath::Abs(pdg)==431) || (TMath::Abs(pdg)==4122) || (TMath::Abs(pdg)==4132) || (TMath::Abs(pdg)==4232) || (TMath::Abs(pdg)==43320)){
995 AliDebug(2, "Charm Mother selected");
998 if(TMath::Abs(pdg) == 11){
999 AliDebug(2, "Mother is electron - look further in hierarchy");
1000 return IsMotherC(imother);
1002 AliDebug(2, "Nothing selected");
1005 AliDebug(2, "Not mother");
1009 //_______________________________________________________________________________________________
1010 Int_t AliHFENonPhotonicElectron::IsMotherB(Int_t tr) const {
1013 // Return the lab of signal mother or -1 if not B
1016 Int_t imother(-1), pdg(-1);
1017 pdg = GetMotherPDG(tr, imother);
1021 if((TMath::Abs(pdg)==511) || (TMath::Abs(pdg)==521) || (TMath::Abs(pdg)==531) || (TMath::Abs(pdg)==5122) || (TMath::Abs(pdg)==5132) || (TMath::Abs(pdg)==5232) || (TMath::Abs(pdg)==53320)){
1022 AliDebug(2, "Bottom Mother selected");
1025 if(TMath::Abs(pdg) == 11){
1026 return IsMotherB(imother);
1027 AliDebug(2, "Mother is electron - look further in hierarchy");
1029 AliDebug(2, "Nothing selected");
1032 AliDebug(2, "Not mother");
1036 //_______________________________________________________________________________________________
1037 Int_t AliHFENonPhotonicElectron::IsMotherEta(Int_t tr) const {
1040 // Return the lab of eta mother or -1 if not eta
1043 Int_t imother(-1), pdg(-1);
1044 pdg = GetMotherPDG(tr, imother);
1048 if(TMath::Abs(pdg) == 221){
1049 AliDebug(2, "Eta mother selected");
1052 if(TMath::Abs(pdg) == 11){
1053 AliDebug(2, "Mother is electron - look further in hierarchy");
1054 return IsMotherEta(imother);
1056 AliDebug(2, "Nothing selected");
1059 AliDebug(2, "Not mother");
1063 //_______________________________________________________________________________________________
1064 Bool_t AliHFENonPhotonicElectron::MakePairDCA(const AliVTrack *inclusive, const AliVTrack *associated, AliVEvent *vEvent, Bool_t isAOD, Double_t &invMass, Double_t &angle) const {
1066 // Make Pairs of electrons using TLorentzVector
1068 Double_t eMass = TDatabasePDG::Instance()->GetParticle(11)->Mass(); //Electron mass in GeV
1069 Double_t bfield = vEvent->GetMagneticField();
1071 AliESDtrack *esdtrack1, *esdtrack2;
1073 // call copy constructor for AODs
1074 esdtrack1 = new AliESDtrack(inclusive);
1075 esdtrack2 = new AliESDtrack(associated);
1077 // call copy constructor for ESDs
1078 esdtrack1 = new AliESDtrack(*(static_cast<const AliESDtrack *>(inclusive)));
1079 esdtrack2 = new AliESDtrack(*(static_cast<const AliESDtrack *>(associated)));
1081 if((!esdtrack1) || (!esdtrack2)){
1087 Double_t xt1 = 0; //radial position track 1 at the DCA point
1088 Double_t xt2 = 0; //radial position track 2 at the DCA point
1089 Double_t dca = esdtrack2->GetDCA(esdtrack1,bfield,xt2,xt1); //DCA track1-track2
1091 // Apply DCA cut already in the function
1097 //Momenta of the track extrapolated to DCA track-track
1098 Double_t p1[3] = {0,0,0};
1099 Double_t p2[3] = {0,0,0};
1100 Bool_t kHasdcaT1 = esdtrack1->GetPxPyPzAt(xt1,bfield,p1); //Track1
1101 Bool_t kHasdcaT2 = esdtrack2->GetPxPyPzAt(xt2,bfield,p2); //Track2
1102 if(!kHasdcaT1 || !kHasdcaT2) AliWarning("It could be a problem in the extrapolation");
1104 TLorentzVector electron1, electron2, mother;
1105 electron1.SetXYZM(p1[0], p1[1], p1[2], eMass);
1106 electron2.SetXYZM(p2[0], p2[1], p2[2], eMass);
1108 mother = electron1 + electron2;
1109 invMass = mother.M();
1110 angle = TVector2::Phi_0_2pi(electron1.Angle(electron2.Vect()));
1117 //_______________________________________________________________________________________________
1118 Bool_t AliHFENonPhotonicElectron::MakePairKF(const AliVTrack *inclusive, const AliVTrack *associated, AliKFVertex &primV, Double_t &invMass, Double_t &angle) const {
1120 // Make pairs of electrons using the AliKF package
1123 //printf("AOD HFE non photonic\n");
1125 Int_t fPDGtrack1 = 11;
1126 if(inclusive->Charge()>0) fPDGtrack1 = -11;
1127 Int_t fPDGtrack2 = 11;
1128 if(associated->Charge()>0) fPDGtrack2 = -11;
1130 AliKFParticle ktrack1(*inclusive, fPDGtrack1);
1131 AliKFParticle ktrack2(*associated, fPDGtrack2);
1132 AliKFParticle recoGamma(ktrack1,ktrack2);
1134 if(recoGamma.GetNDF()<1) return kFALSE; //! Cut on Reconstruction
1136 Double_t chi2OverNDF = recoGamma.GetChi2()/recoGamma.GetNDF();
1137 if(TMath::Sqrt(TMath::Abs(chi2OverNDF))>fChi2OverNDFCut) return kFALSE;
1139 if(fSetMassConstraint){
1143 recoGamma.SetProductionVertex(primV);
1144 recoGamma.SetMassConstraint(0,0.0001);
1148 recoGamma.GetMass(invMass,width);
1149 angle = ktrack1.GetAngle(ktrack2);
1153 //_______________________________________________________________________________________________
1154 Bool_t AliHFENonPhotonicElectron::FilterCategory1Track(const AliVTrack * const track, Bool_t isAOD, Int_t binct){
1156 // Selection of good associated tracks for the pool
1157 // selection is done using strong cuts
1158 // Tracking in the TPC and the ITS is a minimal requirement
1160 Bool_t survivedbackground = kTRUE;
1162 if(!fHFEBackgroundCuts->CheckParticleCuts(AliHFEcuts::kStepRecKineITSTPC + AliHFEcuts::kNcutStepsMCTrack, (TObject *) track)) survivedbackground = kFALSE;
1163 AliDebug(3, Form("First cut: %s\n", survivedbackground == kTRUE ? "yes" : "no"));
1164 if(!fHFEBackgroundCuts->CheckParticleCuts(AliHFEcuts::kStepRecPrim + AliHFEcuts::kNcutStepsMCTrack, (TObject *) track)) survivedbackground = kFALSE;
1165 AliDebug(3, Form("Second cut: %s\n", survivedbackground == kTRUE ? "yes" : "no"));
1167 if(survivedbackground){
1169 AliHFEpidObject hfetrack2;
1171 if(!isAOD) hfetrack2.SetAnalysisType(AliHFEpidObject::kESDanalysis);
1172 else hfetrack2.SetAnalysisType(AliHFEpidObject::kAODanalysis);
1174 hfetrack2.SetRecTrack(track);
1176 hfetrack2.SetCentrality((Int_t)binct);
1177 AliDebug(3,Form("centrality %d and %d",binct,hfetrack2.GetCentrality()));
1178 hfetrack2.SetPbPb();
1181 if(fPIDBackground->IsSelected(&hfetrack2,0x0,"recTrackCont",fPIDBackgroundQA)){
1182 survivedbackground = kTRUE;
1183 } else survivedbackground = kFALSE;
1185 AliDebug(3, Form("PID: %s\n", survivedbackground == kTRUE ? "yes" : "no"));
1186 return survivedbackground;
1189 //_______________________________________________________________________________________________
1190 Bool_t AliHFENonPhotonicElectron::FilterCategory2Track(const AliVTrack * const track, Bool_t isAOD){
1192 // Selection of category 2 tracks: These tracks are exclusively low pt tracks below
1193 // 300 MeV/c which have at least 2 hits in the 4 outer ITS layers and are identified as
1194 // electron candidates by the ITS
1196 if(TMath::Abs(track->Pt()) > 0.3) return kFALSE;
1197 if(TMath::Abs(track->Pt()) < fminPt) return kFALSE;
1198 Int_t nclustersITS(0), nclustersOuter(0);
1200 const AliAODTrack *aodtrack = static_cast<const AliAODTrack *>(track);
1201 if(!(aodtrack->TestFilterBit(AliAODTrack::kTrkGlobalNoDCA) || aodtrack->TestFilterBit(AliAODTrack::kTrkITSsa))) return kFALSE;
1202 if(!aodtrack->IsOn(AliAODTrack::kITSrefit)) return kFALSE;
1203 nclustersITS = aodtrack->GetITSNcls();
1204 for(int ily = 2; ily < 5; ily++)
1205 if(aodtrack->HasPointOnITSLayer(ily)) nclustersOuter++;
1207 const AliESDtrack *esdtrack = static_cast<const AliESDtrack *>(track);
1208 if(esdtrack->GetStatus() & AliESDtrack::kITSpureSA) return kFALSE;
1209 if(esdtrack->GetStatus() & AliESDtrack::kTPCin) return kFALSE;
1210 if(!(esdtrack->GetStatus() & AliESDtrack::kITSrefit)) return kFALSE;
1211 nclustersITS = esdtrack->GetITSclusters(NULL);
1212 for(int ily = 2; ily < 5; ily++)
1213 if(esdtrack->HasPointOnITSLayer(ily)) nclustersOuter++;
1215 if(nclustersITS < 4) return kFALSE;
1216 if(nclustersOuter < 3) return kFALSE;
1219 Double_t nsigmaITS = fPIDBackground->GetPIDResponse()->NumberOfSigmasITS(track, AliPID::kElectron);
1220 fHnsigmaITS->Fill(track->Pt(), nsigmaITS);
1221 if((nsigmaITS - fITSmeanShift) > fITSnSigmaHigh ) return kFALSE;
1222 if((nsigmaITS - fITSmeanShift) < fITSnSigmaLow ) return kFALSE;
1223 // if global track, we apply also TPC PID