PWG1/AliMCInfo.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 The AliMCInfo contains the information about the particles properties
  24 during transportation throuch ALICE Detector
  25
  26 The base Information :
  27 TParticle - fParticle             -  properties of the particle at creation point
  28 AliTrackReference - fXXXRefernces -  TClonesArray of refernces in differnt detectors
  29 fNXXXRef                          -  number of the track refernces in differnt detectors
  30 AliTPCdigitRow    - fTPCRow       -  the map of the hitted rows - (will be repalced by TBits)
  31 fRowsWith*                        - number of rows hitted by particle
  32 fMCtracks         -               - number of turn over of the track inside of the TPC
  33
  34 ++++
  35 some additional information usable for tree draw - TO SPEED UP tree queries
  36 IMPORTANT FOR PROOF FAST PROTOTYPING ANALYSIS
  37
  38
  39
  40
  41 */
  42
  43 #if !defined(__CINT__) || defined(__MAKECINT__)
  44 #include <stdio.h>
  45 //ROOT includes
  46 #include "Rtypes.h"
  47 #include "TClonesArray.h"
  48 //ALIROOT includes
  49 #include "AliTrackReference.h"
  50 #include "AliMCInfo.h"
  51 #endif
  52
  53 //
  54 //
  55
  56 ClassImp(AliTPCdigitRow)
  57 ClassImp(AliMCInfo)
  58
  59
  60
  61 ////////////////////////////////////////////////////////////////////////
  62 AliMCInfo::AliMCInfo():
  63   fTrackRef(),
  64   fTrackRefOut(),
  65   fTRdecay(),
  66   fPrimPart(0),
  67   fParticle(),
  68   fMass(0),
  69   fCharge(0),
  70   fLabel(0),
  71   fEventNr(),
  72   fMCtracks(),
  73   fPdg(0),
  74   fTPCdecay(0),
  75   fRowsWithDigitsInn(0),
  76   fRowsWithDigits(0),
  77   fRowsTrackLength(0),
  78   fPrim(0),
  79   fTPCRow(),
  80   fNTPCRef(0),                    // tpc references counter
  81   fNITSRef(0),                    // ITS references counter
  82   fNTRDRef(0),                    // TRD references counter
  83   fNTOFRef(0),                    // TOF references counter
  84   fTPCReferences(0),
  85   fITSReferences(0),
  86   fTRDReferences(0),
  87   fTOFReferences(0)
  88 {
  89   //
  90   // Default constructor
  91   //
  92   fTPCReferences  = new TClonesArray("AliTrackReference",10);
  93   fITSReferences  = new TClonesArray("AliTrackReference",10);
  94   fTRDReferences  = new TClonesArray("AliTrackReference",10);
  95   fTOFReferences  = new TClonesArray("AliTrackReference",10);
  96   fTRdecay.SetTrack(-1);
  97   fCharge = 0;
  98 }
  99
 100 AliMCInfo::AliMCInfo(const AliMCInfo& info):
 101   TObject(),
 102   fTrackRef(info.fTrackRef),
 103   fTrackRefOut(info.fTrackRefOut),
 104   fTRdecay(info.GetTRdecay()),
 105   fPrimPart(info.fPrimPart),
 106   fParticle(info.fParticle),
 107   fMass(info.GetMass()),
 108   fCharge(info.fCharge),
 109   fLabel(info.fLabel),
 110   fEventNr(info.fEventNr),
 111   fMCtracks(info.fMCtracks),
 112   fPdg(info.fPdg),
 113   fTPCdecay(info.fTPCdecay),
 114   fRowsWithDigitsInn(info.fRowsWithDigitsInn),
 115   fRowsWithDigits(info.fRowsWithDigits),
 116   fRowsTrackLength(info.fRowsTrackLength),
 117   fPrim(info.fPrim),
 118   fTPCRow(info.fTPCRow),
 119   fNTPCRef(info.fNTPCRef),                    // tpc references counter
 120   fNITSRef(info.fNITSRef),                    // ITS references counter
 121   fNTRDRef(info.fNTRDRef),                    // TRD references counter
 122   fNTOFRef(info.fNTOFRef),                    // TOF references counter
 123   fTPCReferences(0),
 124   fITSReferences(0),
 125   fTRDReferences(0),
 126   fTOFReferences(0)
 127 {
 128   //
 129   // copy constructor
 130   //
 131   fTPCReferences = (TClonesArray*)info.fTPCReferences->Clone();
 132   fITSReferences = (TClonesArray*)info.fITSReferences->Clone();
 133   fTRDReferences = (TClonesArray*)info.fTRDReferences->Clone();
 134   fTOFReferences = (TClonesArray*)info.fTOFReferences->Clone();
 135 }
 136
 137
 138 AliMCInfo::~AliMCInfo()
 139 {
 140   //
 141   // Destructor of the class
 142   //
 143   if (fTPCReferences) {
 144     delete fTPCReferences;
 145   }
 146   if (fITSReferences){
 147     delete fITSReferences;
 148   }
 149   if (fTRDReferences){
 150     delete fTRDReferences;
 151   }
 152   if (fTOFReferences){
 153     delete fTOFReferences;
 154   }
 155
 156 }
 157
 158
 159
 160 void AliMCInfo::Update()
 161 {
 162   //
 163   // Update MC info
 164   // Calculates some derived variables
 165   //
 166   //
 167   fMCtracks =1;
 168   if (!fTPCReferences) {
 169     fNTPCRef =0;
 170     return;
 171   }
 172   Float_t direction=1;
 173   //Float_t rlast=0;
 174   fNTPCRef = fTPCReferences->GetEntriesFast();
 175   fNITSRef = fITSReferences->GetEntriesFast();
 176   fNTRDRef = fTRDReferences->GetEntriesFast();
 177   fNTOFRef = fTOFReferences->GetEntriesFast();
 178
 179   for (Int_t iref =0;iref<fTPCReferences->GetEntriesFast();iref++){
 180     AliTrackReference * ref = (AliTrackReference *) fTPCReferences->At(iref);
 181     //Float_t r = (ref->X()*ref->X()+ref->Y()*ref->Y());
 182     Float_t newdirection = ref->X()*ref->Px()+ref->Y()*ref->Py(); //inside or outside
 183     if (iref==0) direction = newdirection;
 184     if ( newdirection*direction<0){
 185       //changed direction
 186       direction = newdirection;
 187       fMCtracks+=1;
 188     }
 189     //rlast=r;
 190   }
 191   //
 192   // decay info
 193   fTPCdecay=kFALSE;
 194   if (fTRdecay.GetTrack()>0){
 195     fDecayCoord[0] = fTRdecay.X();
 196     fDecayCoord[1] = fTRdecay.Y();
 197     fDecayCoord[2] = fTRdecay.Z();
 198     if ( (fTRdecay.R()<250)&&(fTRdecay.R()>85) && (TMath::Abs(fTRdecay.Z())<250) ){
 199       fTPCdecay=kTRUE;
 200     }
 201     else{
 202       fDecayCoord[0] = 0;
 203       fDecayCoord[1] = 0;
 204       fDecayCoord[2] = 0;
 205     }
 206   }
 207   //
 208   //
 209   //digits information update
 210   fRowsWithDigits    = fTPCRow.RowsOn();
 211   fRowsWithDigitsInn = fTPCRow.RowsOn(63); // 63 = number of inner rows
 212   fRowsTrackLength   = fTPCRow.Last() - fTPCRow.First();
 213   //
 214   //
 215   // calculate primary ionization per cm
 216   if (fParticle.GetPDG()){
 217     fMass = fParticle.GetMass();
 218     fCharge = fParticle.GetPDG()->Charge();
 219     if (fTPCReferences->GetEntriesFast()>0){
 220       fTrackRef = *((AliTrackReference*)fTPCReferences->At(0));
 221     }
 222     if (fMass>0){
 223       Float_t p = TMath::Sqrt(fTrackRef.Px()*fTrackRef.Px()+
 224                               fTrackRef.Py()*fTrackRef.Py()+
 225                               fTrackRef.Pz()*fTrackRef.Pz());
 226       if (p>0.001){
 227         Float_t betagama = p /fMass;
 228         fPrim = TPCBetheBloch(betagama);
 229       }else fPrim=0;
 230     }
 231   }else{
 232     fMass =0;
 233     fPrim =0;
 234   }
 235 }
 236
 237
 238
 239
 240 ////////////////////////////////////////////////////////////////////////
 241 AliTPCdigitRow::AliTPCdigitRow()
 242 {
 243   Reset();
 244 }
 245 ////////////////////////////////////////////////////////////////////////
 246 AliTPCdigitRow & AliTPCdigitRow::operator=(const AliTPCdigitRow &digOld)
 247 {
 248   for (Int_t i = 0; i<32; i++) fDig[i] = digOld.fDig[i];
 249   return (*this);
 250 }
 251 ////////////////////////////////////////////////////////////////////////
 252 void AliTPCdigitRow::SetRow(Int_t row)
 253 {
 254   //
 255   // set bit mask for given row
 256   //
 257   if (row >= 8*32) {
 258     //    cerr<<"AliTPCdigitRow::SetRow: index "<<row<<" out of bounds."<<endl;
 259     return;
 260   }
 261   Int_t iC = row/8;
 262   Int_t iB = row%8;
 263   SETBIT(fDig[iC],iB);
 264 }
 265
 266 ////////////////////////////////////////////////////////////////////////
 267 Bool_t AliTPCdigitRow::TestRow(Int_t row) const
 268 {
 269 //
 270 // return kTRUE if row is on
 271 //
 272   Int_t iC = row/8;
 273   Int_t iB = row%8;
 274   return TESTBIT(fDig[iC],iB);
 275 }
 276 ////////////////////////////////////////////////////////////////////////
 277 Int_t AliTPCdigitRow::RowsOn(Int_t upto) const
 278 {
 279 //
 280 // returns number of rows with a digit
 281 // count only rows less equal row number upto
 282 //
 283   Int_t total = 0;
 284   for (Int_t i = 0; i<32; i++) {
 285     for (Int_t j = 0; j < 8; j++) {
 286       if (i*8+j > upto) return total;
 287       if (TESTBIT(fDig[i],j))  total++;
 288     }
 289   }
 290   return total;
 291 }
 292 ////////////////////////////////////////////////////////////////////////
 293 void AliTPCdigitRow::Reset()
 294 {
 295 //
 296 // resets all rows to zero
 297 //
 298   for (Int_t i = 0; i<32; i++) {
 299     fDig[i] <<= 8;
 300   }
 301 }
 302 ////////////////////////////////////////////////////////////////////////
 303 Int_t AliTPCdigitRow::Last() const
 304 {
 305 //
 306 // returns the last row number with a digit
 307 // returns -1 if now digits
 308 //
 309   for (Int_t i = 32-1; i>=0; i--) {
 310     for (Int_t j = 7; j >= 0; j--) {
 311       if TESTBIT(fDig[i],j) return i*8+j;
 312     }
 313   }
 314   return -1;
 315 }
 316 ////////////////////////////////////////////////////////////////////////
 317 Int_t AliTPCdigitRow::First() const
 318 {
 319 //
 320 // returns the first row number with a digit
 321 // returns -1 if now digits
 322 //
 323   for (Int_t i = 0; i<32; i++) {
 324     for (Int_t j = 0; j < 8; j++) {
 325       if (TESTBIT(fDig[i],j)) return i*8+j;
 326     }
 327   }
 328   return -1;
 329 }
 330
 331
 332 //_____________________________________________________________________________
 333 Float_t AliMCInfo::TPCBetheBloch(Float_t bg)
 334 {
 335   //
 336   // Bethe-Bloch energy loss formula
 337   //
 338   const Double_t kp1=0.76176e-1;
 339   const Double_t kp2=10.632;
 340   const Double_t kp3=0.13279e-4;
 341   const Double_t kp4=1.8631;
 342   const Double_t kp5=1.9479;
 343
 344   Double_t dbg = (Double_t) bg;
 345
 346   Double_t beta = dbg/TMath::Sqrt(1.+dbg*dbg);
 347
 348   Double_t aa = TMath::Power(beta,kp4);
 349   Double_t bb = TMath::Power(1./dbg,kp5);
 350
 351   bb=TMath::Log(kp3+bb);
 352
 353   return ((Float_t)((kp2-aa-bb)*kp1/aa));
 354 }
 355