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   fTPCtrackLength(-1),
  79   fPrim(0),
  80   fTPCRow(),
  81   fNTPCRef(0),                    // tpc references counter
  82   fNITSRef(0),                    // ITS references counter
  83   fNTRDRef(0),                    // TRD references counter
  84   fNTOFRef(0),                    // TOF references counter
  85   fTPCReferences(0),
  86   fITSReferences(0),
  87   fTRDReferences(0),
  88   fTOFReferences(0)
  89 {
  90   //
  91   // Default constructor
  92   //
  93   fTPCReferences  = new TClonesArray("AliTrackReference",10);
  94   fITSReferences  = new TClonesArray("AliTrackReference",10);
  95   fTRDReferences  = new TClonesArray("AliTrackReference",10);
  96   fTOFReferences  = new TClonesArray("AliTrackReference",10);
  97   fTRdecay.SetTrack(-1);
  98   fCharge = 0;
  99 }
 100
 101 AliMCInfo::AliMCInfo(const AliMCInfo& info):
 102   TObject(),
 103   fTrackRef(info.fTrackRef),
 104   fTrackRefOut(info.fTrackRefOut),
 105   fTRdecay(info.GetTRdecay()),
 106   fPrimPart(info.fPrimPart),
 107   fParticle(info.fParticle),
 108   fMass(info.GetMass()),
 109   fCharge(info.fCharge),
 110   fLabel(info.fLabel),
 111   fEventNr(info.fEventNr),
 112   fMCtracks(info.fMCtracks),
 113   fPdg(info.fPdg),
 114   fTPCdecay(info.fTPCdecay),
 115   fRowsWithDigitsInn(info.fRowsWithDigitsInn),
 116   fRowsWithDigits(info.fRowsWithDigits),
 117   fRowsTrackLength(info.fRowsTrackLength),
 118   fTPCtrackLength(info.fTPCtrackLength),
 119   fPrim(info.fPrim),
 120   fTPCRow(info.fTPCRow),
 121   fNTPCRef(info.fNTPCRef),                    // tpc references counter
 122   fNITSRef(info.fNITSRef),                    // ITS references counter
 123   fNTRDRef(info.fNTRDRef),                    // TRD references counter
 124   fNTOFRef(info.fNTOFRef),                    // TOF references counter
 125   fTPCReferences(0),
 126   fITSReferences(0),
 127   fTRDReferences(0),
 128   fTOFReferences(0)
 129 {
 130   //
 131   // copy constructor
 132   //
 133   fTPCReferences = (TClonesArray*)info.fTPCReferences->Clone();
 134   fITSReferences = (TClonesArray*)info.fITSReferences->Clone();
 135   fTRDReferences = (TClonesArray*)info.fTRDReferences->Clone();
 136   fTOFReferences = (TClonesArray*)info.fTOFReferences->Clone();
 137 }
 138
 139
 140 AliMCInfo& AliMCInfo::operator=(const AliMCInfo& info) {
 141   //
 142   // Assignment operator
 143   //
 144   this->~AliMCInfo();
 145   new (this) AliMCInfo(info);
 146   return *this;
 147 }
 148
 149
 150 AliMCInfo::~AliMCInfo()
 151 {
 152   //
 153   // Destructor of the class
 154   //
 155   if (fTPCReferences) {
 156     delete fTPCReferences;
 157   }
 158   if (fITSReferences){
 159     delete fITSReferences;
 160   }
 161   if (fTRDReferences){
 162     delete fTRDReferences;
 163   }
 164   if (fTOFReferences){
 165     delete fTOFReferences;
 166   }
 167
 168 }
 169
 170
 171
 172 void AliMCInfo::Update()
 173 {
 174   //
 175   // Update MC info
 176   // Calculates some derived variables
 177   //
 178   //
 179   fMCtracks =1;
 180   //
 181   // decay info
 182   fTPCdecay=kFALSE;
 183   if (fTRdecay.GetTrack()>0){
 184     fDecayCoord[0] = fTRdecay.X();
 185     fDecayCoord[1] = fTRdecay.Y();
 186     fDecayCoord[2] = fTRdecay.Z();
 187     if ( (fTRdecay.R()<250)&&(fTRdecay.R()>85) && (TMath::Abs(fTRdecay.Z())<250) ){
 188       fTPCdecay=kTRUE;
 189     }
 190     else{
 191       fDecayCoord[0] = 0;
 192       fDecayCoord[1] = 0;
 193       fDecayCoord[2] = 0;
 194     }
 195   }
 196   //
 197   //
 198   //digits information update
 199   fRowsWithDigits    = fTPCRow.RowsOn();
 200   fRowsWithDigitsInn = fTPCRow.RowsOn(63); // 63 = number of inner rows
 201   fRowsTrackLength   = fTPCRow.Last() - fTPCRow.First();
 202   //
 203   //
 204   // calculate primary ionization per cm
 205   if (fParticle.GetPDG()){
 206     fMass = fParticle.GetMass();
 207     fCharge = fParticle.GetPDG()->Charge();
 208     if (fTPCReferences->GetEntriesFast()>0){
 209       fTrackRef = *((AliTrackReference*)fTPCReferences->At(0));
 210     }
 211     if (fMass>0){
 212       Float_t p = TMath::Sqrt(fTrackRef.Px()*fTrackRef.Px()+
 213                               fTrackRef.Py()*fTrackRef.Py()+
 214                               fTrackRef.Pz()*fTrackRef.Pz());
 215       if (p>0.001){
 216         Float_t betagama = p /fMass;
 217         fPrim = TPCBetheBloch(betagama);
 218       }else fPrim=0;
 219     }
 220   }else{
 221     fMass =0;
 222     fPrim =0;
 223   }
 224 }
 225
 226
 227
 228
 229 ////////////////////////////////////////////////////////////////////////
 230 AliTPCdigitRow::AliTPCdigitRow()
 231 {
 232   Reset();
 233 }
 234 ////////////////////////////////////////////////////////////////////////
 235 AliTPCdigitRow & AliTPCdigitRow::operator=(const AliTPCdigitRow &digOld)
 236 {
 237   for (Int_t i = 0; i<32; i++) fDig[i] = digOld.fDig[i];
 238   return (*this);
 239 }
 240 ////////////////////////////////////////////////////////////////////////
 241 void AliTPCdigitRow::SetRow(Int_t row)
 242 {
 243   //
 244   // set bit mask for given row
 245   //
 246   if (row >= 8*32) {
 247     //    cerr<<"AliTPCdigitRow::SetRow: index "<<row<<" out of bounds."<<endl;
 248     return;
 249   }
 250   Int_t iC = row/8;
 251   Int_t iB = row%8;
 252   SETBIT(fDig[iC],iB);
 253 }
 254
 255 ////////////////////////////////////////////////////////////////////////
 256 Bool_t AliTPCdigitRow::TestRow(Int_t row) const
 257 {
 258 //
 259 // return kTRUE if row is on
 260 //
 261   Int_t iC = row/8;
 262   Int_t iB = row%8;
 263   return TESTBIT(fDig[iC],iB);
 264 }
 265 ////////////////////////////////////////////////////////////////////////
 266 Int_t AliTPCdigitRow::RowsOn(Int_t upto) const
 267 {
 268 //
 269 // returns number of rows with a digit
 270 // count only rows less equal row number upto
 271 //
 272   Int_t total = 0;
 273   for (Int_t i = 0; i<32; i++) {
 274     for (Int_t j = 0; j < 8; j++) {
 275       if (i*8+j > upto) return total;
 276       if (TESTBIT(fDig[i],j))  total++;
 277     }
 278   }
 279   return total;
 280 }
 281 ////////////////////////////////////////////////////////////////////////
 282 void AliTPCdigitRow::Reset()
 283 {
 284 //
 285 // resets all rows to zero
 286 //
 287   for (Int_t i = 0; i<32; i++) {
 288     fDig[i] <<= 8;
 289   }
 290 }
 291 ////////////////////////////////////////////////////////////////////////
 292 Int_t AliTPCdigitRow::Last() const
 293 {
 294 //
 295 // returns the last row number with a digit
 296 // returns -1 if now digits
 297 //
 298   for (Int_t i = 32-1; i>=0; i--) {
 299     for (Int_t j = 7; j >= 0; j--) {
 300       if TESTBIT(fDig[i],j) return i*8+j;
 301     }
 302   }
 303   return -1;
 304 }
 305 ////////////////////////////////////////////////////////////////////////
 306 Int_t AliTPCdigitRow::First() const
 307 {
 308 //
 309 // returns the first row number with a digit
 310 // returns -1 if now digits
 311 //
 312   for (Int_t i = 0; i<32; i++) {
 313     for (Int_t j = 0; j < 8; j++) {
 314       if (TESTBIT(fDig[i],j)) return i*8+j;
 315     }
 316   }
 317   return -1;
 318 }
 319
 320
 321 //_____________________________________________________________________________
 322 Float_t AliMCInfo::TPCBetheBloch(Float_t bg)
 323 {
 324   //
 325   // Bethe-Bloch energy loss formula
 326   //
 327   const Double_t kp1=0.76176e-1;
 328   const Double_t kp2=10.632;
 329   const Double_t kp3=0.13279e-4;
 330   const Double_t kp4=1.8631;
 331   const Double_t kp5=1.9479;
 332
 333   Double_t dbg = (Double_t) bg;
 334
 335   Double_t beta = dbg/TMath::Sqrt(1.+dbg*dbg);
 336
 337   Double_t aa = TMath::Power(beta,kp4);
 338   Double_t bb = TMath::Power(1./dbg,kp5);
 339
 340   bb=TMath::Log(kp3+bb);
 341
 342   return ((Float_t)((kp2-aa-bb)*kp1/aa));
 343 }
 344
 345
 346 void AliMCInfo::CalcTPCrows(TClonesArray * runArrayTR){
 347   //
 348   // Calculates the numebr of the track references for detectors
 349   // In case of changing direction - curling tracks - the counter is not increasing
 350   //
 351   // Rough calculation
 352   // of the first and last point in the TPC
 353   //
 354   fNTPCRef = 0;
 355   fNITSRef = 0;
 356   fNTRDRef = 0;
 357   fNTOFRef = 0;
 358   Float_t tpcminRadius=250;
 359   Float_t tpcmaxRadius=80;
 360   Float_t dir=0;
 361   Int_t nover=0;
 362
 363   for (Int_t iTrackRef = 0; iTrackRef < runArrayTR->GetEntriesFast(); iTrackRef++) {
 364     //
 365     AliTrackReference *ref = (AliTrackReference*)runArrayTR->At(iTrackRef);
 366     Float_t newdirection = (ref->X()*ref->Px()+ref->Y()*ref->Py()>0)? 1.:-1.; //inside or outside
 367     if (dir*newdirection<0.5) {
 368       nover++;
 369       dir = newdirection;
 370     }
 371     //
 372     if (ref->DetectorId()== AliTrackReference::kTRD){
 373       tpcmaxRadius =250;
 374       fNTRDRef++;
 375     }
 376     if (ref->DetectorId()== AliTrackReference::kITS){
 377       fNITSRef++;
 378       tpcminRadius =90;
 379     }
 380     if (ref->DetectorId()== AliTrackReference::kITS){
 381       fNTOFRef++;
 382       tpcmaxRadius =250;
 383     }
 384     //
 385     if (ref->DetectorId()== AliTrackReference::kTPC){
 386       fNTPCRef++;
 387       if (ref->R()>tpcmaxRadius) tpcmaxRadius = ref->R();
 388       if (ref->R()<tpcminRadius) tpcminRadius = ref->R();
 389     }
 390     if (ref->DetectorId()== AliTrackReference::kDisappeared){
 391       if (TMath::Abs(ref->Z())<250 && TMath::Abs(ref->R()<250))
 392         tpcmaxRadius = ref->R();
 393     }
 394   }
 395   fTPCtrackLength = tpcmaxRadius-tpcminRadius;
 396   fMCtracks=nover;
 397 }