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