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