PWG1/TPC/AliGenKinkInfo.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
  22
  23 */
  24
  25 #if !defined(__CINT__) || defined(__MAKECINT__)
  26 #include <stdio.h>
  27 #include <string.h>
  28 //ROOT includes
  29 #include "TROOT.h"
  30 #include "Rtypes.h"
  31 #include "TFile.h"
  32 #include "TTree.h"
  33 #include "TChain.h"
  34 #include "TCut.h"
  35 #include "TString.h"
  36 #include "TStopwatch.h"
  37 #include "TParticle.h"
  38 #include "TDatabasePDG.h"
  39 #include "TSystem.h"
  40 #include "TCanvas.h"
  41 #include "TPolyLine3D.h"
  42
  43 //ALIROOT includes
  44 #include "AliRun.h"
  45 #include "AliStack.h"
  46 #include "AliSimDigits.h"
  47 #include "AliTPCParam.h"
  48 #include "AliTPC.h"
  49 #include "AliTPCLoader.h"
  50 #include "AliDetector.h"
  51 #include "AliTrackReference.h"
  52 #include "AliTPCParamSR.h"
  53 #include "AliTracker.h"
  54 #include "AliMagF.h"
  55 #include "AliHelix.h"
  56 #include "AliTrackPointArray.h"
  57
  58 #endif
  59 #include "AliGenKinkInfo.h"
  60
  61 //
  62 //
  63
  64 ClassImp(AliGenKinkInfo)
  65
  66
  67
  68 /////////////////////////////////////////////////////////////////////////////////
  69
  70 AliGenKinkInfo::AliGenKinkInfo():
  71   fMCd(),       //info about daughter particle - second particle for V0
  72   fMCm(),       //info about mother particle   - first particle for V0
  73   fMCDist1(0),  //info about closest distance according closest MC - linear DCA
  74   fMCDist2(0),  //info about closest distance parabolic DCA
  75   fMCRr(0),     // rec position of the vertex
  76   fMCR(0)       //exact r position of the vertex
  77 {
  78   //
  79   // default constructor
  80   //
  81   for (Int_t i=0;i<3;i++){
  82     fMCPdr[i]=0;
  83     fMCPd[i]=0;
  84     fMCX[i]=0;
  85     fMCPm[i]=0;
  86     fMCAngle[i]=0;
  87     fMCXr[i] = 0;
  88   }
  89   for (Int_t i=0; i<2; i++) {
  90     fPdg[i]= 0;
  91     fLab[i]=0;
  92   }
  93 }
  94
  95 void AliGenKinkInfo::Update()
  96 {
  97   //
  98   // Update information
  99   //  some redundancy - faster acces to the values in analysis code
 100   //
 101   fMCPd[0] = fMCd.GetParticle().Px();
 102   fMCPd[1] = fMCd.GetParticle().Py();
 103   fMCPd[2] = fMCd.GetParticle().Pz();
 104   fMCPd[3] = fMCd.GetParticle().P();
 105   //
 106   fMCX[0]  = fMCd.GetParticle().Vx();
 107   fMCX[1]  = fMCd.GetParticle().Vy();
 108   fMCX[2]  = fMCd.GetParticle().Vz();
 109   fMCR       = TMath::Sqrt( fMCX[0]*fMCX[0]+fMCX[1]*fMCX[1]);
 110   //
 111   fPdg[0]    = fMCd.GetParticle().GetPdgCode();
 112   fPdg[1]    = fMCm.GetParticle().GetPdgCode();
 113   //
 114   fLab[0]    = fMCd.GetParticle().GetUniqueID();
 115   fLab[1]    = fMCm.GetParticle().GetUniqueID();
 116   //
 117   //
 118   //
 119   Double_t x1[3],p1[3];
 120   TParticle& pdaughter = fMCd.GetParticle();
 121   x1[0] = pdaughter.Vx();
 122   x1[1] = pdaughter.Vy();
 123   x1[2] = pdaughter.Vz();
 124   p1[0] = pdaughter.Px();
 125   p1[1] = pdaughter.Py();
 126   p1[2] = pdaughter.Pz();
 127   Double_t sign = 0;
 128   sign = (pdaughter.GetPDG()->Charge()>0)? -1:1;
 129   AliHelix dhelix1(x1,p1,sign);
 130   //
 131   //
 132   Double_t x2[3],p2[3];
 133   //
 134   TParticle& pmother = fMCm.GetParticle();
 135   x2[0] = pmother.Vx();
 136   x2[1] = pmother.Vy();
 137   x2[2] = pmother.Vz();
 138   p2[0] = pmother.Px();
 139   p2[1] = pmother.Py();
 140   p2[2] = pmother.Pz();
 141   //
 142   const AliTrackReference & pdecay = fMCm.GetTRdecay();
 143   x2[0] = pdecay.X();
 144   x2[1] = pdecay.Y();
 145   x2[2] = pdecay.Z();
 146   p2[0] = pdecay.Px();
 147   p2[1] = pdecay.Py();
 148   p2[2] = pdecay.Pz();
 149   //
 150   sign = (pmother.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   if (points==2){
 169     dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
 170     dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
 171     dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
 172   }
 173   distance1 = TMath::Min(delta1,delta2);
 174   //
 175   //find intersection parabolic
 176   //
 177   points = dhelix1.GetRPHIintersections(mhelix, phase, radius);
 178   delta1=10000,delta2=10000;
 179
 180   if (points>0){
 181     dhelix1.ParabolicDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
 182   }
 183   if (points==2){
 184     dhelix1.ParabolicDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
 185   }
 186
 187   distance2 = TMath::Min(delta1,delta2);
 188   //
 189   if (delta1<delta2){
 190     //get V0 info
 191     dhelix1.Evaluate(phase[0][0],fMCXr);
 192     dhelix1.GetMomentum(phase[0][0],fMCPdr);
 193     mhelix.GetMomentum(phase[0][1],fMCPm);
 194     dhelix1.GetAngle(phase[0][0],mhelix,phase[0][1],fMCAngle);
 195     fMCRr = TMath::Sqrt(radius[0]);
 196   }
 197   else{
 198     dhelix1.Evaluate(phase[1][0],fMCXr);
 199     dhelix1.GetMomentum(phase[1][0],fMCPdr);
 200     mhelix.GetMomentum(phase[1][1],fMCPm);
 201     dhelix1.GetAngle(phase[1][0],mhelix,phase[1][1],fMCAngle);
 202     fMCRr = TMath::Sqrt(radius[1]);
 203   }
 204   //
 205   //
 206   fMCDist1 = TMath::Sqrt(distance1);
 207   fMCDist2 = TMath::Sqrt(distance2);
 208
 209 }
 210
 211
 212 Float_t AliGenKinkInfo::GetQt(){
 213   //
 214   //
 215   Float_t momentumd = TMath::Sqrt(fMCPd[0]*fMCPd[0]+fMCPd[1]*fMCPd[1]+fMCPd[2]*fMCPd[2]);
 216   return TMath::Sin(fMCAngle[2])*momentumd;
 217 }
 218
 219