Fixes for cmake
[u/mrichter/AliRoot.git] / PWG1 / AliGenV0Info.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16
17 ///////////////////////////////////////////////////////////////////////////
18 /*
19
20 Origin: marian.ivanov@cern.ch
21 Container classes with MC infomation for V0 
22
23
24 */
25
26 #if !defined(__CINT__) || defined(__MAKECINT__)
27 #include <stdio.h>
28 #include <string.h>
29 //ROOT includes
30 #include "TROOT.h"
31 #include "Rtypes.h"
32 #include "TFile.h"
33 #include "TTree.h"
34 #include "TChain.h"
35 #include "TCut.h"
36 #include "TString.h"
37 #include "TStopwatch.h"
38 #include "TParticle.h"
39 #include "TDatabasePDG.h"
40 #include "TSystem.h"
41 #include "TCanvas.h"
42 #include "TPolyLine3D.h"
43
44 //ALIROOT includes
45 #include "AliRun.h"
46 #include "AliStack.h"
47 #include "AliSimDigits.h"
48 #include "AliTPCParam.h"
49 #include "AliTPC.h"
50 #include "AliTPCLoader.h"
51 #include "AliDetector.h"
52 #include "AliTrackReference.h"
53 #include "AliTPCParamSR.h"
54 #include "AliTracker.h"
55 #include "AliMagF.h"
56 #include "AliHelix.h"
57 #include "AliTrackPointArray.h"
58
59 #endif
60 #include "AliGenV0Info.h" 
61 //
62 // 
63
64 ClassImp(AliGenV0Info)
65
66
67
68
69
70 /////////////////////////////////////////////////////////////////////////////////
71 AliGenV0Info::AliGenV0Info():
72   fMCd(),      //info about daughter particle - 
73   fMCm(),      //info about mother particle   - first particle for V0
74   fMotherP(),   //particle info about mother particle
75   fMCDist1(0), //info about closest distance according closest MC - linear DCA
76   fMCDist2(0),    //info about closest distance parabolic DCA
77   fMCRr(0),       // rec position of the vertex 
78   fMCR(0),        //exact r position of the vertex
79   fInvMass(0),  //reconstructed invariant mass -
80   fPointAngleFi(0), //point angle fi
81   fPointAngleTh(0), //point angle theta
82   fPointAngle(0)   //point angle full
83 {
84   for (Int_t i=0;i<3; i++){
85    fMCPdr[i]=0;    
86    fMCX[i]=0;     
87    fMCXr[i]=0;    
88    fMCPm[i]=0;    
89    fMCAngle[i]=0; 
90    fMCPd[i]=0;    
91   }
92   fMCPd[3]=0;     
93   for (Int_t i=0; i<2;i++){
94     fPdg[i]=0;   
95     fLab[i]=0;  
96   }
97 }
98
99 void AliGenV0Info::Update(Float_t vertex[3])
100 {
101   //
102   // Update information - calculates derived variables
103   //
104   
105   fMCPd[0] = fMCd.GetParticle().Px();
106   fMCPd[1] = fMCd.GetParticle().Py();
107   fMCPd[2] = fMCd.GetParticle().Pz();
108   fMCPd[3] = fMCd.GetParticle().P();
109   //
110   fMCX[0]  = fMCd.GetParticle().Vx();
111   fMCX[1]  = fMCd.GetParticle().Vy();
112   fMCX[2]  = fMCd.GetParticle().Vz();
113   fMCR       = TMath::Sqrt( fMCX[0]*fMCX[0]+fMCX[1]*fMCX[1]);
114   //
115   fPdg[0]    = fMCd.GetParticle().GetPdgCode();
116   fPdg[1]    = fMCm.GetParticle().GetPdgCode();
117   //
118   fLab[0]    = fMCd.GetParticle().GetUniqueID();
119   fLab[1]    = fMCm.GetParticle().GetUniqueID();
120   //
121   //
122   //
123   Double_t x1[3],p1[3];
124   TParticle& part0 = fMCd.GetParticle();
125   x1[0] = part0.Vx();      
126   x1[1] = part0.Vy();
127   x1[2] = part0.Vz();
128   p1[0] = part0.Px();
129   p1[1] = part0.Py();
130   p1[2] = part0.Pz();
131   if (part0.GetPDG()==0) return;
132
133   Double_t sign = (part0.GetPDG()->Charge()>0)? -1:1;
134   AliHelix dhelix1(x1,p1,sign);
135   //
136   //
137   Double_t x2[3],p2[3];            
138   //
139   TParticle& part1 = fMCm.GetParticle();
140   if (part1.GetPDG()==0) return;
141   x2[0] = part1.Vx();      
142   x2[1] = part1.Vy();      
143   x2[2] = part1.Vz();      
144   p2[0] = part1.Px();
145   p2[1] = part1.Py();
146   p2[2] = part1.Pz();
147   //
148   //
149   if (part1.GetPDG()==0) return;
150   sign = (part1.GetPDG()->Charge()>0) ? -1:1;
151   AliHelix mhelix(x2,p2,sign);
152   
153   //
154   //
155   //
156   //find intersection linear
157   //
158   Double_t distance1, distance2;
159   Double_t phase[2][2],radius[2];
160   Int_t  points = dhelix1.GetRPHIintersections(mhelix, phase, radius);
161   Double_t delta1=10000,delta2=10000;    
162   if (points>0){
163     dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
164     dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
165     dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
166   }
167   else{
168     fInvMass=-1;
169     return;
170   }
171   if (points==2){    
172     dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
173     dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
174     dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
175   }
176   distance1 = TMath::Min(delta1,delta2);
177   //
178   //find intersection parabolic
179   //
180   points = dhelix1.GetRPHIintersections(mhelix, phase, radius);
181   delta1=10000,delta2=10000;  
182   
183   if (points>0){
184     dhelix1.ParabolicDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
185   }
186   if (points==2){    
187     dhelix1.ParabolicDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
188   }
189   
190   distance2 = TMath::Min(delta1,delta2);
191   //
192   if (delta1<delta2){
193     //get V0 info
194     dhelix1.Evaluate(phase[0][0],fMCXr);
195     dhelix1.GetMomentum(phase[0][0],fMCPdr);
196     mhelix.GetMomentum(phase[0][1],fMCPm);
197     dhelix1.GetAngle(phase[0][0],mhelix,phase[0][1],fMCAngle);
198     fMCRr = TMath::Sqrt(radius[0]);
199   }
200   else{
201     dhelix1.Evaluate(phase[1][0],fMCXr);
202     dhelix1.GetMomentum(phase[1][0],fMCPdr);
203     mhelix.GetMomentum(phase[1][1],fMCPm);
204     dhelix1.GetAngle(phase[1][0],mhelix,phase[1][1],fMCAngle);
205     fMCRr = TMath::Sqrt(radius[1]);
206   }     
207   //            
208   //
209   fMCDist1 = TMath::Sqrt(distance1);
210   fMCDist2 = TMath::Sqrt(distance2);      
211   //
212   //
213   // 
214   Float_t v[3] = {fMCXr[0]-vertex[0],fMCXr[1]-vertex[1],fMCXr[2]-vertex[2]};
215   //Float_t v[3] = {fMCXr[0],fMCXr[1],fMCXr[2]};
216   Float_t p[3] = {fMCPdr[0]+fMCPm[0], fMCPdr[1]+fMCPm[1],fMCPdr[2]+fMCPm[2]};
217   Float_t vnorm2 = v[0]*v[0]+v[1]*v[1];
218   Float_t vnorm3 = TMath::Sqrt(v[2]*v[2]+vnorm2);
219   vnorm2 = TMath::Sqrt(vnorm2);
220   Float_t pnorm2 = p[0]*p[0]+p[1]*p[1];
221   Float_t pnorm3 = TMath::Sqrt(p[2]*p[2]+pnorm2);
222   pnorm2 = TMath::Sqrt(pnorm2);
223   //
224   if (vnorm2>0){
225     fPointAngleFi = (v[0]*p[0]+v[1]*p[1])/(vnorm2*pnorm2);
226     fPointAngleTh = (v[2]*p[2]+vnorm2*pnorm2)/(vnorm3*pnorm3);  
227     fPointAngle   = (v[0]*p[0]+v[1]*p[1]+v[2]*p[2])/(vnorm3*pnorm3);
228   }else{
229     fPointAngleFi = 1;
230     fPointAngleTh = 1;  
231     fPointAngle   = 1; 
232   }
233   Double_t mass1 = fMCd.GetMass();
234   Double_t mass2 = fMCm.GetMass();
235
236   
237   Double_t e1    = TMath::Sqrt(mass1*mass1+
238                               fMCPd[0]*fMCPd[0]+
239                               fMCPd[1]*fMCPd[1]+
240                               fMCPd[2]*fMCPd[2]);
241   Double_t e2    = TMath::Sqrt(mass2*mass2+
242                               fMCPm[0]*fMCPm[0]+
243                               fMCPm[1]*fMCPm[1]+
244                               fMCPm[2]*fMCPm[2]);
245   
246   fInvMass =  
247     (fMCPm[0]+fMCPd[0])*(fMCPm[0]+fMCPd[0])+
248     (fMCPm[1]+fMCPd[1])*(fMCPm[1]+fMCPd[1])+
249     (fMCPm[2]+fMCPd[2])*(fMCPm[2]+fMCPd[2]);
250   
251   //  fInvMass = TMath::Sqrt((e1+e2)*(e1+e2)-fInvMass);
252   fInvMass = (e1+e2)*(e1+e2)-fInvMass;
253   if (fInvMass>0) fInvMass = TMath::Sqrt(fInvMass);    
254 }
255
256
257
258
259
260
261
262