]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGHF/hfe/AliAnalysisTaskHFE.cxx
Try to understannd
[u/mrichter/AliRoot.git] / PWGHF / hfe / AliAnalysisTaskHFE.cxx
CommitLineData
01745870 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// The analysis task:
17// Filling an AliCFContainer with the quantities pt, eta and phi
18// for tracks which survivied the particle cuts (MC resp. ESD tracks)
19// Track selection is done using the AliHFE package
20//
21// Author:
22// Raphaelle Bailhache <R.Bailhache@gsi.de>
23// Markus Fasel <M.Fasel@gsi.de>
24// Matus Kalisky <matus.kalisky@cern.ch>
25// MinJung Kweon <minjung@physi.uni-heidelberg.de>
26//
27#include <TAxis.h>
28#include <TBits.h>
29#include <TCanvas.h>
30#include <TChain.h>
31#include <TDirectory.h>
32#include <TFile.h>
33#include <TH3D.h>
34#include <TIterator.h>
35#include <TList.h>
36#include <TLegend.h>
37#include <TMath.h>
38#include <TObjArray.h>
39#include <TObjString.h>
40#include <TParticle.h>
41#include <TProfile.h>
42#include <TString.h>
43#include <TF1.h>
44#include <TTree.h>
45
46#include "AliESDtrackCuts.h"
47#include "AliAnalysisManager.h"
48#include "AliAnalysisUtils.h"
49#include "AliAODInputHandler.h"
50#include "AliAODMCParticle.h"
51#include "AliAODTrack.h"
52#include "AliAODVertex.h"
53#include "AliCentrality.h"
54#include "AliCFContainer.h"
55#include "AliCFManager.h"
56#include "AliESDEvent.h"
57#include "AliESDInputHandler.h"
58#include "AliESDtrack.h"
59#include "AliLog.h"
60#include "AliMCEvent.h"
61#include "AliMCEventHandler.h"
62#include "AliMCParticle.h"
63#include "AliMultiplicity.h"
64#include "AliPID.h"
65#include "AliPIDResponse.h"
66#include "AliOADBContainer.h"
67#include "AliStack.h"
68#include "AliTriggerAnalysis.h"
4437a0d2 69#include "AliTRDTriggerAnalysis.h"
01745870 70#include "AliVVertex.h"
71
72#include "AliHFEcollection.h"
73#include "AliHFEcontainer.h"
74#include "AliHFEcuts.h"
75#include "AliHFEelecbackground.h"
76#include "AliHFENonPhotonicElectron.h"
77#include "AliHFEmcQA.h"
78#include "AliHFEpairs.h"
79#include "AliHFEpid.h"
80#include "AliHFEpidQAmanager.h"
81#include "AliHFEsecVtxs.h"
82#include "AliHFEsecVtx.h"
83#include "AliHFEsignalCuts.h"
84#include "AliHFEtaggedTrackAnalysis.h"
85#include "AliHFEtools.h"
7986e54e 86#include "AliHFEV0taginfo.h"
01745870 87#include "AliHFEvarManager.h"
88#include "AliAnalysisTaskHFE.h"
89#include "AliAODMCHeader.h"
90#include "TClonesArray.h"
91
92ClassImp(AliAnalysisTaskHFE)
93
94//____________________________________________________________
95AliAnalysisTaskHFE::AliAnalysisTaskHFE():
96AliAnalysisTaskSE("PID efficiency Analysis")
97 , fAODMCHeader(NULL)
98 , fAODArrayMCInfo(NULL)
99 , fQAlevel(0)
100 , fPlugins(0)
101 , fCollisionSystem(3)
102 , fFillSignalOnly(kTRUE)
103 , fFillNoCuts(kFALSE)
4437a0d2 104 , fApplyCutAOD(kTRUE)
01745870 105 , fBackGroundFactorApply(kFALSE)
106 , fRemovePileUp(kFALSE)
107 , fIdentifiedAsPileUp(kFALSE)
108 , fIdentifiedAsOutInz(kFALSE)
109 , fPassTheEventCut(kFALSE)
110 , fRejectKinkMother(kTRUE)
111 , fisppMultiBin(kFALSE)
112 , fPbPbUserCentralityBinning(kFALSE)
113 , fRemoveFirstEvent(kFALSE)
114 , fisNonHFEsystematics(kFALSE)
115 , fSpecialTrigger(NULL)
116 , fCentralityF(-1)
117 , fCentralityPercent(-1)
118 , fCentralityEstimator("V0M")
119 , fContributors(0.5)
120 , fWeightBackGround(0.)
121 , fVz(0.0)
122 , fContainer(NULL)
123 , fVarManager(NULL)
124 , fSignalCuts(NULL)
125 , fCFM(NULL)
126 , fTriggerAnalysis(NULL)
127 , fPID(NULL)
128 , fPIDqa(NULL)
618038b6 129 , fTRDTriggerAnalysis(NULL)
01745870 130 , fPIDpreselect(NULL)
131 , fCuts(NULL)
132 , fTaggedTrackCuts(NULL)
133 , fCleanTaggedTrack(kFALSE)
134 , fVariablesTRDTaggedTrack(kFALSE)
135 , fAnalysisUtils(NULL)
136 , fCutspreselect(NULL)
137 , fSecVtx(NULL)
138 , fElecBackGround(NULL)
139 , fMCQA(NULL)
140 , fTaggedTrackAnalysis(NULL)
141 , fExtraCuts(NULL)
142 , fBackgroundSubtraction(NULL)
143 , fTRDTrigger(kFALSE)
144 , fWhichTRDTrigger(0)
7986e54e 145 , fV0Tagger(NULL)
01745870 146 , fQA(NULL)
147 , fOutput(NULL)
148 , fHistMCQA(NULL)
149 , fHistSECVTX(NULL)
150 , fHistELECBACKGROUND(NULL)
151 , fQACollection(NULL)
152{
153 //
154 // Dummy constructor
155 //
156 memset(fElecBackgroundFactor, 0, sizeof(Double_t) * kElecBgSpecies * kBgPtBins * kCentBins * kBgLevels);
157 memset(fkBackGroundFactorArray, 0, sizeof(TF1 *) * 12);
158 memset(fBinLimit, 0, sizeof(Double_t) * (kBgPtBins+1));
159 memset(&fisppMultiBin, kFALSE, sizeof(fisppMultiBin));
160 memset(fCentralityLimits, 0, sizeof(Float_t) * 12);
161
162 SetppAnalysis();
163}
164
165//____________________________________________________________
166AliAnalysisTaskHFE::AliAnalysisTaskHFE(const char * name):
167 AliAnalysisTaskSE(name)
168 , fAODMCHeader(NULL)
169 , fAODArrayMCInfo(NULL)
170 , fQAlevel(0)
171 , fPlugins(0)
172 , fCollisionSystem(3)
173 , fFillSignalOnly(kTRUE)
174 , fFillNoCuts(kFALSE)
4437a0d2 175 , fApplyCutAOD(kTRUE)
01745870 176 , fBackGroundFactorApply(kFALSE)
177 , fRemovePileUp(kFALSE)
178 , fIdentifiedAsPileUp(kFALSE)
179 , fIdentifiedAsOutInz(kFALSE)
180 , fPassTheEventCut(kFALSE)
181 , fRejectKinkMother(kTRUE)
182 , fisppMultiBin(kFALSE)
183 , fPbPbUserCentralityBinning(kFALSE)
184 , fRemoveFirstEvent(kFALSE)
185 , fisNonHFEsystematics(kFALSE)
186 , fSpecialTrigger(NULL)
187 , fCentralityF(-1)
188 , fCentralityPercent(-1)
189 , fCentralityEstimator("V0M")
190 , fContributors(0.5)
191 , fWeightBackGround(0.)
192 , fVz(0.0)
193 , fContainer(NULL)
194 , fVarManager(NULL)
195 , fSignalCuts(NULL)
196 , fCFM(NULL)
197 , fTriggerAnalysis(NULL)
198 , fPID(NULL)
199 , fPIDqa(NULL)
618038b6 200 , fTRDTriggerAnalysis(NULL)
01745870 201 , fPIDpreselect(NULL)
202 , fCuts(NULL)
203 , fTaggedTrackCuts(NULL)
204 , fCleanTaggedTrack(kFALSE)
205 , fVariablesTRDTaggedTrack(kFALSE)
206 , fAnalysisUtils(NULL)
207 , fCutspreselect(NULL)
208 , fSecVtx(NULL)
209 , fElecBackGround(NULL)
210 , fMCQA(NULL)
211 , fTaggedTrackAnalysis(NULL)
212 , fExtraCuts(NULL)
213 , fBackgroundSubtraction(NULL)
214 , fTRDTrigger(kFALSE)
215 , fWhichTRDTrigger(0)
7986e54e 216 , fV0Tagger(NULL)
01745870 217 , fQA(NULL)
218 , fOutput(NULL)
219 , fHistMCQA(NULL)
220 , fHistSECVTX(NULL)
221 , fHistELECBACKGROUND(NULL)
222 , fQACollection(0x0)
223{
224 //
225 // Default constructor
226 //
227 DefineOutput(1, TList::Class());
228 DefineOutput(2, TList::Class());
229
7986e54e 230 fV0Tagger = new AliHFEV0taginfo("Tagger");
01745870 231 fPID = new AliHFEpid("hfePid");
232 fPIDqa = new AliHFEpidQAmanager;
233 fVarManager = new AliHFEvarManager("hfeVarManager");
234 fAnalysisUtils = new AliAnalysisUtils;
618038b6 235 fTRDTriggerAnalysis = new AliTRDTriggerAnalysis();
01745870 236
237 memset(fElecBackgroundFactor, 0, sizeof(Double_t) * kElecBgSpecies * kBgPtBins * kCentBins * kBgLevels);
238 memset(fkBackGroundFactorArray, 0, sizeof(TF1 *) * 12);
239 memset(fBinLimit, 0, sizeof(Double_t) * (kBgPtBins+1));
240 memset(&fisppMultiBin, kFALSE, sizeof(fisppMultiBin));
241 memset(fCentralityLimits, 0, sizeof(Float_t) * 12);
242
243 SetppAnalysis();
244}
245
246//____________________________________________________________
247AliAnalysisTaskHFE::AliAnalysisTaskHFE(const AliAnalysisTaskHFE &ref):
248 AliAnalysisTaskSE(ref)
249 , fAODMCHeader(NULL)
250 , fAODArrayMCInfo(NULL)
251 , fQAlevel(0)
252 , fPlugins(0)
253 , fCollisionSystem(ref.fCollisionSystem)
254 , fFillSignalOnly(ref.fFillSignalOnly)
255 , fFillNoCuts(ref.fFillNoCuts)
256 , fApplyCutAOD(ref.fApplyCutAOD)
257 , fBackGroundFactorApply(ref.fBackGroundFactorApply)
258 , fRemovePileUp(ref.fRemovePileUp)
259 , fIdentifiedAsPileUp(ref.fIdentifiedAsPileUp)
260 , fIdentifiedAsOutInz(ref.fIdentifiedAsOutInz)
261 , fPassTheEventCut(ref.fPassTheEventCut)
262 , fRejectKinkMother(ref.fRejectKinkMother)
263 , fisppMultiBin(ref.fisppMultiBin)
264 , fPbPbUserCentralityBinning(ref.fPbPbUserCentralityBinning)
265 , fRemoveFirstEvent(ref.fRemoveFirstEvent)
266 , fisNonHFEsystematics(ref.fisNonHFEsystematics)
267 , fSpecialTrigger(ref.fSpecialTrigger)
268 , fCentralityF(ref.fCentralityF)
269 , fCentralityPercent(ref.fCentralityPercent)
270 , fCentralityEstimator(ref.fCentralityEstimator)
271 , fContributors(ref.fContributors)
272 , fWeightBackGround(ref.fWeightBackGround)
273 , fVz(ref.fVz)
274 , fContainer(NULL)
275 , fVarManager(NULL)
276 , fSignalCuts(NULL)
277 , fCFM(NULL)
278 , fTriggerAnalysis(NULL)
279 , fPID(NULL)
280 , fPIDqa(NULL)
618038b6 281 , fTRDTriggerAnalysis(NULL)
01745870 282 , fPIDpreselect(NULL)
283 , fCuts(NULL)
284 , fTaggedTrackCuts(NULL)
285 , fCleanTaggedTrack(ref.fCleanTaggedTrack)
286 , fVariablesTRDTaggedTrack(ref.fVariablesTRDTaggedTrack)
287 , fAnalysisUtils(NULL)
288 , fCutspreselect(NULL)
289 , fSecVtx(NULL)
290 , fElecBackGround(NULL)
291 , fMCQA(NULL)
292 , fTaggedTrackAnalysis(NULL)
293 , fExtraCuts(NULL)
294 , fBackgroundSubtraction(NULL)
295 , fTRDTrigger(ref.fTRDTrigger)
296 , fWhichTRDTrigger(ref.fWhichTRDTrigger)
7986e54e 297 , fV0Tagger(NULL)
01745870 298 , fQA(NULL)
299 , fOutput(NULL)
300 , fHistMCQA(NULL)
301 , fHistSECVTX(NULL)
302 , fHistELECBACKGROUND(NULL)
303 , fQACollection(NULL)
304{
305 //
306 // Copy Constructor
307 //
308 ref.Copy(*this);
309}
310
311//____________________________________________________________
312AliAnalysisTaskHFE &AliAnalysisTaskHFE::operator=(const AliAnalysisTaskHFE &ref){
313 //
314 // Assignment operator
315 //
316 if(this == &ref)
317 ref.Copy(*this);
318 return *this;
319}
320
321//____________________________________________________________
322void AliAnalysisTaskHFE::Copy(TObject &o) const {
323 //
324 // Copy into object o
325 //
326 AliAnalysisTaskHFE &target = dynamic_cast<AliAnalysisTaskHFE &>(o);
327 target.fAODMCHeader = fAODMCHeader;
328 target.fAODArrayMCInfo = fAODArrayMCInfo;
329 target.fQAlevel = fQAlevel;
330 target.fPlugins = fPlugins;
331 target.fCollisionSystem = fCollisionSystem;
332 target.fFillSignalOnly = fFillSignalOnly;
333 target.fFillNoCuts = fFillNoCuts;
334 target.fApplyCutAOD = fApplyCutAOD;
335 target.fBackGroundFactorApply = fBackGroundFactorApply;
336 target.fRemovePileUp = fRemovePileUp;
337 target.fIdentifiedAsPileUp = fIdentifiedAsPileUp;
338 target.fIdentifiedAsOutInz = fIdentifiedAsOutInz;
339 target.fPassTheEventCut = fPassTheEventCut;
340 target.fRejectKinkMother = fRejectKinkMother;
341 target.fisppMultiBin = fisppMultiBin;
342 target.fPbPbUserCentralityBinning = fPbPbUserCentralityBinning;
343 target.fRemoveFirstEvent = fRemoveFirstEvent;
344 target.fisNonHFEsystematics = fisNonHFEsystematics;
345 target.fSpecialTrigger = fSpecialTrigger;
346 target.fCentralityF = fCentralityF;
347 target.fCentralityPercent = fCentralityPercent;
348 target.fCentralityEstimator = fCentralityEstimator;
349 target.fContributors = fContributors;
350 target.fWeightBackGround = fWeightBackGround;
351 target.fVz = fVz;
352 target.fContainer = fContainer;
353 target.fVarManager = fVarManager;
354 target.fSignalCuts = fSignalCuts;
355 target.fCFM = fCFM;
356 target.fTriggerAnalysis = fTriggerAnalysis;
357 target.fPID = fPID;
358 target.fPIDqa = fPIDqa;
618038b6 359 target.fTRDTriggerAnalysis = fTRDTriggerAnalysis;
01745870 360 target.fPIDpreselect = fPIDpreselect;
361 target.fCuts = fCuts;
362 target.fTaggedTrackCuts = fTaggedTrackCuts;
363 target.fCleanTaggedTrack = fCleanTaggedTrack;
364 target.fVariablesTRDTaggedTrack = fVariablesTRDTaggedTrack;
365 target.fAnalysisUtils = fAnalysisUtils;
366 target.fCutspreselect = fCutspreselect;
367 target.fSecVtx = fSecVtx;
368 target.fElecBackGround = fElecBackGround;
369 target.fMCQA = fMCQA;
370 target.fTaggedTrackAnalysis = fTaggedTrackAnalysis;
371 target.fExtraCuts = fExtraCuts;
372 target.fBackgroundSubtraction = fBackgroundSubtraction;
373 target.fTRDTrigger = fTRDTrigger;
374 target.fWhichTRDTrigger = fWhichTRDTrigger;
7986e54e 375 target.fV0Tagger = fV0Tagger;
01745870 376 target.fQA = fQA;
377 target.fOutput = fOutput;
378 target.fHistMCQA = fHistMCQA;
379 target.fHistSECVTX = fHistSECVTX;
380 target.fHistELECBACKGROUND = fHistELECBACKGROUND;
381 target.fQACollection = fQACollection;
382}
383
384//____________________________________________________________
385AliAnalysisTaskHFE::~AliAnalysisTaskHFE(){
386 //
387 // Destructor
388 //
389 if(fPID) delete fPID;
390 if(fPIDpreselect) delete fPIDpreselect;
391 if(fVarManager) delete fVarManager;
618038b6 392 if(fTRDTriggerAnalysis) delete fTRDTriggerAnalysis;
01745870 393 if(fCFM) delete fCFM;
394 if(fTriggerAnalysis) delete fTriggerAnalysis;
395 if(fSignalCuts) delete fSignalCuts;
396 if(fSecVtx) delete fSecVtx;
397 if(fMCQA) delete fMCQA;
398 if(fElecBackGround) delete fElecBackGround;
399 if(fBackgroundSubtraction) delete fBackgroundSubtraction;
400 if(fSpecialTrigger) delete fSpecialTrigger;
401 if(fAnalysisUtils) delete fAnalysisUtils;
7986e54e 402 if(fV0Tagger) delete fV0Tagger;
01745870 403 // Delete output objects only if we are not running in PROOF mode because otherwise this produces a crash during merging
404 AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
405 if(mgr && mgr->GetAnalysisType() != AliAnalysisManager::kProofAnalysis){
406 if(fPIDqa) delete fPIDqa;
407 if(fOutput) delete fOutput;
408 if(fQA) delete fQA;
409 }
410}
411
412//____________________________________________________________
413void AliAnalysisTaskHFE::UserCreateOutputObjects(){
414 //
415 // Creating output container and output objects
416 // Here we also Initialize the correction framework container and
417 // the objects for
418 // - PID
419 // - MC QA
420 // - SecVtx
421 // QA histograms are created if requested
422 // Called once per worker
423 //
424 AliDebug(3, "Creating Output Objects");
425
426 // Make lists for Output
427 if(!fQA) fQA = new TList;
428 fQA->SetOwner();
429 if(!fOutput) fOutput = new TList;
430 fOutput->SetOwner();
431
432 // Automatic determination of the analysis mode
433 AliVEventHandler *inputHandler = dynamic_cast<AliVEventHandler *>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
434 if(!TString(inputHandler->IsA()->GetName()).CompareTo("AliAODInputHandler")){
435 SetAODAnalysis();
436 } else {
437 SetESDAnalysis();
438 if(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler())
439 SetHasMCData();
440 }
441 printf("Analysis Mode: %s Analysis\n", IsAODanalysis() ? "AOD" : "ESD");
442 printf("MC Data available %s\n", HasMCData() ? "Yes" : "No");
443
444 // Enable Trigger Analysis
445 fTriggerAnalysis = new AliTriggerAnalysis;
446 fTriggerAnalysis->EnableHistograms();
447 fTriggerAnalysis->SetAnalyzeMC(HasMCData());
448
449 // First Part: Make QA histograms
450 fQACollection = new AliHFEcollection("TaskQA", "QA histos from the Electron Task");
451 fQACollection->CreateTH1F("nElectronTracksEvent", "Number of Electron Candidates", 100, 0, 100);
452 fQACollection->CreateTH1F("nElectron", "Number of electrons", 100, 0, 100);
453 fQACollection->CreateTH2F("radius", "Production Vertex", 100, 0.0, 5.0, 100, 0.0, 5.0);
454 fQACollection->CreateTH1F("nTriggerBit", "Histo Trigger Bit", 22, 0, 22);
e2861c1a 455 fQACollection->CreateTH1F("Filterbegin", "AOD filter of tracks after all cuts", 21, -1, 20);
47faab26 456 fQACollection->CreateTH1F("Filterend", "AOD filter of tracks after all cuts", 21, -1, 20);
01745870 457
458 InitHistoITScluster();
459 InitContaminationQA();
460 fQA->Add(fQACollection);
461
462 // Initialize PID
463 fPID->SetHasMCData(HasMCData());
464 if(!fPID->GetNumberOfPIDdetectors()) fPID->AddDetector("TPC", 0);
465 if(IsQAOn(kPIDqa)){
466 AliInfo("PID QA switched on");
467 fPIDqa->Initialize(fPID);
468 fQA->Add(fPIDqa->MakeList("HFEpidQA"));
469 }
470 fPID->SortDetectors();
471
472 // Background subtraction-------------------------------------------------------------------
473 if (GetPlugin(kNonPhotonicElectron)) {
474 if(!fBackgroundSubtraction) fBackgroundSubtraction = new AliHFENonPhotonicElectron();
475 if(IsAODanalysis()) fBackgroundSubtraction->SetAOD(kTRUE);
476 fBackgroundSubtraction->Init();
477 fOutput->Add(fBackgroundSubtraction->GetListOutput());
478 }
479 //------------------------------------------------------------------------------------------
480
481
482 // Initialize correction Framework and Cuts
483 const Int_t kNcutSteps = AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kNcutStepsRecTrack + AliHFEcuts::kNcutStepsDETrack + AliHFEcuts::kNcutStepsSecvtxTrack;
484 fCFM = new AliCFManager;
485 fCFM->SetNStepParticle(kNcutSteps);
486 MakeParticleContainer();
487 MakeEventContainer();
488 // Temporary fix: Initialize particle cuts with NULL
489 for(Int_t istep = 0; istep < kNcutSteps; istep++)
490 fCFM->SetParticleCutsList(istep, NULL);
491 if(!fCuts){
492 AliWarning("Cuts not available. Default cuts will be used");
493 fCuts = new AliHFEcuts;
494 fCuts->CreateStandardCuts();
495 }
496 if(IsAODanalysis()) fCuts->SetAOD();
497 // Make clone for V0 tagging step
498 fCuts->Initialize(fCFM);
499 if(fCuts->IsQAOn()) fQA->Add(fCuts->GetQAhistograms());
500 fSignalCuts = new AliHFEsignalCuts("HFEsignalCuts", "HFE MC Signal definition");
501 fVarManager->SetSignalCuts(fSignalCuts);
502
503 // add output objects to the List
504 fOutput->AddAt(fContainer, 0);
505 fOutput->AddAt(fCFM->GetEventContainer(), 1);
506
507 // mcQA----------------------------------
508 if (HasMCData() && IsQAOn(kMCqa)) {
509 AliInfo("MC QA on");
510 if(!fMCQA) fMCQA = new AliHFEmcQA;
511 if(!fHistMCQA) fHistMCQA = new TList();
512 fHistMCQA->SetOwner();
513 if(IsPbPb()) fMCQA->SetPbPb();
514 if(fisppMultiBin) fMCQA->SetPPMultiBin();
515 if(TestBit(kTreeStream)){
516 fMCQA->EnableDebugStreamer();
517 }
518 fMCQA->CreatDefaultHistograms(fHistMCQA);
519 fMCQA->SetBackgroundWeightFactor(fElecBackgroundFactor[0][0][0],fBinLimit);
520 fQA->Add(fHistMCQA);
521 }
522
523 // secvtx----------------------------------
524 if (GetPlugin(kSecVtx)) {
525 AliInfo("Secondary Vertex Analysis on");
526 if(!fSecVtx) fSecVtx = new AliHFEsecVtx;
527 fSecVtx->SetHasMCData(HasMCData());
528
529 if(!fHistSECVTX) fHistSECVTX = new TList();
530 fHistSECVTX->SetOwner();
531 fSecVtx->CreateHistograms(fHistSECVTX);
532 fOutput->Add(fHistSECVTX);
533 }
534
535 // background----------------------------------
536 if (GetPlugin(kIsElecBackGround)) {
537 AliInfo("Electron BackGround Analysis on");
538 if(!fElecBackGround){
539 AliWarning("ElecBackGround not available. Default elecbackground will be used");
540 fElecBackGround = new AliHFEelecbackground;
541 }
542 fElecBackGround->SetHasMCData(HasMCData());
543
544 if(!fHistELECBACKGROUND) fHistELECBACKGROUND = new TList();
545 fHistELECBACKGROUND->SetOwner();
546 fElecBackGround->CreateHistograms(fHistELECBACKGROUND);
547 fOutput->Add(fHistELECBACKGROUND);
548 }
549
550 // tagged tracks
551 if(GetPlugin(kTaggedTrackAnalysis)){
552 AliInfo("Analysis on V0-tagged tracks enabled");
553 fTaggedTrackAnalysis = new AliHFEtaggedTrackAnalysis(Form("taggedTrackAnalysis%s", GetName()));
554 fTaggedTrackAnalysis->SetCuts(fTaggedTrackCuts);
555 fTaggedTrackAnalysis->SetClean(fCleanTaggedTrack);
556 AliHFEvarManager *varManager = fTaggedTrackAnalysis->GetVarManager();
557 TObjArray *array = fVarManager->GetVariables();
558 Int_t nvars = array->GetEntriesFast();
559 TString namee;
560 for(Int_t v = 0; v < nvars; v++) {
561 AliHFEvarManager::AliHFEvariable *variable = (AliHFEvarManager::AliHFEvariable *) array->At(v);
562 if(!variable) continue;
563 TString name(((AliHFEvarManager::AliHFEvariable *)variable)->GetName());
564 if(!name.CompareTo("source")) namee = TString("species");
565 else namee = TString(name);
566 Int_t nbins = variable->GetNumberOfBins();
567 if(variable->HasUserDefinedBinning()){
568 varManager->AddVariable(namee, nbins, variable->GetBinning());
569 } else {
570 varManager->AddVariable(namee, nbins, variable->GetMinimum(), variable->GetMaximum());
571 }
572 //printf("For AliTaggedTrackAnalysis, had the variable %s and the one used %s\n",(const char*)variable->GetName(),(const char*) namee);
573 }
574 if(fPIDqa->HasHighResolutionHistos())
575 fTaggedTrackAnalysis->GetPIDqa()->SetHighResolutionHistos();
576 fTaggedTrackAnalysis->SetPID(fPID);
577 fTaggedTrackAnalysis->SetVariablesTRD(fVariablesTRDTaggedTrack);
578 fTaggedTrackAnalysis->InitContainer();
579 fOutput->Add(fTaggedTrackAnalysis->GetContainer());
580 fQA->Add(fTaggedTrackAnalysis->GetPIDQA());
581 fQA->Add(fTaggedTrackAnalysis->GetCutQA());
582 fQA->Add(fTaggedTrackAnalysis->GetQAcollection());
583 }
584
585 //fQA->Print();
586
587 PrintStatus();
588 // Done!!!
589 PostData(1, fOutput);
590 PostData(2, fQA);
591}
592
593//____________________________________________________________
594void AliAnalysisTaskHFE::UserExec(Option_t *){
595 //
596 // Run the analysis
597 //
598
599 //printf("test00\n");
600
601 AliDebug(3, "Starting Single Event Analysis");
602 if(!fInputEvent){
603 AliError("Reconstructed Event not available");
604 //printf("Reconstructed Event not available");
605 return;
606 }
607 if(HasMCData() && IsESDanalysis()){
608 AliDebug(4, Form("MC Event: %p", fMCEvent));
609 if(!fMCEvent){
610 AliError("No MC Event, but MC Data required");
611 //printf("No MC Event, but MC Data required");
612 return;
613 }
614 }
615 if(!fCuts){
616 AliError("HFE cuts not available");
617 //printf("HFE cuts not available");
618 return;
619 }
620 if(!fPID->IsInitialized()){
621 // Initialize PID with the given run number
622 fPID->InitializePID(fInputEvent->GetRunNumber());
623 }
624
625 if(fRemoveFirstEvent){
626 if(fAnalysisUtils->IsFirstEventInChunk(fInputEvent)) return;
627 }
628
629 AliESDEvent *ev = dynamic_cast<AliESDEvent *>(fInputEvent);
4437a0d2 630 if(ev && fTRDTrigger && (fWhichTRDTrigger<6))
01745870 631 {
4437a0d2 632 if(!CheckTRDTriggerESD(ev)) return;
633 }
634 if(fInputEvent && fTRDTrigger && (fWhichTRDTrigger>5))
635 {
636 if(!CheckTRDTrigger(fInputEvent)) return;
01745870 637 }
01745870 638
639 if(IsESDanalysis() && HasMCData()){
640 // Protect against missing MC trees
641 AliMCEventHandler *mcH = dynamic_cast<AliMCEventHandler *>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
642 if(!mcH){
643 AliError("No MC Event Handler available");
644 return;
645 }
646 if(!mcH->InitOk()) return;
647 if(!mcH->TreeK()) return;
648 if(!mcH->TreeTR()) return;
649
650 // Background subtraction-------------------------------------------------------------------
651 if(GetPlugin(kNonPhotonicElectron)) fBackgroundSubtraction->SetMCEvent(fMCEvent);
652 //------------------------------------------------------------------------------------------
653 }
654
655 if(IsAODanalysis() && HasMCData()){
656 // take MC info
657 AliAODEvent *aodE = dynamic_cast<AliAODEvent *>(fInputEvent);
658 if(!aodE){
659 AliError("No AOD Event");
660 return;
661 }
662 fAODMCHeader = dynamic_cast<AliAODMCHeader *>(fInputEvent->FindListObject(AliAODMCHeader::StdBranchName()));
663 if(!fAODMCHeader){
664 AliError("No AliAODMCHeader");
665 //printf("No AliAODMCHeader");
666 return;
667 }
668 fAODArrayMCInfo = dynamic_cast<TClonesArray *>(fInputEvent->FindListObject(AliAODMCParticle::StdBranchName()));
669 if(!fAODArrayMCInfo){
670 AliError("No AOD MC particles");
671 //printf("No AOD MC particles");
672 return;
673 }
674 fSignalCuts->SetMCAODInfo(fAODArrayMCInfo);
675 // Background subtraction-------------------------------------------------------------------
676 if (GetPlugin(kNonPhotonicElectron)) fBackgroundSubtraction->SetAODArrayMCInfo(fAODArrayMCInfo);
677 //------------------------------------------------------------------------------------------
678 }
679
680 //printf("test2\n");
681
682 // need the centrality for everything (MC also)
683 fCentralityF = -1;
684 if(!ReadCentrality()) fCentralityF = -1;
685 //printf("pass centrality\n");
686 //printf("Reading fCentralityF %d\n",fCentralityF);
687
688 // See if pile up and z in the range
689 RejectionPileUpVertexRangeEventCut();
690
691 //printf("test3\n");
692
693 // Protect agains missing
694 if(HasMCData()){
695 //printf("Has MC data\n");
696 fSignalCuts->SetMCEvent(fMCEvent);
697 ProcessMC(); // Run the MC loop + MC QA in case MC Data are available
698 }
699
700 AliPIDResponse *pidResponse = fInputHandler->GetPIDResponse();
701 if(!pidResponse){
702 AliDebug(1, "Using default PID Response");
703 pidResponse = AliHFEtools::GetDefaultPID(HasMCData(), fInputEvent->IsA() == AliESDEvent::Class());
704 }
705 fPID->SetPIDResponse(pidResponse);
706 if(fPIDpreselect) fPIDpreselect->SetPIDResponse(pidResponse);
707
708 // Background subtraction-------------------------------------------------------------------
709 if(GetPlugin(kNonPhotonicElectron)) fBackgroundSubtraction->InitRun(fInputEvent,pidResponse);
710 //------------------------------------------------------------------------------------------
711
712 // Event loop
713 if(IsAODanalysis()){
714 //printf("test4\n");
715 ProcessAOD();
716 } else {
717 const char *specialTrigger = GetSpecialTrigger(fInputEvent->GetRunNumber());
718 // Check Trigger selection
719 if(specialTrigger){
720 AliDebug(2, Form("Special Trigger requested: %s", specialTrigger));
721 if(!(ev && ev->IsTriggerClassFired(specialTrigger))){
722 AliDebug(2, "Event not selected");
723 return;
724 } else AliDebug(2, "Event Selected");
725 } else AliDebug(2, "No Special Trigger requested");
726
727 ProcessESD();
728 }
729 // Done!!!
730 PostData(1, fOutput);
731 PostData(2, fQA);
732}
733
734//____________________________________________________________
735void AliAnalysisTaskHFE::Terminate(Option_t *){
736 //
737 // Terminate not implemented at the moment
738 //
739}
740
741//_______________________________________________________________
742Bool_t AliAnalysisTaskHFE::IsEventInBinZero() {
743 //
744 //
745 //
746
747 //printf("test in IsEventInBinZero\n");
748 if(!fInputEvent){
749 AliError("Reconstructed Event not available");
750 return kFALSE;
751 }
752
753 // check vertex
754 const AliVVertex *vertex = fInputEvent->GetPrimaryVertex();
755 if(!vertex) return kTRUE;
756 //if(vertex) return kTRUE;
757
758 // check tracks
759 if(fInputEvent->GetNumberOfTracks()<=0) return kTRUE;
760 //if(fInputEvent->GetNumberOfTracks()>0) return kTRUE;
761
762
763 return kFALSE;
764
765}
766//____________________________________________________________
767void AliAnalysisTaskHFE::ProcessMC(){
768 //
769 // Runs the MC Loop (filling the container for the MC Cut Steps with the observables pt, eta and phi)
770 // In case MC QA is on also MC QA loop is done
771 //
772 AliDebug(3, "Processing MC Information");
773 Double_t eventContainer [4] = {0., 0., 0., 0.};
774 if(IsESDanalysis()) eventContainer[0] = fMCEvent->GetPrimaryVertex()->GetZ();
775 else eventContainer[0] = fAODMCHeader->GetVtxZ();
776 eventContainer[2] = fCentralityF;
777 eventContainer[3] = fContributors;
778 fVz = eventContainer[0];
779 //printf("z position is %f\n",eventContainer[0]);
780 //if(fCFM->CheckEventCuts(AliHFEcuts::kEventStepGenerated, fMCEvent))
781 fCFM->GetEventContainer()->Fill(eventContainer,AliHFEcuts::kEventStepGenerated);
782 Int_t nElectrons = 0;
783 if(IsESDanalysis()){
784 if(!((fIdentifiedAsPileUp) || (TMath::Abs(fVz) > fCuts->GetVertexRange()) || (fCentralityF < 0))){ //kStepMCGeneratedZOutNoPileUpCentralityFine
785 if (HasMCData() && IsQAOn(kMCqa)) {
786 AliDebug(2, "Running MC QA");
787
788 if(fMCEvent->Stack()){
789 fMCQA->SetMCEvent(fMCEvent);
790 fMCQA->SetGenEventHeader(fMCEvent->GenEventHeader());
791 fMCQA->SetCentrality(fCentralityF);
792 fMCQA->SetPercentrality(static_cast<Int_t>(fCentralityPercent));
793 if(IsPbPb()) { fMCQA->SetPbPb();}
794 else
795 {
796 if(fisppMultiBin) fMCQA->SetPPMultiBin();
797 else fMCQA->SetPP();
798 }
799 fMCQA->Init();
800
801 fMCQA->GetMesonKine();
802
803 // loop over all tracks for decayed electrons
804 for (Int_t igen = 0; igen < fMCEvent->GetNumberOfTracks(); igen++){
805 TParticle* mcpart = fMCEvent->Stack()->Particle(igen);
806 if(!mcpart) continue;
807 fMCQA->GetQuarkKine(mcpart, igen, AliHFEmcQA::kCharm);
808 fMCQA->GetQuarkKine(mcpart, igen, AliHFEmcQA::kBeauty);
809 fMCQA->GetHadronKine(mcpart, AliHFEmcQA::kCharm);
810 fMCQA->GetHadronKine(mcpart, AliHFEmcQA::kBeauty);
811 fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kCharm, AliHFEmcQA::kElectronPDG); // no accept cut
812 fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kBeauty, AliHFEmcQA::kElectronPDG); // no accept cut
813 fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kOthers, AliHFEmcQA::kElectronPDG); // no accept cut
814 }
815 //fMCQA->EndOfEventAna(AliHFEmcQA::kCharm);
816 //fMCQA->EndOfEventAna(AliHFEmcQA::kBeauty);
817 }
818
819 } // end of MC QA loop
820 }
821 // -----------------------------------------------------------------
822 fCFM->SetMCEventInfo(fMCEvent);
823 // fCFM->CheckEventCuts(AliCFManager::kEvtRecCuts, fESD);
824 } else {
825 fCFM->SetMCEventInfo(fInputEvent);
826 }
827 // Run MC loop
828 AliVParticle *mctrack = NULL;
829 Int_t numberofmctracks = 0;
830 if(IsESDanalysis()){
831 numberofmctracks = fMCEvent->GetNumberOfTracks();
832 }
833 else {
834 numberofmctracks = fAODArrayMCInfo->GetEntriesFast();
835 }
836 AliDebug(3, Form("Number of Tracks: %d",numberofmctracks));
837 //printf("Number of MC track %d\n",numberofmctracks);
838 for(Int_t imc = 0; imc <numberofmctracks; imc++){
839 if(IsESDanalysis()) {
840 if(!(mctrack = fMCEvent->GetTrack(imc))) continue;
841 }
842 else {
843 if(!(mctrack = (AliVParticle *) fAODArrayMCInfo->At(imc))) continue;
844 }
845 //printf("Test in ProcessMC\n");
846 AliDebug(4, "Next MC Track");
847 if(ProcessMCtrack(mctrack)) nElectrons++;
848 }
849
850 // fCFM->CheckEventCuts(AliCFManager::kEvtRecCuts, fESD);
851 fQACollection->Fill("nElectron", nElectrons);
852}
853
854//____________________________________________________________
855void AliAnalysisTaskHFE::ProcessESD(){
856 //
857 // Run Analysis of reconstructed event in ESD Mode
858 // Loop over Tracks, filter according cut steps defined in AliHFEcuts
859 //
860 AliDebug(1, Form("Task %s", GetName()));
861 AliDebug(3, "Processing ESD Event");
862 AliESDEvent *fESD = dynamic_cast<AliESDEvent *>(fInputEvent);
863 if(!fESD){
864 AliError("ESD Event required for ESD Analysis");
865 return;
866 }
867
7986e54e 868 // Tag all v0s in current event
869 if(fV0Tagger){
870 fV0Tagger->Reset();
871 fV0Tagger->TagV0Tracks(fESD);
872 }
01745870 873 // Set magnetic field if V0 task on
874 if(fTaggedTrackAnalysis) {
875 fTaggedTrackAnalysis->SetMagneticField(fESD->GetMagneticField());
876 fTaggedTrackAnalysis->SetCentrality(fCentralityF);
4437a0d2 877 if(IsPbPb()) fTaggedTrackAnalysis->SetPbPb();
878 else {
879 if(IspPb()) fTaggedTrackAnalysis->SetpPb();
880 else fTaggedTrackAnalysis->SetPP();
881 }
882
01745870 883 }
884
885 // Do event Normalization
886 Double_t eventContainer[4];
887 eventContainer[0] = 0.;
888 if(HasMCData()) eventContainer[0] = fVz;
889 else {
890 const AliESDVertex* vtxESD = fESD->GetPrimaryVertex();
891 if(vtxESD) eventContainer[0] = vtxESD->GetZ();
892 }
893 eventContainer[1] = 0.;
894 eventContainer[2] = fCentralityF;
895 eventContainer[3] = fContributors;
896 if(fTriggerAnalysis->IsOfflineTriggerFired(fESD, AliTriggerAnalysis::kV0AND))
897 eventContainer[1] = 1.;
898
899 //
900 fCFM->GetEventContainer()->Fill(eventContainer, AliHFEcuts::kEventStepRecNoCut);
901
902 //
903 if(fIdentifiedAsPileUp) return;
904 fCFM->GetEventContainer()->Fill(eventContainer, AliHFEcuts::kEventStepRecNoPileUp);
905
906 //
907 if(TMath::Abs(fCentralityF) < 0) return;
908 fCFM->GetEventContainer()->Fill(eventContainer, AliHFEcuts::kEventStepRecCentralityOk);
909 //printf("In ProcessESD %f\n",fCentralityF);
910
911 //
912 if(fIdentifiedAsOutInz) return;
913 fCFM->GetEventContainer()->Fill(eventContainer, AliHFEcuts::kEventStepZRange);
914
915 //
916 if(!fPassTheEventCut) return;
917 fCFM->GetEventContainer()->Fill(eventContainer, AliHFEcuts::kEventStepReconstructed);
918
919
920 fContainer->NewEvent();
921
922 if (GetPlugin(kIsElecBackGround)) {
923 fElecBackGround->SetEvent(fESD);
924 }
925 if (GetPlugin(kSecVtx)) {
926 fSecVtx->SetEvent(fESD);
927 fSecVtx->GetPrimaryCondition();
928 }
929
930 if(HasMCData()){
931 if (GetPlugin(kSecVtx)) {
932 fSecVtx->SetMCEvent(fMCEvent);
933 fSecVtx->SetMCQA(fMCQA);
934 }
935 if (GetPlugin(kIsElecBackGround)) {
936 fElecBackGround->SetMCEvent(fMCEvent);
937 }
938 }
939
940 Double_t container[10];
941 memset(container, 0, sizeof(Double_t) * 10);
942 // container for the output THnSparse
943 Double_t dataDca[6]; // [source, pT, dca, centrality]
944 Int_t nElectronCandidates = 0;
945 AliESDtrack *track = NULL, *htrack = NULL;
946 AliMCParticle *mctrack = NULL;
947 AliMCParticle *mctrackmother = NULL;
948
949 Bool_t signal = kTRUE;
950
951 fCFM->SetRecEventInfo(fESD);
952
953 // Get Number of contributors to the primary vertex for multiplicity-dependent correction
954 Int_t ncontribVtx = 0;
955 const AliESDVertex *priVtx = fESD->GetPrimaryVertexTracks();
956 if(priVtx){
957 ncontribVtx = priVtx->GetNContributors();
958 }
959
960 // minjung for IP QA(temporary ~ 2weeks)
961 if(!fExtraCuts){
962 fExtraCuts = new AliHFEextraCuts("hfetmpCuts","HFE tmp Cuts");
963 }
964 fExtraCuts->SetRecEventInfo(fESD);
965
966 // Electron background analysis
967 if (GetPlugin(kIsElecBackGround)) {
968
969 AliDebug(2, "Running BackGround Analysis");
970
971 fElecBackGround->Reset();
972
973 } // end of electron background analysis
974
975
976 // Background subtraction-------------------------------------------------------------------
977 if (GetPlugin(kNonPhotonicElectron)) fBackgroundSubtraction->FillPoolAssociatedTracks(fInputEvent, fCentralityF);
978 //------------------------------------------------------------------------------------------
979
980 //
981 // Loop ESD
982 //
983 AliDebug(3, Form("Number of Tracks: %d", fESD->GetNumberOfTracks()));
984 for(Int_t itrack = 0; itrack < fESD->GetNumberOfTracks(); itrack++){
985 AliDebug(4, "New ESD track");
986 track = fESD->GetTrack(itrack);
987 track->SetESDEvent(fESD);
4437a0d2 988
01745870 989 // fill counts of v0-identified particles
dab4f95d 990 AliPID::EParticleType v0pid = fV0Tagger ? fV0Tagger->GetV0Info(track->GetID()) : AliPID::kUnknown;
01745870 991 // here the tagged track analysis will run
7986e54e 992 if(fTaggedTrackAnalysis && v0pid != AliPID::kUnknown){
01745870 993 AliDebug(1, Form("Track identified as %s", AliPID::ParticleName(v0pid)));
994 fTaggedTrackAnalysis->ProcessTrack(track, v0pid);
995 AliDebug(1, "V0 PID done");
996 }
997
998
999 //Fill non-HFE source containers at reconstructed events cut step
1000 AliDebug(3, Form("Doing track %d, %p", itrack, track));
1001
1002
1003 //////////////////////////////////////
1004 // preselect
1005 //////////////////////////////////////
1006 if(fPIDpreselect || fCutspreselect) {
1007 if(!PreSelectTrack(track)) continue;
1008 }
1009
1010 signal = kTRUE;
1011
1012 // Fill step without any cut
1013
1014 if(HasMCData()){
1015 // Check if it is electrons near the vertex
1016 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(track->GetLabel()))))) continue;
1017
1018 if(fFillSignalOnly && !fCFM->CheckParticleCuts(AliHFEcuts::kStepMCGenerated, mctrack)) signal = kFALSE;
1019 else AliDebug(3, "Signal Electron");
1020
1021 // Fill K pt for Ke3 contributions
1022 if(mctrack && (TMath::Abs(mctrack->Particle()->GetPdgCode())==321)) fQACollection->Fill("Kptspectra",mctrack->Pt());
1023 else if(mctrack && (TMath::Abs(mctrack->Particle()->GetPdgCode())==130)) fQACollection->Fill("K0Lptspectra",mctrack->Pt());
1024 }
1025 // Cache new Track information inside the var manager
1026 fVarManager->NewTrack(track, mctrack, fCentralityF, -1, signal);
1027
1028 if(fFillNoCuts) {
1029 if(signal || !fFillSignalOnly){
1030 fVarManager->FillContainer(fContainer, "recTrackContReco", AliHFEcuts::kStepRecNoCut, kFALSE);
1031 fVarManager->FillContainer(fContainer, "recTrackContMC", AliHFEcuts::kStepRecNoCut, kTRUE);
1032 }
1033 }
1034
1035 // RecKine: ITSTPC cuts
1036 if(!ProcessCutStep(AliHFEcuts::kStepRecKineITSTPC, track)) continue;
1037
1038
1039 // RecPrim
1040 if(fRejectKinkMother) {
1041 if(track->GetKinkIndex(0) != 0) continue; } // Quick and dirty fix to reject both kink mothers and daughters
1042 if(!ProcessCutStep(AliHFEcuts::kStepRecPrim, track)) continue;
1043
1044 // HFEcuts: ITS layers cuts
1045 if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsITS, track)) continue;
1046
1047 // HFE cuts: TOF PID and mismatch flag
1048 if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTOF, track)) continue;
1049
1050 // HFE cuts: TPC PID cleanup
1051 if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTPC, track)) continue;
1052
1053 // HFEcuts: Nb of tracklets TRD0
1054 if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTRD, track)) continue;
1055
1056 // Fill correlation maps before PID
1057 if(signal && fContainer->GetCorrelationMatrix("correlationstepbeforePID")) {
1058 //printf("Fill correlation maps before PID\n");
1059 fVarManager->FillCorrelationMatrix(fContainer->GetCorrelationMatrix("correlationstepbeforePID"));
1060 }
1061
1062 if(HasMCData()){
1063 FillProductionVertex(track);
1064
1065 if(fMCQA && signal){
1066 fMCQA->SetCentrality(fCentralityF);
1067 if(mctrack && (TMath::Abs(mctrack->Particle()->GetPdgCode()) == 11)){
1068 Double_t weightElecBgV0[kBgLevels] = {0.,0.,0.};
1069 Double_t hfeimpactRtmp=0., hfeimpactnsigmaRtmp=0.;
1070 fExtraCuts->GetHFEImpactParameters(track, hfeimpactRtmp, hfeimpactnsigmaRtmp);
1071 UChar_t itsPixel = track->GetITSClusterMap();
1072 Double_t ilyrhit=0, ilyrstat=0;
1073 for(Int_t ilyr=0; ilyr<6; ilyr++){
1074 if(TESTBIT(itsPixel, ilyr)) ilyrhit += TMath::Power(2,ilyr);
1075 if(fExtraCuts->CheckITSstatus(fExtraCuts->GetITSstatus(track,ilyr))) ilyrstat += TMath::Power(2,ilyr);
1076 }
1077 fMCQA->SetITSInfo(ilyrhit,ilyrstat);
1078 fMCQA->SetHFEImpactParameters(hfeimpactRtmp, hfeimpactnsigmaRtmp);
1079 fMCQA->SetTrkKine(track->Pt(),track->Eta(), track->Phi());
1080 fMCQA->SetContainerStep(3);
1081 for(Int_t iLevel = 0; iLevel < kBgLevels; iLevel++){
1082 weightElecBgV0[iLevel] = fMCQA->GetWeightFactor(mctrack, iLevel); // positive:conversion e, negative: nonHFE
1083 if(!fisNonHFEsystematics || IsPbPb())break;
1084 }
1085
1086 if(fisNonHFEsystematics){
1087 //Fill additional containers for electron source distinction
1088 Int_t elecSource = 0;
1089 elecSource = fMCQA->GetElecSource(mctrack->Particle());
1090 const Char_t *sourceName[kElecBgSpecies]={"Pion","Eta","Omega","Phi","EtaPrime","Rho"};
1091 const Char_t *levelName[kBgLevels]={"Best","Lower","Upper"};
1092 Int_t iName = 0;
1093 for(Int_t iSource = AliHFEmcQA::kPi0; iSource <= AliHFEmcQA::kGammaRho0; iSource++){
1094 if((iSource == AliHFEmcQA::kElse)||(iSource == AliHFEmcQA::kMisID)) continue;
1095 if(elecSource == iSource){
1096 for(Int_t iLevel = 0; iLevel < kBgLevels; iLevel++){
1097 if(weightElecBgV0[iLevel]>0){
1098 fVarManager->FillContainer(fContainer, Form("conversionElecs%s%s",sourceName[iName], levelName[iLevel]), 3, kFALSE, weightElecBgV0[iLevel]);
1099 }
1100 else if(weightElecBgV0[iLevel]<0){
1101 fVarManager->FillContainer(fContainer, Form("mesonElecs%s%s",sourceName[iName], levelName[iLevel]), 3, kFALSE, -1*weightElecBgV0[iLevel]);
1102 }
1103 if(IsPbPb())break;
1104 }
1105 break;
1106 }
1107 iName++;
1108 if(iName == kElecBgSpecies)iName = 0;
1109 }
1110 }
1111 //else{
1112 if(weightElecBgV0[0]>0) {
1113 fVarManager->FillContainer(fContainer, "conversionElecs", 3, kFALSE, weightElecBgV0[0]);
1114 fVarManager->FillContainer(fContainer, "conversionElecs", 4, kTRUE, weightElecBgV0[0]);
1115 }
1116 else if(weightElecBgV0[0]<0) {
1117 fVarManager->FillContainer(fContainer, "mesonElecs", 3, kFALSE, -1*weightElecBgV0[0]);
1118 fVarManager->FillContainer(fContainer, "mesonElecs", 4, kTRUE, -1*weightElecBgV0[0]);
1119 }
1120 //}
1121 }
1122 }
1123
1124 Double_t hfeimpactR4all=0., hfeimpactnsigmaR4all=0.;
1125 Int_t sourceDca =-1;
1126 if(mctrack && (TMath::Abs(mctrack->Particle()->GetPdgCode()) == 211)){
1127 if(track->Pt()>4.){
1128 fExtraCuts->GetHFEImpactParameters(track, hfeimpactR4all, hfeimpactnsigmaR4all);
1129 dataDca[0]=0; //pion
1130 dataDca[1]=track->Pt();
1131 dataDca[2]=hfeimpactR4all;
1132 dataDca[3]=fCentralityF;
1133 dataDca[4] = v0pid;
1134 dataDca[5] = double(track->Charge());
1135 fQACollection->Fill("Dca", dataDca);
1136 }
1137 }
1138 else if(mctrack && (TMath::Abs(mctrack->Particle()->GetPdgCode()) == 11)){ // to increas statistics for Martin
1139 if(signal){
1140 fExtraCuts->GetHFEImpactParameters(track, hfeimpactR4all, hfeimpactnsigmaR4all);
1141 if(fSignalCuts->IsCharmElectron(track)){
1142 sourceDca=1;
1143 }
1144 else if(fSignalCuts->IsBeautyElectron(track)){
1145 sourceDca=2;
1146 }
1147 else if(fSignalCuts->IsGammaElectron(track)){
1148 sourceDca=3;
1149 }
1150 else if(fSignalCuts->IsNonHFElectron(track)){
1151 sourceDca=4;
1152 }
1153 else if(fSignalCuts->IsJpsiElectron(track)){
1154 sourceDca=5;
1155 }
1156 else {
1157 sourceDca=6;
1158 }
1159 dataDca[0]=sourceDca;
1160 dataDca[1]=track->Pt();
1161 dataDca[2]=hfeimpactR4all;
1162 dataDca[3]=fCentralityF;
1163 dataDca[4] = v0pid;
1164 dataDca[5] = double(track->Charge());
1165 if(signal) fQACollection->Fill("Dca", dataDca);
1166 }
1167 }
1168 }
1169
1170 AliHFEpidObject hfetrack;
1171 hfetrack.SetAnalysisType(AliHFEpidObject::kESDanalysis);
1172 hfetrack.SetRecTrack(track);
1173 if(HasMCData()) hfetrack.SetMCTrack(mctrack);
1174 hfetrack.SetCentrality(fCentralityF);
1175 hfetrack.SetMulitplicity(ncontribVtx);
4437a0d2 1176 if(IsPbPb()) hfetrack.SetPbPb();
1177 else {
1178 if(IspPb()) hfetrack.SetpPb();
1179 else hfetrack.SetPP();
1180 }
01745870 1181 fPID->SetVarManager(fVarManager);
1182 if(!fPID->IsSelected(&hfetrack, fContainer, "recTrackCont", fPIDqa)) continue;
1183 nElectronCandidates++;
1184
1185 // Background subtraction------------------------------------------------------------------------------------------
1186 if (GetPlugin(kNonPhotonicElectron)) {
1187 Int_t indexmother = -1;
1188 Int_t mcsource = -1;
1189 if(HasMCData()){
1190 mcsource = fBackgroundSubtraction->FindMother(mctrack->GetLabel(),indexmother);
1191 }
1192 fBackgroundSubtraction->LookAtNonHFE(itrack, track, fInputEvent, 1, fCentralityF, -1, mcsource, indexmother);
1193 }
1194 //-----------------------------------------------------------------------------------------------------------------
1195
1196 // Temporary histogram for chi2/ITS cluster
1197 if(IsPbPb()) {
1198 TBits shared = track->GetTPCSharedMap();
1199 Int_t sharebit=0;
1200 if(shared.CountBits() >= 2) sharebit=1;
1201
1202 Double_t itschi2percluster = 0.0;
1203 Double_t itsnbcls = static_cast<Double_t>(track->GetNcls(0));
1204 if(itsnbcls > 0) itschi2percluster = track->GetITSchi2()/itsnbcls;
1205
1206 Double_t itsChi2[7] = {track->Pt(),track->Eta(), track->Phi(),
1207 static_cast<Double_t>(fCentralityF),static_cast<Double_t>(track->GetTPCsignalN()), static_cast<Double_t>(sharebit),itschi2percluster};
1208 fQACollection->Fill("fChi2perITScluster", itsChi2);
1209 }
1210 else{
1211
1212 Double_t itschi2percluster = 0.0;
1213 Double_t itsnbcls = static_cast<Double_t>(track->GetNcls(0));
1214 if(itsnbcls > 0) itschi2percluster = track->GetITSchi2()/itsnbcls;
1215
1216 Double_t itsChi2[3] = {track->Pt(), static_cast<Double_t>(fCentralityF), itschi2percluster};
1217 fQACollection->Fill("fChi2perITScluster", itsChi2);
1218 }
1219
1220 // Fill Histogram for Hadronic Background
1221 if(HasMCData()){
1222 if(mctrack && (TMath::Abs(mctrack->Particle()->GetPdgCode()) != 11))
1223 fVarManager->FillContainer(fContainer, "hadronicBackground", UInt_t(0), kFALSE);
1224 else if(mctrack){
1225 // Fill Ke3 contributions
1226 Int_t glabel=TMath::Abs(mctrack->GetMother());
1227 if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
1228 if(TMath::Abs(mctrackmother->Particle()->GetPdgCode())==321)
1229 fQACollection->Fill("Ke3Kecorr",mctrack->Pt(),mctrackmother->Pt());
1230 else if(TMath::Abs(mctrackmother->Particle()->GetPdgCode())==130)
1231 fQACollection->Fill("Ke3K0Lecorr",mctrack->Pt(),mctrackmother->Pt());
1232 }
1233 }
1234 }
1235
1236 // Fill Containers
1237 if(signal) {
1238 // Apply weight for background contamination
1239 if(fBackGroundFactorApply) {
4437a0d2 1240 if(IsPbPb() && fCentralityF >= 0) fWeightBackGround = fkBackGroundFactorArray[fCentralityF >= 0 ? fCentralityF : 0]->Eval(TMath::Abs(track->P()));
1241 else fWeightBackGround = fkBackGroundFactorArray[0]->Eval(TMath::Abs(track->P())); // for pp and pPb
01745870 1242
1243 if(fWeightBackGround < 0.0) fWeightBackGround = 0.0;
1244 else if(fWeightBackGround > 1.0) fWeightBackGround = 1.0;
1245 // weightBackGround as special weight
1246 fVarManager->FillContainer(fContainer, "hadronicBackground", 1, kFALSE, fWeightBackGround);
1247 }
1248 fVarManager->FillCorrelationMatrix(fContainer->GetCorrelationMatrix("correlationstepafterPID"));
1249 }
1250
1251 Bool_t bTagged=kFALSE;
1252 if(GetPlugin(kSecVtx)) {
1253 AliDebug(2, "Running Secondary Vertex Analysis");
1254 if(fSecVtx->Process(track) && signal) {
1255 fVarManager->FillContainer(fContainer, "recTrackContSecvtxReco", AliHFEcuts::kStepHFEcutsSecvtx, kFALSE);
1256 fVarManager->FillContainer(fContainer, "recTrackContSecvtxMC", AliHFEcuts::kStepHFEcutsSecvtx, kTRUE);
1257 bTagged=kTRUE;
1258 }
1259 }
1260
1261 // Electron background analysis
1262 if (GetPlugin(kIsElecBackGround)) {
1263
1264 AliDebug(2, "Running BackGround Analysis");
1265
1266 for(Int_t jtrack = 0; jtrack < fESD->GetNumberOfTracks(); jtrack++){
1267 htrack = fESD->GetTrack(jtrack);
1268 if ( itrack == jtrack ) continue;
1269 fElecBackGround->PairAnalysis(track, htrack);
1270 }
1271 } // end of electron background analysis
1272
1273 if (GetPlugin(kDEstep)) {
1274 Double_t weightElecBgV0[kBgLevels] = {0.,0.,0.,};
1275 Int_t elecSource = 0;
1276 Double_t hfeimpactR=0., hfeimpactnsigmaR=0.;
1277 fExtraCuts->GetHFEImpactParameters(track, hfeimpactR, hfeimpactnsigmaR);
1278 if(HasMCData())
1279 {
1280 if(mctrack && (TMath::Abs(mctrack->Particle()->GetPdgCode()) != 11)){
1281 fQACollection->Fill("hadronsBeforeIPcut",track->Pt());
1282 }
1283 if(fMCQA && signal) {
1284
1285 fMCQA->SetContainerStep(0);
1286 for(Int_t iLevel = 0; iLevel < kBgLevels; iLevel++){
1287 weightElecBgV0[iLevel] = fMCQA->GetWeightFactor(mctrack, iLevel); // positive:conversion e, negative: nonHFE
1288 if(!fisNonHFEsystematics || IsPbPb())break;
1289 }
1290
1291 if(fisNonHFEsystematics){
1292 //Fill additional containers for electron source distinction
1293 elecSource = fMCQA->GetElecSource(mctrack->Particle());
1294 const Char_t *sourceName[kElecBgSpecies]={"Pion","Eta","Omega","Phi","EtaPrime","Rho"};
1295 const Char_t *levelName[kBgLevels]={"Best","Lower","Upper"};
1296 Int_t iName = 0;
1297 for(Int_t iSource = AliHFEmcQA::kPi0; iSource <= AliHFEmcQA::kGammaRho0; iSource++){
1298 if((iSource == AliHFEmcQA::kElse)||(iSource == AliHFEmcQA::kMisID)) continue;
1299 if(elecSource == iSource){
1300 for(Int_t iLevel = 0; iLevel < kBgLevels; iLevel++){
1301 if(weightElecBgV0[iLevel]>0) fVarManager->FillContainer(fContainer, Form("conversionElecs%s%s",sourceName[iName], levelName[iLevel]), 0, kFALSE, weightElecBgV0[iLevel]);
1302 else if(weightElecBgV0[iLevel]<0) fVarManager->FillContainer(fContainer, Form("mesonElecs%s%s",sourceName[iName], levelName[iLevel]), 0, kFALSE, -1*weightElecBgV0[iLevel]);
1303 if(IsPbPb())break;
1304 }
1305 break;
1306 }
1307 iName++;
1308 if(iName == kElecBgSpecies)iName = 0;
1309 }
1310 }
1311 //else{
1312 if(weightElecBgV0[0]>0) {
1313 fVarManager->FillContainer(fContainer, "conversionElecs", 0, kFALSE, weightElecBgV0[0]);
1314 fVarManager->FillContainer(fContainer, "conversionElecs", 5, kTRUE, weightElecBgV0[0]);
1315 }
1316 else if(weightElecBgV0[0]<0) {
1317 fVarManager->FillContainer(fContainer, "mesonElecs", 0, kFALSE, -1*weightElecBgV0[0]);
1318 fVarManager->FillContainer(fContainer, "mesonElecs", 5, kTRUE, -1*weightElecBgV0[0]);
1319 }
1320 //}
1321 if(bTagged){ // bg estimation for the secondary vertex tagged signals
1322 if(weightElecBgV0[0]>0) fVarManager->FillContainer(fContainer, "conversionElecs", 2, kFALSE, weightElecBgV0[0]);
1323 else if(weightElecBgV0[0]<0) fVarManager->FillContainer(fContainer, "mesonElecs", 2, kFALSE, -1*weightElecBgV0[0]);
1324 }
1325 }
1326 } // end of MC
1327
1328 dataDca[0]=-1; //for data, don't know the origin
1329 dataDca[1]=track->Pt();
1330 dataDca[2]=hfeimpactR;
1331 dataDca[3]=fCentralityF;
1332 dataDca[4] = v0pid;
1333 dataDca[5] = double(track->Charge());
1334 if (!HasMCData()) fQACollection->Fill("Dca", dataDca);
1335
1336 // Fill Containers for impact parameter analysis
1337 if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepHFEcutsDca + AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kNcutStepsRecTrack,track)) continue;
1338 if(signal) {
1339 // Apply weight for background contamination after ip cut
1340 if(fBackGroundFactorApply) {
1341 fWeightBackGround = fkBackGroundFactorArray[0]->Eval(TMath::Abs(track->P())); // pp case
1342 if(fWeightBackGround < 0.0) fWeightBackGround = 0.0;
1343 else if(fWeightBackGround > 1.0) fWeightBackGround = 1.0;
1344 // weightBackGround as special weight
1345 fVarManager->FillContainer(fContainer, "hadronicBackground", 2, kFALSE, fWeightBackGround);
1346 }
1347 }
1348
1349 if(HasMCData()){
1350 if(fMCQA && signal) {
1351 fMCQA->SetContainerStep(1);
1352 for(Int_t iLevel = 0; iLevel < kBgLevels; iLevel++){
1353 weightElecBgV0[iLevel] = fMCQA->GetWeightFactor(mctrack, iLevel); // positive:conversion e, negative: nonHFE
1354 if(!fisNonHFEsystematics || IsPbPb())break;
1355 }
1356 if(fisNonHFEsystematics){
1357 //Fill additional containers for electron source distinction
1358 elecSource = fMCQA->GetElecSource(mctrack->Particle());
1359 const Char_t *sourceName[kElecBgSpecies]={"Pion","Eta","Omega","Phi","EtaPrime","Rho"};
1360 const Char_t *levelName[kBgLevels]={"Best","Lower","Upper"};
1361 Int_t iName = 0;
1362 for(Int_t iSource = AliHFEmcQA::kPi0; iSource <= AliHFEmcQA::kGammaRho0; iSource++){
1363 if((iSource == AliHFEmcQA::kElse)||(iSource == AliHFEmcQA::kMisID)) continue;
1364 if(elecSource == iSource){
1365 for(Int_t iLevel = 0; iLevel < kBgLevels; iLevel++){
1366 if(weightElecBgV0[iLevel]>0) fVarManager->FillContainer(fContainer, Form("conversionElecs%s%s",sourceName[iName], levelName[iLevel]), 1, kFALSE, weightElecBgV0[iLevel]);
1367 else if(weightElecBgV0[iLevel]<0) fVarManager->FillContainer(fContainer, Form("mesonElecs%s%s",sourceName[iName], levelName[iLevel]), 1, kFALSE, -1*weightElecBgV0[iLevel]);
1368 if(IsPbPb())break;
1369 }
1370 break;
1371 }
1372 iName++;
1373 if(iName == kElecBgSpecies)iName = 0;
1374 }
1375 }
1376 // else{
1377 if(weightElecBgV0[0]>0) {
1378 fVarManager->FillContainer(fContainer, "conversionElecs", 1, kFALSE, weightElecBgV0[0]);
1379 fVarManager->FillContainer(fContainer, "conversionElecs", 6, kTRUE, weightElecBgV0[0]);
1380 }
1381 else if(weightElecBgV0[0]<0) {
1382 fVarManager->FillContainer(fContainer, "mesonElecs", 1, kFALSE, -1*weightElecBgV0[0]);
1383 fVarManager->FillContainer(fContainer, "mesonElecs", 6, kTRUE, -1*weightElecBgV0[0]);
1384 }
1385 //}
1386 }
1387 }
1388 if(signal) {
1389 fVarManager->FillContainer(fContainer, "recTrackContDEReco", AliHFEcuts::kStepHFEcutsDca, kFALSE);
1390 fVarManager->FillContainer(fContainer, "recTrackContDEMC", AliHFEcuts::kStepHFEcutsDca, kTRUE);
1391 fVarManager->FillCorrelationMatrix(fContainer->GetCorrelationMatrix("correlationstepafterDE"));
1392 }
1393 if(HasMCData()){
1394 if(mctrack && (TMath::Abs(mctrack->Particle()->GetPdgCode()) != 11)){
1395 fQACollection->Fill("hadronsAfterIPcut",track->Pt());
1396 }
1397 }
1398 }
1399
1400 }
1401
1402 // Background subtraction-------------------------------------------------------------------
1403 if (GetPlugin(kNonPhotonicElectron)) fBackgroundSubtraction->CountPoolAssociated(fInputEvent, fCentralityF);
1404 //------------------------------------------------------------------------------------------
1405
1406 fQACollection->Fill("nElectronTracksEvent", nElectronCandidates);
1407}
1408
1409//____________________________________________________________
1410void AliAnalysisTaskHFE::ProcessAOD(){
1411 //
1412 // Run Analysis in AOD Mode
1413 // Function is still in development
1414 //
1415 //printf("Process AOD\n");
1416 AliDebug(3, "Processing AOD Event");
1417 Double_t eventContainer[4];
1418 eventContainer[0] = 0.0;
1419 if(HasMCData()) eventContainer[0] = fVz;
1420 else {
1421 if(fInputEvent->GetPrimaryVertex()) eventContainer[0] = fInputEvent->GetPrimaryVertex()->GetZ();
1422 }
1423 eventContainer[1] = 1.; // No Information available in AOD analysis, assume all events have V0AND
1424 eventContainer[2] = fCentralityF;
1425 eventContainer[3] = fContributors;
1426
1427 //printf("value event container %f, %f, %f, %f\n",eventContainer[0],eventContainer[1],eventContainer[2],eventContainer[3]);
1428
1429 AliAODEvent *fAOD = dynamic_cast<AliAODEvent *>(fInputEvent);
1430 if(!fAOD){
1431 AliError("AOD Event required for AOD Analysis");
1432 return;
1433 }
e2861c1a 1434
1435 // Tag all v0s in current event
1436 if(fV0Tagger){
1437 fV0Tagger->Reset();
1438 fV0Tagger->TagV0Tracks(fAOD);
1439 }
1440 // Set magnetic field if V0 task on
1441 if(fTaggedTrackAnalysis) {
1442 fTaggedTrackAnalysis->SetMagneticField(fAOD->GetMagneticField());
1443 fTaggedTrackAnalysis->SetCentrality(fCentralityF);
1444 if(IsPbPb()) fTaggedTrackAnalysis->SetPbPb();
1445 else {
dab4f95d 1446 if(IspPb()) fTaggedTrackAnalysis->SetpPb();
1447 else fTaggedTrackAnalysis->SetPP();
e2861c1a 1448 }
1449 }
01745870 1450
1451 //printf("Will fill\n");
1452 //
1453 fCFM->GetEventContainer()->Fill(eventContainer, AliHFEcuts::kEventStepRecNoCut);
1454 //printf("Fill\n");
1455 //
1456 if(fIdentifiedAsPileUp) return;
1457 fCFM->GetEventContainer()->Fill(eventContainer, AliHFEcuts::kEventStepRecNoPileUp);
1458
1459 //
1460 if(fIdentifiedAsOutInz) return;
1461 fCFM->GetEventContainer()->Fill(eventContainer, AliHFEcuts::kEventStepZRange);
1462
1463 //
1464 if(!fPassTheEventCut) return;
1465 fCFM->GetEventContainer()->Fill(eventContainer, AliHFEcuts::kEventStepReconstructed);
1466 //printf("pass\n");
1467
1468 fContainer->NewEvent();
1469
1470 fCFM->SetRecEventInfo(fAOD);
1471
1472 if(!fExtraCuts){
1473 fExtraCuts = new AliHFEextraCuts("hfeExtraCuts","HFE Extra Cuts");
1474 }
1475 fExtraCuts->SetRecEventInfo(fAOD);
1476
1477 // Get Number of contributors to the primary vertex for multiplicity-dependent correction
1478 Int_t ncontribVtx = 0;
1479 AliAODVertex *priVtx = fAOD->GetPrimaryVertex();
1480 if(priVtx){
1481 ncontribVtx = priVtx->GetNContributors();
1482 }
1483
1484 // Look for kink mother
1485 Int_t numberofvertices = fAOD->GetNumberOfVertices();
1486 Double_t listofmotherkink[numberofvertices];
1487 Int_t numberofmotherkink = 0;
1488 for(Int_t ivertex=0; ivertex < numberofvertices; ivertex++) {
1489 AliAODVertex *aodvertex = fAOD->GetVertex(ivertex);
1490 if(!aodvertex) continue;
1491 if(aodvertex->GetType()==AliAODVertex::kKink) {
1492 AliAODTrack *mother = (AliAODTrack *) aodvertex->GetParent();
1493 if(!mother) continue;
1494 Int_t idmother = mother->GetID();
1495 listofmotherkink[numberofmotherkink] = idmother;
1496 //printf("ID %d\n",idmother);
1497 numberofmotherkink++;
1498 }
1499 }
1500 //printf("Number of kink mother in the events %d\n",numberofmotherkink);
1501
1502 // Background subtraction-------------------------------------------------------------------
1503 if (GetPlugin(kNonPhotonicElectron)) fBackgroundSubtraction->FillPoolAssociatedTracks(fInputEvent, fCentralityF);
1504 //------------------------------------------------------------------------------------------
1505
1506 // Loop over tracks
1507 AliAODTrack *track = NULL;
1508 AliAODMCParticle *mctrack = NULL;
1509 Double_t dataDca[6]; // [source, pT, dca, centrality]
1510 Int_t nElectronCandidates = 0;
1511 Bool_t signal;
1512
1513 //printf("Number of track %d\n",(Int_t) fAOD->GetNumberOfTracks());
1514 for(Int_t itrack = 0; itrack < fAOD->GetNumberOfTracks(); itrack++){
1515 track = fAOD->GetTrack(itrack); mctrack = NULL;
1516 if(!track) continue;
e2861c1a 1517
1518 // fill counts of v0-identified particles
dab4f95d 1519 AliPID::EParticleType v0pid = fV0Tagger ? fV0Tagger->GetV0Info(track->GetID()) : AliPID::kUnknown;
e2861c1a 1520 // here the tagged track analysis will run
1521 if(fTaggedTrackAnalysis && v0pid != AliPID::kUnknown){
1522 AliDebug(1, Form("Track identified as %s", AliPID::ParticleName(v0pid)));
1523 fTaggedTrackAnalysis->ProcessTrack(track, v0pid);
1524 AliDebug(1, "V0 PID done");
1525 }
01745870 1526
1527 signal = kTRUE;
1528 if(HasMCData()){
01745870 1529 Int_t label = TMath::Abs(track->GetLabel());
1530 if(label && label < fAODArrayMCInfo->GetEntriesFast())
1531 mctrack = dynamic_cast<AliAODMCParticle *>(fAODArrayMCInfo->At(label));
1532 if(fFillSignalOnly && !fCFM->CheckParticleCuts(AliHFEcuts::kStepMCGenerated, mctrack)) signal = kFALSE;
1533 }
1534
1535 fVarManager->NewTrack(track, mctrack, fCentralityF, -1, signal);
1536
1537 if(fFillNoCuts) {
1538 if(signal || !fFillSignalOnly){
1539 fVarManager->FillContainer(fContainer, "recTrackContReco", AliHFEcuts::kStepRecNoCut, kFALSE);
1540 fVarManager->FillContainer(fContainer, "recTrackContMC", AliHFEcuts::kStepRecNoCut, kTRUE);
1541 }
1542 }
1543
e2861c1a 1544 // begin AOD QA
1545 fQACollection->Fill("Filterbegin", -1);
1546 for(Int_t k=0; k<20; k++) {
1547 Int_t u = 1<<k;
1548 if((track->TestFilterBit(u))) {
1549 fQACollection->Fill("Filterbegin", k);
1550 }
1551 }
1552
01745870 1553 if(fApplyCutAOD) {
1554 //printf("Apply cuts\n");
1555 // RecKine: ITSTPC cuts
1556 if(!ProcessCutStep(AliHFEcuts::kStepRecKineITSTPC, track)) continue;
1557
1558 // Reject kink mother
1559 if(fRejectKinkMother) {
4437a0d2 1560 Bool_t kinkmotherpass = kTRUE;
1561 for(Int_t kinkmother = 0; kinkmother < numberofmotherkink; kinkmother++) {
1562 if(track->GetID() == listofmotherkink[kinkmother]) {
1563 kinkmotherpass = kFALSE;
1564 continue;
1565 }
1566 }
1567 if(!kinkmotherpass) continue;
01745870 1568 }
1569
1570 // RecPrim
1571 if(!ProcessCutStep(AliHFEcuts::kStepRecPrim, track)) continue;
1572
1573 // HFEcuts: ITS layers cuts
1574 if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsITS, track)) continue;
1575
1576 // HFE cuts: TOF PID and mismatch flag
1577 if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTOF, track)) continue;
1578
1579 // HFE cuts: TPC PID cleanup
1580 if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTPC, track)) continue;
1581
1582 // HFEcuts: Nb of tracklets TRD0
1583 if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTRD, track)) continue;
1584 }
1585
1586 // Fill correlation maps before PID
1587 if(signal && fContainer->GetCorrelationMatrix("correlationstepbeforePID")) {
1588 //printf("Fill correlation maps before PID\n");
1589 fVarManager->FillCorrelationMatrix(fContainer->GetCorrelationMatrix("correlationstepbeforePID"));
1590 }
1591
1592 if(HasMCData()){
1593 Double_t hfeimpactR4all=0., hfeimpactnsigmaR4all=0.;
1594 Int_t sourceDca =-1;
1595 if(mctrack && (TMath::Abs(mctrack->GetPdgCode()) == 211)){
1596 if(track->Pt()>4.){
1597 fExtraCuts->GetHFEImpactParameters(track, hfeimpactR4all, hfeimpactnsigmaR4all);
1598 dataDca[0]=0; //pion
1599 dataDca[1]=track->Pt();
1600 dataDca[2]=hfeimpactR4all;
1601 dataDca[3]=fCentralityF;
1602 dataDca[4] = -1; // not store V0 for the moment
1603 dataDca[5] = double(track->Charge());
1604 fQACollection->Fill("Dca", dataDca);
1605 }
1606 }
1607 else if(mctrack && (TMath::Abs(mctrack->GetPdgCode()) == 11)){ // to increas statistics for Martin
1608 if(signal){
1609 fExtraCuts->GetHFEImpactParameters(track, hfeimpactR4all, hfeimpactnsigmaR4all);
1610 if(fSignalCuts->IsCharmElectron(track)){
1611 sourceDca=1;
1612 }
1613 else if(fSignalCuts->IsBeautyElectron(track)){
1614 sourceDca=2;
1615 }
1616 else if(fSignalCuts->IsGammaElectron(track)){
1617 sourceDca=3;
1618 }
1619 else if(fSignalCuts->IsNonHFElectron(track)){
1620 sourceDca=4;
1621 }
1622 else if(fSignalCuts->IsJpsiElectron(track)){
1623 sourceDca=5;
1624 }
1625 else {
1626 sourceDca=6;
1627 }
1628 dataDca[0]=sourceDca;
1629 dataDca[1]=track->Pt();
1630 dataDca[2]=hfeimpactR4all;
1631 dataDca[3]=fCentralityF;
1632 dataDca[4] = -1; // not store V0 for the moment
1633 dataDca[5] = double(track->Charge());
1634 if(signal) fQACollection->Fill("Dca", dataDca);
1635 }
1636 }
1637 }
1638
1639 //printf("Will process to PID\n");
1640
1641 // track accepted, do PID
1642 AliHFEpidObject hfetrack;
1643 hfetrack.SetAnalysisType(AliHFEpidObject::kAODanalysis);
1644 hfetrack.SetRecTrack(track);
1645 if(HasMCData()) hfetrack.SetMCTrack(mctrack);
1646 hfetrack.SetCentrality(fCentralityF);
1647 hfetrack.SetMulitplicity(ncontribVtx); // for correction
4437a0d2 1648 if(IsPbPb()) hfetrack.SetPbPb();
1649 else{
1650 if(IspPb()) hfetrack.SetpPb();
1651 else hfetrack.SetPP();
1652 }
01745870 1653 fPID->SetVarManager(fVarManager);
1654 if(!fPID->IsSelected(&hfetrack, fContainer, "recTrackCont", fPIDqa)) continue;
1655 // we will do PID here as soon as possible
1656
1657 // Background subtraction----------------------------------------------------------------------------------------------
1658 if (GetPlugin(kNonPhotonicElectron)) {
1659 Int_t indexmother = -1;
1660 Int_t mcsource = -1;
39427568 1661 if(HasMCData() && mctrack) mcsource = fBackgroundSubtraction->FindMother(TMath::Abs(track->GetLabel()),indexmother);
01745870 1662 fBackgroundSubtraction->LookAtNonHFE(itrack, track, fInputEvent, 1, fCentralityF, -1,mcsource, indexmother);
1663 }
1664 //---------------------------------------------------------------------------------------------------------------------
1665
1666 // end AOD QA
1667 fQACollection->Fill("Filterend", -1);
1668 for(Int_t k=0; k<20; k++) {
1669 Int_t u = 1<<k;
1670 if((track->TestFilterBit(u))) {
1671 fQACollection->Fill("Filterend", k);
1672 }
1673 }
1674
1675 // Apply weight for background contamination
1676 //Double_t weightBackGround = 1.0;
1677 if(signal) {
1678 // Apply weight for background contamination
1679 if(fBackGroundFactorApply) {
4437a0d2 1680 if(IsPbPb() && fCentralityF >= 0) fWeightBackGround = fkBackGroundFactorArray[fCentralityF >= 0 ? fCentralityF : 0]->Eval(TMath::Abs(track->P()));
01745870 1681 else fWeightBackGround = fkBackGroundFactorArray[0]->Eval(TMath::Abs(track->P())); // pp case
1682
1683 if(fWeightBackGround < 0.0) fWeightBackGround = 0.0;
1684 else if(fWeightBackGround > 1.0) fWeightBackGround = 1.0;
1685 // weightBackGround as special weight
1686 fVarManager->FillContainer(fContainer, "hadronicBackground", 1, kFALSE, fWeightBackGround);
1687 }
1688 fVarManager->FillCorrelationMatrix(fContainer->GetCorrelationMatrix("correlationstepafterPID"));
1689 }
1690
1691 nElectronCandidates++;
1692
1693 if (GetPlugin(kDEstep)) {
1694 if (!HasMCData()){
1695 Double_t hfeimpactR=0., hfeimpactnsigmaR=0.;
1696 fExtraCuts->GetHFEImpactParameters(track, hfeimpactR, hfeimpactnsigmaR);
1697 dataDca[0]=-1; //for data, don't know the origin
1698 dataDca[1]=track->Pt();
1699 dataDca[2]=hfeimpactR;
1700 dataDca[3]=fCentralityF;
1701 dataDca[4] = -1; // not store V0 for the moment
1702 dataDca[5] = double(track->Charge());
1703 fQACollection->Fill("Dca", dataDca);
1704 }
1705
1706 // Fill Containers for impact parameter analysis
1707 if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepHFEcutsDca + AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kNcutStepsRecTrack,track)) continue;
1708 if(signal) {
1709 // Apply weight for background contamination after ip cut
1710 if(fBackGroundFactorApply) {
1711 fWeightBackGround = fkBackGroundFactorArray[0]->Eval(TMath::Abs(track->P())); // pp case
1712 if(fWeightBackGround < 0.0) fWeightBackGround = 0.0;
1713 else if(fWeightBackGround > 1.0) fWeightBackGround = 1.0;
1714 // weightBackGround as special weight
1715 fVarManager->FillContainer(fContainer, "hadronicBackground", 2, kFALSE, fWeightBackGround);
1716 }
1717 }
1718 }
1719 }
1720
1721 // Background subtraction-------------------------------------------------------------------
1722 if (GetPlugin(kNonPhotonicElectron)) fBackgroundSubtraction->CountPoolAssociated(fInputEvent, fCentralityF);
1723 //------------------------------------------------------------------------------------------
1724
1725 fQACollection->Fill("nElectronTracksEvent", nElectronCandidates);
1726}
1727
1728//____________________________________________________________
1729Bool_t AliAnalysisTaskHFE::ProcessMCtrack(AliVParticle *track){
1730 //
1731 // Filter the Monte Carlo Track
1732 // Additionally Fill a THnSparse for Signal To Background Studies
1733 // Works for AOD and MC analysis Type
1734 //
1735 fVarManager->NewTrack(track, NULL, fCentralityF, -1, kTRUE);
1736
1737
1738 Double_t vertex[3] = {0.,0.,0.}; // Production vertex cut to mask gammas which are NOT supposed to have hits in the first ITS layer(s)
1739 if(IsESDanalysis()){
1740 AliMCParticle *mctrack = dynamic_cast<AliMCParticle *>(track);
1741 if(mctrack){
1742 vertex[0] = mctrack->Particle()->Vx();
1743 vertex[1] = mctrack->Particle()->Vy();
1744 }
1745 } else {
1746 AliAODMCParticle *aodmctrack = dynamic_cast<AliAODMCParticle *>(track);
1747 if(aodmctrack) aodmctrack->XvYvZv(vertex);
1748 }
1749
1750 //printf("MC Generated\n");
1751 if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepMCGenerated, track)) return kFALSE;
1752 //printf("MC Generated pass\n");
1753 fVarManager->FillContainer(fContainer, "MCTrackCont", AliHFEcuts::kStepMCGenerated, kFALSE);
1754
1755 // Step GeneratedZOutNoPileUp
1756 if((fIdentifiedAsPileUp) || (TMath::Abs(fVz) > fCuts->GetVertexRange()) || (fCentralityF < 0)) return kFALSE;
1757 fVarManager->FillContainer(fContainer, "MCTrackCont", AliHFEcuts::kStepMCGeneratedZOutNoPileUpCentralityFine, kFALSE);
1758 //printf("In ProcessMCtrack %f\n",fCentralityF);
1759
1760 // Step Generated Event Cut
1761 if(!fPassTheEventCut) return kFALSE;
1762 fVarManager->FillContainer(fContainer, "MCTrackCont", AliHFEcuts::kStepMCGeneratedEventCut, kFALSE);
1763
e2861c1a 1764 if(IsESDanalysis()){
1765 if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepMCInAcceptance, track)) return kFALSE;
1766 fVarManager->FillContainer(fContainer, "MCTrackCont", AliHFEcuts::kStepMCInAcceptance, kFALSE);
1767 }
01745870 1768 return kTRUE;
1769}
1770
1771//____________________________________________________________
1772Bool_t AliAnalysisTaskHFE::PreSelectTrack(AliESDtrack *track) const {
1773 //
1774 // Preselect tracks
1775 //
1776
1777
1778 Bool_t survived = kTRUE;
1779
1780 if(fCutspreselect) {
1781 //printf("test preselect\n");
1782 if(!fCutspreselect->IsSelected(track)) survived=kFALSE;
1783 }
1784 //printf("survived %d\n",(Int_t)survived);
1785
1786 if(survived && fPIDpreselect){
1787 // Apply PID
1788 AliHFEpidObject hfetrack;
1789 hfetrack.SetAnalysisType(AliHFEpidObject::kESDanalysis);
1790 hfetrack.SetRecTrack(track);
1791 if(!fPIDpreselect->IsSelected(&hfetrack)) {
1792 //printf("Did not pass AliHFEcuts::kPID\n");
1793 survived = kFALSE;
1794 }
1795 //else printf("Pass AliHFEcuts::kPID\n");
1796 }
1797
1798 return survived;
1799
1800}
1801//____________________________________________________________
1802void AliAnalysisTaskHFE::MakeEventContainer(){
1803 //
1804 // Create the event container for the correction framework and link it
1805 // 1st bin: Vertex z-position
1806 // 2nd bin: V0AND decision (normalization to sigma_inel)
1807 // 3rd bin: Centrality class (for pp defined as number of contributors in vertex.)
1808 // 4th bin: Number of contributors > 0
1809 //
1810
1811 const Int_t kNvar = 4; // number of variables on the grid:
1812 Int_t nBins[kNvar] = {120, 2, 11, 2};
1813 Double_t binMin[kNvar] = {-30. , 0., 0.0, 0.};
1814 Double_t binMax[kNvar] = {30., 2., 11.0, 2.};
1815
1816 AliCFContainer *evCont = new AliCFContainer("eventContainer", "Container for events", AliHFEcuts::kNcutStepsEvent, kNvar, nBins);
1817
1818 Double_t *vertexBins = AliHFEtools::MakeLinearBinning(nBins[0], binMin[0], binMax[0]);
1819 Double_t *v0andBins = AliHFEtools::MakeLinearBinning(nBins[1], binMin[1], binMax[1]);
1820 Double_t *centralityBins = AliHFEtools::MakeLinearBinning(nBins[2], binMin[2], binMax[2]);
1821 Double_t *contributorsBins = AliHFEtools::MakeLinearBinning(nBins[3], binMin[3], binMax[3]);
1822 evCont->SetBinLimits(0, vertexBins);
1823 evCont->SetBinLimits(1, v0andBins);
1824 evCont->SetBinLimits(2, centralityBins);
1825 evCont->SetBinLimits(3, contributorsBins);
1826 delete[] vertexBins; delete[] v0andBins; delete[] centralityBins; delete[] contributorsBins;
1827
1828 fCFM->SetEventContainer(evCont);
1829}
1830
1831//____________________________________________________________
1832void AliAnalysisTaskHFE::MakeParticleContainer(){
1833 //
1834 // Create the particle container for the correction framework manager and
1835 // link it
1836 //
1837 if(!fContainer) fContainer = new AliHFEcontainer("trackContainer");
1838 fVarManager->DefineVariables(fContainer);
1839
1840 // Create Correction Framework containers
1841 fContainer->CreateContainer("MCTrackCont", "Track Container filled with MC information", AliHFEcuts::kNcutStepsMCTrack);
1842 fContainer->CreateContainer("recTrackContReco", "Track Container filled with MC information", AliHFEcuts::kNcutStepsRecTrack + fPID->GetNumberOfPIDdetectors());
1843 fContainer->CreateContainer("recTrackContMC", "Track Container filled with MC information", AliHFEcuts::kNcutStepsRecTrack + fPID->GetNumberOfPIDdetectors());
1844
1845 fContainer->CreateContainer("hadronicBackground", "Container for Hadronic Background", 3);
1846 fContainer->CreateContainer("recTrackContDEReco", "Container for displaced electron analysis with Reco information", 1);
1847 fContainer->CreateContainer("recTrackContDEMC", "Container for displaced electron analysis with MC information", 1);
1848 fContainer->CreateContainer("recTrackContSecvtxReco", "Container for secondary vertexing analysis with Reco information", 1);
1849 fContainer->CreateContainer("recTrackContSecvtxMC", "Container for secondary vertexing analysis with MC information", 1);
1850
1851 if(HasMCData()){
1852 fContainer->CreateContainer("conversionElecs", "Container for weighted conversion electrons",7);
1853 fContainer->CreateContainer("mesonElecs", "Container for weighted electrons from meson decays",7);
1854 fContainer->Sumw2("conversionElecs");
1855 fContainer->Sumw2("mesonElecs");
1856
1857 if(fisNonHFEsystematics){
1858 const Char_t *sourceName[kElecBgSpecies]={"Pion","Eta","Omega","Phi","EtaPrime","Rho"};
1859 const Char_t *levelName[kBgLevels]={"Best","Lower","Upper"};
1860 for(Int_t iSource = 0; iSource < kElecBgSpecies; iSource++){
1861 for(Int_t iLevel = 0; iLevel < kBgLevels; iLevel++){
1862 fContainer->CreateContainer(Form("conversionElecs%s%s",sourceName[iSource],levelName[iLevel]), Form("Container for weighted conversion electrons from %s grandm., %s level",sourceName[iSource],levelName[iLevel]),5);
1863 fContainer->CreateContainer(Form("mesonElecs%s%s",sourceName[iSource],levelName[iLevel]), Form("Container for weighted electrons from %s decays, %s level",sourceName[iSource],levelName[iLevel]),5);
1864 fContainer->Sumw2(Form("conversionElecs%s%s",sourceName[iSource],levelName[iLevel]));
1865 fContainer->Sumw2(Form("mesonElecs%s%s",sourceName[iSource],levelName[iLevel]));
1866 if(IsPbPb())break;
1867 }
1868 }
1869 }
1870 //fContainer->CreateContainer("charmElecs", "Container for weighted charm electrons",2);
1871 }
1872
1873 fContainer->CreateCorrelationMatrix("correlationstepafterPID","THnSparse with correlations");
1874 fContainer->CreateCorrelationMatrix("correlationstepafterDE","THnSparse with correlations");
1875 if(!fVarManager->IsVariableDefined("centrality")) {
1876 //printf("Create the two other correlation maps\n");
1877 fContainer->CreateCorrelationMatrix("correlationstepbeforePID","THnSparse with correlations");
1878 fContainer->CreateCorrelationMatrix("correlationstepafterTOF","THnSparse with correlations");
1879 }
1880
1881 // Define the step names
1882 for(UInt_t istep = 0; istep < AliHFEcuts::kNcutStepsMCTrack; istep++){
1883 fContainer->SetStepTitle("MCTrackCont", AliHFEcuts::MCCutName(istep), istep);
1884 }
1885 for(UInt_t istep = 0; istep < AliHFEcuts::kNcutStepsRecTrack; istep++){
1886 fContainer->SetStepTitle("recTrackContReco", AliHFEcuts::RecoCutName(istep), istep);
1887 fContainer->SetStepTitle("recTrackContMC", AliHFEcuts::RecoCutName(istep), istep);
1888 }
1889 for(UInt_t ipid = 0; ipid < fPID->GetNumberOfPIDdetectors(); ipid++){
1890 fContainer->SetStepTitle("recTrackContReco", fPID->SortedDetectorName(ipid), AliHFEcuts::kNcutStepsRecTrack + ipid);
1891 fContainer->SetStepTitle("recTrackContMC", fPID->SortedDetectorName(ipid), AliHFEcuts::kNcutStepsRecTrack + ipid);
1892 }
1893}
1894//____________________________________________________________
1895void AliAnalysisTaskHFE::InitContaminationQA(){
1896 //
1897 // Add QA for Impact Parameter cut
1898 //
1899
1900 TObjArray *array = fVarManager->GetVariables();
1901 Int_t nvars = array->GetEntriesFast();
1902 for(Int_t v = 0; v < nvars; v++) {
1903 AliHFEvarManager::AliHFEvariable *variable = (AliHFEvarManager::AliHFEvariable *) array->At(v);
1904 if(!variable) continue;
1905 TString name(((AliHFEvarManager::AliHFEvariable *)variable)->GetName());
1906 if(!name.CompareTo("pt")) {
1907 const Int_t nBinPt = variable->GetNumberOfBins();
1908 const Double_t *kPtRange = variable->GetBinning();
1909
1910 fQACollection->CreateTH1Farray("hadronsBeforeIPcut", "Hadrons before IP cut", nBinPt, kPtRange);
1911 fQACollection->CreateTH1Farray("hadronsAfterIPcut", "Hadrons after IP cut", nBinPt, kPtRange);
1912
1913 fQACollection->CreateTH2Farray("Ke3Kecorr", "Ke3 decay e and K correlation; Ke3K p_{t}; Ke3e p_{t}; ", nBinPt, kPtRange, 20,0.,20.);
1914 fQACollection->CreateTH2Farray("Ke3K0Lecorr", "Ke3 decay e and K0L correlation; Ke3K0L p_{t}; Ke3e p_{t}; ", nBinPt, kPtRange, 20,0.,20.);
1915 fQACollection->CreateTH1Farray("Kptspectra", "Charged Kaons: MC p_{t} ", nBinPt, kPtRange);
1916 fQACollection->CreateTH1Farray("K0Lptspectra", "K0L: MC p_{t} ", nBinPt, kPtRange);
1917
1918 const Double_t kDCAbound[2] = {-0.2, 0.2};
1919
1920 const Int_t nDimDca=6;
1921 const Int_t nBinDca[nDimDca] = { 8, nBinPt, 800, 12, 6, 2};
1922 Double_t minimaDca[nDimDca] = { -1., 0., kDCAbound[0], -1., -1, -1.1};
1923 Double_t maximaDca[nDimDca] = { 7., 20., kDCAbound[1], 11., 5, 1.1};
1924
1925 Double_t *sourceBins = AliHFEtools::MakeLinearBinning(nBinDca[0], minimaDca[0], maximaDca[0]);
1926 Double_t *dcaBins = AliHFEtools::MakeLinearBinning(nBinDca[2], minimaDca[2], maximaDca[2]);
1927 Double_t *centralityBins = AliHFEtools::MakeLinearBinning(nBinDca[3], minimaDca[3], maximaDca[3]);
1928 Double_t *v0PIDBins = AliHFEtools::MakeLinearBinning(nBinDca[4], minimaDca[4], maximaDca[4]);
1929 Double_t *chargeBins = AliHFEtools::MakeLinearBinning(nBinDca[5], minimaDca[5], maximaDca[5]);
1930
1931 fQACollection->CreateTHnSparseNoLimits("Dca", "Dca; source (0-all, 1-charm,etc); pT [GeV/c]; dca; centrality bin; v0pid; charge", nDimDca, nBinDca);
1932 ((THnSparse*)(fQACollection->Get("Dca")))->SetBinEdges(0, sourceBins);
1933 ((THnSparse*)(fQACollection->Get("Dca")))->SetBinEdges(1, kPtRange);
1934 ((THnSparse*)(fQACollection->Get("Dca")))->SetBinEdges(2, dcaBins);
1935 ((THnSparse*)(fQACollection->Get("Dca")))->SetBinEdges(3, centralityBins);
1936 ((THnSparse*)(fQACollection->Get("Dca")))->SetBinEdges(4, v0PIDBins);
1937 ((THnSparse*)(fQACollection->Get("Dca")))->SetBinEdges(5, chargeBins);
1938
1939 break;
1940 }
1941 }
1942
1943}
1944
1945//____________________________________________________________
1946void AliAnalysisTaskHFE::InitHistoITScluster(){
1947 //
1948 // Initialize a temporary histogram to monitor the chi2/ITS cluster
1949 if(IsPbPb()) {
1950 const Int_t kNDim = 7;
1951 const Int_t kNBins[kNDim] = {88, 20,90,11, 160, 2, 1000};
1952 const Double_t kMin[kNDim] = {0.1, -1,0, 0.,0., 0, 0.};
1953 const Double_t kMax[kNDim] = {20., 1, 2.*TMath::Pi(), 11.,160, 2, 100.};
1954 fQACollection->CreateTHnSparse("fChi2perITScluster", "chi2/ITS cluster; p_{T} (GeV/c);eta;phi; centrality class;nclus;sharebit; #chi^{2}/ITS cluster", kNDim, kNBins, kMin, kMax);
1955 fQACollection->BinLogAxis("fChi2perITScluster", 0);
1956 }
1957 else
1958 {
1959 const Int_t kNDim = 3;
1960 const Int_t kNBins[kNDim] = {44, 11, 1000};
1961 const Double_t kMin[kNDim] = {0.1, 0., 0.};
1962 const Double_t kMax[kNDim] = {20., 11., 100.};
1963 fQACollection->CreateTHnSparse("fChi2perITScluster", "chi2/ITS cluster; p_{T} (GeV/c); centrality class; #chi^{2}/ITS cluster", kNDim, kNBins, kMin, kMax);
1964 fQACollection->BinLogAxis("fChi2perITScluster", 0);
1965 }
1966}
1967
1968//____________________________________________________________
1969void AliAnalysisTaskHFE::SelectSpecialTrigger(const Char_t *trgclust, Int_t runMin, Int_t runMax){
1970 //
1971 // Select only events triggered by a special trigeer cluster
1972 //
1973 if(!fSpecialTrigger) fSpecialTrigger = new AliOADBContainer("SpecialTrigger");
1974 fSpecialTrigger->AppendObject(new TObjString(trgclust), runMin, runMax);
1975}
1976
1977//____________________________________________________________
1978const Char_t * AliAnalysisTaskHFE::GetSpecialTrigger(Int_t run){
1979 //
1980 // Derive selected trigger string for given run
1981 //
1982 if(!fSpecialTrigger) return NULL;
1983 TObjString *trg = dynamic_cast<TObjString *>(fSpecialTrigger->GetObject(run));
1984 if(!trg) return NULL;
1985 return trg->String().Data();
1986}
1987
1988//____________________________________________________________
1989void AliAnalysisTaskHFE::PrintStatus() const {
1990 //
1991 // Print Analysis status
1992 //
1993 printf("\n\tAnalysis Settings\n\t========================================\n\n");
1994 printf("\tSecondary Vertex finding: %s\n", GetPlugin(kSecVtx) ? "YES" : "NO");
1995 printf("\tPrimary Vertex resolution: %s\n", GetPlugin(kPriVtx) ? "YES" : "NO");
1996 printf("\tDisplaced electron analysis step: %s\n", GetPlugin(kDEstep) ? "YES" : "NO");
1997 printf("\tTagged Track Analysis: %s\n", GetPlugin(kTaggedTrackAnalysis) ? "YES" : "NO");
1998 printf("\n");
1999 printf("\tParticle Identification Detectors:\n");
2000 fPID->PrintStatus();
2001 printf("\n");
2002 printf("\tQA: \n");
2003 printf("\t\tPID: %s\n", IsQAOn(kPIDqa) ? "YES" : "NO");
2004 printf("\t\tCUTS: %s\n", (fCuts != NULL && fCuts->IsQAOn()) ? "YES" : "NO");
2005 printf("\t\tMC: %s\n", IsQAOn(kMCqa) ? "YES" : "NO");
2006 printf("\n");
2007}
2008
2009//____________________________________________________________
2010Bool_t AliAnalysisTaskHFE::FillProductionVertex(const AliVParticle * const track) const{
2011 //
2012 // Find the production vertex of the associated MC track
2013 //
2014 if(!fMCEvent) return kFALSE;
2015 const AliVParticle *mctrack = NULL;
2016 TString objectType = track->IsA()->GetName();
2017 if(objectType.CompareTo("AliESDtrack") == 0 || objectType.CompareTo("AliAODTrack") == 0){
2018 // Reconstructed track
2019 mctrack = fMCEvent->GetTrack(TMath::Abs(track->GetLabel()));
2020 } else {
2021 // MCParticle
2022 mctrack = track;
2023 }
2024
2025 if(!mctrack) return kFALSE;
2026
2027 Double_t xv = 0.0;
2028 Double_t yv = 0.0;
2029
2030 if(TString(mctrack->IsA()->GetName()).CompareTo("AliMCParticle") == 0){
2031 // case MCParticle
2032 const AliMCParticle *mcpart = dynamic_cast<const AliMCParticle *>(mctrack);
2033 if(mcpart){
2034 xv = mcpart->Xv();
2035 yv = mcpart->Yv();
2036 }
2037 } else {
2038 // case AODMCParticle
2039 const AliAODMCParticle *mcpart = dynamic_cast<const AliAODMCParticle *>(mctrack);
2040 if(mcpart){
2041 xv = mcpart->Xv();
2042 yv = mcpart->Yv();
2043 }
2044 }
2045
2046 //printf("xv %f, yv %f\n",xv,yv);
2047 fQACollection->Fill("radius", TMath::Abs(xv),TMath::Abs(yv));
2048
2049 return kTRUE;
2050
2051}
2052//__________________________________________
2053void AliAnalysisTaskHFE::SwitchOnPlugin(Int_t plug){
2054 //
2055 // Switch on Plugin
2056 // Available:
2057 // - Primary vertex studies
2058 // - Secondary vertex Studies
2059 // - Post Processing
2060 //
2061 switch(plug){
2062 case kPriVtx: SETBIT(fPlugins, plug); break;
2063 case kSecVtx: SETBIT(fPlugins, plug); break;
2064 case kIsElecBackGround: SETBIT(fPlugins, plug); break;
2065 case kPostProcess: SETBIT(fPlugins, plug); break;
2066 case kDEstep: SETBIT(fPlugins, plug); break;
2067 case kTaggedTrackAnalysis: SETBIT(fPlugins, plug); break;
2068 case kNonPhotonicElectron: SETBIT(fPlugins, plug); break;
2069 default: AliError("Unknown Plugin");
2070 };
2071}
2072//__________________________________________
2073Bool_t AliAnalysisTaskHFE::ProcessCutStep(Int_t cutStep, AliVParticle *track){
2074 //
2075 // Check single track cuts for a given cut step
2076 // Fill the particle container
2077 //
2078 const Int_t kMCOffset = AliHFEcuts::kNcutStepsMCTrack;
2079 if(!fCFM->CheckParticleCuts(cutStep + kMCOffset, track)) return kFALSE;
2080 if(fVarManager->IsSignalTrack()) {
2081 fVarManager->FillContainer(fContainer, "recTrackContReco", cutStep, kFALSE);
2082 fVarManager->FillContainer(fContainer, "recTrackContMC", cutStep, kTRUE);
2083 }
2084 return kTRUE;
2085}
2086//___________________________________________________
2087Bool_t AliAnalysisTaskHFE::ReadCentrality() {
2088 //
2089 // Recover the centrality of the event from ESD or AOD
2090 //
2091
2092 Float_t fCentralityLimitstemp[12];
2093 Float_t fCentralityLimitsdefault[12]= {0.,5.,10., 20., 30., 40., 50., 60.,70.,80., 90., 100.};
2094 if(!fPbPbUserCentralityBinning) memcpy(fCentralityLimitstemp,fCentralityLimitsdefault,sizeof(fCentralityLimitsdefault));
2095 else memcpy(fCentralityLimitstemp,fCentralityLimits,sizeof(fCentralityLimitsdefault));
2096
2097
2098 Int_t bin = -1;
4437a0d2 2099 if(IsPbPb()||IspPb()) {
01745870 2100 // Centrality
2101 AliCentrality *centrality = fInputEvent->GetCentrality();
2102 fCentralityPercent = centrality->GetCentralityPercentile(fCentralityEstimator.Data());
2103 //printf("centrality %f\n",fCentralityPercent);
2104
2105 for(Int_t ibin = 0; ibin < 11; ibin++){
2106 if(fCentralityPercent >= fCentralityLimitstemp[ibin] && fCentralityPercent < fCentralityLimitstemp[ibin+1]){
2107 bin = ibin;
2108 //printf("test bin %f, low %f, high %f, %d\n",fCentralityPercent,fCentralityLimitstemp[ibin],fCentralityLimitstemp[ibin+1],ibin);
2109 break;
2110 }
2111 }
2112
2113 if(bin == -1) bin = 11; // Overflow
2114 } else {
2115 // PP: Tracklet multiplicity, use common definition
2116 Int_t itsMultiplicity = GetITSMultiplicity(fInputEvent);
2117 Int_t multiplicityLimits[8] = {0, 1, 9, 17, 25, 36, 60, 500};
2118 for(Int_t ibin = 0; ibin < 7; ibin++){
2119 if(itsMultiplicity >= multiplicityLimits[ibin] && itsMultiplicity < multiplicityLimits[ibin + 1]){
2120 bin = ibin;
2121 break;
2122 }
2123 }
2124 if(bin == -1) bin = 7; // Overflow
2125 }
2126 fCentralityF = bin;
2127 AliDebug(2, Form("Centrality class %d\n", fCentralityF));
2128
2129
2130 // contributors, to be outsourced
2131 const AliVVertex *vtx;
2132 if(IsAODanalysis()){
2133 AliAODEvent *fAOD = dynamic_cast<AliAODEvent *>(fInputEvent);
2134 if(!fAOD){
2135 AliError("AOD Event required for AOD Analysis");
2136 return kFALSE;
2137 }
2138 vtx = fAOD->GetPrimaryVertex();
2139 } else {
2140 AliESDEvent *fESD = dynamic_cast<AliESDEvent *>(fInputEvent);
2141 if(!fESD){
2142 AliError("ESD Event required for ESD Analysis");
2143 return kFALSE;
2144 }
2145 vtx = fESD->GetPrimaryVertex() ;
2146 }
2147 if(!vtx){
2148 fContributors = 0.5;
2149 return kFALSE;
2150 }
2151 else {
2152 Int_t contributorstemp = vtx->GetNContributors();
2153 if( contributorstemp <= 0) {
2154 fContributors = 0.5;
2155 //printf("Number of contributors %d and vz %f\n",contributorstemp,vtx->GetZ());
2156 }
2157 else fContributors = 1.5;
2158 //printf("Number of contributors %d\n",contributorstemp);
2159 }
2160 return kTRUE;
2161}
2162
2163//___________________________________________________
2164Int_t AliAnalysisTaskHFE::GetITSMultiplicity(AliVEvent *ev){
2165 //
2166 // Definition of the Multiplicity according to the JPSI group (F. Kramer)
2167 //
2168 Int_t nTracklets = 0;
2169 Int_t nAcc = 0;
2170 Double_t etaRange = 1.6;
2171
2172 if (ev->IsA() == AliAODEvent::Class()) {
2173 AliAODTracklets *tracklets = ((AliAODEvent*)ev)->GetTracklets();
2174 nTracklets = tracklets->GetNumberOfTracklets();
2175 for (Int_t nn = 0; nn < nTracklets; nn++) {
2176 Double_t theta = tracklets->GetTheta(nn);
2177 Double_t eta = -TMath::Log(TMath::Tan(theta/2.0));
2178 if (TMath::Abs(eta) < etaRange) nAcc++;
2179 }
2180 } else if (ev->IsA() == AliESDEvent::Class()) {
2181 nTracklets = ((AliESDEvent*)ev)->GetMultiplicity()->GetNumberOfTracklets();
2182 for (Int_t nn = 0; nn < nTracklets; nn++) {
2183 Double_t eta = ((AliESDEvent*)ev)->GetMultiplicity()->GetEta(nn);
2184 if (TMath::Abs(eta) < etaRange) nAcc++;
2185 }
2186 } else return -1;
2187
2188 return nAcc;
2189}
2190
2191//___________________________________________________
2192void AliAnalysisTaskHFE::RejectionPileUpVertexRangeEventCut() {
2193 //
2194 // Recover the centrality of the event from ESD or AOD
2195 //
2196 if(IsAODanalysis()){
2197
2198 AliAODEvent *fAOD = dynamic_cast<AliAODEvent *>(fInputEvent);
2199 if(!fAOD){
2200 AliError("AOD Event required for AOD Analysis");
2201 return;
2202 }
2203 // PileUp
2204 fIdentifiedAsPileUp = kFALSE;
2205 if(fRemovePileUp && fAOD->IsPileupFromSPD()) fIdentifiedAsPileUp = kTRUE;
2206 // Z vertex
2207 fIdentifiedAsOutInz = kFALSE;
2208 //printf("Z vertex %f and out %f\n",fAOD->GetPrimaryVertex()->GetZ(),fCuts->GetVertexRange());
2209 if(TMath::Abs(fAOD->GetPrimaryVertex()->GetZ()) > fCuts->GetVertexRange()) fIdentifiedAsOutInz = kTRUE;
2210 // Event Cut
2211 fPassTheEventCut = kTRUE;
2212 if(!fCFM->CheckEventCuts(AliHFEcuts::kEventStepReconstructed, fAOD)) fPassTheEventCut = kFALSE;
2213
2214
2215 } else {
2216
2217 AliDebug(3, "Processing ESD Centrality");
2218 AliESDEvent *fESD = dynamic_cast<AliESDEvent *>(fInputEvent);
2219 if(!fESD){
2220 AliError("ESD Event required for ESD Analysis");
2221 return;
2222 }
2223 // PileUp
2224 fIdentifiedAsPileUp = kFALSE;
2225 if(fRemovePileUp && fESD->IsPileupFromSPD()) fIdentifiedAsPileUp = kTRUE;
2226
2227
2228
2229 // Z vertex
2230 fIdentifiedAsOutInz = kFALSE;
2231 Bool_t findvertex = kTRUE;
2232 const AliESDVertex* vtxESD = fESD->GetPrimaryVertex();
2233 if((!vtxESD) || (vtxESD->GetNContributors() <= 0)) findvertex = kFALSE;
2234 if(findvertex) {
2235 if(TMath::Abs(vtxESD->GetZ()) > fCuts->GetVertexRange()) fIdentifiedAsOutInz = kTRUE;
2236 }
2237
2238 //Event Cut
2239 fPassTheEventCut = kTRUE;
2240 if(!fCFM->CheckEventCuts(AliHFEcuts::kEventStepReconstructed, fESD)) fPassTheEventCut = kFALSE;
2241
2242 }
2243
2244}
8c07fb0d 2245
2246//___________________________________________________
4437a0d2 2247Bool_t AliAnalysisTaskHFE::CheckTRDTriggerESD(AliESDEvent *ev) {
8c07fb0d 2248//
2249// Check TRD trigger; pPb settings
2250//
2251 Bool_t cint8=kFALSE;
2252 Bool_t cint7=kFALSE;
2253 Bool_t cint5=kFALSE;
2254 Bool_t trdtrgevent=kFALSE;
2255
4437a0d2 2256
2257 // mb selection of WU events
8c07fb0d 2258 if(fWhichTRDTrigger==1)
2259 {
4437a0d2 2260 if(ev->IsTriggerClassFired("CINT7WU-B-NOPF-ALL"))
8c07fb0d 2261 {
2262 DrawTRDTrigger(ev);
2263 return kTRUE;
2264 }
4437a0d2 2265 else return kFALSE;
8c07fb0d 2266 }
618038b6 2267
4437a0d2 2268
2269 // HSE no cleanup
8c07fb0d 2270 if(fWhichTRDTrigger==2)
2271 {
4437a0d2 2272 cint8= ev->IsTriggerClassFired("CINT8WUHSE-B-NOPF-CENT");
2273 cint7= ev->IsTriggerClassFired("CINT7WUHSE-B-NOPF-CENT");
2274 cint5= (ev->IsTriggerClassFired("CINT5WU-B-NOPF-ALL")) &&
2275 (ev->GetHeader()->GetL1TriggerInputs() & (1 << 10));
2276 if((cint7==kFALSE)&&(cint8==kFALSE)&&(cint5==kFALSE)) return kFALSE;
8c07fb0d 2277 else
2278 {
2279 DrawTRDTrigger(ev);
2280 return kTRUE;
2281 }
2282 }
2283
4437a0d2 2284
2285
618038b6 2286 //HQU no cleanup
8c07fb0d 2287 if(fWhichTRDTrigger==3)
2288 {
8c07fb0d 2289 cint8= ev->IsTriggerClassFired("CINT8WUHQU-B-NOPF-CENT");
2290 cint7= ev->IsTriggerClassFired("CINT7WUHQU-B-NOPF-CENT");
2291 cint5= (ev->IsTriggerClassFired("CINT5WU-B-NOPF-ALL")) &&
2292 (ev->GetHeader()->GetL1TriggerInputs() & (1 << 12));
8c07fb0d 2293 if((cint7==kFALSE)&&(cint8==kFALSE)&&(cint5==kFALSE)) return kFALSE;
2294 else
2295 {
2296 DrawTRDTrigger(ev);
2297 return kTRUE;
2298 }
8c07fb0d 2299 }
618038b6 2300
4437a0d2 2301
2302
2303 return trdtrgevent;
2304
2305}
2306
2307
2308//___________________________________________________
2309Bool_t AliAnalysisTaskHFE::CheckTRDTrigger(AliVEvent *ev) {
2310//
2311// Check TRD trigger; pPb settings
2312//
2313
2314 fTRDTriggerAnalysis->CalcTriggers(ev);
2315
2316 // HSE cleanup
2317 if(fWhichTRDTrigger==6)
2318 {
2319 if(!fTRDTriggerAnalysis->IsFired(AliTRDTriggerAnalysis::kHSE))
2320 {
2321 return kFALSE;
2322 }
2323 else
2324 {
2325 // DrawTRDTrigger(ev);
2326 return kTRUE;
2327 }
2328 }
2329
2330
2331
618038b6 2332 // HQU cleanup
4437a0d2 2333 if(fWhichTRDTrigger==7)
8c07fb0d 2334 {
618038b6 2335
2336 if(!fTRDTriggerAnalysis->IsFired(AliTRDTriggerAnalysis::kHQU))
2337 {
2338 return kFALSE;
2339 }
8c07fb0d 2340 else
2341 {
4437a0d2 2342 // DrawTRDTrigger(ev);
8c07fb0d 2343 return kTRUE;
2344 }
2345 }
4437a0d2 2346
2347 // HSE or HQU cleanup
2348 if(fWhichTRDTrigger==8)
2349 {
2350 if((fTRDTriggerAnalysis->IsFired(AliTRDTriggerAnalysis::kHSE))||(fTRDTriggerAnalysis->IsFired(AliTRDTriggerAnalysis::kHQU)))
2351 {
2352 // DrawTRDTrigger(ev);
2353 return kTRUE;
2354 }
2355 else
2356 {
2357 return kFALSE;
2358 }
2359 }
2360
2361 return kFALSE;
8c07fb0d 2362
2363}
2364
2365void AliAnalysisTaskHFE::DrawTRDTrigger(AliESDEvent *ev) {
2366
2367 Int_t ntriggerbit=0;
2368 fQACollection->Fill("nTriggerBit",ntriggerbit);
2369 if(ev->IsTriggerClassFired("CINT7-B-NOPF-ALLNOTRD"))
2370 {
2371 ntriggerbit=2;
2372 fQACollection->Fill("nTriggerBit",ntriggerbit);
2373 }
2374 if(ev->IsTriggerClassFired("CINT7WU-B-NOPF-ALL"))
2375 {
2376 ntriggerbit=3;
2377 fQACollection->Fill("nTriggerBit",ntriggerbit);
2378 if(ev->IsTriggerClassFired("CINT7WUHSE-B-NOPF-CENT")) {
2379 ntriggerbit=18;
2380 fQACollection->Fill("nTriggerBit",ntriggerbit);
2381 }
2382 if(ev->IsTriggerClassFired("CINT7WUHQU-B-NOPF-CENT")) {
2383 ntriggerbit=19;
2384 fQACollection->Fill("nTriggerBit",ntriggerbit);
2385 }
2386 }
2387 if(ev->IsTriggerClassFired("CINT7WUHJT-B-NOPF-CENT"))
2388 {
2389 ntriggerbit=4;
2390 fQACollection->Fill("nTriggerBit",ntriggerbit);
2391
2392 if(ev->IsTriggerClassFired("CINT7WUHSE-B-NOPF-CENT")) {
2393 ntriggerbit=13;
2394 fQACollection->Fill("nTriggerBit",ntriggerbit);
2395 }
2396 if(ev->IsTriggerClassFired("CINT7WUHQU-B-NOPF-CENT")) {
2397 ntriggerbit=14;
2398 fQACollection->Fill("nTriggerBit",ntriggerbit);
2399 if(ev->IsTriggerClassFired("CINT7WUHSE-B-NOPF-CENT")) {
2400 ntriggerbit=17;
2401 fQACollection->Fill("nTriggerBit",ntriggerbit);
2402 }
2403 }
2404 }
2405 if(ev->IsTriggerClassFired("CINT7WUHQU-B-NOPF-CENT"))
2406 {
2407 ntriggerbit=5;
2408 fQACollection->Fill("nTriggerBit",ntriggerbit);
2409 if(ev->IsTriggerClassFired("CINT7WUHSE-B-NOPF-CENT")) {
2410 ntriggerbit=11;
2411 fQACollection->Fill("nTriggerBit",ntriggerbit);
2412 }
2413 if((!(ev->IsTriggerClassFired("CINT7WUHSE-B-NOPF-CENT")))&&(!(ev->IsTriggerClassFired("CINT7WUHJT-B-NOPF-CENT")))) {
2414 ntriggerbit=21;
2415 fQACollection->Fill("nTriggerBit",ntriggerbit);
2416
2417 /*
2418 Int_t nTrdTracks = ev->GetNumberOfTrdTracks();
2419 for (Int_t iTrack = 0; iTrack < nTrdTracks; ++iTrack) {
2420 AliESDTrdTrack* trdTrack = ev->GetTrdTrack(iTrack);
2421 printf("GTU track %3i: pt = %5.1f, PID = %3i\n", iTrack, trdTrack->Pt(), trdTrack->GetPID());
2422 }*/
2423
2424
2425 }
2426
2427 }
2428 if(ev->IsTriggerClassFired("CINT7WUHSE-B-NOPF-CENT"))
2429 {
2430 ntriggerbit=6;
2431 fQACollection->Fill("nTriggerBit",ntriggerbit);
2432 if(ev->IsTriggerClassFired("CINT7WUHQU-B-NOPF-CENT")) {
2433 ntriggerbit=12;
2434 fQACollection->Fill("nTriggerBit",ntriggerbit);
2435 }
2436 if(ev->IsTriggerClassFired("CINT7WUHSE-B-NOPF-FAST")){
2437 ntriggerbit=15;
2438 fQACollection->Fill("nTriggerBit",ntriggerbit);
2439
2440 if((!(ev->IsTriggerClassFired("CINT7WUHQU-B-NOPF-CENT")))&&(!(ev->IsTriggerClassFired("CINT7WUHJT-B-NOPF-CENT")))) {
2441 ntriggerbit=20;
2442 fQACollection->Fill("nTriggerBit",ntriggerbit);
2443 /*
2444 Int_t nTrdTracks = ev->GetNumberOfTrdTracks();
2445 for (Int_t iTrack = 0; iTrack < nTrdTracks; ++iTrack) {
2446 AliESDTrdTrack* trdTrack = ev->GetTrdTrack(iTrack);
2447 printf("HSE GTU track %3i: pt = %5.1f, PID = %3i\n", iTrack, trdTrack->Pt(), trdTrack->GetPID());
2448 } */
2449
2450 }
2451
2452 }
2453
2454 }
2455 if(ev->IsTriggerClassFired("CEMC7WUHEE-B-NOPF-CENT")) {
2456 ntriggerbit=7;
2457 fQACollection->Fill("nTriggerBit",ntriggerbit);
2458 }
2459 if(ev->IsTriggerClassFired("CINT7WUHJT-B-NOPF-FAST")){
2460 ntriggerbit=8;
2461 fQACollection->Fill("nTriggerBit",ntriggerbit);
2462 }
2463 if(ev->IsTriggerClassFired("CINT7WUHQU-B-NOPF-FAST")){
2464 ntriggerbit=9;
2465 fQACollection->Fill("nTriggerBit",ntriggerbit);
2466 }
2467 if(ev->IsTriggerClassFired("CINT7WUHSE-B-NOPF-FAST")){
2468 ntriggerbit=10;
2469 fQACollection->Fill("nTriggerBit",ntriggerbit);
2470 if(ev->IsTriggerClassFired("CINT7WUHSE-B-NOPF-CENT")) {
2471 ntriggerbit=16;
2472 fQACollection->Fill("nTriggerBit",ntriggerbit);
2473 }
2474 }
2475 if(ntriggerbit==0) fQACollection->Fill("nTriggerBit",1);
2476
01745870 2477}