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::UpdateKF(const AliESDVertex &vertex, Int_t pdg0, Int_t pdg1, Float_t mass){
 266   //
 267   // Calculate properties of V0 vertex using different type of constraints
 268   //
 269   fKFrecChi2NC=0;     //  ONLINE V0 finder non constrained chi2
 270   fKFrecChi2C=0;      //  ONLINE V0 finder   constrained chi2 - prim vertex
 271   fKFrecChi2CM=0;     //  ONLINE V0 finder   constrained chi2 - prim vertex+mass
 272   if (fKFRecNC) {delete fKFRecNC; fKFRecNC=0;}
 273   if (fKFRecC)  {delete fKFRecC;  fKFRecC=0;}
 274   if (fKFRecCM) {delete fKFRecCM; fKFRecCM=0;}
 275   //
 276   fKFrecOffChi2NC=0;  // OFFLINE V0 finder - non constrained chi2
 277   fKFrecOffChi2C=0;   // OFFLINE V0 finder -     constrained chi2 - prim vertex
 278   fKFrecOffChi2CM=0;  // OFFLINE V0 finder -     constrained chi2 - prim vertex+mass
 279   if (fKFOffRecNC) {delete fKFOffRecNC; fKFOffRecNC=0;}
 280   if (fKFOffRecC)  {delete fKFOffRecC;  fKFOffRecC=0;}
 281   if (fKFOffRecCM) {delete fKFOffRecCM; fKFOffRecCM=0;}
 282   if (fV0Status==0) return; //
 283   //
 284   AliKFVertex primVtx(vertex);
 285   //
 286   if (fV0rec &&
 287       TMath::Abs(fV0rec->GetParamN()->GetSigmaY2())>0.000000001&&
 288       TMath::Abs(fV0rec->GetParamP()->GetSigmaY2())>0.000000001
 289       ){
 290     //
 291     Double_t x, y, z;
 292     AliKFParticle p1( *(fV0rec->GetParamN()), pdg0 );
 293     AliKFParticle p2( *(fV0rec->GetParamP()), pdg1 );
 294     //
 295     fKFRecNC  = new AliKFParticle;
 296     fV0rec->GetXYZ(x,y,z);
 297     fKFRecNC->SetVtxGuess(x,y,z);
 298     *(fKFRecNC)+=p1;
 299     *(fKFRecNC)+=p2;
 300     fKFrecChi2NC =fKFRecNC->GetChi2() ;
 301     //
 302     fKFRecC  = new AliKFParticle;
 303     fV0rec->GetXYZ(x,y,z);
 304     fKFRecC->SetVtxGuess(x,y,z);
 305     *(fKFRecC)+=p1;
 306     *(fKFRecC)+=p2;
 307     fKFRecC->SetProductionVertex(primVtx);
 308     fKFrecChi2C =fKFRecC->GetChi2();
 309     //
 310     fKFRecCM  = new AliKFParticle;
 311     fV0rec->GetXYZ(x,y,z);
 312     fKFRecCM->SetVtxGuess(x,y,z);
 313     *(fKFRecCM)+=p1;
 314     *(fKFRecCM)+=p2;
 315     fKFRecCM->SetProductionVertex(primVtx);
 316     fKFRecCM->SetMassConstraint(mass);
 317     fKFrecChi2CM =fKFRecCM->GetChi2();
 318   }
 319
 320   if (fV0recOff &&
 321       TMath::Abs(fV0recOff->GetParamN()->GetSigmaY2())>0.000000001&&
 322       TMath::Abs(fV0recOff->GetParamP()->GetSigmaY2())>0.000000001
 323       ){
 324     //
 325     Double_t x, y, z;
 326     AliKFParticle p1( *(fV0recOff->GetParamN()), pdg0 );
 327     AliKFParticle p2( *(fV0recOff->GetParamP()), pdg1 );
 328     //
 329     fKFOffRecNC  = new AliKFParticle;
 330     fV0recOff->GetXYZ(x,y,z);
 331     fKFOffRecNC->SetVtxGuess(x,y,z);
 332     *(fKFOffRecNC)+=p1;
 333     *(fKFOffRecNC)+=p2;
 334     fKFrecOffChi2NC =fKFOffRecNC->GetChi2() ;
 335     //
 336     fKFOffRecC  = new AliKFParticle;
 337     fV0recOff->GetXYZ(x,y,z);
 338     fKFOffRecC->SetVtxGuess(x,y,z);
 339     *(fKFOffRecC)+=p1;
 340     *(fKFOffRecC)+=p2;
 341     fKFOffRecC->SetProductionVertex(primVtx);
 342     fKFrecOffChi2C =fKFOffRecC->GetChi2();
 343     //
 344     fKFOffRecCM  = new AliKFParticle;
 345     fV0recOff->GetXYZ(x,y,z);
 346     fKFOffRecCM->SetVtxGuess(x,y,z);
 347     *(fKFOffRecCM)+=p1;
 348     *(fKFOffRecCM)+=p2;
 349     fKFOffRecCM->SetProductionVertex(primVtx);
 350     fKFOffRecCM->SetMassConstraint(mass);
 351     fKFrecOffChi2CM =fKFOffRecCM->GetChi2();
 352   }
 353
 354
 355 }
 356
 357