]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGLF/STRANGENESS/Hypernuclei/AliAnalysisTaskHelium3PiMC.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
.so cleanup: more gSystem->Load()
[u/mrichter/AliRoot.git] / PWGLF / STRANGENESS / Hypernuclei / AliAnalysisTaskHelium3PiMC.cxx
CommitLineData
c683985a 1/**************************************************************************
2 * Contributors are not mentioned at all. *
3 * *
4 * Permission to use, copy, modify and distribute this software and its *
5 * documentation strictly for non-commercial purposes is hereby granted *
6 * without fee, provided that the above copyright notice appears in all *
7 * copies and that both the copyright notice and this permission notice *
8 * appear in the supporting documentation. The authors make no claims *
9 * about the suitability of this software for any purpose. It is *
10 * provided "as is" without express or implied warranty. *
11 **************************************************************************/
12//
13//-----------------------------------------------------------------
14// AliAnalysisTaskHelium3Pi class
15//-----------------------------------------------------------------
16
17
18class TTree;
19class TParticle;
20class TVector3;
21
22#include "AliAnalysisManager.h"
23#include <AliMCEventHandler.h>
24#include <AliMCEvent.h>
25#include <AliStack.h>
26
27class AliESDVertex;
28class AliAODVertex;
29class AliESDv0;
30class AliAODv0;
31class AliCascadeVertexer;
32
33#include <iostream>
34#include "AliAnalysisTaskSE.h"
35#include "TList.h"
36#include "TH1.h"
37#include "TH2.h"
38#include "TH3.h"
39#include "TNtuple.h"
40#include "TGraph.h"
41#include "TCutG.h"
42#include "TF1.h"
43#include "TCanvas.h"
44#include "TMath.h"
45#include "TChain.h"
46#include "Riostream.h"
47#include "AliLog.h"
48#include "AliCascadeVertexer.h"
49#include "AliESDEvent.h"
50#include "AliESDtrack.h"
51#include "AliExternalTrackParam.h"
52#include "AliAODEvent.h"
53#include "AliInputEventHandler.h"
54#include "AliESDcascade.h"
55#include "AliAODcascade.h"
56#include "AliAnalysisTaskHelium3PiMC.h"
57#include "AliESDtrackCuts.h"
58#include "AliCentrality.h"
59#include "TString.h"
60#include <TDatime.h>
61#include <TRandom3.h>
62using std::endl;
63using std::cout;
64
65const Int_t AliAnalysisTaskHelium3PiMC::fgNrot = 15;
66
67
68ClassImp(AliAnalysisTaskHelium3PiMC)
69
70//________________________________________________________________________
71AliAnalysisTaskHelium3PiMC::AliAnalysisTaskHelium3PiMC()
72: AliAnalysisTaskSE(),
73 fAnalysisType("ESD"),
74 fCollidingSystems(0),
75 fDataType("SIM"),
76 fListHistCascade(0),
77 fHistEventMultiplicity(0),
78 fHistTrackMultiplicity(0),
79 fHistMCMultiplicityTracks(0),
80 fHistMCEta(0),
81 fHistMCPt(0),
82 fHistMCTheta(0),
83 fHistMCDecayPosition(0),
84 fHistMCDecayRho(0),
85 fhRigidityHevsMomPiMC(0),
86 fhRigidityHevsMomPiRec(0),
87 fhInvMassMC(0),
88 fhInvMassMum(0),
89 fhInvMassRec(0),
90 fhInvMassRec1(0),
91 fhInvMassRec2(0),
92 fhInvMassRec3(0),
93 fhInvMassRec4(0),
94 fhInvMassRec5(0),
95 fhInvMassRec6(0),
96 fhInvMassRec7(0),
97 fhHeMCRigidity(0),
98 fhPioneMC(0),
99 hBBTPCnoCuts(0),
100 fhBBTPC(0),
101 fhBBTPCNegativePions(0),
102 fhBBTPCPositivePions(0),
103 fhBBTPCHe3(0),
104 fHistProvaDCA(0),
105 fHistPercentileVsTrackNumber(0),
106 fhHeDCAXY(0),
107 fhHeDCAZ(0),
108 fhPiDCAXY(0),
109 fhPiDCAZ(0),
110 hITSClusterMap(0),
111 fNtuple1(0),
112 fNtuple2(0)
113
114{
115 // Dummy Constructor(0);
116}
117
118//________________________________________________________________________
119AliAnalysisTaskHelium3PiMC::AliAnalysisTaskHelium3PiMC(const char *name)
120 : AliAnalysisTaskSE(name),
121 fAnalysisType("ESD"),
122 fCollidingSystems(0),
123 fDataType("SIM"),
124 fListHistCascade(0),
125 fHistEventMultiplicity(0),
126 fHistTrackMultiplicity(0),
127 fHistMCMultiplicityTracks(0),
128 fHistMCEta(0),
129 fHistMCPt(0),
130 fHistMCTheta(0),
131 fHistMCDecayPosition(0),
132 fHistMCDecayRho(0),
133 fhRigidityHevsMomPiMC(0),
134 fhRigidityHevsMomPiRec(0),
135 fhInvMassMC(0),
136 fhInvMassMum(0),
137 fhInvMassRec(0),
138 fhInvMassRec1(0),
139 fhInvMassRec2(0),
140 fhInvMassRec3(0),
141 fhInvMassRec4(0),
142 fhInvMassRec5(0),
143 fhInvMassRec6(0),
144 fhInvMassRec7(0),
145 fhHeMCRigidity(0),
146 fhPioneMC(0),
147 hBBTPCnoCuts(0),
148 fhBBTPC(0),
149 fhBBTPCNegativePions(0),
150 fhBBTPCPositivePions(0),
151 fhBBTPCHe3(0),
152 fHistProvaDCA(0),
153 fHistPercentileVsTrackNumber(0),
154 fhHeDCAXY(0),
155 fhHeDCAZ(0),
156 fhPiDCAXY(0),
157 fhPiDCAZ(0),
158 hITSClusterMap(0),
159 fNtuple1(0),
160 fNtuple2(0)
161
162
163{
164 // Define input and output slots here
165 // Input slot #0 works with a TChain
166 //DefineInput(0, TChain::Class());
167 // Output slot #0 writes into a TList container (Cascade)
168 DefineOutput(1, TList::Class());
169}
170//_______________________________________________________
171AliAnalysisTaskHelium3PiMC::~AliAnalysisTaskHelium3PiMC()
172{
173 // Destructor
174 if (fListHistCascade) {
175 delete fListHistCascade;
176 fListHistCascade = 0;
177 }
178
179}
180//=================DEFINITION BETHE BLOCH==============================
181
182Double_t AliAnalysisTaskHelium3PiMC::BetheBloch(Double_t betaGamma,Double_t charge,Bool_t isPbPb) {
183
184 Double_t kp1, kp2, kp3, kp4, kp5;
185
186 if(isPbPb){
187
188 //pass2 //to be checked
189 kp1 = 5.2*charge*charge;
190 kp2 = 8.98482806165147636e+00;
191 kp3 = 1.54000000000000005e-05;
192 kp4 = 2.30445734159456084e+00;
193 kp5 = 2.25624744086878559e+00;
194
195
196 }
197
198 else{
199
200 //pass2 // to be defined
201 kp1 = 5.2*charge*charge;
202 kp2 = 8.98482806165147636e+00;
203 kp3 = 1.54000000000000005e-05;
204 kp4 = 2.30445734159456084e+00;
205 kp5 = 2.25624744086878559e+00;
206
207 }
208
209 Double_t beta = betaGamma / TMath::Sqrt(1.0 + betaGamma * betaGamma);
210
211 Double_t aa = TMath::Power(beta, kp4);
212 Double_t bb = TMath::Power(1.0 / betaGamma, kp5);
213
214 bb = TMath::Log(kp3 + bb);
215
216 Double_t out = (kp2 - aa - bb) * kp1 / aa;
217
218 return out;
219
220}
221
222//==================DEFINITION OF OUTPUT OBJECTS==============================
223
224void AliAnalysisTaskHelium3PiMC::UserCreateOutputObjects()
225{
226 fListHistCascade = new TList();
227 fListHistCascade->SetOwner(); // IMPORTANT!
228
229 if(! fHistEventMultiplicity ){
230 fHistEventMultiplicity = new TH1F( "fHistEventMultiplicity" , "Nb of Events" , 6 , -1, 5 );
231 fHistEventMultiplicity->GetXaxis()->SetTitle("Event Type");
232 fListHistCascade->Add(fHistEventMultiplicity);
233 }
234
235 if(! fHistTrackMultiplicity ){
236
237 fHistTrackMultiplicity = new TH1F( "fHistTrackMultiplicity" , "Nb of Tracks" , 25000,0, 25000 );
238 fHistTrackMultiplicity->GetXaxis()->SetTitle("Number of tracks");
239 fListHistCascade->Add(fHistTrackMultiplicity);
240 }
241
242 if(! fHistMCMultiplicityTracks){
243 fHistMCMultiplicityTracks =new TH1F("fHistMCMultiplicityTracks","MC Multiplicity Tracks",1000,0,1000);
244 fHistMCMultiplicityTracks->GetXaxis()->SetTitle("MC Number of tracks");
245 fListHistCascade->Add(fHistMCMultiplicityTracks);
246 }
247 if(!fHistMCEta ){
248 fHistMCEta=new TH1F("fHistMCEta","MC eta",1000,-3,3);
249 fHistMCEta->GetXaxis()->SetTitle("Injected Eta");
250 fListHistCascade->Add(fHistMCEta);
251 }
252 if(! fHistMCPt){
253 fHistMCPt =new TH1F("fHistMCPt","MC pt",1000,0,20);
254 fHistMCPt->GetXaxis()->SetTitle("Injected Pt");
255 fListHistCascade->Add(fHistMCPt);
256 }
257 if(!fHistMCTheta ){
258 fHistMCTheta=new TH1F("fHistMCTheta","MC theta",1000,-6,6);
259 fHistMCTheta->GetXaxis()->SetTitle("Injected Theta");
260 fListHistCascade->Add(fHistMCTheta);
261 }
262 if(!fHistMCDecayPosition){
263 fHistMCDecayPosition =new TH1F("fHistMCDecayPosition","MC Decay Position",10000,0,1000);
264 fHistMCDecayPosition->GetXaxis()->SetTitle("Decay Position");
265 fListHistCascade->Add(fHistMCDecayPosition);
266 }
267 if(!fHistMCDecayRho ){
268 fHistMCDecayRho =new TH1F("fHistMCDecayRho","MC decay position 3d",10000,0,1000);
269 fHistMCDecayRho->GetXaxis()->SetTitle("Decay rho");
270 fListHistCascade->Add(fHistMCDecayRho);
271 }
272
273 if(!fhRigidityHevsMomPiMC ){
274 fhRigidityHevsMomPiMC=new TH2F("fhRigidityHevsMomPiMC","Rigidity He vs Mom Pi MC",20,0,10,300,0,30);
275 fhRigidityHevsMomPiMC->GetXaxis()->SetTitle("He3 Rigidity");
276 fhRigidityHevsMomPiMC->GetYaxis()->SetTitle("Pi momentum");
277 fListHistCascade->Add(fhRigidityHevsMomPiMC);
278 }
279
280 if(! fhRigidityHevsMomPiRec){
281 fhRigidityHevsMomPiRec=new TH2F("fhRigidityHevsMomPiRec","Rigidity He vs Mom Pi Rec",20,0,10,300,0,30);
282 fhRigidityHevsMomPiRec->GetXaxis()->SetTitle("He3 Rigidity");
283 fhRigidityHevsMomPiRec->GetYaxis()->SetTitle("Pi momentum");
284 fListHistCascade->Add(fhRigidityHevsMomPiRec);
285 }
286
287 if(!fhInvMassMC){
288 fhInvMassMC=new TH1F("fhInvMassMC","fhInvMassMC",800,2.,6.);
289 fhInvMassMC->GetXaxis()->SetTitle("(He3,#pi) InvMass");
290 fListHistCascade->Add(fhInvMassMC);
291 }
292
293 if(!fhInvMassMum){
294 fhInvMassMum=new TH1F("fhInvMassMum","fhInvMassMum",800,2.,6.);
295 fhInvMassMum->GetXaxis()->SetTitle("(He3,#pi) InvMass");
296 fListHistCascade->Add(fhInvMassMum);
297 }
298
299 if(!fhInvMassRec){
300 fhInvMassRec=new TH1F("fhInvMassRec","fhInvMassRec",800,2.,6.);
301 fhInvMassRec->GetXaxis()->SetTitle("(He3,#pi) InvMass");
302 fListHistCascade->Add(fhInvMassRec);
303 }
304
305 if(!fhInvMassRec1){
306 fhInvMassRec1=new TH1F("fhInvMassRec1","No Altri tagli",800,2.,6.);
307 fhInvMassRec1->GetXaxis()->SetTitle("(He3,#pi) InvMass");
308 fListHistCascade->Add(fhInvMassRec1);
309 }
310 if(!fhInvMassRec2){
311 fhInvMassRec2=new TH1F("fhInvMassRec2","DCA pi > 0.1",800,2.,6.);
312 fhInvMassRec2->GetXaxis()->SetTitle("(He3,#pi) InvMass");
313 fListHistCascade->Add(fhInvMassRec2);
314 }
315 if(!fhInvMassRec3){
316 fhInvMassRec3=new TH1F("fhInvMassRec3","DCA He > 0.05",800,2.,6.);
317 fhInvMassRec3->GetXaxis()->SetTitle("(He3,#pi) InvMass");
318 fListHistCascade->Add(fhInvMassRec3);
319 }
320 if(!fhInvMassRec4){
321 fhInvMassRec4=new TH1F("fhInvMassRec4","DCA tracks < 1 cm",800,2.,6.);
322 fhInvMassRec4->GetXaxis()->SetTitle("(He3,#pi) InvMass");
323 fListHistCascade->Add(fhInvMassRec4);
324 }
325 if(!fhInvMassRec5){
326 fhInvMassRec5=new TH1F("fhInvMassRec5","Condizione xn+xp",800,2.,6.);
327 fhInvMassRec5->GetXaxis()->SetTitle("(He3,#pi) InvMass");
328 fListHistCascade->Add(fhInvMassRec5);
329 }
330 if(!fhInvMassRec6){
331 fhInvMassRec6=new TH1F("fhInvMassRec6","Ho fatto V0 ",800,2.,6.);
332 fhInvMassRec6->GetXaxis()->SetTitle("(He3,#pi) InvMass");
333 fListHistCascade->Add(fhInvMassRec6);
334 }
335 if(!fhInvMassRec7){
336 fhInvMassRec7=new TH1F("fhInvMassRec7","V0+Taglio CPA",800,2.,6.);
337 fhInvMassRec7->GetXaxis()->SetTitle("(He3,#pi) InvMass");
338 fListHistCascade->Add(fhInvMassRec7);
339 }
340
341 if(!fhHeMCRigidity ){
342 fhHeMCRigidity=new TH1F("fhHeMCRigidity","He3 rigidity distribution",200,0,20);
343 fhHeMCRigidity->GetXaxis()->SetTitle("He3 rigidity");
344 fListHistCascade->Add(fhHeMCRigidity);
345 }
346 if(!fhPioneMC ){
347 fhPioneMC=new TH1F("hPioneMC","Pion mom distribution",200,0,50);
348 fhPioneMC->GetXaxis()->SetTitle("Pion momentum");
349 fListHistCascade->Add(fhPioneMC);
350 }
351
352 if(!hBBTPCnoCuts ){
353 hBBTPCnoCuts=new TH2F("hBBTPCnoCuts","scatterPlot TPC no cuts",2000,-10,10,1000,0,3000);
354 hBBTPCnoCuts->GetXaxis()->SetTitle("p/Z (GeV/#it{c})");
355 hBBTPCnoCuts->GetYaxis()->SetTitle("TPC Signal (a.u)");
356 fListHistCascade->Add(hBBTPCnoCuts);
357 }
358 if(!fhBBTPC ){
359 fhBBTPC=new TH2F("fhBBTPC","scatterPlot TPC",2000,-10,10,1000,0,3000);
360 fhBBTPC->GetXaxis()->SetTitle("p/Z (GeV/#it{c})");
361 fhBBTPC->GetYaxis()->SetTitle("TPC Signal (a.u)");
362 fListHistCascade->Add(fhBBTPC);
363 }
364 if(!fhBBTPCNegativePions ){
365 fhBBTPCNegativePions=new TH2F("fhBBTPCNegativePions","scatterPlot Neg Pions",2000,-10,10,1000,0,3000);
366 fhBBTPCNegativePions->GetXaxis()->SetTitle("p/Z (GeV/#it{c})");
367 fhBBTPCNegativePions->GetYaxis()->SetTitle("TPC Signal (a.u)");
368 fListHistCascade->Add(fhBBTPCNegativePions);
369 }
370 if(!fhBBTPCPositivePions ){
371 fhBBTPCPositivePions=new TH2F("fhBBTPCPositivePions","scatterPlot Pos Pions",2000,-10,10,1000,0,3000);
372 fhBBTPCPositivePions->GetXaxis()->SetTitle("p/Z (GeV/#it{c})");
373 fhBBTPCPositivePions->GetYaxis()->SetTitle("TPC Signal (a.u)");
374 fListHistCascade->Add(fhBBTPCPositivePions);
375 }
376 if(!fhBBTPCHe3 ){
377 fhBBTPCHe3=new TH2F("fhBBTPCHe3","scatterPlot TPC - He3",2000,-10,10,1000,0,3000);
378 fhBBTPCHe3->GetXaxis()->SetTitle("p/Z (GeV/#it{c})");
379 fhBBTPCHe3->GetYaxis()->SetTitle("TPC Signal (a.u)");
380 fListHistCascade->Add(fhBBTPCHe3);
381 }
382 if(!fHistProvaDCA ){
383 fHistProvaDCA=new TH2F("fHistProvaDCA","fHistProvaDCA",1000,-50,50,1000,0,100);
384 fHistProvaDCA->GetXaxis()->SetTitle("xn+xp");
385 fHistProvaDCA->GetYaxis()->SetTitle("dca tracks");
386 fListHistCascade->Add(fHistProvaDCA);
387 }
388
389 if(!hITSClusterMap ){
390 hITSClusterMap=new TH1F("hITSClusterMap","hITSClusterMap",65,-1,64);
391 fListHistCascade->Add(hITSClusterMap);
392 }
393
394 if(!fHistPercentileVsTrackNumber){
395 fHistPercentileVsTrackNumber=new TH2F("fHistPercentileVsTrackNumber","fHistPercentileVsTrackNumber",120,-3,117,2500,0,25000);
396 fHistPercentileVsTrackNumber->GetXaxis()->SetTitle("Percentile");
397 fHistPercentileVsTrackNumber->GetYaxis()->SetTitle("Tracks Number");
398 fListHistCascade->Add(fHistPercentileVsTrackNumber);
399 }
400
401 if(!fhHeDCAXY){
402 fhHeDCAXY=new TH1F("fhHeDCAXY","fhHeDCAXY",800,-4,4);
403 fListHistCascade->Add(fhHeDCAXY);
404 }
405 if(!fhHeDCAZ){
406 fhHeDCAZ=new TH1F("fhHeDCAZ","fhHeDCAZ",800,-30,30);
407 fListHistCascade->Add(fhHeDCAZ);
408 }
409 if(!fhPiDCAXY){
410 fhPiDCAXY=new TH1F("fhPiDCAXY","fhPiDCAXY",800,-4,4);
411 fListHistCascade->Add(fhPiDCAXY);
412 }
413 if(!fhPiDCAZ){
414 fhPiDCAZ=new TH1F("fhPiDCAZ","fhPiDCAZ",800,-30,30);
415 fListHistCascade->Add(fhPiDCAZ);
416 }
417
418 if(! fNtuple1 ) {
419 fNtuple1 = new TNtuple("fNtuple1","Ntuple1","runNumber:evNumber:TrackNumber:percentile:xPrimaryVertex:yPrimaryVertex:zPrimaryVertex:xSecondaryVertex:ySecondaryVertex:zSecondaryVertex:dcaTracks:CosPointingAngle:DCAV0toPrimaryVertex:HeSign:HepInTPC:HeTPCsignal:DcaHeToPrimVertex:HeEta:momHex:momHey:momHez:momHeAtSVx:momHeAtSVy:momHeAtSVz:HeTPCNcls:HeimpactXY:HeimpactZ:isTOFHe:HeBeta:HeITSClusterMap:IsHeITSRefit:PionSign:PionpInTPC:PionTPCsignal:DcaPionToPrimVertex:PionEta:momPionx:momPiony:momPionz:momNegPionAtSVx:momNegPionAtSVy:momNegPionAtSVz:PionTPCNcls:PionimpactXY:PionimpactZ:isTOFPion:PionBeta:PionITSClusterMap:IsPiITSRefit:PDGCodeNeg:PDCCodePos:motherPDGNeg:motherPDGPos:labelPi:labelHe:mumidNeg:mumidPos");
420
421 fListHistCascade->Add(fNtuple1);
422 }
423
424 if(! fNtuple2 ) {
425
426 fNtuple2 = new TNtuple("fNtuple2","Ntuple2","run:event:iMC:Centrality:PVx:PVy:PVz:PDGcodeMum:MotherIndex:SVxD0:SVyD0:SVzD0:SVxD1:SVyD1:SVzD1:SV3d:EtaMum:YMum:ThetaMum:PhiMum:PxMum:PyMum:PzMum:PdgDaughter0:PdgDaughter1:PxD0:PyD0:PzD0:PxD1:PyD1:PzD1");
427
428 fListHistCascade->Add(fNtuple2);
429 }
430
431 PostData(1,fListHistCascade);
432
433}// end UserCreateOutputObjects
434
435
436
437//====================== USER EXEC ========================
438
439void AliAnalysisTaskHelium3PiMC::UserExec(Option_t *)
440{
441 //_______________________________________________________________________
442
443 //!*********************!//
444 //! Define variables !//
445 //!*********************!//
446 Float_t vett1[60];
447 for(Int_t i=0;i<60;i++) vett1[i]=0;
448
449 Float_t vett2[40];
450 for(Int_t i=0;i<40;i++) vett2[i]=0;
451
c683985a 452
f26537d8 453 Double_t pinTPC=0.,TPCSignal=0.;
c683985a 454 Double_t xPrimaryVertex=0.,yPrimaryVertex=0.,zPrimaryVertex=0.;
455
456 ULong_t status=0;
457 ULong_t statusT=0;
458 ULong_t statusPi=0;
459
f26537d8 460 Bool_t isTPC=kFALSE,isTOFHe3=kFALSE,isTOFPi=kFALSE;
c683985a 461
462 Double_t fPos[3]={0.,0.,0.};
463 Double_t runNumber=0.;
464 Double_t evNumber=0.;
465
466 Int_t id0 = 0, id1 = 0;
467 Double_t mcDecayPosXD0 = 0, mcDecayPosYD0 = 0, mcDecayPosR = 0, mcDecayPosZD0 =0, mcDecayPosRho=0.;
468 Double_t mcDecayPosXD1 = 0, mcDecayPosYD1 = 0, mcDecayPosZD1 =0;
469
470 Double_t lEtaCurrentPart =0., lPtCurrentPart = 0.,lThetaCurrentPart = 0., lPhiCurrentPart = 0.;
471 Int_t iCurrentMother = 0;
f26537d8 472 Double_t mcPosX = 0., mcPosY = 0.,mcPosZ = 0.;
c683985a 473
f26537d8 474 Double_t lPdgCurrentDaughter0 = 0, lPdgCurrentDaughter1= 0., /*lPdgCurrentMother=0.,*/lPdgCurrentDaughter =0;
c683985a 475
476 Double_t PxD0 = 0, PyD0 = 0,PzD0 = 0;
477 Double_t PxD1 = 0, PyD1 = 0,PzD1 = 0;
478
c683985a 479 Int_t lNbMCPart = 0;
480 Int_t lPdgcodeCurrentPart = 0;
481 //!----------------------------------------------------------------
482
483 //! A set of very loose parameters for cuts
484
485 Double_t fgChi2max=33.; //! max chi2
486 Double_t fgDNmin=0.05; //! min imp parameter for the 1st daughter = 500um
487 Double_t fgDCAmax=1.; //! max DCA between the daughter tracks in cm
488 Double_t fgCPAmin=0.9; //! min cosine of V0's pointing angle
489 Double_t fgRmin=0.1; //! min radius of the fiducial volume = 1 mm
490 Double_t fgRmax=200.; //! max radius of the fiducial volume = 2 m
491
492 //------------------------------------------
493 // create pointer to event
494
495 AliVEvent *event = InputEvent();
496 if (!event) { Printf("ERROR: Could not retrieve event"); return; }
497
498
499
500// AliVEvent *lESDevent = InputEvent();
501// if (!lESDevent) {
502// Printf("ERROR: Could not retrieve event");
503// return;
504// }
505
506 Info("AliAnalysisTaskHelium3PiMC","Starting UserExec");
507
508 SetDataType("SIM");
f26537d8 509 AliStack *stack = 0;
c683985a 510 if(fDataType == "SIM") {
511
512 // Main loop
513 // Called for EACH event
514 AliMCEvent *mcEvent = MCEvent();
515 if (!mcEvent) {
516 Printf("ERROR: Could not retrieve MC event");
517 return;
518 }
519
520 Printf("MC particles: %d", mcEvent->GetNumberOfTracks());
521
522 // set up a stack for use in check for primary/stable particles
523 stack = mcEvent->Stack();
524 if( !stack ) { Printf( "Stack not available"); return; }
525 }
526
f26537d8 527 else{
528 Printf( "This Task Works Only on Simulation");
529 return;
530 }
531 AliESDEvent *lESDevent = 0;
c683985a 532
533 //********************************** Connect to the InputEvent ******//
534
535 //Int_t TrackNumber = 0;
536 if(fAnalysisType == "ESD"){
537 lESDevent = dynamic_cast<AliESDEvent*>(event);
538 if (!lESDevent) {
539 Printf("ERROR: lESDevent not available \n");
540 return;
541 }
542 }
543
f26537d8 544 else{
545 Printf("This Analysis Works Only for ESD\n");
546 return;
547 }
548
c683985a 549 //*****************//
550 //* Centrality *//
551 //*****************//
552
f26537d8 553 AliCentrality *centrality = (AliCentrality*)lESDevent->GetCentrality();
c683985a 554
555 Float_t percentile=centrality->GetCentralityPercentile("V0M");
556
557 //------------------------------
558
559 runNumber = lESDevent->GetRunNumber();
560 evNumber =lESDevent->GetEventNumberInFile();
561
562 //---------------------
563
564 // Int_t primary = stack->GetNprimary();
565 Int_t label =-1;
566
c683985a 567 lNbMCPart = stack->GetNtrack();
568
569 fHistMCMultiplicityTracks->Fill(lNbMCPart); //histo
570
571 TArrayD MomPionsMC(lNbMCPart); //Neg pions
572 Int_t nPionsMC=0;
573 TArrayD MomHeMC(lNbMCPart); //helium3
574 Int_t nHeMC=0;
575
576 //------ Trimomento pion
577 TArrayD PxPionsMC(lNbMCPart);
578 Int_t nPxPionsMC=0;
579 TArrayD PyPionsMC(lNbMCPart);
580 Int_t nPyPionsMC=0;
581 TArrayD PzPionsMC(lNbMCPart);
582 Int_t nPzPionsMC=0;
583 //------ Trimomento He
584 TArrayD PxHeMC(lNbMCPart);
585 Int_t nPxHeMC=0;
586 TArrayD PyHeMC(lNbMCPart);
587 Int_t nPyHeMC=0;
588 TArrayD PzHeMC(lNbMCPart);
589 Int_t nPzHeMC=0;
590
591 //Mass Definition
592
593 Double_t Helium3Mass = 2.80839;
594 Double_t PionMass = 0.13957;
595
596 TLorentzVector vPionMC,vHeliumMC,vSumMC;
597 TLorentzVector vPionMum,vHeliumMum,vSumMum;
598 TLorentzVector vPionRec,vHeliumRec,vSumRec;
599 Bool_t isTwoBody=kFALSE;
600
601 for (Int_t iMC=0; iMC<stack->GetNtrack(); iMC++)
602 {
603 TParticle *p0 = stack->Particle(iMC);
604
605 if (!p0) {
606 //Printf("ERROR: particle with label %d not found in stack (mc loop)", iMc);
607 continue;
608 }
609
610 lPdgcodeCurrentPart = p0->GetPdgCode();
611
612 if(lPdgcodeCurrentPart == 1000020030 || lPdgcodeCurrentPart == -1000020030 ){
613
614 MomHeMC[nHeMC++]=p0->P();
615
616 PxHeMC[nPxHeMC++]=p0->Px();
617 PyHeMC[nPyHeMC++]=p0->Py();
618 PzHeMC[nPzHeMC++]=p0->Pz();
619
620 fhHeMCRigidity->Fill(p0->P()/2);
621 }
622
623 if(lPdgcodeCurrentPart == 211 || lPdgcodeCurrentPart == -211 ){
624
625 MomPionsMC[nPionsMC++]=p0->P();
626
627 PxPionsMC[nPxPionsMC++]=p0->Px();
628 PyPionsMC[nPyPionsMC++]=p0->Py();
629 PzPionsMC[nPzPionsMC++]=p0->Pz();
630
631 fhPioneMC->Fill(p0->P());
632 }
633
634 if ( lPdgcodeCurrentPart == 1010010030 || lPdgcodeCurrentPart == -1010010030 ){
635
636 lEtaCurrentPart = p0->Eta();
637 lPtCurrentPart = p0->Pt();
638 lThetaCurrentPart = p0->Theta();
639 lPhiCurrentPart = p0->Phi();
640 iCurrentMother = p0->GetFirstMother();
641
642 fHistMCEta->Fill(lEtaCurrentPart);
643 fHistMCPt->Fill(lPtCurrentPart);
644 fHistMCTheta->Fill(lThetaCurrentPart);
645
f26537d8 646 // if (iCurrentMother == -1){lPdgCurrentMother=0; } else {lPdgCurrentMother = stack->Particle(iCurrentMother)->GetPdgCode();}
c683985a 647
648 mcPosX = p0->Vx();
649 mcPosY = p0->Vy();
650 mcPosZ = p0->Vz();
c683985a 651
652 isTwoBody=kFALSE;
653
654 for(Int_t i=p0->GetFirstDaughter(); i<= p0->GetLastDaughter(); i++){
655 TParticle *pDaughter = stack->Particle(i);
656 lPdgCurrentDaughter= pDaughter->GetPdgCode();
657 cout<<lPdgCurrentDaughter<<endl;
658 if(lPdgCurrentDaughter == 1000020030 || lPdgCurrentDaughter ==-1000020030 ){
659 isTwoBody=kTRUE;
660
661 }
662 }
663
664 if(isTwoBody){
665
666 for(Int_t i=p0->GetFirstDaughter(); i<= p0->GetLastDaughter(); i++){
667
668 TParticle *pDaughter = stack->Particle(i);
669
670 lPdgCurrentDaughter= pDaughter->GetPdgCode();
671
672 if(lPdgCurrentDaughter == 211 || lPdgCurrentDaughter == -211 ){
673 id0 = i;
674 }
675
676 if(lPdgCurrentDaughter == 1000020030 || lPdgCurrentDaughter == -1000020030 ){
677 id1 = i;
678 }
679 }
680
681 TParticle *pDaughter0 = stack->Particle(id0);
682 TParticle *pDaughter1 = stack->Particle(id1);
683 lPdgCurrentDaughter0 = pDaughter0->GetPdgCode();
684 lPdgCurrentDaughter1 = pDaughter1->GetPdgCode();
685
686 // Decay Radius and Production Radius
687
688 if ( id0 <= lNbMCPart && id0 > 0 && id1 <= lNbMCPart && id1 > 0) {
689
690 lPdgCurrentDaughter0 = pDaughter0->GetPdgCode();
691 lPdgCurrentDaughter1 = pDaughter1->GetPdgCode();
692
693 PxD0 = pDaughter0->Px();
694 PyD0 = pDaughter0->Py();
695 PzD0 = pDaughter0->Pz();
696
697 PxD1 = pDaughter1->Px();
698 PyD1 = pDaughter1->Py();
699 PzD1 = pDaughter1->Pz();
700
701 mcDecayPosXD0 = pDaughter0->Vx();
702 mcDecayPosYD0 = pDaughter0->Vy();
703 mcDecayPosZD0 = pDaughter0->Vz();
704
705 mcDecayPosXD1 = pDaughter0->Vx();
706 mcDecayPosYD1 = pDaughter0->Vy();
707 mcDecayPosZD1 = pDaughter0->Vz();
708
709 mcDecayPosR = TMath::Sqrt(mcDecayPosXD0*mcDecayPosXD0+mcDecayPosYD0*mcDecayPosYD0);
710 fHistMCDecayPosition->Fill(mcDecayPosR);
711
712 mcDecayPosRho = TMath::Sqrt(mcDecayPosXD0*mcDecayPosXD0+mcDecayPosYD0*mcDecayPosYD0+mcDecayPosZD0*mcDecayPosZD0);
713 fHistMCDecayRho->Fill(mcDecayPosRho);
714
715 //---- Initial mass Test
716
717 vHeliumMum.SetXYZM(PxD1,PyD1,PzD1,Helium3Mass);
718 vPionMum.SetXYZM(PxD0,PyD0,PzD0,PionMass);
719 vSumMum=vHeliumMum+vPionMum;
720
721 fhInvMassMum->Fill(vSumMum.M());
722
723 //Ntupla hyper triton
724
725 vett2[0]=(Float_t)lESDevent->GetRunNumber();
726 vett2[1]=(Float_t)lESDevent->GetEventNumberInFile();
727 vett2[2]=(Float_t)iMC;
728 vett2[3]=(Float_t)percentile;
729 vett2[4]=(Float_t)mcPosX;
730 vett2[5]=(Float_t)mcPosY;
731 vett2[6]=(Float_t)mcPosZ;
732 vett2[7]=(Float_t)lPdgcodeCurrentPart;
733 vett2[8]=(Float_t)iCurrentMother;
734 vett2[9]=(Float_t)mcDecayPosXD0;
735 vett2[10]=(Float_t)mcDecayPosYD0;
736 vett2[11]=(Float_t)mcDecayPosZD0;
737 vett2[12]=(Float_t)mcDecayPosXD1;
738 vett2[13]=(Float_t)mcDecayPosYD1;
739 vett2[14]=(Float_t)mcDecayPosZD1;
740 vett2[15]=(Float_t)mcDecayPosRho;
741 vett2[16]=(Float_t)lEtaCurrentPart;
742 vett2[17]=(Float_t)p0->Y();
743 vett2[18]=(Float_t)lThetaCurrentPart;
744 vett2[19]=(Float_t)lPhiCurrentPart;
745 vett2[20]=(Float_t)p0->Px();
746 vett2[21]=(Float_t)p0->Py();
747 vett2[22]=(Float_t)p0->Pz();
748 vett2[23]=(Float_t)lPdgCurrentDaughter0;
749 vett2[24]=(Float_t)lPdgCurrentDaughter1;
750 vett2[25]=(Float_t)PxD0; //pion
751 vett2[26]=(Float_t)PyD0;
752 vett2[27]=(Float_t)PzD0;
753 vett2[28]=(Float_t)PxD1; //He3
754 vett2[29]=(Float_t)PyD1;
755 vett2[30]=(Float_t)PzD1;
756
757 fNtuple2->Fill(vett2);
758
759 }//if check daughters index
760 }//is twobody
761 } // if mother
762 } // Kinetic Track loop ends here
763
764 // Loop phase - space
765
766 Double_t HeMomMC =0;
767 Double_t PionMomMC=0;
768 Double_t PxHeMc=0, PyHeMc=0,PzHeMc=0;
769 Double_t PxPionMc=0, PyPionMc=0,PzPionMc=0;
770
771 for(Int_t l=0; l < nHeMC; l++){
772
773 HeMomMC=MomHeMC[l];
774
775 PxHeMc=PxHeMC[l];
776 PyHeMc=PyHeMC[l];
777 PzHeMc=PzHeMC[l];
778
779 for(Int_t k=0; k < nPionsMC; k++){
780
781 PionMomMC=MomPionsMC[k];
782
783 PxPionMc=PxPionsMC[k];
784 PyPionMc=PyPionsMC[k];
785 PzPionMc=PzPionsMC[k];
786
787 fhRigidityHevsMomPiMC->Fill(HeMomMC/2,PionMomMC);
788
789 vHeliumMC.SetXYZM(PxHeMc,PyHeMc,PzHeMc,Helium3Mass);
790 vPionMC.SetXYZM(PxPionMc,PyPionMc,PzPionMc,PionMass);
791 vSumMC=vHeliumMC+vPionMC;
792
793 fhInvMassMC->Fill(vSumMC.M());
794
795 }
796
797 } // end loop phase space
798
799 //-------------- RECONSTRUCTION -------------------
800
801 fHistEventMultiplicity->Fill(0);
802
803 Double_t lMagneticField=lESDevent->GetMagneticField();
804
805 Int_t TrackNumber = -1;
806
807 // ANALISYS reconstructed tracks
808
809 // Primary vertex cut
810
811 const AliESDVertex *vtx = lESDevent->GetPrimaryVertexTracks();
812
813 if(vtx->GetNContributors()<1) {
814
815 // SPD vertex cut
816 vtx = lESDevent->GetPrimaryVertexSPD();
817
818 if(vtx->GetNContributors()<1) {
819 Info("AliAnalysisTaskHelium3PiMC","No good vertex, skip event");
820 return; // NO GOOD VERTEX, SKIP EVENT
821 }
822 }
823
824 fHistEventMultiplicity->Fill(1); // analyzed events with PV
825
e690d4d0 826 xPrimaryVertex=vtx->GetX();
827 yPrimaryVertex=vtx->GetY();
828 zPrimaryVertex=vtx->GetZ();
c683985a 829
830 TrackNumber = lESDevent->GetNumberOfTracks();
831 fHistTrackMultiplicity->Fill(TrackNumber); //tracce per evento
832
833 fHistPercentileVsTrackNumber->Fill(percentile,TrackNumber);
834
835 if (TrackNumber<2) return;
836
837 fHistEventMultiplicity->Fill(2);
838
839 //Find Pair candidates
840
841 TArrayI PionsTPC(TrackNumber); //Neg pions
842 Int_t nPionsTPC=0;
843
844 TArrayI HeTPC(TrackNumber); //helium3
845 Int_t nHeTPC=0;
846
847 // find pairs candidates phase daughter tracks rec
848
849 TArrayD MomPionsRec(TrackNumber); //Neg pions
850 Int_t nPionsRec=0;
851
852 TArrayD MomHeRec(TrackNumber); //helium3
853 Int_t nHeRec=0;
854
855 //------ Trimomento pion
856 TArrayD PxPionsRec(TrackNumber);
857 Int_t nPxPionsRec=0;
858 TArrayD PyPionsRec(TrackNumber);
859 Int_t nPyPionsRec=0;
860 TArrayD PzPionsRec(TrackNumber);
861 Int_t nPzPionsRec=0;
862
863 //------ Trimomento He
864 TArrayD PxHeRec(TrackNumber);
865 Int_t nPxHeRec=0;
866 TArrayD PyHeRec(TrackNumber);
867 Int_t nPyHeRec=0;
868 TArrayD PzHeRec(TrackNumber);
869 Int_t nPzHeRec=0;
870
871 Float_t impactXY=-999, impactZ=-999;
872 Float_t impactXYpi=-999, impactZpi=-999;
873
874 Int_t PDGCodePos=0;
875 Int_t PDGCodeNeg=0;
876 Int_t motherPDGNeg=0;
877 Int_t motherPDGPos=0;
878
c683985a 879 //! SELECTIONS:
880 //! - No ITSpureSA
881 //! - ITSrefit
882 //! - TPCrefit
883 //! - ITSmap !=0
884
885 // ******************* Track Cuts Definitions ********************//
886
887 //! ITS
888
889 AliESDtrackCuts* esdtrackCutsITS = new AliESDtrackCuts("esdtrackCutsITS");
890 esdtrackCutsITS->SetRequireITSStandAlone(kFALSE);
891 esdtrackCutsITS->SetRequireITSPureStandAlone(kFALSE);
892
893 //! TPC
894
895 Int_t minclsTPC=60;
896 Double_t maxchi2perTPCcl=5.;
897
898 AliESDtrackCuts* esdtrackCutsTPC = new AliESDtrackCuts("esdtrackCutsTPC");
899 esdtrackCutsTPC->SetRequireTPCRefit(kTRUE);
900 esdtrackCutsTPC->SetAcceptKinkDaughters(kFALSE);
901 esdtrackCutsTPC->SetMinNClustersTPC(minclsTPC);
902 esdtrackCutsTPC->SetMaxChi2PerClusterTPC(maxchi2perTPCcl);
903
904 //*************************************************************
905
906 for (Int_t j=0; j<TrackNumber; j++) { //loop on tracks
907
908 AliESDtrack *esdtrack=lESDevent->GetTrack(j);
909
910 if(!esdtrack) {
911 AliError(Form("ERROR: Could not retrieve esdtrack %d",j));
912 continue;
913 }
914
915 hBBTPCnoCuts->Fill(esdtrack->GetSign()*esdtrack->P(),esdtrack->GetTPCsignal());
916
917 // ************** Track cuts ****************
918
919 status = (ULong_t)esdtrack->GetStatus();
920
921 isTPC = (((status) & (AliESDtrack::kTPCin)) != 0);
c683985a 922
923 Bool_t IsTrackAcceptedTPC = esdtrackCutsTPC->AcceptTrack(esdtrack);
924 Bool_t IsTrackAcceptedITS = esdtrackCutsITS->AcceptTrack(esdtrack);
925
926 if (!(IsTrackAcceptedTPC && IsTrackAcceptedITS)) continue;
927
928 //----------------------------------------------
929
930 //****** Cuts from AliV0Vertex.cxx *************
931
932 Double_t d=esdtrack->GetD(xPrimaryVertex,yPrimaryVertex,lMagneticField);
933 // if (TMath::Abs(d)<fgDPmin) continue;
934 if (TMath::Abs(d)>fgRmax) continue;
935
936 //---- (Usefull) Stuff
937
938 TPCSignal=esdtrack->GetTPCsignal();
939
940 if (TPCSignal<10)continue;
941
942 if(!isTPC)continue;
943
944 if(!esdtrack->GetTPCInnerParam())continue;
945
946 AliExternalTrackParam trackIn(*esdtrack->GetInnerParam());
947 pinTPC = trackIn.GetP();
948
c683985a 949 fhBBTPC->Fill(pinTPC*esdtrack->GetSign(),TPCSignal);
950
951 d=esdtrack->GetD(xPrimaryVertex,yPrimaryVertex,lMagneticField);
952 // if (TMath::Abs(d)<fgDPmin) continue;
953 if (TMath::Abs(d)>fgRmax) continue;
954
955 label = TMath::Abs(esdtrack->GetLabel());
956
957 if (label>=10000000) {
958 // Underlying event. 10000000 is the
959 // value of fkMASKSTEP in AliRunDigitizer
960 cout <<"Strange, there should be no underlying event"<<endl;
961 }
962
963 else {
964 if (label>=lNbMCPart) {
965 cout <<"Strange, label outside the range"<< endl;
966 continue;
967 }
968 }
969
970 TParticle * part = stack->Particle(label);
971
972 Int_t PDGCode=part->GetPdgCode();
c683985a 973
974 if(PDGCode==-211)
975 fhBBTPCNegativePions->Fill(esdtrack->GetSign()*esdtrack->P(),esdtrack->GetTPCsignal());
976
977 if(PDGCode==+211)
978 fhBBTPCPositivePions->Fill(esdtrack->GetSign()*esdtrack->P(),esdtrack->GetTPCsignal());
979
980
981 // if(PDGCode == 211){
982
983 if(PDGCode==-211 || PDGCode==+211){
984
985 PionsTPC[nPionsTPC++]=j;
986
987 esdtrack->GetImpactParameters(impactXY, impactZ);
988 fhPiDCAXY->Fill(impactXY);
989 fhPiDCAZ->Fill(impactZ);
990
991 MomPionsRec[nPionsRec++]=esdtrack->P();
992
993 PxPionsRec[nPxPionsRec++]=esdtrack->Px();
994 PyPionsRec[nPyPionsRec++]=esdtrack->Py();
995 PzPionsRec[nPzPionsRec++]=esdtrack->Pz();
996
997 }
998
999 if(PDGCode==1000020030 ||PDGCode==-1000020030 ){
1000
1001
1002 HeTPC[nHeTPC++]=j;
1003
1004 fhBBTPCHe3->Fill(esdtrack->GetSign()*esdtrack->P(),esdtrack->GetTPCsignal());
1005
1006 esdtrack->GetImpactParameters(impactXY, impactZ);
1007 fhHeDCAXY->Fill(impactXY);
1008 fhHeDCAZ->Fill(impactZ);
1009
1010 MomHeRec[nHeRec++]=esdtrack->P();
1011
1012 PxHeRec[nPxHeRec++]=esdtrack->Px();
1013 PyHeRec[nPyHeRec++]=esdtrack->Py();
1014 PzHeRec[nPzHeRec++]=esdtrack->Pz();
1015
1016 }
1017
1018 } // ! track
1019
1020 //-------------- LOOP pairs 1 -------------
1021 // Fill phase space and inva mass before cuts
1022
1023 Double_t HeMomRec =0;
1024 Double_t PionMomRec=0;
1025 Double_t PxHeReco=0, PyHeReco=0,PzHeReco=0;
1026 Double_t PxPionReco=0, PyPionReco=0,PzPionReco=0;
1027
1028 for(Int_t l=0; l < nHeRec; l++){
1029
1030 HeMomRec=MomHeRec[l];
1031
1032 PxHeReco=PxHeRec[l];
1033 PyHeReco=PyHeRec[l];
1034 PzHeReco=PzHeRec[l];
1035
1036 for(Int_t k=0; k < nPionsRec; k++){
1037
1038 PionMomRec=MomPionsRec[k];
1039
1040 PxPionReco=PxPionsRec[k];
1041 PyPionReco=PyPionsRec[k];
1042 PzPionReco=PzPionsRec[k];
1043
1044 fhRigidityHevsMomPiRec->Fill(HeMomRec,PionMomRec);
1045
1046 vHeliumRec.SetXYZM(2*PxHeReco,2*PyHeReco,2*PzHeReco,Helium3Mass);
1047 vPionRec.SetXYZM(PxPionReco,PyPionReco,PzPionReco,PionMass);
1048 vSumRec=vHeliumRec+vPionRec;
1049
1050 fhInvMassRec->Fill(vSumRec.M());
1051 }
1052
1053 } // fine loop phase space
1054
1055 //--------------- LOOP PAIRS ----------------------//
1056
1057 Double_t DcaHeToPrimVertex=0;
1058 Double_t DcaPionToPrimVertex=0;
1059
1060 impactXY=-999, impactZ=-999;
1061 impactXYpi=-999, impactZpi=-999;
1062
c683985a 1063
1064 // Vettori per il PxPyPz
1065
1066 Double_t momPionVett[3];
1067 for(Int_t i=0;i<3;i++)momPionVett[i]=0;
1068
1069 Double_t momHeVett[3];
1070 for(Int_t i=0;i<3;i++)momHeVett[i]=0;
1071
1072 //At SV
1073
1074 Double_t momPionVettAt[3];
1075 for(Int_t i=0;i<3;i++)momPionVettAt[i]=0;
1076
1077 Double_t momHeVettAt[3];
1078 for(Int_t i=0;i<3;i++)momHeVettAt[i]=0;
1079
1080 Bool_t IsHeITSRefit,IsPiITSRefit ;
1081
1082 //----------- My 2nd Vertex Finder
1083
1084 for (Int_t k=0; k < nPionsTPC; k++) { //! Pions Loop
1085
1086 DcaPionToPrimVertex=0.;
1087 DcaHeToPrimVertex=0;
1088
1089 Int_t PionIdx=PionsTPC[k];
1090
f26537d8 1091 AliESDtrack *PionTrack=lESDevent->GetTrack(PionIdx);
c683985a 1092
1093 statusPi = (ULong_t)PionTrack->GetStatus();
1094 IsPiITSRefit = ((statusPi) & (AliESDtrack::kITSrefit));
1095
1096 Int_t labelPi = TMath::Abs(PionTrack->GetLabel());
1097 TParticle * partNeg = stack->Particle(labelPi);
1098 PDGCodeNeg=partNeg->GetPdgCode();
1099
1100 Int_t mumidNeg = partNeg->GetFirstMother();
1101 if(mumidNeg>-1){
1102 TParticle *motherNeg=(TParticle*)stack->Particle(mumidNeg);
1103 motherPDGNeg = motherNeg->GetPdgCode();
1104 }
1105
1106 if (PionTrack)
1107 DcaPionToPrimVertex = TMath::Abs(PionTrack->GetD(xPrimaryVertex, yPrimaryVertex,lMagneticField)); //OK
1108
1109 if(DcaPionToPrimVertex<0.1)continue;
1110
1111 AliExternalTrackParam trackInPion(*PionTrack);
1112
1113 for (Int_t i=0; i<nHeTPC; i++){ //! Helium Loop
1114
1115 Int_t HeIdx=HeTPC[i];
1116
f26537d8 1117 AliESDtrack *HeTrack=lESDevent->GetTrack(HeIdx);
c683985a 1118
1119 statusT= (ULong_t)HeTrack->GetStatus();
1120 IsHeITSRefit = ((statusT) & (AliESDtrack::kITSrefit));
1121
1122 Int_t labelHe = TMath::Abs(HeTrack->GetLabel());
1123 TParticle * partPos = stack->Particle(labelHe);
1124 PDGCodePos=partPos->GetPdgCode();
1125
1126 Int_t mumidPos = partPos->GetFirstMother();
1127 if(mumidPos>-1){
1128 TParticle *motherPos=(TParticle*)stack->Particle(mumidPos);
1129 motherPDGPos = motherPos->GetPdgCode();
1130 }
1131
1132 if (HeTrack)
1133 DcaHeToPrimVertex = TMath::Abs(HeTrack->GetD(xPrimaryVertex, yPrimaryVertex,lMagneticField)); //OK
1134
1135 AliExternalTrackParam trackInHe(*HeTrack);
1136
1137 HeTrack->PxPyPz(momHeVett);
1138 PionTrack->PxPyPz(momPionVett);
1139
1140 vHeliumRec.SetXYZM(2*momHeVett[0],2*momHeVett[1],2*momHeVett[2],Helium3Mass);
1141 vPionRec.SetXYZM(momPionVett[0],momPionVett[1],momPionVett[2],PionMass);
1142 vSumRec=vHeliumRec+vPionRec;
1143
1144 fhInvMassRec1->Fill(vSumRec.M());
1145
1146 fhInvMassRec2->Fill(vSumRec.M());
1147
1148 if ( DcaPionToPrimVertex < fgDNmin) //OK
1149 if ( DcaHeToPrimVertex < fgDNmin) continue; //OK
1150
1151 fhInvMassRec3->Fill(vSumRec.M());
1152
1153 Double_t xn, xp;
1154 Double_t dca=0.;
1155
1156 dca= PionTrack->GetDCA(HeTrack,lMagneticField,xn,xp); //!distance between two tracks (Neg to Pos)
1157 fHistProvaDCA->Fill(xn-xp,dca);
1158 if (dca > fgDCAmax) continue;
1159
1160 fhInvMassRec4->Fill(vSumRec.M());
1161
1162 if ((xn+xp) > 2*fgRmax) continue;
1163 if ((xn+xp) < 2*fgRmin) continue;
1164 fhInvMassRec5->Fill(vSumRec.M());
1165
1166 //CORREZIONE da AliV0Vertex
1167
1168 Bool_t corrected=kFALSE;
1169 if ((trackInPion.GetX() > 3.) && (xn < 3.)) {
1170 //correct for the beam pipe material
1171 corrected=kTRUE;
1172 }
1173 if ((trackInHe.GetX() > 3.) && (xp < 3.)) {
1174 //correct for the beam pipe material
1175 corrected=kTRUE;
1176 }
1177 if (corrected) {
1178 dca=trackInPion.GetDCA(&trackInHe,lMagneticField,xn,xp);
1179 if (dca > fgDCAmax) continue;
1180 if ((xn+xp) > 2*fgRmax) continue;
1181 if ((xn+xp) < 2*fgRmin) continue;
1182 }
1183
1184 //=============================================//
1185 // Make a "V0" with Tracks //
1186 //=============================================//
1187
1188 trackInPion.PropagateTo(xn,lMagneticField);
1189 trackInHe.PropagateTo(xp,lMagneticField);
1190
1191 AliESDv0 vertex(trackInPion,PionIdx,trackInHe,HeIdx);
1192 if (vertex.GetChi2V0() > fgChi2max) continue;
1193 fhInvMassRec6->Fill(vSumRec.M());
1194
1195 Float_t CosPointingAngle=vertex.GetV0CosineOfPointingAngle(xPrimaryVertex,yPrimaryVertex,zPrimaryVertex); //PointingAngle
1196 if (CosPointingAngle < fgCPAmin) continue;
1197
1198 fhInvMassRec7->Fill(vSumRec.M());
1199
1200 vertex.SetDcaV0Daughters(dca);
1201 vertex.SetV0CosineOfPointingAngle(CosPointingAngle);
1202
1203 fPos[0]=vertex.Xv();
1204 fPos[1]=vertex.Yv();
1205 fPos[2]=vertex.Zv();
1206
1207
1208
1209 Double_t raggio=TMath::Sqrt(fPos[0]*fPos[0]+fPos[1]*fPos[1]+fPos[2]*fPos[2]);
1210 HeTrack->GetPxPyPzAt(raggio,lMagneticField,momHeVettAt);
1211 PionTrack->GetPxPyPzAt(raggio,lMagneticField,momPionVettAt);
1212
1213 //------------------------------------------------------------------------//
1214
1215 HeTrack->GetImpactParameters(impactXY, impactZ);
1216
1217 PionTrack->GetImpactParameters(impactXYpi, impactZpi);
1218
1219 Float_t timeTOFHe= HeTrack->GetTOFsignal(); // ps
1220 Float_t trackLenghtTOFHe= HeTrack->GetIntegratedLength(); // cm
1221
1222 Float_t timeTOFPi= PionTrack->GetTOFsignal(); // ps
1223 Float_t trackLenghtTOFPi= PionTrack->GetIntegratedLength(); // cm
1224
1225 //----------------------------------------------------------------------//
1226
1227 vett1[0]=(Float_t)runNumber;
1228 vett1[1]=(Float_t)evNumber;
1229 vett1[2]=(Float_t)lNbMCPart;
1230 vett1[3]=(Float_t)percentile;
1231 vett1[4]=(Float_t)xPrimaryVertex; //PRIMARY
1232 vett1[5]=(Float_t)yPrimaryVertex;
1233 vett1[6]=(Float_t)zPrimaryVertex;
1234 vett1[7]=(Float_t)fPos[0]; //SECONDARY
1235 vett1[8]=(Float_t)fPos[1];
1236 vett1[9]=(Float_t)fPos[2];
1237 vett1[10]=(Float_t)dca; //between 2 tracks
1238 vett1[11]=(Float_t)CosPointingAngle; //cosPointingAngle da V0
1239 vett1[12]=(Float_t)vertex.GetD(xPrimaryVertex,yPrimaryVertex,zPrimaryVertex);
1240 vett1[13]=(Float_t)HeTrack->GetSign(); //He
1241 vett1[14]=(Float_t)trackInHe.GetP();
1242 vett1[15]=(Float_t)HeTrack->GetTPCsignal();
1243 vett1[16]=(Float_t)DcaHeToPrimVertex;
1244 vett1[17]=(Float_t)HeTrack->Eta();
1245 vett1[18]=(Float_t)momHeVett[0];
1246 vett1[19]=(Float_t)momHeVett[1];
1247 vett1[20]=(Float_t)momHeVett[2];
1248 vett1[21]=(Float_t)momHeVettAt[0];
1249 vett1[22]=(Float_t)momHeVettAt[1];
1250 vett1[23]=(Float_t)momHeVettAt[2];
1251 vett1[24]=(Float_t)HeTrack->GetTPCNcls();
1252 vett1[25]=(Float_t)impactXY;
1253 vett1[26]=(Float_t)impactZ;
1254 vett1[27]=(Float_t)isTOFHe3;
1255 vett1[28]=(Float_t)(trackLenghtTOFHe/timeTOFHe)/2.99792458e-2;
1256 vett1[29]=(Float_t)HeTrack->GetITSClusterMap();
1257 vett1[30]=(Float_t)IsHeITSRefit;
1258 vett1[31]=(Float_t)PionTrack->GetSign(); //Pion
1259 vett1[32]=(Float_t)trackInPion.GetP();
1260 vett1[33]=(Float_t)PionTrack->GetTPCsignal();
1261 vett1[34]=(Float_t)DcaPionToPrimVertex;
1262 vett1[35]=(Float_t)PionTrack->Eta();
1263 vett1[36]=(Float_t)momPionVett[0];
1264 vett1[37]=(Float_t)momPionVett[1];
1265 vett1[38]=(Float_t)momPionVett[2];
1266 vett1[39]=(Float_t)momPionVettAt[0];
1267 vett1[40]=(Float_t)momPionVettAt[1];
1268 vett1[41]=(Float_t)momPionVettAt[2];
1269 vett1[42]=(Float_t)PionTrack->GetTPCNcls();
1270 vett1[43]=(Float_t)impactXYpi;
1271 vett1[44]=(Float_t)impactZpi;
1272 vett1[45]=(Float_t)isTOFPi;
1273 vett1[46]=(Float_t)(trackLenghtTOFPi/timeTOFPi)/2.99792458e-2;
1274 vett1[47]=(Float_t)PionTrack->GetITSClusterMap();
1275 vett1[48]=(Float_t)IsPiITSRefit;
1276 vett1[49]=(Float_t)PDGCodeNeg;
1277 vett1[50]=(Float_t)PDGCodePos;
1278 vett1[51]=(Float_t)motherPDGNeg;
1279 vett1[52]=(Float_t)motherPDGPos;
1280 vett1[53]=(Float_t)labelPi;
1281 vett1[54]=(Float_t)labelHe;
1282 vett1[55]=(Float_t)mumidNeg;
1283 vett1[56]=(Float_t)mumidPos;
1284
1285 fNtuple1->Fill(vett1);
1286
1287 }// positive TPC
1288
1289 } //negative tpc
1290
1291 PostData(1,fListHistCascade);
1292
1293}// end userexec
1294
1295
1296//________________________________________________________________________
1297//
1298void AliAnalysisTaskHelium3PiMC::Terminate(Option_t *)
1299{
1300 // Draw result to the screen
1301 // Called once at the end of the query
1302}
1303