remove obselete clean QA macros command
[u/mrichter/AliRoot.git] / PWG1 / AliESDRecV0Info.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 //  Comparison class for V0 information                                      //
21 //  responsible: 
22 //  marian.ivanov@cern.ch                                                    //
23 //
24 //
25
26  
27
28
29
30 #include <stdio.h>
31 #include <string.h>
32 //ROOT includes
33 #include "Rtypes.h"
34 //
35 //ALIROOT includes
36 //
37 #include "AliESDtrack.h"
38 #include "AliTPCParam.h"
39 #include "AliTrackReference.h"
40 #include "AliTPCParamSR.h"
41 #include "AliESD.h"
42 #include "AliESDfriend.h"
43 #include "AliESDtrack.h"
44 #include "AliTPCseed.h"
45 #include "AliITStrackMI.h"
46 #include "AliTRDtrack.h"
47 #include "AliHelix.h"
48 #include "AliESDVertex.h"
49 #include "AliExternalTrackParam.h"
50 #include "AliESDkink.h"
51 #include "AliESDv0.h"
52 #include "AliV0.h"
53 #include "AliKFParticle.h"
54 #include "AliKFVertex.h"
55 //
56 #include "AliTreeDraw.h"
57 #include "AliMCInfo.h"
58 #include "AliGenKinkInfo.h"
59 #include "AliGenV0Info.h"
60
61
62 #include "AliESDRecV0Info.h"
63
64
65
66 ClassImp(AliESDRecV0Info)
67
68
69 AliESDRecV0Info:: AliESDRecV0Info():
70   TObject(),
71   fT1(),               //track1
72   fT2(),               //track2  
73   fDist1(0),           //info about closest distance according closest MC - linear DCA
74   fDist2(0),           //info about closest distance parabolic DCA
75   fInvMass(0),         //reconstructed invariant mass -
76   //
77   fDistMinR(0),        // distance at minimal radius
78   fRr(0),              // rec position of the vertex 
79   fPointAngleFi(0),    //point angle fi
80   fPointAngleTh(0),    //point angle theta
81   fPointAngle(0),      //point angle full
82   fV0Status(0),        // status of the kink
83   fV0tpc(0),           // Vo information from reconsturction according TPC
84   fV0its(0),           // Vo information from reconsturction according ITS
85   fV0rec(0),           // V0 information form the reconstruction
86   fV0recOff(0),        // V0 information form the reconstruction - OFFLINE
87   fMultiple(0),        // how man times V0 was recostructed 
88   fRecStatus(0),       // status form the reconstuction
89   fV0MultipleOn(0),    // how man times was V0 reconstucted
90   fV0MultipleOff(0),   // how man times was V0 reconstucted
91   //
92   fKFrecChi2NC(0),     //  ONLINE V0 finder non constrained chi2  
93   fKFrecChi2C(0),      //  ONLINE V0 finder   constrained chi2 - prim vertex  
94   fKFrecChi2CM(0),     //  ONLINE V0 finder   constrained chi2 - prim vertex+mass 
95   fKFRecNC(0),         //  non constrained  
96   fKFRecC(0),          //  constrained vertex
97   fKFRecCM(0),         //  constrained vertex+mass
98   fKFrecOffChi2NC(0),  // OFFLINE V0 finder - non constrained chi2  
99   fKFrecOffChi2C(0),   // OFFLINE V0 finder -     constrained chi2 - prim vertex  
100   fKFrecOffChi2CM(0),  // OFFLINE V0 finder -     constrained chi2 - prim vertex+mass
101   fKFOffRecNC(0),      //  non constrained  
102   fKFOffRecC(0),       //  constrained vertex
103   fKFOffRecCM(0)       //  constrained vertex+mass
104 {
105   //
106   // default constructor
107   //
108   fV0tpc    = new AliV0();
109   fV0its    = new AliV0();
110   fV0rec    = new AliV0();
111   fV0recOff = new AliV0();
112 }
113
114
115 void  AliESDRecV0Info::Update(Float_t vertex[3])
116
117
118   if ( (fT1.fStatus[1]>0)&& (fT2.fStatus[1]>0)){
119     Float_t distance1,distance2;
120     Double_t xx[3],pp[3];
121     //
122     Double_t xd[3],pd[3],signd;
123     Double_t xm[3],pm[3],signm;
124     //
125     //
126     if (fT1.fITSOn&&fT2.fITSOn){
127       for (Int_t i=0;i<3;i++){
128         xd[i] = fT2.fITSinR1[i];
129         pd[i] = fT2.fITSinP1[i];
130         xm[i] = fT1.fITSinR1[i];
131         pm[i] = fT1.fITSinP1[i];
132       }
133     }
134     else{
135       
136       for (Int_t i=0;i<3;i++){
137         xd[i] = fT2.fTPCinR1[i];
138         pd[i] = fT2.fTPCinP1[i];
139         xm[i] = fT1.fTPCinR1[i];
140         pm[i] = fT1.fTPCinP1[i];
141       }
142     }
143     //
144     //
145     signd =  fT2.fSign<0 ? -1:1;
146     signm =  fT1.fSign<0 ? -1:1;
147
148     AliHelix dhelix1(xd,pd,signd);
149     dhelix1.GetMomentum(0,pp,0);
150     dhelix1.Evaluate(0,xx);      
151     // 
152     //  Double_t x2[3],p2[3];
153     //            
154     AliHelix mhelix(xm,pm,signm);    
155     //
156     //find intersection linear
157     //
158     Double_t phase[2][2],radius[2];
159     Int_t  points = dhelix1.GetRPHIintersections(mhelix, phase, radius,200);
160     Double_t delta1=10000,delta2=10000;  
161
162     if (points==1){
163       fRs[0] = TMath::Sqrt(radius[0]);
164       fRs[1] = TMath::Sqrt(radius[0]);
165     }
166     if (points==2){
167       fRs[0] =TMath::Min(TMath::Sqrt(radius[0]),TMath::Sqrt(radius[1]));
168       fRs[1] =TMath::Max(TMath::Sqrt(radius[0]),TMath::Sqrt(radius[1]));
169     }
170     
171     if (points>0){
172       dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
173       dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
174       dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
175     }
176     if (points==2){    
177       dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
178       dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
179       dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
180     }
181     if (points==1){
182       fRs[0] = TMath::Sqrt(radius[0]);
183       fRs[1] = TMath::Sqrt(radius[0]);
184       fDistMinR = delta1;
185     }
186     if (points==2){
187       if (radius[0]<radius[1]){
188         fRs[0] = TMath::Sqrt(radius[0]);
189         fRs[1] = TMath::Sqrt(radius[1]);
190         fDistMinR = delta1;
191       }
192       else{
193         fRs[0] = TMath::Sqrt(radius[1]);
194         fRs[1] = TMath::Sqrt(radius[0]);
195         fDistMinR = delta2;
196       }
197     }
198     //
199     //
200     distance1 = TMath::Min(delta1,delta2);
201     //
202     //find intersection parabolic
203     //
204     points = dhelix1.GetRPHIintersections(mhelix, phase, radius);
205     delta1=10000,delta2=10000;  
206     
207     if (points>0){
208       dhelix1.ParabolicDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
209     }
210     if (points==2){    
211       dhelix1.ParabolicDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
212     }
213     
214     distance2 = TMath::Min(delta1,delta2);
215     if (distance2>100) fDist2 =100;
216     return;
217     if (delta1<delta2){
218       //get V0 info
219       dhelix1.Evaluate(phase[0][0],fXr);
220       dhelix1.GetMomentum(phase[0][0],fPdr);
221       mhelix.GetMomentum(phase[0][1],fPm);
222       dhelix1.GetAngle(phase[0][0],mhelix,phase[0][1],fAngle);
223       fRr = TMath::Sqrt(radius[0]);
224     }
225     else{
226       dhelix1.Evaluate(phase[1][0],fXr);
227       dhelix1.GetMomentum(phase[1][0], fPdr);
228       mhelix.GetMomentum(phase[1][1], fPm);
229       dhelix1.GetAngle(phase[1][0],mhelix,phase[1][1],fAngle);
230       fRr = TMath::Sqrt(radius[1]);
231     }
232     fDist1 = TMath::Sqrt(distance1);
233     fDist2 = TMath::Sqrt(distance2);      
234     
235     if (fDist2<10.5){
236       Double_t x,alpha,param[5],cov[15];
237       //
238       fT1.GetESDtrack()->GetInnerExternalParameters(alpha,x,param);
239       fT1.GetESDtrack()->GetInnerExternalCovariance(cov);
240       AliExternalTrackParam paramm(x,alpha,param,cov);
241       //
242       fT2.GetESDtrack()->GetInnerExternalParameters(alpha,x,param);
243       fT2.GetESDtrack()->GetInnerExternalCovariance(cov);
244       AliExternalTrackParam paramd(x,alpha,param,cov);
245     }    
246     //            
247     //   
248     
249     Float_t v[3] = {fXr[0]-vertex[0],fXr[1]-vertex[1],fXr[2]-vertex[2]};
250     Float_t p[3] = {fPdr[0]+fPm[0], fPdr[1]+fPm[1],fPdr[2]+fPm[2]};
251     
252     Float_t vnorm2 = v[0]*v[0]+v[1]*v[1];
253     Float_t vnorm3 = TMath::Sqrt(v[2]*v[2]+vnorm2);
254     vnorm2 = TMath::Sqrt(vnorm2);
255     Float_t pnorm2 = p[0]*p[0]+p[1]*p[1];
256     Float_t pnorm3 = TMath::Sqrt(p[2]*p[2]+pnorm2);
257     pnorm2 = TMath::Sqrt(pnorm2);
258     
259     fPointAngleFi = (v[0]*p[0]+v[1]*p[1])/(vnorm2*pnorm2);
260     fPointAngleTh = (v[2]*p[2]+vnorm2*pnorm2)/(vnorm3*pnorm3);  
261     fPointAngle   = (v[0]*p[0]+v[1]*p[1]+v[2]*p[2])/(vnorm3*pnorm3);
262   }
263 }
264
265 void  AliESDRecV0Info::Reset(){
266   //
267   // Reset status and all counters
268   //
269   fDist1=-1;    //info about closest distance according closest MC - linear DCA
270   fDist2=-1;    //info about closest distance parabolic DCA
271   fInvMass=-1;  //reconstructed invariant mass -
272   //
273   fDistMinR=-1; // distance at minimal radius
274   fRr=-1;       // rec position of the vertex 
275   fLab[0]=-20;   //MC label of the partecle
276   fLab[1]=-10;   //MC label of the partecle
277   fPointAngleFi=0; //point angle fi
278   fPointAngleTh=0; //point angle theta
279   fPointAngle=0;   //point angle full
280   //
281   fV0Status= -100;       // status of the V0
282   fMultiple=0;           // how man times V0 was recostructed 
283   fRecStatus=0;          // status form the reconstuction - 1 reconstructed - -1 fake
284   fV0MultipleOn=0;       // how man times was V0 reconstucted - onfly
285   fV0MultipleOff=0;      // how man times was V0 reconstucted - offline
286   //
287   // AliKF variables - variables to make a selection + resoluton study
288   //
289   fKFrecChi2NC=0;     //  ONLINE V0 finder non constrained chi2  
290   fKFrecChi2C=0;      //  ONLINE V0 finder   constrained chi2 - prim vertex  
291   fKFrecChi2CM=0;     //  ONLINE V0 finder   constrained chi2 - prim vertex+mass 
292   //
293   fKFrecOffChi2NC=0;  // OFFLINE V0 finder - non constrained chi2  
294   fKFrecOffChi2C=0;   // OFFLINE V0 finder -     constrained chi2 - prim vertex  
295   fKFrecOffChi2CM=0;  // OFFLINE V0 finder -     constrained chi2 - prim vertex+mass
296 }
297
298
299
300 void  AliESDRecV0Info::UpdateKF(const AliESDVertex &vertex, Int_t pdg0, Int_t pdg1, Float_t mass){
301   //
302   // Calculate properties of V0 vertex using different type of constraints 
303   //
304   fKFrecChi2NC=0;     //  ONLINE V0 finder non constrained chi2  
305   fKFrecChi2C=0;      //  ONLINE V0 finder   constrained chi2 - prim vertex  
306   fKFrecChi2CM=0;     //  ONLINE V0 finder   constrained chi2 - prim vertex+mass 
307   if (fKFRecNC) {delete fKFRecNC; fKFRecNC=0;}
308   if (fKFRecC)  {delete fKFRecC;  fKFRecC=0;}
309   if (fKFRecCM) {delete fKFRecCM; fKFRecCM=0;}
310   //
311   fKFrecOffChi2NC=0;  // OFFLINE V0 finder - non constrained chi2  
312   fKFrecOffChi2C=0;   // OFFLINE V0 finder -     constrained chi2 - prim vertex  
313   fKFrecOffChi2CM=0;  // OFFLINE V0 finder -     constrained chi2 - prim vertex+mass
314   if (fKFOffRecNC) {delete fKFOffRecNC; fKFOffRecNC=0;}
315   if (fKFOffRecC)  {delete fKFOffRecC;  fKFOffRecC=0;}
316   if (fKFOffRecCM) {delete fKFOffRecCM; fKFOffRecCM=0;}
317   if (fV0Status==0) return; //
318   //
319   AliKFVertex primVtx(vertex);
320   //
321   if (fV0rec && 
322       TMath::Abs(fV0rec->GetParamN()->GetSigmaY2())>0.000000001&&
323       TMath::Abs(fV0rec->GetParamP()->GetSigmaY2())>0.000000001
324       ){
325     //
326     Double_t x, y, z;
327     AliKFParticle p1( *(fV0rec->GetParamN()), pdg0 );
328     AliKFParticle p2( *(fV0rec->GetParamP()), pdg1 );
329     //
330     fKFRecNC  = new AliKFParticle;
331     fV0rec->GetXYZ(x,y,z);
332     fKFRecNC->SetVtxGuess(x,y,z);
333     *(fKFRecNC)+=p1;
334     *(fKFRecNC)+=p2;
335     fKFrecChi2NC =fKFRecNC->GetChi2() ;
336     //
337     fKFRecC  = new AliKFParticle;
338     fV0rec->GetXYZ(x,y,z);
339     fKFRecC->SetVtxGuess(x,y,z);
340     *(fKFRecC)+=p1;
341     *(fKFRecC)+=p2;
342     fKFRecC->SetProductionVertex(primVtx);
343     fKFrecChi2C =fKFRecC->GetChi2();
344     //
345     fKFRecCM  = new AliKFParticle;
346     fV0rec->GetXYZ(x,y,z);
347     fKFRecCM->SetVtxGuess(x,y,z);
348     *(fKFRecCM)+=p1;
349     *(fKFRecCM)+=p2;
350     fKFRecCM->SetProductionVertex(primVtx);
351     fKFRecCM->SetMassConstraint(mass);
352     fKFrecChi2CM =fKFRecCM->GetChi2();    
353   }
354
355   if (fV0recOff && 
356       TMath::Abs(fV0recOff->GetParamN()->GetSigmaY2())>0.000000001&&
357       TMath::Abs(fV0recOff->GetParamP()->GetSigmaY2())>0.000000001
358       ){
359     //
360     Double_t x, y, z;
361     AliKFParticle p1( *(fV0recOff->GetParamN()), pdg0 );
362     AliKFParticle p2( *(fV0recOff->GetParamP()), pdg1 );
363     //
364     fKFOffRecNC  = new AliKFParticle;
365     fV0recOff->GetXYZ(x,y,z);
366     fKFOffRecNC->SetVtxGuess(x,y,z);
367     *(fKFOffRecNC)+=p1;
368     *(fKFOffRecNC)+=p2;
369     fKFrecOffChi2NC =fKFOffRecNC->GetChi2() ;
370     //
371     fKFOffRecC  = new AliKFParticle;
372     fV0recOff->GetXYZ(x,y,z);
373     fKFOffRecC->SetVtxGuess(x,y,z);
374     *(fKFOffRecC)+=p1;
375     *(fKFOffRecC)+=p2;
376     fKFOffRecC->SetProductionVertex(primVtx);
377     fKFrecOffChi2C =fKFOffRecC->GetChi2();
378     //
379     fKFOffRecCM  = new AliKFParticle;
380     fV0recOff->GetXYZ(x,y,z);
381     fKFOffRecCM->SetVtxGuess(x,y,z);
382     *(fKFOffRecCM)+=p1;
383     *(fKFOffRecCM)+=p2;
384     fKFOffRecCM->SetProductionVertex(primVtx);
385     fKFOffRecCM->SetMassConstraint(mass);
386     fKFrecOffChi2CM =fKFOffRecCM->GetChi2();    
387   }
388
389
390 }
391
392