PWGPP/TPC/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] = {{0,0},{0,0}};
 160   Double_t radius[2] = {0};
 161   Int_t  points = dhelix1.GetRPHIintersections(mhelix, phase, radius);
 162   Double_t delta1=10000,delta2=10000;
 163   if (points>0){
 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     dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
 167   }
 168   else{
 169     fInvMass=-1;
 170     return;
 171   }
 172   if (points==2){
 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     dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
 176   }
 177   distance1 = TMath::Min(delta1,delta2);
 178   //
 179   //find intersection parabolic
 180   //
 181   points = dhelix1.GetRPHIintersections(mhelix, phase, radius);
 182   delta1=10000,delta2=10000;
 183
 184   if (points>0){
 185     dhelix1.ParabolicDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
 186   }
 187   if (points==2){
 188     dhelix1.ParabolicDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
 189   }
 190
 191   distance2 = TMath::Min(delta1,delta2);
 192   //
 193   if (delta1<delta2){
 194     //get V0 info
 195     dhelix1.Evaluate(phase[0][0],fMCXr);
 196     dhelix1.GetMomentum(phase[0][0],fMCPdr);
 197     mhelix.GetMomentum(phase[0][1],fMCPm);
 198     dhelix1.GetAngle(phase[0][0],mhelix,phase[0][1],fMCAngle);
 199     fMCRr = TMath::Sqrt(radius[0]);
 200   }
 201   else{
 202     dhelix1.Evaluate(phase[1][0],fMCXr);
 203     dhelix1.GetMomentum(phase[1][0],fMCPdr);
 204     mhelix.GetMomentum(phase[1][1],fMCPm);
 205     dhelix1.GetAngle(phase[1][0],mhelix,phase[1][1],fMCAngle);
 206     fMCRr = TMath::Sqrt(radius[1]);
 207   }
 208   //
 209   //
 210   fMCDist1 = TMath::Sqrt(distance1);
 211   fMCDist2 = TMath::Sqrt(distance2);
 212   //
 213   //
 214   //
 215   Float_t v[3] = {fMCXr[0]-vertex[0],fMCXr[1]-vertex[1],fMCXr[2]-vertex[2]};
 216   //Float_t v[3] = {fMCXr[0],fMCXr[1],fMCXr[2]};
 217   Float_t p[3] = {fMCPdr[0]+fMCPm[0], fMCPdr[1]+fMCPm[1],fMCPdr[2]+fMCPm[2]};
 218   Float_t vnorm2 = v[0]*v[0]+v[1]*v[1];
 219   Float_t vnorm3 = TMath::Sqrt(v[2]*v[2]+vnorm2);
 220   vnorm2 = TMath::Sqrt(vnorm2);
 221   Float_t pnorm2 = p[0]*p[0]+p[1]*p[1];
 222   Float_t pnorm3 = TMath::Sqrt(p[2]*p[2]+pnorm2);
 223   pnorm2 = TMath::Sqrt(pnorm2);
 224   //
 225   if (vnorm2>0){
 226     fPointAngleFi = (v[0]*p[0]+v[1]*p[1])/(vnorm2*pnorm2);
 227     fPointAngleTh = (v[2]*p[2]+vnorm2*pnorm2)/(vnorm3*pnorm3);
 228     fPointAngle   = (v[0]*p[0]+v[1]*p[1]+v[2]*p[2])/(vnorm3*pnorm3);
 229   }else{
 230     fPointAngleFi = 1;
 231     fPointAngleTh = 1;
 232     fPointAngle   = 1;
 233   }
 234   Double_t mass1 = fMCd.GetMass();
 235   Double_t mass2 = fMCm.GetMass();
 236
 237
 238   Double_t e1    = TMath::Sqrt(mass1*mass1+
 239                               fMCPd[0]*fMCPd[0]+
 240                               fMCPd[1]*fMCPd[1]+
 241                               fMCPd[2]*fMCPd[2]);
 242   Double_t e2    = TMath::Sqrt(mass2*mass2+
 243                               fMCPm[0]*fMCPm[0]+
 244                               fMCPm[1]*fMCPm[1]+
 245                               fMCPm[2]*fMCPm[2]);
 246
 247   fInvMass =
 248     (fMCPm[0]+fMCPd[0])*(fMCPm[0]+fMCPd[0])+
 249     (fMCPm[1]+fMCPd[1])*(fMCPm[1]+fMCPd[1])+
 250     (fMCPm[2]+fMCPd[2])*(fMCPm[2]+fMCPd[2]);
 251
 252   //  fInvMass = TMath::Sqrt((e1+e2)*(e1+e2)-fInvMass);
 253   fInvMass = (e1+e2)*(e1+e2)-fInvMass;
 254   if (fInvMass>0) fInvMass = TMath::Sqrt(fInvMass);
 255 }
 256
 257
 258
 259
 260
 261
 262
 263