TENDER becomes Tender, removing .so
[u/mrichter/AliRoot.git] / PWGGA / GammaConv / AliAODConversionMother.cxx
1 #include "AliAODConversionMother.h"
2 #include "AliKFConversionMother.h"
3 #include "AliKFParticle.h"
4
5
6 // Author D. Lohner (Daniel.Lohner@cern.ch)
7
8 using namespace std;
9
10 ClassImp(AliAODConversionMother)
11
12 AliAODConversionMother::AliAODConversionMother() :
13 AliAODConversionParticle(),
14     fMCLabel(-1),
15     fChi2(-1),
16     fOpeningAngle(-1),
17     fAlpha(-1),
18     fWeight(1),
19     fdcaBetweenPhotons(1),
20     fdcaZPrimVtx(100),
21     fdcaRPrimVtx(100),
22     fQuality(0),
23     fTrueMeson(0)
24    
25 {
26         fLabel[0] = -1;
27         fLabel[1] = -1;
28         fLabel[2] = 0;
29    
30    fProductionVtx[0]=0;
31    fProductionVtx[1]=0;
32    fProductionVtx[2]=0;
33   
34 }
35
36 AliAODConversionMother::AliAODConversionMother(AliKFConversionMother *kf):
37 AliAODConversionParticle(),
38 fMCLabel(kf->GetMCLabel()),
39 fChi2(kf->GetChi2()),
40 fOpeningAngle(kf->GetOpeningAngle()),
41 fAlpha(kf->GetAlpha()),
42 fWeight(1),
43 fdcaBetweenPhotons(100),
44 fdcaZPrimVtx(100),
45 fdcaRPrimVtx(100),
46 fQuality(0),
47 fTrueMeson(0)
48 {
49     // Set 4momentu
50     SetPxPyPzE(kf->GetPx(),kf->GetPy(),kf->GetPz(),kf->GetE());
51
52    fProductionVtx[0]=0;
53    fProductionVtx[1]=0;
54    fProductionVtx[2]=0;
55   
56      //Set Decay Photon Labels
57     fLabel[0]=kf->GetGammaLabel(0);
58     fLabel[1]=kf->GetGammaLabel(1);
59     fLabel[2]=0;
60 }
61
62 AliAODConversionMother::AliAODConversionMother(AliAODConversionPhoton *y1,AliAODConversionPhoton *y2):
63 AliAODConversionParticle(),
64 fMCLabel(-1),
65 fChi2(-1),
66 fOpeningAngle(-1),
67 fAlpha(-1),
68 fWeight(1),
69 fdcaBetweenPhotons(1),
70 fdcaZPrimVtx(100),
71 fdcaRPrimVtx(100),
72 fQuality(0),
73 fTrueMeson(0)
74 {
75     // Set 4momentum
76     SetPxPyPzE(y1->Px()+y2->Px(),y1->Py()+y2->Py(),y1->Pz()+y2->Pz(),y1->E()+y2->E());
77
78     // Calculate Opening Angle
79     TVector3 v1(y1->Px(),y1->Py(),y1->Pz());
80     TVector3 v2(y2->Px(),y2->Py(),y2->Pz());
81     fOpeningAngle=v1.Angle(v2);
82     fdcaBetweenPhotons = CalculateDistanceBetweenPhotons(y1,y2,fProductionVtx);
83     DetermineMesonQuality(y1,y2);
84     // Calculate Alpha
85     if((y1->E()+y2->E()) != 0){
86                 fAlpha=TMath::Abs((y1->E()-y2->E())/(y1->E()+y2->E()));
87     }
88
89     // Set Chi2 to the mean chi2 of gammas
90  //   fChi2=0.5*(y1->GetChi2perNDF()+y2->GetChi2perNDF());
91
92     //Set Decay Photon Labels
93     fLabel[0]=-1;
94     fLabel[1]=-1;
95     fLabel[2]=0;
96 }
97
98 AliAODConversionMother::AliAODConversionMother(AliAODConversionMother *meson,AliAODConversionPhoton *gamma):
99 AliAODConversionParticle(),
100 fMCLabel(-1),
101 fChi2(-1),
102 fOpeningAngle(-1),
103 fAlpha(-1),
104 fWeight(1),
105 fdcaBetweenPhotons(1),
106 fdcaZPrimVtx(100),
107 fdcaRPrimVtx(100),
108 fQuality(0),
109 fTrueMeson(0)
110 {
111     // Set 4momentum
112     SetPxPyPzE(meson->Px()+gamma->Px(),meson->Py()+gamma->Py(),meson->Pz()+gamma->Pz(),meson->E()+gamma->E());
113
114     // Calculate Opening Angle
115     TVector3 v1(meson->Px(),meson->Py(),meson->Pz());
116     TVector3 v2(gamma->Px(),gamma->Py(),gamma->Pz());
117     fOpeningAngle=v1.Angle(v2);
118      
119         fProductionVtx[0]=0;
120         fProductionVtx[1]=0;
121         fProductionVtx[2]=0;
122
123     // Calculate Alpha
124     if((meson->E()+gamma->E()) != 0){
125                 fAlpha=TMath::Abs((meson->E()-gamma->E())/(meson->E()+gamma->E()));
126     }
127
128     // Set Chi2 to the mean chi2 of gammas
129         // fChi2=0.5*(y1->GetChi2perNDF()+y2->GetChi2perNDF());
130
131     //Set Decay Photon Labels
132     fLabel[0]=-1;
133     fLabel[1]=-1;
134     fLabel[2]=0;
135 }
136
137
138 AliAODConversionMother::~AliAODConversionMother() {
139     // empty standard destructor
140 }
141
142 TParticle *AliAODConversionMother::GetMCParticle(AliStack *fMCStack){
143     if(!fMCStack){AliError("MC Stack not defined");return 0x0;}
144
145     if(fMCLabel>-1){
146                 return fMCStack->Particle(fMCLabel);
147     }
148     return 0x0;
149 }
150
151 Bool_t AliAODConversionMother::IsTrueMeson(AliStack *fMCStack,Int_t pdgcode){
152     TParticle *part=GetMCParticle(fMCStack);
153
154     if(part){
155         // Check if it is a true photon
156         if(part->GetPdgCode()==pdgcode){
157             return kTRUE;
158         }
159     }
160     return kFALSE;
161 }
162
163 Float_t AliAODConversionMother::CalculateDistanceBetweenPhotons(AliAODConversionPhoton* y1, AliAODConversionPhoton* y2 , Double_t prodPoint[3]){
164
165    TVector3 a(y1->GetConversionX(),y1->GetConversionY(),y1->GetConversionZ());
166    TVector3 b(y1->GetPx(),y1->GetPy(),y1->GetPz());
167    TVector3 c(y2->GetConversionX(),y2->GetConversionY(),y2->GetConversionZ());
168    TVector3 d(y2->GetPx(),y2->GetPy(),y2->GetPz());
169    
170    TVector3 n = b.Cross(d);
171    TVector3 nn = n.Unit();
172    
173    Double_t dist = 0;
174    if (n.Mag() == 0){
175       TVector3 e = a-c;
176       if (d.Mag() != 0){
177          dist = TMath::Abs((e.Cross(d)).Mag())/TMath::Abs(d.Mag());
178       }
179       prodPoint[0] = 0;
180       prodPoint[1] = 0;
181       prodPoint[2] = 0;
182    } else {   
183       dist = TMath::Abs(n.Dot(c-a))/TMath::Abs(n.Mag());
184       Double_t lambda = (b.Dot(d) * (a-c).Dot(d) - d.Dot(d) * (a-c).Dot(b))/(b.Dot(b) * d.Dot(d) - TMath::Power(b.Dot(d),2));
185       Double_t mu = ((a-c).Dot(d) * b.Dot(b) - (a-c).Dot(b) * b.Dot(d) )/(b.Dot(b) * d.Dot(d) - TMath::Power(b.Dot(d),2));
186       
187       TVector3 S1 = a + lambda* b;
188       TVector3 S2 = c + mu* d;
189       TVector3 Prod = S1 + 0.5*dist*(S2-S1).Unit();
190       prodPoint[0] = Prod(0);
191       prodPoint[1] = Prod(1);
192       prodPoint[2] = Prod(2);
193       
194    }
195    if (dist > 1000) dist = 999.;
196    return dist;
197 }
198
199 ///________________________________________________________________________
200 void AliAODConversionMother::CalculateDistanceOfClossetApproachToPrimVtx(const AliVVertex* primVertex){
201
202    Double_t primCo[3] = {primVertex->GetX(),primVertex->GetY(),primVertex->GetZ()};
203
204    Double_t absoluteP = TMath::Sqrt(TMath::Power(Px(),2) + TMath::Power(Py(),2) + TMath::Power(Pz(),2));   
205    Double_t p[3] = {Px()/absoluteP,Py()/absoluteP,Pz()/absoluteP};
206    Double_t CP[3];
207    
208    CP[0] =  fProductionVtx[0] - primCo[0];
209    CP[1] =  fProductionVtx[1] - primCo[1];
210    CP[2] =  fProductionVtx[2] - primCo[2];
211    
212    Double_t Lambda = - (CP[0]*p[0]+CP[1]*p[1]+CP[2]*p[2])/(p[0]*p[0]+p[1]*p[1]+p[2]*p[2]);
213    
214    Double_t S[3];
215    S[0] = fProductionVtx[0] + p[0]*Lambda;
216    S[1] = fProductionVtx[1] + p[1]*Lambda;
217    S[2] = fProductionVtx[2] + p[2]*Lambda;
218    
219    fdcaRPrimVtx = TMath::Sqrt( TMath::Power(primCo[0]-S[0],2) + TMath::Power(primCo[1]-S[1],2));
220    fdcaZPrimVtx = primCo[2]-S[2];
221    
222    
223 //    cout << "DCA z: " << dca[1] << "\t DCA r: " << dca[0] << "\t DCA 3d: " << TMath::Sqrt(dca[1]*dca[1] + dca[0]*dca[0]) << endl;
224    
225    
226    return;
227 }
228
229 ///________________________________________________________________________
230 void AliAODConversionMother::DetermineMesonQuality(AliAODConversionPhoton* y1, AliAODConversionPhoton* y2){
231    UChar_t photonQA1 = y1->GetPhotonQuality();
232    UChar_t photonQA2 = y2->GetPhotonQuality();
233    
234    if (photonQA1 == 0 || photonQA2 == 0){
235       fQuality = 0;
236       return;
237    } 
238    if (photonQA1 == 1 && photonQA2 == 1){
239       fQuality = 1;
240       return;
241    }               
242    if (photonQA1 == 2 && photonQA2 == 2){
243       fQuality = 4;
244       return;
245    }               
246    if (photonQA1 == 3 && photonQA2 == 3){
247       fQuality = 6;
248       return;
249    }         
250    if (photonQA1 == 1){
251        if (photonQA2 == 2){
252         fQuality = 2;  
253         return;
254        }
255        if (photonQA2 == 3){
256         fQuality = 3;  
257         return;
258        }
259    }
260    if (photonQA2 == 1){
261        if (photonQA1 == 2){
262         fQuality = 2;  
263         return;
264        }
265        if (photonQA1 == 3){
266         fQuality = 3;  
267         return;
268        }
269    }
270    if ((photonQA1 == 2 && photonQA2 == 3)|| (photonQA1 == 3 && photonQA2 == 2)){
271         fQuality = 5;  
272         return;
273    }
274    
275 }