1 /**************************************************************************
2 * Author : Nicole Alice Martin (nicole.alice.martin@cern.ch) *
4 * Contributors are mentioned in the code where appropriate. *
6 * Permission to use, copy, modify and distribute this software and its *
7 * documentation strictly for non-commercial purposes is hereby granted *
8 * without fee, provided that the above copyright notice appears in all *
9 * copies and that both the copyright notice and this permission notice *
10 * appear in the supporting documentation. The authors make no claims *
11 * about the suitability of this software for any purpose. It is *
12 * provided "as is" without express or implied warranty. *
13 **************************************************************************/
15 //-----------------------------------------------------------------
16 // AliAnalysisTaskLambdaNAOD class
17 // task for the investigation of (anti-)lambda-n bound state
18 // uses the V0 finders, based on AODs or ESDS
19 //-----------------------------------------------------------------
33 #include <Riostream.h>
41 #include "THnSparse.h"
48 #include <TDatabasePDG.h>
50 #include "AliAnalysisManager.h"
51 #include "AliAnalysisTask.h"
52 #include "AliInputEventHandler.h"
54 #include "AliCentrality.h"
55 #include "AliV0vertexer.h"
56 #include "AliPIDResponse.h"
57 #include "AliMultiplicity.h"
58 #include "AliVertexerTracks.h"
60 #include "AliVEvent.h"
61 #include "AliVTrack.h"
63 #include "AliESDInputHandler.h"
64 #include "AliESDEvent.h"
65 #include "AliESDtrackCuts.h"
66 #include "AliESDpid.h"
69 #include "AliAODInputHandler.h"
70 #include "AliAODEvent.h"
71 #include "AliAODTrack.h"
74 #include "AliCFContainer.h"
75 #include "AliKFParticle.h"
76 #include "AliKFVertex.h"
78 #include "AliMCEventHandler.h"
79 #include "AliMCEvent.h"
82 #include "AliAnalysisTaskLambdaNAOD.h"
86 ClassImp(AliAnalysisTaskLambdaNAOD)
89 //________________________________________________________________________
90 AliAnalysisTaskLambdaNAOD::AliAnalysisTaskLambdaNAOD()
91 : AliAnalysisTaskSE(),
98 fHistCentralityClass10(0),
99 fHistCentralityPercentile(0),
101 fHistTriggerStatAfterEventSelection(0),
102 fHistLambdaNeutronPtGen(0),
103 fHistLambdaNeutronPtGenMinBias(0),
104 fHistLambdaNeutronPtGenCentral(0),
105 fHistLambdaNeutronPtGenSemiCentral(0),
106 fHistAntiLambdaNeutronPtGen(0),
107 fHistAntiLambdaNeutronPtGenMinBias(0),
108 fHistAntiLambdaNeutronPtGenCentral(0),
109 fHistAntiLambdaNeutronPtGenSemiCentral(0),
110 fHistLambdaNeutronInvaMassGen(0),
111 fHistAntiLambdaNeutronInvaMassGen(0),
112 fHistLambdaNeutronDecayLengthGen(0),
113 fHistAntiLambdaNeutronDecayLengthGen(0),
114 fHistLambdaNeutronPtAso(0),
115 fHistLambdaNeutronPtAsoCuts(0),
116 fHistAntiLambdaNeutronPtAso(0),
117 fHistAntiLambdaNeutronPtAsoCuts(0),
118 fHistLambdaNeutronInvaMassAso(0),
119 fHistAntiLambdaNeutronInvaMassAso(0),
120 fHistLambdaNeutronDecayLengthAso(0),
121 fHistAntiLambdaNeutronDecayLengthAso(0),
123 fHistArmenterosPodolanskiDeuteronPion(0),
124 fHistArmenterosPodolanskiAntiDeuteronPion(0),
133 fOutputContainer(NULL)
135 // default Constructor
137 // Define input and output slots here
140 //________________________________________________________________________
141 AliAnalysisTaskLambdaNAOD::AliAnalysisTaskLambdaNAOD(const char *name)
142 : AliAnalysisTaskSE(name),
143 fAnalysisType("ESD"),
149 fHistCentralityClass10(0),
150 fHistCentralityPercentile(0),
152 fHistTriggerStatAfterEventSelection(0),
153 fHistLambdaNeutronPtGen(0),
154 fHistLambdaNeutronPtGenMinBias(0),
155 fHistLambdaNeutronPtGenCentral(0),
156 fHistLambdaNeutronPtGenSemiCentral(0),
157 fHistAntiLambdaNeutronPtGen(0),
158 fHistAntiLambdaNeutronPtGenMinBias(0),
159 fHistAntiLambdaNeutronPtGenCentral(0),
160 fHistAntiLambdaNeutronPtGenSemiCentral(0),
161 fHistLambdaNeutronInvaMassGen(0),
162 fHistAntiLambdaNeutronInvaMassGen(0),
163 fHistLambdaNeutronDecayLengthGen(0),
164 fHistAntiLambdaNeutronDecayLengthGen(0),
165 fHistLambdaNeutronPtAso(0),
166 fHistLambdaNeutronPtAsoCuts(0),
167 fHistAntiLambdaNeutronPtAso(0),
168 fHistAntiLambdaNeutronPtAsoCuts(0),
169 fHistLambdaNeutronInvaMassAso(0),
170 fHistAntiLambdaNeutronInvaMassAso(0),
171 fHistLambdaNeutronDecayLengthAso(0),
172 fHistAntiLambdaNeutronDecayLengthAso(0),
174 fHistArmenterosPodolanskiDeuteronPion(0),
175 fHistArmenterosPodolanskiAntiDeuteronPion(0),
184 fOutputContainer(NULL)
188 // Define input and output slots here
189 // Input slot #0 works with a TChain
190 DefineInput(0, TChain::Class());
191 // Output slot #0 writes into a TH1 container
192 DefineOutput(1, TObjArray::Class());
193 DefineOutput(2, TTree::Class());
195 //ESD Track Cuts for v0's
196 fESDtrackCutsV0 = new AliESDtrackCuts("AliESDtrackCutsV0","AliESDtrackCutsV0");
197 fESDtrackCutsV0->SetAcceptKinkDaughters(kFALSE);
198 fESDtrackCutsV0->SetMinNClustersTPC(60);
199 fESDtrackCutsV0->SetMaxChi2PerClusterTPC(6); //set to 6 for sytematics
200 //fESDtrackCutsV0->SetMaxChi2PerClusterTPC(5);
201 fESDtrackCutsV0->SetRequireTPCRefit(kTRUE);
202 // fESDtrackCutsV0->SetMinNClustersITS(1); // to be tested if this cut is not too strong
203 fESDtrackCutsV0->SetEtaRange(-0.9,0.9);
211 //____________________________________________________________
212 AliAnalysisTaskLambdaNAOD::~AliAnalysisTaskLambdaNAOD(){
216 //if (fESD) delete fESD;
217 if (fESDtrackCutsV0) delete fESDtrackCutsV0;
221 //____________________________________________________________
222 const Int_t AliAnalysisTaskLambdaNAOD::fgkPdgCode[] = {
227 1000010020, //Deuteron
228 -1000010020, //Anti-Deuteron
229 1000020030, //Helium3
230 -1000020030, //Anti-Helium3
233 1010000020, //LambdaNeutron
234 -1010000020, //Anti-Lambda-Neutron
235 1010010030, //Hypertriton
236 -1010010030 //Anti-Hypertriton
239 //____________________________________________________________
240 const Double_t AliAnalysisTaskLambdaNAOD::fgkMass[] = {
245 2.80941 //Helium3 Quelle: Wolfram Alpha
248 //________________________________________________________________________
249 void AliAnalysisTaskLambdaNAOD::UserCreateOutputObjects(){
252 AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
253 AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
254 fPIDResponse = inputHandler->GetPIDResponse();
256 // Create histograms for MC
257 //generated histogramms
258 fHistLambdaNeutronPtGen = new TH1F("fHistLambdaNeutronPtGen", "Generated Lambdan", 201, 0., 10.1);
259 fHistLambdaNeutronPtGen->GetYaxis()->SetTitle("Counts");
260 fHistLambdaNeutronPtGen->GetXaxis()->SetTitle("#it{p}_{T} (GeV/#it{c})");
262 fHistLambdaNeutronPtGenMinBias = new TH1F("fHistLambdaNeutronPtGenMinBias", "Generated Lambdan", 201, 0., 10.1);
263 fHistLambdaNeutronPtGenMinBias->GetYaxis()->SetTitle("Counts");
264 fHistLambdaNeutronPtGenMinBias->GetXaxis()->SetTitle("#it{p}_{T} (GeV/#it{c})");
266 fHistLambdaNeutronPtGenCentral = new TH1F("fHistLambdaNeutronPtGenCentral", "Generated Lambdan", 201, 0., 10.1);
267 fHistLambdaNeutronPtGenCentral->GetYaxis()->SetTitle("Counts");
268 fHistLambdaNeutronPtGenCentral->GetXaxis()->SetTitle("#it{p}_{T} (GeV/#it{c})");
270 fHistLambdaNeutronPtGenSemiCentral = new TH1F("fHistLambdaNeutronPtGenSemiCentral", "Generated Lambdan", 201, 0., 10.1);
271 fHistLambdaNeutronPtGenSemiCentral->GetYaxis()->SetTitle("Counts");
272 fHistLambdaNeutronPtGenSemiCentral->GetXaxis()->SetTitle("#it{p}_{T} (GeV/#it{c})");
274 fHistAntiLambdaNeutronPtGen = new TH1F("fHistAntiLambdaNeutronPtGen", "Generated #bar{#Lambdan}", 201, 0., 10.1);
275 fHistAntiLambdaNeutronPtGen->GetYaxis()->SetTitle("Counts");
276 fHistAntiLambdaNeutronPtGen->GetXaxis()->SetTitle("#it{p}_{T} (GeV/#it{c})");
278 fHistAntiLambdaNeutronPtGenMinBias = new TH1F("fHistAntiLambdaNeutronPtGenMinBias", "Generated #bar{#Lambdan}", 201, 0., 10.1);
279 fHistAntiLambdaNeutronPtGenMinBias->GetYaxis()->SetTitle("Counts");
280 fHistAntiLambdaNeutronPtGenMinBias->GetXaxis()->SetTitle("#it{p}_{T} (GeV/#it{c})");
282 fHistAntiLambdaNeutronPtGenCentral = new TH1F("fHistAntiLambdaNeutronPtGenCentral", "Generated #bar{#Lambdan}", 201, 0., 10.1);
283 fHistAntiLambdaNeutronPtGenCentral->GetYaxis()->SetTitle("Counts");
284 fHistAntiLambdaNeutronPtGenCentral->GetXaxis()->SetTitle("#it{p}_{T} (GeV/#it{c})");
286 fHistAntiLambdaNeutronPtGenSemiCentral = new TH1F("fHistAntiLambdaNeutronPtGenSemiCentral", "Generated #bar{#Lambdan}", 201, 0., 10.1);
287 fHistAntiLambdaNeutronPtGenSemiCentral->GetYaxis()->SetTitle("Counts");
288 fHistAntiLambdaNeutronPtGenSemiCentral->GetXaxis()->SetTitle("#it{p}_{T} (GeV/#it{c})");
290 fHistLambdaNeutronInvaMassGen = new TH1F("fHistLambdaNeutronInvaMassGen", "Generated #Lambdan ", 100, 2.0, 2.1);
291 fHistLambdaNeutronInvaMassGen->GetYaxis()->SetTitle("Counts");
292 fHistLambdaNeutronInvaMassGen->GetXaxis()->SetTitle("Invariant mass (d #pi^{-}) (GeV/#it{c}^{2})");
294 fHistAntiLambdaNeutronInvaMassGen = new TH1F("fHistAntiLambdaNeutronInvaMassGen", "Generated #bar{#Lambdan}", 100, 2.0, 2.1);
295 fHistAntiLambdaNeutronInvaMassGen->GetYaxis()->SetTitle("Counts");
296 fHistAntiLambdaNeutronInvaMassGen->GetXaxis()->SetTitle("Invariant mass (#bar{d} #pi^{+}) (GeV/#it{c}^{2})");
298 fHistLambdaNeutronDecayLengthGen = new TH1F("fHistLambdaNeutronDecayLengthGen", "Generated #Lambdan", 401, 0., 400.1);
299 fHistLambdaNeutronDecayLengthGen->GetYaxis()->SetTitle("Counts");
300 fHistLambdaNeutronDecayLengthGen->GetXaxis()->SetTitle("#it{decay length} (cm)");
302 fHistAntiLambdaNeutronDecayLengthGen = new TH1F("fHistAntiLambdaNeutronDecayLengthGen", "Generated #bar{#Lambdan}", 401, 0., 400.1);
303 fHistAntiLambdaNeutronDecayLengthGen->GetYaxis()->SetTitle("Counts");
304 fHistAntiLambdaNeutronDecayLengthGen->GetXaxis()->SetTitle("#it{decay length} (cm)");
306 //assoziated (reconstracted) histogramms
307 fHistLambdaNeutronPtAso = new TH1F("fHistLambdaNeutronPtAso", "Associated #Lambdan", 201, 0., 10.1);
308 fHistLambdaNeutronPtAso->GetYaxis()->SetTitle("Counts");
309 fHistLambdaNeutronPtAso->GetXaxis()->SetTitle("#it{p}_{T} (GeV/#it{c})");
311 fHistLambdaNeutronPtAsoCuts = new TH1F("fHistLambdaNeutronPtAsoCuts", "Associated #Lambdan Cuts", 201, 0., 10.1);
312 fHistLambdaNeutronPtAsoCuts->GetYaxis()->SetTitle("Counts");
313 fHistLambdaNeutronPtAsoCuts->GetXaxis()->SetTitle("#it{p}_{T} (GeV/#it{c})");
315 fHistAntiLambdaNeutronPtAsoCuts = new TH1F("fHistAntiLambdaNeutronPtAsoCuts", " associated #bar{#Lambdan} Cuts", 201, 0., 10.1);
316 fHistAntiLambdaNeutronPtAsoCuts->GetYaxis()->SetTitle("Counts");
317 fHistAntiLambdaNeutronPtAsoCuts->GetXaxis()->SetTitle("#it{p}_{T} (GeV/#it{c})");
319 fHistAntiLambdaNeutronPtAso = new TH1F("fHistAntiLambdaNeutronPtAso", " associated #bar{#Lambdan}", 201, 0., 10.1);
320 fHistAntiLambdaNeutronPtAso->GetYaxis()->SetTitle("Counts");
321 fHistAntiLambdaNeutronPtAso->GetXaxis()->SetTitle("#it{p}_{T} (GeV/#it{c})");
323 fHistLambdaNeutronInvaMassAso = new TH1F("fHistLambdaNeutronInvaMassAso", "Associated #Lambdan", 100, 2.0, 2.1);
324 fHistLambdaNeutronInvaMassAso->GetYaxis()->SetTitle("Counts");
325 fHistLambdaNeutronInvaMassAso->GetXaxis()->SetTitle("Invariant mass (d #pi^{-}) (GeV/#it{c}^{2})");
327 fHistAntiLambdaNeutronInvaMassAso = new TH1F("fHistAntiLambdaNeutronInvaMassAso", " Associated #bar{#Lambdan}", 100, 2.0, 2.1);
328 fHistAntiLambdaNeutronInvaMassAso->GetYaxis()->SetTitle("Counts");
329 fHistAntiLambdaNeutronInvaMassAso->GetXaxis()->SetTitle("Invariant mass (#bar{d} #pi^{+}) (GeV/#it{c}^{2})");
331 fHistLambdaNeutronDecayLengthAso = new TH1F("fHistLambdaNeutronDecayLengthAso", "Associated #Lambdan", 401, 0., 400.1);
332 fHistLambdaNeutronDecayLengthAso->GetYaxis()->SetTitle("Counts");
333 fHistLambdaNeutronDecayLengthAso->GetXaxis()->SetTitle("#it{decay length} (cm)");
335 fHistAntiLambdaNeutronDecayLengthAso = new TH1F("fHistAntiLambdaNeutronDecayLengthAso", "Associated #bar{#Lambdan}", 401, 0., 400.1);
336 fHistAntiLambdaNeutronDecayLengthAso->GetYaxis()->SetTitle("Counts");
337 fHistAntiLambdaNeutronDecayLengthAso->GetXaxis()->SetTitle("#it{decay length} (cm)");
339 //fHistLambdaNeutronPtAso = new TH1F("fHistLambdaNeutronPtAso", "Associated #Lambdan", 201, 0., 10.1);
340 //fHistLambdaNeutronPtAso->GetYaxis()->SetTitle("Counts");
341 //Fhistlambdaneutronptaso->GetXaxis()->SetTitle("#it{p}_{T} (GeV/#it{c})");
344 //------------ Tree and branch definitions ----------------//
345 fTreeV0 = new TTree("tree", "fTreeV0");
346 fTreeV0->Branch("fItrk", &fItrk, "fItrk/I");
348 //fTreeV0->Branch("fV0object",&fV0object,"fV0object[fItrk]");
349 fTreeV0->Branch("fV0finder",fV0finder,"fV0finder[fItrk]/I");
350 fTreeV0->Branch("fkMB",fkMB,"fkMB[fItrk]/I");
351 fTreeV0->Branch("fkCentral",fkCentral,"fkCentral[fItrk]/I");
352 fTreeV0->Branch("fkSemiCentral",fkSemiCentral,"fkSemiCentral[fItrk]/I");
353 fTreeV0->Branch("fkEMCEJE",fkEMCEJE,"fkEMCEJE[fItrk]/I");
354 fTreeV0->Branch("fkEMCEGA",fkEMCEGA,"fkEMCEGA[fItrk]/I");
356 fTreeV0->Branch("fCentralityClass10",fCentralityClass10,"fCentralityClass10[fItrk]/I");
357 fTreeV0->Branch("fCentralityPercentile",fCentralityPercentile,"fCentralityPercentile[fItrk]/I");
358 fTreeV0->Branch("fMultiplicity",fMultiplicity,"fMultipicity[fItrk]/I");
359 fTreeV0->Branch("fRefMultiplicity",fRefMultiplicity,"fRefMultipicity[fItrk]/I");
361 fTreeV0->Branch("fPtotN",fPtotN,"fPtotN[fItrk]/D");
362 fTreeV0->Branch("fPtotP",fPtotP,"fPtotP[fItrk]/D");
363 fTreeV0->Branch("fMotherPt",fMotherPt,"fMotherPt[fItrk]/D");
364 fTreeV0->Branch("fdEdxN",fdEdxN,"fdEdxN[fItrk]/D");
365 fTreeV0->Branch("fdEdxP",fdEdxP,"fdEdxP[fItrk]/D");
366 fTreeV0->Branch("fSignN",fSignN,"fSignN[fItrk]/D");
367 fTreeV0->Branch("fSignP",fSignP,"fSignP[fItrk]/D");
369 fTreeV0->Branch("fSigmadEdxPionPos",fSigmadEdxPionPos,"fSigmadEdxPionPos[fItrk]/D");
370 fTreeV0->Branch("fSigmadEdxPionNeg",fSigmadEdxPionNeg,"fSigmadEdxPionNeg[fItrk]/D");
372 fTreeV0->Branch("fDCAv0",fDCAv0,"fDCAv0[fItrk]/F"); //Dca v0 Daughters
373 fTreeV0->Branch("fCosinePAv0",fCosinePAv0,"fCosinePAv0[fItrk]/F"); //Cosine of Pionting Angle
374 fTreeV0->Branch("fDecayRadiusTree",fDecayRadiusTree,"fDecayRadiusTree[fItrk]/F"); //decay radius
376 fTreeV0->Branch("fChargeComboDeuteronPionTree",fChargeComboDeuteronPionTree,"fChargeComboDeuteronPionTree[fItrk]/I");
377 fTreeV0->Branch("fInvaMassDeuteronPionTree",fInvaMassDeuteronPionTree,"fInvaMassDeuteronPionTree[fItrk]/D"); //invariant mass
379 fTreeV0->Branch("fIsCorrectlyAssociated",fIsCorrectlyAssociated,"fIsCorrectlyAssociated[fItrk]/O"); //associated hypertriton
381 fTreeV0->Branch("fAmenterosAlphaTree",fAmenterosAlphaTree,"fAmenterosAlphaTree[fItrk]/D");
382 fTreeV0->Branch("fAmenterosQtTree",fAmenterosQtTree,"fAmenterosQtTree[fItrk]/D");
383 fTreeV0->Branch("fRotationTree",fRotationTree,"fRotationTree[fItrk]/I");
385 fTreeV0->Branch("fImpactParameterDeuteronPos",fImpactParameterDeuteronPos,"fImpactParameterDeuteronPos[fItrk]/D");
386 fTreeV0->Branch("fImpactParameterDeuteronNeg",fImpactParameterDeuteronNeg,"fImpactParameterDeuteronNeg[fItrk]/D");
387 fTreeV0->Branch("fImpactParameterPionPos",fImpactParameterPionPos,"fImpactParameterPionPos[fItrk]/D");
388 fTreeV0->Branch("fImpactParameterPionNeg",fImpactParameterPionNeg,"fImpactParameterPionNeg[fItrk]/D");
390 fTreeV0->Branch("fImpactParameterDeuteronPosAliKF",fImpactParameterDeuteronPosAliKF,"fImpactParameterDeuteronPosAliKF[fItrk]/D");
391 fTreeV0->Branch("fImpactParameterDeuteronNegAliKF",fImpactParameterDeuteronNegAliKF,"fImpactParameterDeuteronNegAliKF[fItrk]/D");
392 fTreeV0->Branch("fImpactParameterPionPosAliKF",fImpactParameterPionPosAliKF,"fImpactParameterPionPosAliKF[fItrk]/D");
393 fTreeV0->Branch("fImpactParameterPionNegAliKF",fImpactParameterPionNegAliKF,"fImpactParameterPionNegAliKF[fItrk]/D");
395 fTreeV0->Branch("fMinNClustersTPCPos",fMinNClustersTPCPos,"fMinNClustersTPCPos[fItrk]/s");
396 fTreeV0->Branch("fMinNClustersTPCNeg",fMinNClustersTPCNeg,"fMinNClustersTPCNeg[fItrk]/s");
398 fTreeV0->Branch("fMaxChi2PerClusterTPCPos",fMaxChi2PerClusterTPCPos,"fMaxChi2PerClusterTPCPos[fItrk]/D");
399 fTreeV0->Branch("fMaxChi2PerClusterTPCNeg",fMaxChi2PerClusterTPCNeg,"fMaxChi2PerClusterTPCNeg[fItrk]/D");
401 //Armenteros-Podolanski
402 fHistArmenterosPodolanskiDeuteronPion= new TH2F("fHistArmenterosPodolanskiDeuteronPion", "Armenteros-Podolanski d #pi^{-}", 200,-1.0,1.0, 500,0,1);
403 fHistArmenterosPodolanskiDeuteronPion->GetXaxis()->SetTitle("#alpha");
404 fHistArmenterosPodolanskiDeuteronPion->GetYaxis()->SetTitle("q_{t}");
405 fHistArmenterosPodolanskiDeuteronPion->SetMarkerStyle(kFullCircle);
407 fHistArmenterosPodolanskiAntiDeuteronPion= new TH2F("fHistArmenterosPodolanskiAntiDeuteronPion", "Armenteros-Podolanski #bar{d} #pi^{+}", 200,-1.0,1.0, 500,0,1);
408 fHistArmenterosPodolanskiAntiDeuteronPion->GetXaxis()->SetTitle("#alpha");
409 fHistArmenterosPodolanskiAntiDeuteronPion->GetYaxis()->SetTitle("q_{t}");
410 fHistArmenterosPodolanskiAntiDeuteronPion->SetMarkerStyle(kFullCircle);
412 // control histograms
413 fof = new TH2F("fof", "OnTheFlyStatus ",5,0.5,5.5,2,-0.5,1.5);
414 fof->GetYaxis()->SetBinLabel(1,"offline");
415 fof->GetYaxis()->SetBinLabel(2,"onTheFly");
416 fof->GetXaxis()->SetBinLabel(1,"total");
417 fof->GetXaxis()->SetBinLabel(2,"dcaCut");
418 fof->GetXaxis()->SetBinLabel(3,"cosCut");
419 fof->GetXaxis()->SetBinLabel(4,"nucleonPID");
420 fof->GetXaxis()->SetBinLabel(5,"pionPID");
421 //fof->GetXaxis()->SetBinLabel(6,"decayRadiusCut");
422 //fof->SetMarkerStyle(kFullCircle);
424 //histogram to count number of events
425 fHistCentralityClass10 = new TH1F("fHistCentralityClass10", "centrality with class10", 11, -0.5, 10.5);
426 fHistCentralityClass10->GetXaxis()->SetTitle("Centrality");
427 fHistCentralityClass10->GetYaxis()->SetTitle("Entries");
429 fHistCentralityPercentile = new TH1F("fHistCentralityPercentile", "centrality with percentile", 101, -0.1, 100.1);
430 fHistCentralityPercentile->GetXaxis()->SetTitle("Centrality");
431 fHistCentralityPercentile->GetYaxis()->SetTitle("Entries");
433 //trigger statitics histogram
434 fHistTriggerStat = new TH1F("fHistTriggerStat","Trigger statistics", fNTriggers,-0.5,fNTriggers-0.5);
435 const Char_t* aTriggerNames[] = { "kMB", "kCentral", "kSemiCentral", "kEMCEJE", "kEMCEGA" };
436 for ( Int_t ii=0; ii < fNTriggers; ii++ )
437 fHistTriggerStat->GetXaxis()->SetBinLabel(ii+1, aTriggerNames[ii]);
439 fHistTriggerStatAfterEventSelection = new TH1F("fHistTriggerStatAfterEventSelection","Trigger statistics after event selection", fNTriggers,-0.5,fNTriggers-0.5);
440 for ( Int_t ii=0; ii < fNTriggers; ii++ )
441 fHistTriggerStatAfterEventSelection->GetXaxis()->SetBinLabel(ii+1, aTriggerNames[ii]);
444 fHistDeDxQA = new TH3F("fHistDeDxQA", "QA dE/dx", 400, -6.0, 6.0, 500, 0.0, 2000, 15, -0.5, 14.5);
445 fHistDeDxQA->GetYaxis()->SetTitle("TPC Signal");
446 fHistDeDxQA->GetXaxis()->SetTitle("p (GeV/c)");
447 fHistDeDxQA->GetZaxis()->SetBinLabel(0,"all pos v0 tracks");
448 fHistDeDxQA->GetZaxis()->SetBinLabel(1,"all neg v0 tracks");
449 fHistDeDxQA->GetZaxis()->SetBinLabel(2,"all neg deuteron");
450 fHistDeDxQA->GetZaxis()->SetBinLabel(3,"all pos deuteron");
451 fHistDeDxQA->GetZaxis()->SetBinLabel(6,"all selected tracks");
452 fHistDeDxQA->GetZaxis()->SetBinLabel(7,"neg deuteron for deuteron+pion");
453 fHistDeDxQA->GetZaxis()->SetBinLabel(8,"pos pion for deuteron+pion");
454 fHistDeDxQA->GetZaxis()->SetBinLabel(9,"pos deuteron for deuteron+pion");
455 fHistDeDxQA->GetZaxis()->SetBinLabel(10,"neg pion for deuteron+pion");
457 //Define and fill the OutputContainer
458 fOutputContainer = new TObjArray(1);
459 fOutputContainer->SetOwner(kTRUE);
460 fOutputContainer->SetName(GetName()) ;
461 fOutputContainer->Add(fHistArmenterosPodolanskiDeuteronPion);
462 fOutputContainer->Add(fHistArmenterosPodolanskiAntiDeuteronPion);
463 fOutputContainer->Add(fof);
464 fOutputContainer->Add(fHistDeDxQA);
465 fOutputContainer->Add(fHistCentralityClass10);
466 fOutputContainer->Add(fHistCentralityPercentile);
467 fOutputContainer->Add(fHistTriggerStat);
468 fOutputContainer->Add(fHistTriggerStatAfterEventSelection);
469 fOutputContainer->Add(fHistLambdaNeutronPtGen);
470 fOutputContainer->Add(fHistLambdaNeutronPtGenMinBias);
471 fOutputContainer->Add(fHistLambdaNeutronPtGenCentral);
472 fOutputContainer->Add(fHistLambdaNeutronPtGenSemiCentral);
473 fOutputContainer->Add(fHistAntiLambdaNeutronPtGen);
474 fOutputContainer->Add(fHistAntiLambdaNeutronPtGenMinBias);
475 fOutputContainer->Add(fHistAntiLambdaNeutronPtGenCentral);
476 fOutputContainer->Add(fHistAntiLambdaNeutronPtGenSemiCentral);
477 fOutputContainer->Add(fHistLambdaNeutronInvaMassGen);
478 fOutputContainer->Add(fHistAntiLambdaNeutronInvaMassGen);
479 fOutputContainer->Add(fHistLambdaNeutronDecayLengthGen);
480 fOutputContainer->Add(fHistAntiLambdaNeutronDecayLengthGen);
481 fOutputContainer->Add(fHistLambdaNeutronPtAso);
482 fOutputContainer->Add(fHistLambdaNeutronPtAsoCuts);
483 fOutputContainer->Add(fHistAntiLambdaNeutronPtAso);
484 fOutputContainer->Add(fHistAntiLambdaNeutronPtAsoCuts);
485 fOutputContainer->Add(fHistLambdaNeutronInvaMassAso);
486 fOutputContainer->Add(fHistAntiLambdaNeutronInvaMassAso);
487 fOutputContainer->Add(fHistLambdaNeutronDecayLengthAso);
488 fOutputContainer->Add(fHistAntiLambdaNeutronDecayLengthAso);
491 //________________________________________________________________________
492 void AliAnalysisTaskLambdaNAOD::UserExec(Option_t *){
495 // Called for each event
497 //cout << "katze1" << endl;
500 //get Event-Handler for the trigger information
501 fEventHandler= dynamic_cast<AliInputEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
502 if (!fEventHandler) {
503 AliError("Could not get InputHandler");
509 AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
511 //Printf("ERROR: Could not retrieve MC event handler");
515 AliMCEvent* mcEvent = 0x0;
516 AliStack* stack = 0x0;
517 if (eventHandler) mcEvent = eventHandler->MCEvent();
519 //Printf("ERROR: Could not retrieve MC event");
524 stack = mcEvent->Stack();
528 //look for the generated particles
529 //MCGenerated(stack);
531 //cout << "katze2" << endl;
532 if (SetupEvent() < 0) {
533 PostData(1, fOutputContainer);
538 AliESDEvent *fESDevent = 0x0;
539 AliAODEvent *fAODevent = 0x0;
544 AliError("Cannot get pid response");
549 Int_t centralityClass10 = -1;
550 Int_t centralityPercentile = -1;
551 Double_t vertex[3] = {-100.0, -100.0, -100.0};
553 //Initialisation of the event and basic event cuts:
555 if (fAnalysisType == "ESD") {
557 fESDevent=dynamic_cast<AliESDEvent*>(InputEvent());
559 AliWarning("ERROR: fESDevent not available \n");
564 //1.) vertex existence
565 /*const AliESDVertex *vertexESD = fESDevent->GetPrimaryVertexTracks();
566 if (vertexESD->GetNContributors()<1)
569 vertexESD = fESDevent->GetPrimaryVertexSPD();
570 if(vertexESD->GetNContributors()<1) {
571 PostData(1, fOutputContainer);
576 vertexESD->GetXYZ(vertex);
578 //2. vertex position within 10 cm
579 if (TMath::Abs(vertexESD->GetZv()) > 10) return;*/
581 const AliESDVertex *vertexTracks = fESDevent->GetPrimaryVertexTracks();
582 if (vertexTracks->GetNContributors()<1) vertexTracks = 0x0;
584 const AliESDVertex *vertexSPD = fESDevent->GetPrimaryVertexSPD();
585 if (vertexSPD->GetNContributors()<1) vertexSPD = 0x0;
587 //cout << "before" << endl;
589 if(!vertexTracks || !vertexSPD) return;
590 //cout << "after" <<endl;
592 //if (vertexTracks && vertexSPD){
593 //cout << "Vertex: " << TMath::Abs(vertexTracks->GetZv()) << endl;
594 if (TMath::Abs(vertexTracks->GetZv()) > 10 || TMath::Abs(vertexSPD->GetZv()) > 10) return;
598 //centrality selection in PbPb
599 if (fESDevent->GetEventSpecie() == 4) //species == 4 == PbPb
601 AliCentrality *esdCentrality = fESDevent->GetCentrality();
602 centralityClass10 = esdCentrality->GetCentralityClass10("V0M"); // centrality percentile determined with V0
603 centralityPercentile = esdCentrality->GetCentralityPercentile("V0M"); // centrality percentile determined with V0
604 //cout << "************************************EventSpecie " << fESDevent->GetEventSpecie() << " centrality "<< centrality << endl;
606 if (centralityPercentile < 0. || centralityPercentile > 80. ) return; //select only events with centralities between 0 and 80 %
610 //cout << "EventSpecie " << fESDevent->GetEventSpecie() << " centrality "<< centrality << endl;
611 //count centrality classes
612 fHistCentralityClass10->Fill(centralityClass10);
613 fHistCentralityPercentile->Fill(centralityPercentile);
615 if(fTriggerFired[0] == kTRUE)fHistTriggerStatAfterEventSelection->Fill(0);
616 if(fTriggerFired[1] == kTRUE)fHistTriggerStatAfterEventSelection->Fill(1);
617 if(fTriggerFired[2] == kTRUE)fHistTriggerStatAfterEventSelection->Fill(2);
618 if(fTriggerFired[3] == kTRUE)fHistTriggerStatAfterEventSelection->Fill(3);
619 if(fTriggerFired[4] == kTRUE)fHistTriggerStatAfterEventSelection->Fill(4);
624 else if (fAnalysisType == "AOD") {
626 fAODevent = dynamic_cast<AliAODEvent*>( InputEvent() );
628 AliWarning("ERROR: lAODevent not available \n");
632 const AliAODVertex *vertexAOD = fAODevent->GetPrimaryVertex();
633 if (vertexAOD->GetNContributors()<1)
636 vertexAOD = fAODevent->GetPrimaryVertexSPD();
637 if(vertexAOD->GetNContributors()<1) {
638 PostData(1, fOutputContainer);
642 vertexAOD->GetXYZ(vertex);
644 //2. vertex position within 10 cm
645 if (TMath::Abs(vertex[2]) > 10) return;
647 //centrality selection
648 AliCentrality *aodCentrality = fAODevent->GetCentrality();
649 centralityClass10 = aodCentrality->GetCentralityClass10("V0M"); // centrality percentile determined with V0
650 centralityPercentile = aodCentrality->GetCentralityPercentile("V0M");
651 if (centralityPercentile < 0. || centralityPercentile > 80.) return; //select only events with centralities between 0 and 80 %
653 // count number of events
654 fHistCentralityClass10->Fill(centralityClass10);
658 Printf("Analysis type (ESD or AOD) not specified \n");
670 runNumber = (InputEvent())->GetRunNumber();
672 Int_t nTrackMultiplicity = -1;
673 nTrackMultiplicity = (InputEvent())->GetNumberOfTracks();
674 //const AliMulitiplicity *spdMultiplicity = 0x0;
675 //if (fAnalysisType == "ESD")spdMultiplicity= fESDevent->GetMultiplicity();
676 Int_t refMultTpc = -1;
677 //const AliMulitiplicity *spdMultiplicity = 0x0;
682 // EstimateMultiplicity(&multi1,&multi2,&multi3, 0.9,kTRUE, kTRUE);
684 if (fAnalysisType == "ESD"){
685 fESDevent->EstimateMultiplicity(multi1,multi2,multi3,0.9,kTRUE,kTRUE);
686 refMultTpc = AliESDtrackCuts::GetReferenceMultiplicity(fESDevent, kTRUE);
687 // const AliMultiplicity *spdMultiplicity = fESDevent->GetMultiplicity();
688 //cout << "nTrackMultiplicity: " << nTrackMultiplicity << " multi1: " << multi1 << " multi2: " << multi2 << " multi3: " << multi3 << " refMultTpc: " << refMultTpc << endl;
691 //look for the generated particles
692 MCGenerated(stack,refMultTpc);
694 //AliKFParticle settings
695 AliKFParticle::SetField(fInputEvent->GetMagneticField());
696 AliKFVertex primVtx(*(fInputEvent->GetPrimaryVertex()));
698 for (Int_t ivertex=0; ivertex<(InputEvent()->GetNumberOfV0s()); ivertex++) { //beginn v0 loop
700 AliESDv0 * v0ESD = 0x0;
701 AliAODv0 * v0AOD = 0x0;
703 Bool_t v0ChargesCorrect = kTRUE;
704 Bool_t testTrackCuts = kFALSE;
705 //Bool_t testFilterBit = kFALSE;
708 Int_t onFlyStatus = 5;
711 Float_t cosPointing= 2;
712 Float_t decayRadius= -1;
714 AliVTrack * trackN = 0x0;
715 AliVTrack * trackP = 0x0;
717 Double_t ptotN = -1000;
718 Double_t ptotP = -1000;
720 Int_t chargeComboDeuteronPion = -1;
722 Double_t momentumPion[3]={0,0,0};
723 Double_t momentumPionRot[3]={0,0,0};
724 Double_t momentumDeuteron[3]={0,0,0};
725 Double_t momentumDeuteronRot[3]={0,0,0};
726 Double_t momentumMother[3] = {0,0,0};
728 Double_t transversMomentumPion = 0;
729 Double_t transversMomentumDeuteron = 0;
731 Double_t transversMomentumMother = 0;
732 Double_t longitudinalMomentumMother = 0;
734 Double_t totalEnergyMother = 0;
735 Double_t energyPion = 0;
736 Double_t energyDeuteron = 0;
738 Double_t rapidity = 2;
740 TVector3 vecPion(0,0,0);
741 TVector3 vecDeuteron(0,0,0);
742 TVector3 vecMother(0,0,0);
747 Double_t thetaPion = 0;
748 Double_t thetaDeuteron = 0;
751 Double_t invaMassDeuteronPion = 0;
759 fV0finder[fItrk] = -1;
761 fkCentral[fItrk] = -1;
762 fkSemiCentral[fItrk] = -1;
763 fkEMCEJE[fItrk] = -1;
764 fkEMCEGA[fItrk] = -1;
766 fCentralityClass10[fItrk] = -1;
767 fCentralityPercentile[fItrk] = -1;
768 fMultiplicity[fItrk] = -1;
769 fRefMultiplicity[fItrk] = -1;
771 fPtotN[fItrk] = -1000;
772 fPtotP[fItrk] = -1000;
773 fMotherPt[fItrk] = -1000;
780 fSigmadEdxPionPos[fItrk] = -1;
781 fSigmadEdxPionNeg[fItrk] = -1;
784 fCosinePAv0[fItrk] = -2;
785 fDecayRadiusTree[fItrk] = -1;
787 fInvaMassDeuteronPionTree[fItrk] = 0;
788 fChargeComboDeuteronPionTree[fItrk] = -1;
790 fAmenterosAlphaTree[fItrk] = 2;
791 fAmenterosQtTree[fItrk] = -1;
793 fImpactParameterDeuteronPos[fItrk] = -1;
794 fImpactParameterDeuteronNeg[fItrk] = -1;
796 fImpactParameterPionPos[fItrk] = -1;
797 fImpactParameterPionNeg[fItrk] = -1;
799 fImpactParameterDeuteronPosAliKF[fItrk] = -1;
800 fImpactParameterDeuteronNegAliKF[fItrk] = -1;
802 fImpactParameterPionPosAliKF[fItrk] = -1;
803 fImpactParameterPionNegAliKF[fItrk] = -1;
805 fMinNClustersTPCPos[fItrk] = 0;
806 fMinNClustersTPCNeg[fItrk] = 0;
808 fMaxChi2PerClusterTPCPos[fItrk] = -1;
809 fMaxChi2PerClusterTPCNeg[fItrk] = -1;
812 if(fAnalysisType == "ESD")v0ESD = fESDevent->GetV0(ivertex);
813 if(fAnalysisType == "AOD")v0AOD = fAODevent->GetV0(ivertex);
816 //distinguish between the two V0 finders: 0 offline V0 finder / 1 online V0 finder
818 if(fAnalysisType == "ESD") of = v0ESD->GetOnFlyStatus();
819 if(fAnalysisType == "AOD") of = v0AOD->GetOnFlyStatus();
821 if(of) onFlyStatus= 1;
822 if(!of) onFlyStatus= 0;
824 if(onFlyStatus==0)fof->Fill(1,0);
825 if(onFlyStatus==1)fof->Fill(1,1);
827 //Get dca, cos of pointing angle and decay radius
830 if(fAnalysisType == "ESD")
832 dcaV0 = v0ESD->GetDcaV0Daughters();
833 cosPointing = v0ESD->GetV0CosineOfPointingAngle();
834 decayRadius = v0ESD->GetRr();
837 if(fAnalysisType == "AOD")
839 dcaV0 = v0AOD->DcaV0Daughters();
840 cosPointing = v0AOD->CosPointingAngle(vertex);
841 decayRadius = v0AOD->DecayLengthV0(vertex);
845 // select coresponding tracks
846 if(fAnalysisType == "ESD")
848 trackN = fESDevent->GetTrack(v0ESD->GetIndex(0));
849 trackP = fESDevent->GetTrack(v0ESD->GetIndex(1));
851 if (trackN->Charge() > 0) // switch because of bug in V0 interface
853 trackN = fESDevent->GetTrack(v0ESD->GetIndex(1));
854 trackP = fESDevent->GetTrack(v0ESD->GetIndex(0));
855 v0ChargesCorrect = kFALSE;
858 testTrackCuts = TrackCuts(trackN,testTrackCuts);
859 if(testTrackCuts == kFALSE) continue;
860 testTrackCuts = TrackCuts(trackP,testTrackCuts);
861 if(testTrackCuts == kFALSE) continue;
864 if(fAnalysisType == "AOD")
866 trackN = dynamic_cast<AliVTrack*>(v0AOD->GetDaughter(0));
867 trackP = dynamic_cast<AliVTrack*>(v0AOD->GetDaughter(1));
869 if (trackN->Charge() > 0) // switch because of bug in V0 interface
871 trackN = dynamic_cast<AliVTrack*>(v0AOD->GetDaughter(1));
872 trackP = dynamic_cast<AliVTrack*>(v0AOD->GetDaughter(0));
873 v0ChargesCorrect = kFALSE;
875 //Test-Filterbit - only use for track base analysis -- not for V0 candidates -- the AOD V0 candidates sould be a copy of the ESD candidates, BUT NOT for AOD0086 -> here there wqas a wider cos-cutto have more candidates in the AODs
876 //testFilterBit = FilterBit(trackN,testFilterBit);
877 //if(testFilterBit == kFALSE) continue;
878 //testFilterBit = FilterBit(trackP,testFilterBit);
879 //if(testFilterBit == kFALSE) continue;
881 /*testTrackCuts = TrackCuts(trackN,testTrackCuts);
882 if(testTrackCuts == kFALSE) continue;
883 testTrackCuts = TrackCuts(trackP,testTrackCuts);
884 if(testTrackCuts == kFALSE) continue;*/
887 //Get the total momentum for each track ---> at the inner readout of the TPC???? // momenta a always positive
888 if(fAnalysisType == "AOD") {
889 ptotN = trackN->P(); //GetInnerParam()->GetP(); ******************FIXME: InnerParam
890 ptotP = trackP->P(); //GetInnerParam()->GetP(); ******************FIXME: InnerParam
893 if(fAnalysisType == "ESD") {
894 ptotN = MomentumInnerParam(trackN,ptotN);
895 ptotP = MomentumInnerParam(trackP,ptotP);
898 //fill QA dEdx with all V0 candidates
899 if(trackP) fHistDeDxQA->Fill(ptotP*trackP->Charge(), trackP->GetTPCsignal(),0);
900 if(trackN) fHistDeDxQA->Fill(ptotN*trackN->Charge(), trackN->GetTPCsignal(),1);
902 if (dcaV0 > 3 || dcaV0 < 1e-20) continue;
904 //Check how much the dca cut reduces the statistic (background) for the different VO-finder
905 if(onFlyStatus==0)fof->Fill(2,0);
906 if(onFlyStatus==1)fof->Fill(2,1);
908 if (cosPointing < 0.9) continue;
910 //Check how much the cos-of-the-pointing-angle-cut reduces the statistic (background) for the different VO-finder
911 if(onFlyStatus==0)fof->Fill(3,0);
912 if(onFlyStatus==1)fof->Fill(3,1);
914 //deuteron PID via specific energy loss in the TPC
915 Bool_t isDeuteron[3] = {kFALSE,kFALSE,kFALSE}; //0 = posDeuteron, 1 = negDeuteron, 2 = trackN is deuteron
917 DeuteronPID(trackP,trackN,ptotP,ptotN,runNumber,isDeuteron);
919 //if(!isDeuteron[0] && !isDeuteron[1]) continue;
920 if((isDeuteron[0]==kFALSE) && (isDeuteron[1]==kFALSE)) continue;
923 //Check how much the nuclei PID cuts reduce the statistics (background) for the two V0-finders
924 if(onFlyStatus==0)fof->Fill(4,0);
925 if(onFlyStatus==1)fof->Fill(4,1);
927 //Fill the QA dEdx with deuterons and helium3 after the nuclei PID cut
928 if(isDeuteron[1]) fHistDeDxQA->Fill(ptotN*trackN->Charge(), trackN->GetTPCsignal(),2);
929 if(isDeuteron[0]) fHistDeDxQA->Fill(ptotP*trackP->Charge(), trackP->GetTPCsignal(),3);
931 //deuteron PID via specific energy loss in the TPC
932 Bool_t isPion[2] = {kFALSE,kFALSE}; //0 = posPion, 1 = negPion
934 PionPID(trackP,trackN,ptotP,ptotN,runNumber,isPion);
936 //if(isDeuteron[0] && !isPion[1]) continue; //pos deuteron and neg Pion
937 //if(isDeuteron[1] && !isPion[0]) continue; //neg deuteron and pos Pion
939 //Check how much the pion PID cuts reduce the statistics (background) for the two V0-finders
940 if(onFlyStatus==0)fof->Fill(5,0);
941 if(onFlyStatus==1)fof->Fill(5,1);
944 //Save the different charge combinations to differentiat between particles and anit-particles:
945 // -/+ = Anti-deuteron + positive Pion
946 // -/- = Anti-deuteron + negative Pion
947 // +/- = deuteron + negative Pion
948 // +/+ = deuteron + positive Pion
951 if (trackN->Charge()<0 && trackP->Charge()<0 && isDeuteron[1]==kTRUE && isPion[1]==kTRUE) chargeComboDeuteronPion = 1; // -/-
952 if (trackN->Charge()>0 && trackP->Charge()>0 && isDeuteron[0]==kTRUE && isPion[0]==kTRUE ) chargeComboDeuteronPion = 3; // +/+
955 if (trackN->Charge()<0 && trackP->Charge()>0 && isDeuteron[1]==kTRUE && isPion[0]==kTRUE ) chargeComboDeuteronPion = 0; // -/+
956 //if (trackN->Charge()>0 && trackP->Charge()<0 && isDeuteron[1]==kTRUE && isPion[0]==kTRUE ) chargeComboDeuteronPion = 0; // -/+ //should not exist because of charge correction
957 //if (trackN->Charge()>0 && trackP->Charge()<0 && isDeuteron[0]==kTRUE && isPion[1]==kTRUE) chargeComboDeuteronPion = 2; // +/- //should not exist because of charge correction
958 if (trackN->Charge()<0 && trackP->Charge()>0 && isDeuteron[0]==kTRUE && isPion[1]==kTRUE) chargeComboDeuteronPion = 2; // +/-
960 //if(chargeComboDeuteronPion==0) cout << "chargeN: " << trackN->Charge() << " chargeP: " << trackP->Charge() << endl;
963 //Fill the QA dEdx with all selected tracks after the PID cuts
964 fHistDeDxQA->Fill(ptotP*trackP->Charge(), trackP->GetTPCsignal(),6);
965 fHistDeDxQA->Fill(ptotN*trackN->Charge(), trackN->GetTPCsignal(),6);
969 //get the momenta of the daughters
971 if(fAnalysisType == "ESD")
973 if(chargeComboDeuteronPion==0){ //anti-deuteron (isDeuteron[1]==kTRUE), positives pion (isPion[0]==kTRUE), trackN corresponds to anti-deuteron, tackP corresponds to pion
975 v0ESD->GetPPxPyPz(momentumPion[0], momentumPion[1], momentumPion[2]);
976 v0ESD->GetNPxPyPz(momentumDeuteron[0], momentumDeuteron[1], momentumDeuteron[2]);
977 if (!v0ChargesCorrect) {
979 v0ESD->GetNPxPyPz(momentumPion[0], momentumPion[1], momentumPion[2]);
980 v0ESD->GetPPxPyPz(momentumDeuteron[0], momentumDeuteron[1], momentumDeuteron[2]);
984 if(chargeComboDeuteronPion==2){ //deuteron (isDeuteron[0]==kTRUE), negative pion (isPion[1]==kTRUE), trackP corresponds to deuteron, tackN corresponds to pion
986 v0ESD->GetNPxPyPz(momentumPion[0], momentumPion[1], momentumPion[2]);
987 v0ESD->GetPPxPyPz(momentumDeuteron[0], momentumDeuteron[1], momentumDeuteron[2]);
988 if (!v0ChargesCorrect){
990 v0ESD->GetPPxPyPz(momentumPion[0], momentumPion[1], momentumPion[2]);
991 v0ESD->GetNPxPyPz(momentumDeuteron[0], momentumDeuteron[1], momentumDeuteron[2]);
995 if(chargeComboDeuteronPion==1 || chargeComboDeuteronPion==3){ //trackN corresponds to deuteron or pion, trackP corresponds to deuteron or pion
996 if(isDeuteron[2]==kTRUE){ //trackN corresponds to deuteron, trackP corresponds to pion
998 v0ESD->GetPPxPyPz(momentumPion[0], momentumPion[1], momentumPion[2]);
999 v0ESD->GetNPxPyPz(momentumDeuteron[0], momentumDeuteron[1], momentumDeuteron[2]);
1000 if (!v0ChargesCorrect) {
1002 v0ESD->GetNPxPyPz(momentumPion[0], momentumPion[1], momentumPion[2]);
1003 v0ESD->GetPPxPyPz(momentumDeuteron[0], momentumDeuteron[1], momentumDeuteron[2]);
1006 if(isDeuteron[2]==kFALSE){ //trackP corresponds tp deuteron, trackN corresponds to pion
1008 v0ESD->GetNPxPyPz(momentumPion[0], momentumPion[1], momentumPion[2]);
1009 v0ESD->GetPPxPyPz(momentumDeuteron[0], momentumDeuteron[1], momentumDeuteron[2]);
1010 if (!v0ChargesCorrect) {
1012 v0ESD->GetPPxPyPz(momentumPion[0], momentumPion[1], momentumPion[2]);
1013 v0ESD->GetNPxPyPz(momentumDeuteron[0], momentumDeuteron[1], momentumDeuteron[2]);
1019 //get the momenta of the mother
1020 v0ESD->GetPxPyPz(momentumMother[0],momentumMother[1],momentumMother[2]);
1024 if(fAnalysisType == "AOD")
1026 if(chargeComboDeuteronPion==0){ //anti-deuteron (isDeuteron[1]==kTRUE), positives pion (isPion[0]==kTRUE), trackN corresponds to anti-deuteron, tackP corresponds to pion
1028 momentumPion[0] = v0AOD->MomPosX();
1029 momentumPion[1] = v0AOD->MomPosY();
1030 momentumPion[2] = v0AOD->MomPosZ();
1032 momentumDeuteron[0] = v0AOD->MomNegX();
1033 momentumDeuteron[1] = v0AOD->MomNegY();
1034 momentumDeuteron[2] = v0AOD->MomNegZ();
1036 if (!v0ChargesCorrect){
1038 momentumPion[0] = v0AOD->MomNegX();
1039 momentumPion[1] = v0AOD->MomNegY();
1040 momentumPion[2] = v0AOD->MomNegZ();
1042 momentumDeuteron[0] = v0AOD->MomPosX();
1043 momentumDeuteron[1] = v0AOD->MomPosY();
1044 momentumDeuteron[2] = v0AOD->MomPosZ();
1048 if (chargeComboDeuteronPion==2){ //deuteron (isDeuteron[0]==kTRUE), negative pion (isPion[1]==kTRUE), trackP corresponds to deuteron, tackN corresponds to pion
1050 momentumPion[0] = v0AOD->MomNegX();
1051 momentumPion[1] = v0AOD->MomNegY();
1052 momentumPion[2] = v0AOD->MomNegZ();
1054 momentumDeuteron[0] = v0AOD->MomPosX();
1055 momentumDeuteron[1] = v0AOD->MomPosY();
1056 momentumDeuteron[2] = v0AOD->MomPosZ();
1058 if (!v0ChargesCorrect){
1060 momentumPion[0] = v0AOD->MomPosX();
1061 momentumPion[1] = v0AOD->MomPosY();
1062 momentumPion[2] = v0AOD->MomPosZ();
1064 momentumDeuteron[0] = v0AOD->MomNegX();
1065 momentumDeuteron[1] = v0AOD->MomNegY();
1066 momentumDeuteron[2] = v0AOD->MomNegZ();
1070 if(chargeComboDeuteronPion==1 || chargeComboDeuteronPion==3){ //trackN correponds to deuteron or pion, trackP corresponds to deuteron or pion
1071 if(isDeuteron[2]==kTRUE){ //trackN correponds to deuteron, trackP corresponds to pion
1073 momentumPion[0] = v0AOD->MomPosX();
1074 momentumPion[1] = v0AOD->MomPosY();
1075 momentumPion[2] = v0AOD->MomPosZ();
1077 momentumDeuteron[0] = v0AOD->MomNegX();
1078 momentumDeuteron[1] = v0AOD->MomNegY();
1079 momentumDeuteron[2] = v0AOD->MomNegZ();
1081 if (!v0ChargesCorrect){
1083 momentumPion[0] = v0AOD->MomNegX();
1084 momentumPion[1] = v0AOD->MomNegY();
1085 momentumPion[2] = v0AOD->MomNegZ();
1087 momentumDeuteron[0] = v0AOD->MomPosX();
1088 momentumDeuteron[1] = v0AOD->MomPosY();
1089 momentumDeuteron[2] = v0AOD->MomPosZ();
1093 if(isDeuteron[2]==kFALSE){ //trackP corresponds to deuteron, trackN corresponds to pion
1095 momentumPion[0] = v0AOD->MomNegX();
1096 momentumPion[1] = v0AOD->MomNegY();
1097 momentumPion[2] = v0AOD->MomNegZ();
1099 momentumDeuteron[0] = v0AOD->MomPosX();
1100 momentumDeuteron[1] = v0AOD->MomPosY();
1101 momentumDeuteron[2] = v0AOD->MomPosZ();
1103 if (!v0ChargesCorrect){
1105 momentumPion[0] = v0AOD->MomPosX();
1106 momentumPion[1] = v0AOD->MomPosY();
1107 momentumPion[2] = v0AOD->MomPosZ();
1109 momentumDeuteron[0] = v0AOD->MomNegX();
1110 momentumDeuteron[1] = v0AOD->MomNegY();
1111 momentumDeuteron[2] = v0AOD->MomNegZ();
1115 //get the momenta of the mother
1116 momentumMother[0] = v0AOD->MomV0X();
1117 momentumMother[1] = v0AOD->MomV0Y();
1118 momentumMother[2] = v0AOD->MomV0Z();
1123 transversMomentumPion = TMath::Sqrt(momentumPion[0]*momentumPion[0]+momentumPion[1]*momentumPion[1]);
1124 transversMomentumDeuteron = TMath::Sqrt(momentumDeuteron[0]*momentumDeuteron[0]+momentumDeuteron[1]*momentumDeuteron[1]);
1126 transversMomentumMother=TMath::Sqrt((momentumDeuteron[0]+momentumPion[0])*(momentumDeuteron[0]+momentumPion[0])+(momentumDeuteron[1]+momentumPion[1])*(momentumDeuteron[1]+momentumPion[1]));
1128 longitudinalMomentumMother = TMath::Sqrt((momentumDeuteron[2]+momentumPion[2])*(momentumDeuteron[2]+momentumPion[2]));
1130 energyDeuteron = TMath::Sqrt(momentumDeuteron[0]*momentumDeuteron[0]+momentumDeuteron[1]*momentumDeuteron[1]+momentumDeuteron[2]*momentumDeuteron[2]+fgkMass[kMassDeuteron]*fgkMass[kMassDeuteron]);
1132 energyPion = TMath::Sqrt(momentumPion[0]*momentumPion[0]+momentumPion[1]*momentumPion[1]+momentumPion[2]*momentumPion[2]+fgkMass[kMassPion]*fgkMass[kMassPion]);
1134 totalEnergyMother = energyPion + energyDeuteron;
1137 rapidity = 0.5 * TMath::Log((totalEnergyMother+longitudinalMomentumMother)/( totalEnergyMother-longitudinalMomentumMother));
1139 if(rapidity > 1.0 || rapidity < -1) continue;
1142 //Armanteros-Podolanski
1143 vecPion.SetXYZ(momentumPion[0],momentumPion[1],momentumPion[2]);
1144 vecDeuteron.SetXYZ(momentumDeuteron[0],momentumDeuteron[1],momentumDeuteron[2]);
1145 vecMother.SetXYZ(momentumMother[0],momentumMother[1],momentumMother[2]);
1147 thetaPion = acos((vecPion * vecMother)/(vecPion.Mag() * vecMother.Mag()));
1148 thetaDeuteron = acos((vecDeuteron * vecMother)/(vecDeuteron.Mag() * vecMother.Mag()));
1150 alpha = ((vecPion.Mag())*cos(thetaPion)-(vecDeuteron.Mag())*cos(thetaDeuteron))/
1151 ((vecDeuteron.Mag())*cos(thetaDeuteron)+(vecPion.Mag())*cos(thetaPion)) ;
1152 qt = vecDeuteron.Mag()*sin(thetaDeuteron);
1155 //Rotation for background calculation
1156 //Int_t rotation=1; // =1 signal, =2 Rotation of the pion , =3 Rotation of the deuteron
1158 //Double_t fStartAnglePhi=TMath::Pi();
1159 //Double_t fConeAnglePhi=TMath::Pi(); //-0.174;
1160 //phi = fStartAnglePhi+(2*gRandom->Rndm()-1)*fConeAnglePhi;
1162 for(Int_t rotation=1;rotation<4;rotation++){ //loop for rotation
1164 Double_t fStartAnglePhi=TMath::Pi();
1165 Double_t fConeAnglePhi=TMath::Pi(); //-0.174;
1166 phi = fStartAnglePhi+(2*gRandom->Rndm()-1)*fConeAnglePhi;
1168 //calculate new rotated momenta
1169 momentumPionRot[0]=TMath::Cos(phi)*momentumPion[0]-TMath::Sin(phi)*momentumPion[1];
1170 momentumPionRot[1]=TMath::Sin(phi)*momentumPion[0]+TMath::Cos(phi)*momentumPion[1];
1172 momentumDeuteronRot[0]=TMath::Cos(phi)*momentumDeuteron[0]-TMath::Sin(phi)*momentumDeuteron[1];
1173 momentumDeuteronRot[1]=TMath::Sin(phi)*momentumDeuteron[0]+TMath::Cos(phi)*momentumDeuteron[1];
1175 //invariant mass calculations
1176 fIsCorrectlyAssociated[fItrk] = kFALSE;
1178 //factor for the invariant mass calculation, which only include the pion
1179 e12 = fgkMass[kMassPion]*fgkMass[kMassPion]
1180 +momentumPion[0]*momentumPion[0]
1181 +momentumPion[1]*momentumPion[1]
1182 +momentumPion[2]*momentumPion[2];
1184 r12 = fgkMass[kMassPion]*fgkMass[kMassPion]
1185 +momentumPionRot[0]*momentumPionRot[0]
1186 +momentumPionRot[1]*momentumPionRot[1]
1187 +momentumPion[2]*momentumPion[2];
1189 //factor for the invariant mass calculation, which only include the deuterons
1190 d22 = fgkMass[kMassDeuteron]*fgkMass[kMassDeuteron]
1191 +momentumDeuteron[0]*momentumDeuteron[0]
1192 +momentumDeuteron[1]*momentumDeuteron[1]
1193 +momentumDeuteron[2]*momentumDeuteron[2];
1194 dr22 = fgkMass[kMassDeuteron]*fgkMass[kMassDeuteron]
1195 +momentumDeuteronRot[0]*momentumDeuteronRot[0]
1196 +momentumDeuteronRot[1]*momentumDeuteronRot[1]
1197 +momentumDeuteron[2]*momentumDeuteron[2];
1199 if(rotation == 1){ //signal
1200 invaMassDeuteronPion = TMath::Sqrt(TMath::Max(fgkMass[kMassPion]*fgkMass[kMassPion]
1201 +fgkMass[kMassDeuteron]*fgkMass[kMassDeuteron]
1202 + 2.*(TMath::Sqrt(e12*d22)
1203 -momentumPion[0]*momentumDeuteron[0]
1204 -momentumPion[1]*momentumDeuteron[1]
1205 -momentumPion[2]*momentumDeuteron[2]), 0.));
1211 labelN = trackN->GetLabel();
1212 labelP = trackP->GetLabel();
1214 TParticle *tparticleDaughterN = stack->Particle(TMath::Abs(labelN));
1215 TParticle *tparticleDaughterP = stack->Particle(TMath::Abs(labelP));
1217 Int_t labelMotherN = tparticleDaughterN->GetFirstMother();
1218 Int_t labelMotherP = tparticleDaughterP->GetFirstMother();
1220 TParticle *tparticleMotherN = stack->Particle(TMath::Abs(labelMotherN));
1221 TParticle *tparticleMotherP = stack->Particle(TMath::Abs(labelMotherP));
1224 if(tparticleMotherN->GetPdgCode() == fgkPdgCode[kPDGLambdaNeutron] && tparticleMotherP->GetPdgCode() == fgkPdgCode[kPDGLambdaNeutron] && onFlyStatus ==1 && labelMotherN == labelMotherP ){//check mother PDG and fill the histogramms only for the only V0 finder
1226 //cout << "number of daughters: " << tparticleMotherN->GetNDaughters() << endl;;
1227 if(tparticleMotherN->GetNDaughters() > 2.) continue;
1229 Int_t labelSecondDaughter = tparticleMotherP->GetDaughter(1);
1230 Int_t labelFirstDaughter = labelSecondDaughter-1;
1232 TParticle *tparticleFirstDaughter = stack->Particle(TMath::Abs(labelFirstDaughter));
1233 TParticle *tparticleSecondDaughter = stack->Particle(TMath::Abs(labelSecondDaughter));
1235 if(tparticleFirstDaughter->GetPdgCode() == fgkPdgCode[kPDGDeuteron]){//check first daughter PDG
1237 if(tparticleSecondDaughter->GetPdgCode() == fgkPdgCode[kPDGPionMinus]){//check second daughter PDG
1239 if(invaMassDeuteronPion < 2.02) continue;
1241 fHistLambdaNeutronPtAso->Fill(transversMomentumMother);
1242 fHistLambdaNeutronInvaMassAso->Fill(invaMassDeuteronPion);
1243 fIsCorrectlyAssociated[fItrk] = kTRUE;
1244 fHistLambdaNeutronDecayLengthAso->Fill((decayRadius*2.054)/(tparticleMotherP->P()));
1248 (cosPointing > 0.999) &&
1249 (decayRadius < 50.0) &&
1250 (ptotP > 0.2)){//trackP->GetTPCsignal()>180 &&&& decayRadius > 1.5 && decayRadius< 50
1252 fHistLambdaNeutronPtAsoCuts->Fill(transversMomentumMother);
1254 }//end check second daughter PDG
1255 }//end check first daughter PDG
1256 }//end LambdaNeutron
1258 //Anti-LambdaNeutron
1259 if(tparticleMotherN->GetPdgCode() == fgkPdgCode[kPDGAntiLambdaNeutron] && tparticleMotherP->GetPdgCode() == fgkPdgCode[kPDGAntiLambdaNeutron] && onFlyStatus ==1 && labelMotherN == labelMotherP){//check mother PDG and fill the histogramms only for the only V0 finder
1261 Int_t labelSecondDaughter = tparticleMotherN->GetDaughter(1);
1262 Int_t labelFirstDaughter = labelSecondDaughter-1;
1264 TParticle *tparticleFirstDaughter = stack->Particle(TMath::Abs(labelFirstDaughter));
1265 TParticle *tparticleSecondDaughter = stack->Particle(TMath::Abs(labelSecondDaughter));
1267 if(tparticleFirstDaughter->GetPdgCode() == fgkPdgCode[kPDGAntiDeuteron]) {//check first daughter PDG
1269 if(tparticleSecondDaughter->GetPdgCode() == fgkPdgCode[kPDGPionPlus]) {//check second daughter PDG
1271 if(invaMassDeuteronPion < 2.02) continue;
1273 fIsCorrectlyAssociated[fItrk] = kTRUE;
1274 fHistAntiLambdaNeutronPtAso->Fill(transversMomentumMother);
1275 fHistAntiLambdaNeutronInvaMassAso->Fill(invaMassDeuteronPion);
1276 fHistAntiLambdaNeutronDecayLengthAso->Fill((decayRadius*2.054)/(tparticleMotherP->P()));
1279 if(dcaV0 < 1. && cosPointing > 0.99 && ptotN > 0.2)//&& trackN->GetTPCsignal()>180 && decayRadius > 1.5 && decayRadius< 50
1281 fHistAntiLambdaNeutronPtAsoCuts->Fill(transversMomentumMother);
1283 }//end check second daughter PDG
1284 }//end check first daughter PDG
1285 }//end Anti-LambdaNeutron
1287 }//end rotation == 1, signal
1289 if(rotation == 2){ // rotation of the pion
1290 invaMassDeuteronPion = TMath::Sqrt(TMath::Max(fgkMass[kMassPion]*fgkMass[kMassPion]
1291 +fgkMass[kMassDeuteron]*fgkMass[kMassDeuteron]
1292 + 2.*(TMath::Sqrt(r12*d22)
1293 -momentumPionRot[0]*momentumDeuteron[0]
1294 -momentumPionRot[1]*momentumDeuteron[1]
1295 -momentumPion[2]*momentumDeuteron[2]), 0.));
1296 }//end rotation == 2, rotation of the pion
1298 if(rotation == 3){// Rotation of the deuteron
1299 invaMassDeuteronPion = TMath::Sqrt(TMath::Max(fgkMass[kMassPion]*fgkMass[kMassPion]
1300 +fgkMass[kMassDeuteron]*fgkMass[kMassDeuteron]
1301 + 2.*(TMath::Sqrt(e12*dr22)
1302 -momentumPion[0]*momentumDeuteronRot[0]
1303 -momentumPion[1]*momentumDeuteronRot[1]
1304 -momentumPion[2]*momentumDeuteron[2]), 0.));
1305 }//end rotation == 3, rotation of the deuteron
1307 //fill the THnSparse and the tree variables
1309 //tree variables which are independent of the particle-species
1310 /*fV0finder[fItrk] = onFlyStatus;
1311 fkMB[fItrk] = fTriggerFired[0];
1312 fkCentral[fItrk] = fTriggerFired[1];
1313 fkSemiCentral[fItrk] = fTriggerFired[2];
1314 fkEMCEJE[fItrk] = fTriggerFired[3];
1315 fkEMCEGA[fItrk] = fTriggerFired[4];
1317 fPtotN[fItrk] = ptotN;
1318 fPtotP[fItrk] = ptotP;
1319 fMotherPt[fItrk] = transversMomentumMother;
1321 fdEdxN[fItrk] = trackN->GetTPCsignal();
1322 fdEdxP[fItrk] = trackP->GetTPCsignal();
1323 fSignN[fItrk] = trackN->Charge();
1324 fSignP[fItrk] = trackP->Charge();
1326 fDCAv0[fItrk] = dcaV0;
1327 fCosinePAv0[fItrk] = cosPointing;
1328 fDecayRadiusTree[fItrk] = decayRadius;
1330 fAmenterosAlphaTree[fItrk] = alpha;
1331 fAmenterosQtTree[fItrk] = qt;
1332 fRotationTree[fItrk] = rotation;*/
1335 if (isDeuteron[0] == kTRUE) //pos deuteron
1337 /* fInvaMassDeuteronPionTree[fItrk] = invaMassDeuteronPion;
1338 fChargeComboDeuteronPionTree[fItrk] = chargeComboDeuteronPion;
1342 if(invaMassDeuteronPion < 2.1)
1344 if(chargeComboDeuteronPion == 2)fHistDeDxQA->Fill(ptotP*(trackP->Charge()), trackP->GetTPCsignal(),9);
1345 if(chargeComboDeuteronPion == 2)fHistDeDxQA->Fill(ptotN*(trackN->Charge()), trackN->GetTPCsignal(),10);
1347 fHistArmenterosPodolanskiDeuteronPion->Fill(alpha*(-1),qt);
1350 if (isDeuteron[1] == kTRUE)
1352 /*fInvaMassDeuteronPionTree[fItrk] = invaMassDeuteronPion;
1353 fChargeComboDeuteronPionTree[fItrk] = chargeComboDeuteronPion;
1357 if(invaMassDeuteronPion < 2.1)
1359 if(chargeComboDeuteronPion == 0)fHistDeDxQA->Fill(ptotN*trackN->Charge(), trackN->GetTPCsignal(),7);
1360 if(chargeComboDeuteronPion == 0)fHistDeDxQA->Fill(ptotP*trackP->Charge(), trackP->GetTPCsignal(),8);
1362 fHistArmenterosPodolanskiAntiDeuteronPion->Fill(alpha,qt);
1365 UShort_t numberOfTPCclustersPos = 0;
1366 UShort_t numberOfTPCclustersNeg = 0;
1368 numberOfTPCclustersPos = TPCclusters(trackP,numberOfTPCclustersPos);
1369 numberOfTPCclustersNeg = TPCclusters(trackN,numberOfTPCclustersNeg);
1371 Double_t numberOfChi2clustersTPCPos = 10.;
1372 Double_t numberOfChi2clustersTPCNeg = 10.;
1374 numberOfChi2clustersTPCPos = TPCchi2(trackP,numberOfChi2clustersTPCPos,numberOfTPCclustersPos);
1375 numberOfChi2clustersTPCNeg = TPCchi2(trackN,numberOfChi2clustersTPCNeg,numberOfTPCclustersNeg);
1378 Double_t dcaPosToVertex = 0;
1379 Double_t dcaNegToVertex = 0;
1381 dcaPosToVertex = ImpactParameter(trackP,dcaPosToVertex);
1382 dcaNegToVertex = ImpactParameter(trackN,dcaNegToVertex);
1384 //ImpactParameter with AliKF
1385 //Get the primary vertex from the ITS for the following calculations, because it is the most precise determination of the primary vertex
1386 Double_t itsVertex[3]={0.,0.,0.};
1388 if(fAnalysisType == "ESD"){
1389 itsVertex[0]=fESDevent->GetPrimaryVertexSPD()->GetX();
1390 itsVertex[1]=fESDevent->GetPrimaryVertexSPD()->GetY();
1391 itsVertex[2]=fESDevent->GetPrimaryVertexSPD()->GetZ();
1394 Double_t impactParameterPionPos = -1;
1395 Double_t impactParameterPionNeg = -1;
1397 Double_t impactParameterDeuteronPos = -1;
1398 Double_t impactParameterDeuteronNeg = -1;
1400 Double_t impactParameterPionPosAliKF = -1;
1401 Double_t impactParameterPionNegAliKF = -1;
1403 Double_t impactParameterDeuteronPosAliKF = -1;
1404 Double_t impactParameterDeuteronNegAliKF = -1;
1406 AliKFParticle* negPionKF;
1407 AliKFParticle* posPionKF;
1408 AliKFParticle* negDeuteronKF;
1409 AliKFParticle* posDeuteronKF;
1411 if(chargeComboDeuteronPion==0){ //anti-deuteron (isDeuteron[1]==kTRUE), positives pion (isPion[0]==kTRUE), trackN corresponds to anti-deuteron, tackP corresponds to pion
1413 posPionKF = new AliKFParticle (*(trackP),fgkPdgCode[kPDGPionPlus]);
1414 impactParameterPionPosAliKF = posPionKF->GetDistanceFromVertex(itsVertex);
1415 impactParameterPionPos = dcaPosToVertex;
1416 negDeuteronKF = new AliKFParticle (*(trackN),fgkPdgCode[kPDGPionMinus]);
1417 impactParameterDeuteronNegAliKF = negDeuteronKF->GetDistanceFromVertex(itsVertex);
1418 impactParameterDeuteronNeg = dcaNegToVertex;
1421 if(chargeComboDeuteronPion==2){ //deuteron (isDeuteron[0]==kTRUE), negative pion (isPion[1]==kTRUE), trackP corresponds to deuteron, tackN corresponds to pion
1423 negPionKF = new AliKFParticle (*(trackN),fgkPdgCode[kPDGPionMinus]);
1424 impactParameterPionNegAliKF = negPionKF->GetDistanceFromVertex(itsVertex);
1425 impactParameterPionNeg = dcaNegToVertex;
1426 posDeuteronKF = new AliKFParticle (*(trackP),fgkPdgCode[kPDGPionPlus]);
1427 impactParameterDeuteronPosAliKF = posDeuteronKF->GetDistanceFromVertex(itsVertex);
1428 impactParameterDeuteronPos = dcaPosToVertex;
1431 if(chargeComboDeuteronPion==1){ //trackN corresponds to deuteron or pion, trackP corresponds to deuteron or pion
1432 if(isDeuteron[2]==kTRUE){ //trackN corresponds to deuteron, trackP corresponds to pion
1434 negPionKF = new AliKFParticle (*(trackP),fgkPdgCode[kPDGPionPlus]);
1435 impactParameterPionNegAliKF = negPionKF->GetDistanceFromVertex(itsVertex);
1436 impactParameterPionNeg = dcaPosToVertex;
1437 negDeuteronKF = new AliKFParticle (*(trackN),fgkPdgCode[kPDGPionMinus]);
1438 impactParameterDeuteronNegAliKF = negDeuteronKF->GetDistanceFromVertex(itsVertex);
1439 impactParameterDeuteronNeg = dcaNegToVertex;
1441 if(isDeuteron[2]==kFALSE){ //trackP corresponds tp deuteron, trackN corresponds to pion
1443 negPionKF = new AliKFParticle (*(trackN),fgkPdgCode[kPDGPionPlus]);
1444 impactParameterPionNegAliKF = negPionKF->GetDistanceFromVertex(itsVertex);
1445 impactParameterPionNeg = dcaNegToVertex;
1446 negDeuteronKF = new AliKFParticle (*(trackP),fgkPdgCode[kPDGPionPlus]);
1447 impactParameterDeuteronNegAliKF = negDeuteronKF->GetDistanceFromVertex(itsVertex);
1448 impactParameterDeuteronNeg = dcaPosToVertex;
1452 if(chargeComboDeuteronPion==3){ //trackN corresponds to deuteron or pion, trackP corresponds to deuteron or pion
1453 if(isDeuteron[2]==kTRUE){ //trackN corresponds to deuteron, trackP corresponds to pion
1455 posPionKF = new AliKFParticle (*(trackP),fgkPdgCode[kPDGPionPlus]);
1456 impactParameterPionPosAliKF = posPionKF->GetDistanceFromVertex(itsVertex);
1457 impactParameterPionPos = dcaPosToVertex;
1458 posDeuteronKF = new AliKFParticle (*(trackN),fgkPdgCode[kPDGPionMinus]);
1459 impactParameterDeuteronPosAliKF = posDeuteronKF->GetDistanceFromVertex(itsVertex);
1460 impactParameterDeuteronPos = dcaNegToVertex;
1462 if(isDeuteron[2]==kFALSE){ //trackP corresponds tp deuteron, trackN corresponds to pion
1464 posPionKF = new AliKFParticle (*(trackN),fgkPdgCode[kPDGPionPlus]);
1465 impactParameterPionPosAliKF = posPionKF->GetDistanceFromVertex(itsVertex);
1466 impactParameterPionPos = dcaNegToVertex;
1467 posDeuteronKF = new AliKFParticle (*(trackP),fgkPdgCode[kPDGPionPlus]);
1468 impactParameterDeuteronPosAliKF = posDeuteronKF->GetDistanceFromVertex(itsVertex);
1469 impactParameterDeuteronPos = dcaPosToVertex;
1473 if((isDeuteron[0] == kTRUE) || (isDeuteron[1] == kTRUE)){
1474 //tree variables which are independent of the particle-species
1475 fV0finder[fItrk] = onFlyStatus;
1476 fkMB[fItrk] = fTriggerFired[0];
1477 fkCentral[fItrk] = fTriggerFired[1];
1478 fkSemiCentral[fItrk] = fTriggerFired[2];
1479 fkEMCEJE[fItrk] = fTriggerFired[3];
1480 fkEMCEGA[fItrk] = fTriggerFired[4];
1482 fCentralityClass10[fItrk]= centralityClass10;
1483 fCentralityPercentile[fItrk]= centralityPercentile;
1484 fMultiplicity[fItrk] = multi1;
1485 fRefMultiplicity[fItrk] = refMultTpc;
1487 fPtotN[fItrk] = ptotN;
1488 fPtotP[fItrk] = ptotP;
1489 fMotherPt[fItrk] = transversMomentumMother;
1491 fdEdxN[fItrk] = trackN->GetTPCsignal();
1492 fdEdxP[fItrk] = trackP->GetTPCsignal();
1493 fSignN[fItrk] = trackN->Charge();
1494 fSignP[fItrk] = trackP->Charge();
1496 fSigmadEdxPionPos[fItrk] = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(trackP, AliPID::kPion));
1497 fSigmadEdxPionNeg[fItrk] = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(trackN, AliPID::kPion));
1499 fDCAv0[fItrk] = dcaV0;
1500 fCosinePAv0[fItrk] = cosPointing;
1501 fDecayRadiusTree[fItrk] = decayRadius;
1503 fAmenterosAlphaTree[fItrk] = alpha;
1504 fAmenterosQtTree[fItrk] = qt;
1505 fRotationTree[fItrk] = rotation;
1507 fInvaMassDeuteronPionTree[fItrk] = invaMassDeuteronPion;
1508 fChargeComboDeuteronPionTree[fItrk] = chargeComboDeuteronPion;
1510 fImpactParameterDeuteronPos[fItrk] = impactParameterDeuteronPos;
1511 fImpactParameterDeuteronNeg[fItrk] = impactParameterDeuteronNeg;
1513 fImpactParameterPionPos[fItrk] = impactParameterPionPos;
1514 fImpactParameterPionNeg[fItrk] = impactParameterPionNeg;
1516 fImpactParameterDeuteronPosAliKF[fItrk] = impactParameterDeuteronPosAliKF;
1517 fImpactParameterDeuteronNegAliKF[fItrk] = impactParameterDeuteronNegAliKF;
1519 fImpactParameterPionPosAliKF[fItrk] = impactParameterPionPosAliKF;
1520 fImpactParameterPionNegAliKF[fItrk] = impactParameterPionNegAliKF;
1522 fMinNClustersTPCPos[fItrk] = numberOfTPCclustersPos;
1523 fMinNClustersTPCNeg[fItrk] = numberOfTPCclustersNeg;
1525 // cout << "numberOfTPCclustersPos: " << numberOfTPCclustersPos << " fMinNClustersTPCPos[fItrk]: " << fMinNClustersTPCPos[fItrk] << endl;
1527 fMaxChi2PerClusterTPCPos[fItrk] = numberOfChi2clustersTPCPos;
1528 fMaxChi2PerClusterTPCNeg[fItrk] = numberOfChi2clustersTPCNeg;
1532 //cout << "fItrk: " << fItrk-1 << " fChargeComboDeuteronPionTree[fItrk]: " << fChargeComboDeuteronPionTree[fItrk-1] << " fV0finder[fItrk]: " << fV0finder[fItrk-1] << " fInvaMassDeuteronPionTree[fItrk]: " << fInvaMassDeuteronPionTree[fItrk-1] << " fMaxChi2PerClusterTPCPos[fItrk]: " << fMaxChi2PerClusterTPCPos[fItrk-1] << " numberOfTPCclustersPos: " << numberOfTPCclustersPos << " fMinNClustersTPCPos[fItrk]: " << fMinNClustersTPCPos[fItrk-1] << " fCentralityPercentile[fItrk]: " << fCentralityPercentile[fItrk-1] << " fMotherPt[fItrk]: " << fMotherPt[fItrk-1] << endl;
1534 }//end rotation loop
1541 // Post output data.
1542 PostData(1, fOutputContainer);
1543 PostData(2, fTreeV0);
1546 //________________________________________________________________________
1547 void AliAnalysisTaskLambdaNAOD::Terminate(const Option_t *)
1549 // Draw result to the screen
1550 // Called once at the end of the query
1552 //get output data and darw 'fHistPt'
1553 if (!GetOutputData(0)) return;
1554 //TH1F *hist=(TH1F*)(((TObjArray*)GetOutputData(0))->FindObject("fHistPt"));
1555 //if (hist) hist->Draw();
1558 //_____________________________________________________________________________
1559 Int_t AliAnalysisTaskLambdaNAOD::Initialize() {
1568 //________________________________________________________________________
1569 Int_t AliAnalysisTaskLambdaNAOD::SetupEvent() {
1570 // Setup Reading of event
1574 // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
1576 // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
1578 //Bool_t isTriggered = IsTriggered();
1582 //_____________________________________________________________________________
1583 void AliAnalysisTaskLambdaNAOD::ResetEvent() {
1588 for (Int_t ii = 0; ii < fNTriggers; ++ii)
1589 fTriggerFired[ii] = kFALSE;
1593 //_______________________________________________________________________
1594 void AliAnalysisTaskLambdaNAOD::BinLogAxis(const THnSparse *h, Int_t axisNumber) {
1596 // Method for the correct logarithmic binning of histograms
1598 TAxis *axis = h->GetAxis(axisNumber);
1599 int bins = axis->GetNbins();
1601 Double_t from = axis->GetXmin();
1602 Double_t to = axis->GetXmax();
1603 Double_t *newBins = new Double_t[bins + 1];
1606 Double_t factor = pow(to/from, 1./bins);
1608 for (int i = 1; i <= bins; i++) {
1609 newBins[i] = factor * newBins[i-1];
1611 axis->Set(bins, newBins);
1615 //________________________________________________________________________
1616 Bool_t AliAnalysisTaskLambdaNAOD::IsTriggered() {
1617 // Check if Event is triggered and fill Trigger Histogram
1619 if ((fEventHandler->IsEventSelected() & AliVEvent::kMB)) fTriggerFired[0] = kTRUE;
1620 if ((fEventHandler->IsEventSelected() & AliVEvent::kCentral)) fTriggerFired[1] = kTRUE;
1621 if ((fEventHandler->IsEventSelected() & AliVEvent::kSemiCentral)) fTriggerFired[2] = kTRUE;
1622 if ((fEventHandler->IsEventSelected() & AliVEvent::kEMCEJE)) fTriggerFired[3] = kTRUE;
1623 if ((fEventHandler->IsEventSelected() & AliVEvent::kEMCEGA)) fTriggerFired[4] = kTRUE;
1625 Bool_t isTriggered = kFALSE;
1627 for (Int_t ii=0; ii<fNTriggers; ++ii) {
1628 if(fTriggerFired[ii]) {
1629 isTriggered = kTRUE;
1630 fHistTriggerStat->Fill(ii);
1636 //______________________________________________________________________________
1637 Bool_t AliAnalysisTaskLambdaNAOD::DeuteronPID(AliVTrack *trackP, AliVTrack *trackN, Double_t ptotP, Double_t ptotN, Int_t runNumber, Bool_t isDeuteron[3]) {
1639 //define the arrays for the Bethe-Bloch-Parameters
1640 Double_t parDeuteron[5] = {0,0,0,0,0};
1642 if(runNumber < 166500) //LHC10h
1644 parDeuteron[0] = 1.45802; // ALEPH parameters for deuterons (pass2)
1645 parDeuteron[1] = 27.4992;
1646 parDeuteron[2] = 4.00313e-15;
1647 parDeuteron[3] = 2.48485;
1648 parDeuteron[4] = 8.31768;
1651 if(runNumber > 166500) //LHC11h
1653 parDeuteron[0] = 1.17; // ALEPH parameters for deuterons (pass2)
1654 parDeuteron[1] = 26.1144;
1655 parDeuteron[2] = 4.00313e-15;
1656 parDeuteron[3] = 2.72969 ;
1657 parDeuteron[4] = 9.15038;
1660 if(fMCtrue && runNumber > 166500){ //LHC11h MC (local production)
1661 parDeuteron[0] = 4.449; // ALEPH parameters for deuterons (pass2)
1662 parDeuteron[1] = 6.91865;
1663 parDeuteron[2] = 0.0183501;
1664 parDeuteron[3] = 2.49437;
1665 parDeuteron[4] = 2.62616;
1669 //define expected signals for the various species
1670 Double_t expSignalDeuteronN = 0;
1671 Double_t expSignalDeuteronP = 0;
1673 isDeuteron[0] = kFALSE;
1674 isDeuteron[1] = kFALSE;
1675 isDeuteron[2] = kFALSE;
1680 expSignalDeuteronN = AliExternalTrackParam::BetheBlochAleph(ptotN/(fgkMass[kMassDeuteron]),parDeuteron[0],parDeuteron[1],parDeuteron[2],parDeuteron[3],parDeuteron[4]);
1681 expSignalDeuteronP = AliExternalTrackParam::BetheBlochAleph(ptotP/(fgkMass[kMassDeuteron]),parDeuteron[0],parDeuteron[1],parDeuteron[2],parDeuteron[3],parDeuteron[4]);
1683 if(trackP->GetTPCsignal() >= 100 && //needed to reduce the size of the tree -- below the deuterons merge with the other particle spezies; standard value 110, variation for systematics
1684 trackP->GetTPCsignal() < 1200 && //needed to reduce the size of the tree -- anyway above the TPC does not work properly (?)
1685 (TMath::Abs(trackP->GetTPCsignal() - expSignalDeuteronP)/expSignalDeuteronP) < 0.32 && //for systemactics -- normal value 0.24
1688 if(trackP->Charge() >0) isDeuteron[0] = kTRUE; //pos deuteron
1689 if(trackP->Charge() <0) isDeuteron[1] = kTRUE; //neg deuteron
1692 if(trackN->GetTPCsignal() >= 100 && //needed to reduce the size of the tree -- below the deuterons merge with the other particle spezies
1693 trackN->GetTPCsignal() < 1200 && //needed to reduce the size of the tree -- anyway above the TPC does not work properly (?)
1694 (TMath::Abs(trackN->GetTPCsignal() - expSignalDeuteronN)/expSignalDeuteronN) < 0.32 && //for systemactics -- normal value 0.24
1697 isDeuteron[2] = kTRUE;
1699 if(trackN->Charge() >0) isDeuteron[0] = kTRUE; //pos deuteron
1700 if(trackN->Charge() <0) isDeuteron[1] = kTRUE; //neg deuteron
1707 if(runNumber < 166500) //2010
1709 expSignalDeuteronN = 0.65*AliExternalTrackParam::BetheBlochAleph(ptotN/(fgkMass[kMassDeuteron]),parDeuteron[0],parDeuteron[1],parDeuteron[2],parDeuteron[3],parDeuteron[4]);
1710 expSignalDeuteronP = 0.65*AliExternalTrackParam::BetheBlochAleph(ptotP/(fgkMass[kMassDeuteron]),parDeuteron[0],parDeuteron[1],parDeuteron[2],parDeuteron[3],parDeuteron[4]);
1712 if(runNumber > 166500) //2011
1714 expSignalDeuteronN = AliExternalTrackParam::BetheBlochAleph(ptotN/(fgkMass[kMassDeuteron]),parDeuteron[0],parDeuteron[1],parDeuteron[2],parDeuteron[3],parDeuteron[4]);
1715 expSignalDeuteronP = AliExternalTrackParam::BetheBlochAleph(ptotP/(fgkMass[kMassDeuteron]),parDeuteron[0],parDeuteron[1],parDeuteron[2],parDeuteron[3],parDeuteron[4]);
1718 if(trackP->GetTPCsignal() < 1200 &&
1719 (TMath::Abs(trackP->GetTPCsignal() - expSignalDeuteronP)/expSignalDeuteronP) < 0.2 &&
1722 if(trackP->Charge() >0) isDeuteron[0] = kTRUE; //pos deuteron
1723 if(trackP->Charge() <0) isDeuteron[1] = kTRUE; //neg deuteron
1726 if(trackN->GetTPCsignal() < 1200 &&
1727 (TMath::Abs(trackN->GetTPCsignal() - expSignalDeuteronN)/expSignalDeuteronN) < 0.2 &&
1730 isDeuteron[2] = kTRUE;
1732 if(trackN->Charge() >0) isDeuteron[0] = kTRUE; //pos deuteron
1733 if(trackN->Charge() <0) isDeuteron[1] = kTRUE; //neg deuteron
1743 //______________________________________________________________________________
1744 Bool_t AliAnalysisTaskLambdaNAOD::PionPID(AliVTrack *trackP, AliVTrack *trackN, Double_t ptotP, Double_t ptotN, Int_t runNumber, Bool_t isPion[2]) {
1746 //Pion PID via specific energy loss in the TPC
1747 //define the array for the Bethe-Bloch-Parameters (only necessary for local MC)
1748 Double_t parPion[5] = {0,0,0,0,0};
1750 if(runNumber < 166500) { //LHC10h
1751 parPion[0] = 1.45802; // ALEPH parameters for pions (pass2)
1752 parPion[1] = 27.4992;
1753 parPion[2] = 4.00313e-15;
1754 parPion[3] = 2.48485;
1755 parPion[4] = 8.31768;
1758 if(runNumber > 166500) { //LHC11h
1759 parPion[0] = 1.11243; // ALEPH parameters for pions (pass2)
1760 parPion[1] = 26.1144;
1761 parPion[2] = 4.00313e-15;
1762 parPion[3] = 2.72969;
1763 parPion[4] = 9.15038;
1766 if(fMCtrue && runNumber > 166500){ //LHC11h MC (local production)
1767 parPion[0] = 1.064; // ALEPH parameters for deuterons (pass2)
1769 parPion[2] = 4.00313e-15;
1770 parPion[3] = 2.703 ;
1774 Double_t expSignalPionP = 0;
1775 Double_t expSignalPionN = 0;
1779 if(runNumber < 166500){ //2010
1780 expSignalPionP = 0.7*AliExternalTrackParam::BetheBlochAleph(ptotP/(fgkMass[kMassPion]),parPion[0],parPion[1],parPion[2],parPion[3],parPion[4]);
1781 expSignalPionN = 0.7*AliExternalTrackParam::BetheBlochAleph(ptotN/(fgkMass[kMassPion]),parPion[0],parPion[1],parPion[2],parPion[3],parPion[4]);
1783 if(runNumber > 166500){ //2011
1784 expSignalPionP = 1.1*AliExternalTrackParam::BetheBlochAleph(ptotP/(fgkMass[kMassPion]),parPion[0],parPion[1],parPion[2],parPion[3],parPion[4]);
1785 expSignalPionN = 1.1*AliExternalTrackParam::BetheBlochAleph(ptotN/(fgkMass[kMassPion]),parPion[0],parPion[1],parPion[2],parPion[3],parPion[4]);
1794 if(TMath::Abs(fPIDResponse->NumberOfSigmasTPC(trackP, AliPID::kPion))<4){ //for systematics -- normal value 3
1796 if(trackP->Charge()>0 ) isPion[0] = kTRUE; //pos pion
1797 if(trackP->Charge()<0 ) isPion[1] = kTRUE; //neg pion
1799 if(TMath::Abs(fPIDResponse->NumberOfSigmasTPC(trackN, AliPID::kPion))<4){ //for systematics -- normal value 3
1801 if(trackN->Charge()>0 ) isPion[0] = kTRUE; //pos pion
1802 if(trackN->Charge()<0 ) isPion[1] = kTRUE; //neg pion
1808 if(TMath::Abs(trackP->GetTPCsignal() - expSignalPionP)/expSignalPionP < 0.2
1810 if(trackP->Charge()>0) isPion[0] = kTRUE; //pos pion
1811 if(trackP->Charge()<0) isPion[1] = kTRUE; //neg pion
1813 if(TMath::Abs(trackN->GetTPCsignal() - expSignalPionN)/expSignalPionN < 0.2
1815 if(trackN->Charge()>0) isPion[0] = kTRUE; //pos pion
1816 if(trackN->Charge()<0) isPion[1] = kTRUE; //neg pion
1823 //______________________________________________________________________________
1824 Bool_t AliAnalysisTaskLambdaNAOD::TrackCuts(AliVTrack *track, Bool_t testTrackCuts) {
1826 testTrackCuts = kFALSE;
1828 AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
1829 //if(!esdtrack) testTrackCuts = kFALSE;
1830 if (fESDtrackCutsV0->AcceptTrack(esdtrack)) testTrackCuts = kTRUE;
1831 //if( testTrackCuts == kTRUE) cout << "testTrackCuts im test: " << testTrackCuts << endl;
1834 return testTrackCuts;
1836 //______________________________________________________________________________
1837 /*Bool_t AliAnalysisTaskLambdaNAOD::FilterBit(AliVTrack *track, Bool_t testFilterBit) {
1839 testFilterBit = kFALSE;
1841 AliAODTrack *aodtrack = static_cast<AliAODTrack *>(track);
1842 //if(!aodtrack) testFilterBit = kFALSE;
1843 if (aodtrack->TestFilterBit(7)) testFilterBit = kTRUE;
1845 //if(testFilterBit == kTRUE) cout << "testFilterBit im test: " << testFilterBit<< endl;
1847 return testFilterBit;
1849 //______________________________________________________________________
1850 void AliAnalysisTaskLambdaNAOD::MCGenerated(AliStack* stack, Int_t multiplicity)
1853 // Monte Carlo for genenerated particles
1854 if (fMCtrue) //MC loop
1859 for(stackN = 0; stackN < stack->GetNtrack(); stackN++) //loop over stack
1862 const TParticle *tparticleMother = stack->Particle(stackN);
1863 Long_t pdgCodeMother = tparticleMother->GetPdgCode();
1865 //check which particle the mother is
1868 if(pdgCodeMother == fgkPdgCode[kPDGLambdaNeutron]) //check mother PDG
1870 MCTwoBodyDecay(stack,tparticleMother,pdgCodeMother,fgkPdgCode[kPDGDeuteron],fgkPdgCode[kPDGPionMinus],fgkMass[kMassDeuteron],fgkMass[kMassPion],multiplicity);
1873 //Anti-LambdaNeutron
1874 if(pdgCodeMother == fgkPdgCode[kPDGAntiLambdaNeutron]) //check mother PDG
1876 MCTwoBodyDecay(stack,tparticleMother,pdgCodeMother,fgkPdgCode[kPDGAntiDeuteron],fgkPdgCode[kPDGPionPlus],fgkMass[kMassDeuteron],fgkMass[kMassPion],multiplicity);
1880 }//end loop over stack
1885 //_____________________________________________
1886 void AliAnalysisTaskLambdaNAOD::MCTwoBodyDecay(AliStack* stack, const TParticle *tparticleMother, Long_t PDGMother, Long_t PDGFirstDaughter, Long_t PDGSecondDaughter, Double_t massFirstDaughter, Double_t massSecondDaughter, Int_t multiplicity){ //function that calculates the invariant mass and the transverse momentum for MC
1888 Double_t momentumFirstDaughterGen[3]={0,0,0};
1889 Double_t momentumSecondDaughterGen[3]={0,0,0};
1891 Double_t energyFirstDaughterGen = 0;
1892 Double_t energySecondDaughterGen = 0;
1894 Double_t transversMomentumMotherGen = 0;
1895 Double_t longitudinalMomentumMotherGen = 0;
1896 Double_t totalEnergyMotherGen = 0;
1898 Double_t rapidityGen = 2;
1900 //Int_t labelFirstDaughter = tparticleMother->GetDaughter(0);
1901 Int_t labelSecondDaughter = tparticleMother->GetDaughter(1);
1902 Int_t labelFirstDaughter = labelSecondDaughter -1;
1904 TParticle *tparticleFirstDaughter = stack->Particle(TMath::Abs(labelFirstDaughter));
1905 TParticle *tparticleSecondDaughter = stack->Particle(TMath::Abs(labelSecondDaughter));
1907 if(tparticleFirstDaughter->GetPdgCode() == PDGFirstDaughter) //check first daughter PDG
1909 if(tparticleSecondDaughter->GetPdgCode() == PDGSecondDaughter) //check second daughter PDG
1912 momentumFirstDaughterGen[0] = tparticleFirstDaughter->Px();
1913 momentumFirstDaughterGen[1] = tparticleFirstDaughter->Py();
1914 momentumFirstDaughterGen[2] = tparticleFirstDaughter->Pz();
1916 momentumSecondDaughterGen[0] = tparticleSecondDaughter->Px();
1917 momentumSecondDaughterGen[1] = tparticleSecondDaughter->Py();
1918 momentumSecondDaughterGen[2] = tparticleSecondDaughter->Pz();
1920 energyFirstDaughterGen = tparticleFirstDaughter->Energy();
1921 energySecondDaughterGen = tparticleSecondDaughter->Energy();
1923 transversMomentumMotherGen=TMath::Sqrt((momentumFirstDaughterGen[0]+momentumSecondDaughterGen[0])*(momentumFirstDaughterGen[0]+momentumSecondDaughterGen[0])+(momentumFirstDaughterGen[1]+momentumSecondDaughterGen[1])*(momentumFirstDaughterGen[1]+momentumSecondDaughterGen[1]));
1926 //longitudinal momentum of mother
1927 longitudinalMomentumMotherGen = TMath::Sqrt((momentumFirstDaughterGen[2]+momentumSecondDaughterGen[2])*(momentumFirstDaughterGen[2]+momentumSecondDaughterGen[2]));
1929 //total energy of mother
1930 totalEnergyMotherGen = energyFirstDaughterGen + energySecondDaughterGen;
1933 rapidityGen = 0.5 * TMath::Log( (totalEnergyMotherGen+longitudinalMomentumMotherGen)/(totalEnergyMotherGen-longitudinalMomentumMotherGen));
1935 if(rapidityGen > 1.0 || rapidityGen < -1 ) return;
1937 //cout << "decay lenght first daughter :" << tparticleFirstDaughter->Rho() << " decay lenght second daughter :" << tparticleSecondDaughter->Rho()<< endl;
1938 //fHistLambdaNeutronDecayLengthGen->Fill(tparticleFirstDaughter->Rho());
1940 //calculate the invariant mass
1941 Double_t firstDaughterPart = 0;
1942 Double_t secondDaughterPart = 0;
1943 Double_t invaMass = 0;
1945 firstDaughterPart = massFirstDaughter*massFirstDaughter
1946 +momentumFirstDaughterGen[0]*momentumFirstDaughterGen[0]
1947 +momentumFirstDaughterGen[1]*momentumFirstDaughterGen[1]
1948 +momentumFirstDaughterGen[2]*momentumFirstDaughterGen[2];
1949 secondDaughterPart = massSecondDaughter*massSecondDaughter
1950 +momentumSecondDaughterGen[0]*momentumSecondDaughterGen[0]
1951 +momentumSecondDaughterGen[1]*momentumSecondDaughterGen[1]
1952 +momentumSecondDaughterGen[2]*momentumSecondDaughterGen[2];
1954 invaMass = TMath::Sqrt(TMath::Max(massFirstDaughter*massFirstDaughter
1955 +massSecondDaughter*massSecondDaughter
1956 + 2.*(TMath::Sqrt(firstDaughterPart*secondDaughterPart)
1957 -momentumFirstDaughterGen[0]*momentumSecondDaughterGen[0]
1958 -momentumFirstDaughterGen[1]*momentumSecondDaughterGen[1]
1959 -momentumFirstDaughterGen[2]*momentumSecondDaughterGen[2]), 0.));
1964 case fgkPdgCode[kPDGLambdaNeutron] :
1965 fHistLambdaNeutronPtGen->Fill(transversMomentumMotherGen);
1966 if(multiplicity < 2500)fHistLambdaNeutronPtGenMinBias->Fill(transversMomentumMotherGen);
1967 if(multiplicity > 1500 && multiplicity < 2750)fHistLambdaNeutronPtGenCentral->Fill(transversMomentumMotherGen);
1968 if(multiplicity > 300 && multiplicity < 2000)fHistLambdaNeutronPtGenSemiCentral->Fill(transversMomentumMotherGen);
1969 fHistLambdaNeutronInvaMassGen->Fill(invaMass);
1971 case fgkPdgCode[kPDGAntiLambdaNeutron] :
1972 fHistAntiLambdaNeutronPtGen->Fill(transversMomentumMotherGen);
1973 if(multiplicity < 2500)fHistAntiLambdaNeutronPtGenMinBias->Fill(transversMomentumMotherGen);
1974 if(multiplicity > 1500 && multiplicity < 2750)fHistAntiLambdaNeutronPtGenCentral->Fill(transversMomentumMotherGen);
1975 if(multiplicity > 300 && multiplicity < 2000)fHistAntiLambdaNeutronPtGenSemiCentral->Fill(transversMomentumMotherGen);
1976 fHistAntiLambdaNeutronInvaMassGen->Fill(invaMass);
1979 printf("should not happen!!!! \n");
1985 if(PDGMother == fgkPdgCode[kPDGLambdaNeutron])
1987 fHistLambdaNeutronPtGen->Fill(transversMomentumMotherGen);
1988 if(multiplicity < 2500)fHistLambdaNeutronPtGenMinBias->Fill(transversMomentumMotherGen);
1989 if(multiplicity > 1500 && multiplicity < 2750)fHistLambdaNeutronPtGenCentral->Fill(transversMomentumMotherGen);
1990 if(multiplicity > 300 && multiplicity < 2000)fHistLambdaNeutronPtGenSemiCentral->Fill(transversMomentumMotherGen);
1991 fHistLambdaNeutronInvaMassGen->Fill(invaMass);
1992 fHistLambdaNeutronDecayLengthGen->Fill(((tparticleFirstDaughter->Rho())*2.054)/(tparticleMother->P()));
1994 //Anti-LambdaNeutron
1995 if(PDGMother == fgkPdgCode[kPDGAntiLambdaNeutron])
1997 fHistAntiLambdaNeutronPtGen->Fill(transversMomentumMotherGen);
1998 if(multiplicity < 2500)fHistAntiLambdaNeutronPtGenMinBias->Fill(transversMomentumMotherGen);
1999 if(multiplicity > 1500 && multiplicity < 2750)fHistAntiLambdaNeutronPtGenCentral->Fill(transversMomentumMotherGen);
2000 if(multiplicity > 300 && multiplicity < 2000)fHistAntiLambdaNeutronPtGenSemiCentral->Fill(transversMomentumMotherGen);
2001 fHistAntiLambdaNeutronInvaMassGen->Fill(invaMass);
2002 fHistAntiLambdaNeutronDecayLengthGen->Fill(((tparticleFirstDaughter->Rho())*2.054)/(tparticleMother->P()));
2006 }//end of check second daughter PDG
2007 }//end of check first daughter PDG
2010 //_____________________________________________
2011 Double_t AliAnalysisTaskLambdaNAOD::MomentumInnerParam(AliVTrack *track, Double_t ptot){ //function to get the momentum at the inner wall of the TPC
2013 AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
2014 ptot = esdtrack->GetInnerParam()->GetP();
2018 //_____________________________________________
2019 UShort_t AliAnalysisTaskLambdaNAOD::TPCclusters(AliVTrack *track, UShort_t numberOfTPCclusters){ //function to get the number of clusters used for each track
2021 AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
2022 numberOfTPCclusters = esdtrack->GetTPCNcls();
2024 return numberOfTPCclusters;
2026 //_____________________________________________
2027 Double_t AliAnalysisTaskLambdaNAOD::TPCchi2(AliVTrack *track, Double_t numberOfChi2clustersTPC, UShort_t numberOfTPCclusters){ //function to get the chi2 per clusters used for each track
2029 AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
2030 numberOfChi2clustersTPC = esdtrack->GetTPCchi2()/numberOfTPCclusters;
2032 return numberOfChi2clustersTPC;
2034 //_____________________________________________
2035 Double_t AliAnalysisTaskLambdaNAOD::ImpactParameter(AliVTrack *track, Double_t dcaToVertex){ //function to get the number of clusters used for each track
2037 Float_t tdcaToVertex[2] = {-1,-1};
2039 AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
2040 esdtrack->GetImpactParameters(tdcaToVertex[0],tdcaToVertex[1]);
2041 dcaToVertex = TMath::Sqrt(tdcaToVertex[0]*tdcaToVertex[0]+tdcaToVertex[1]*tdcaToVertex[1]);
2045 //________________________________________________________________________
2046 /*void AliAnalysisTaskLambdaNAOD::RotateKFParticle(AliKFParticle * kfParticle,Double_t angle, const AliVEvent * const ev){
2048 // Before rotate needs to be moved to position 0,0,0, ; move back after rotation
2049 if (!kfParticle) return;
2055 dx = ev->GetPrimaryVertex()->GetX()-0.;
2056 dy = ev->GetPrimaryVertex()->GetY()-0.;
2057 dz = ev->GetPrimaryVertex()->GetZ()-0.;
2060 kfParticle->X() = kfParticle->GetX() - dx;
2061 kfParticle->Y() = kfParticle->GetY() - dy;
2062 kfParticle->Z() = kfParticle->GetZ() - dz;
2065 // Rotate the kf particle
2066 Double_t c = cos(angle);
2067 Double_t s = sin(angle);
2070 for( Int_t i=0; i<8; i++ ){
2071 for( Int_t j=0; j<8; j++){
2075 for( int i=0; i<8; i++ ){
2078 mA[0][0] = c; mA[0][1] = s;
2079 mA[1][0] = -s; mA[1][1] = c;
2080 mA[3][3] = c; mA[3][4] = s;
2081 mA[4][3] = -s; mA[4][4] = c;
2086 for( Int_t i=0; i<8; i++ ){
2088 for( Int_t k=0; k<8; k++){
2089 mAp[i]+=mA[i][k] * kfParticle->GetParameter(k);
2093 for( Int_t i=0; i<8; i++){
2094 kfParticle->Parameter(i) = mAp[i];
2097 for( Int_t i=0; i<8; i++ ){
2098 for( Int_t j=0; j<8; j++ ){
2100 for( Int_t k=0; k<8; k++ ){
2101 mAC[i][j]+= mA[i][k] * kfParticle->GetCovariance(k,j);
2106 for( Int_t i=0; i<8; i++ ){
2107 for( Int_t j=0; j<=i; j++ ){
2109 for( Int_t k=0; k<8; k++){
2110 xx+= mAC[i][k]*mA[j][k];
2112 kfParticle->Covariance(i,j) = xx;
2116 kfParticle->X() = kfParticle->GetX() + dx;
2117 kfParticle->Y() = kfParticle->GetY() + dy;
2118 kfParticle->Z() = kfParticle->GetZ() + dz;