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