PWG1/AliGenInfo.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 "TSystem.h"
  39 #include "TCanvas.h"
  40 #include "TGeometry.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 "AliGenInfo.h"
  60 //
  61 //
  62
  63 ClassImp(AliTPCdigitRow);
  64 ClassImp(AliMCInfo);
  65 ClassImp(AliGenV0Info)
  66 ClassImp(AliGenKinkInfo)
  67
  68
  69
  70 ////////////////////////////////////////////////////////////////////////
  71 AliMCInfo::AliMCInfo():
  72   fTrackRef(),
  73   fTrackRefOut(),
  74   fTRdecay(),
  75   fPrimPart(0),
  76   fParticle(),
  77   fMass(0),
  78   fCharge(0),
  79   fLabel(0),
  80   fEventNr(),
  81   fMCtracks(),
  82   fPdg(0),
  83   fTPCdecay(0),
  84   fRowsWithDigitsInn(0),
  85   fRowsWithDigits(0),
  86   fRowsTrackLength(0),
  87   fPrim(0),
  88   fTPCRow(),
  89   fNTPCRef(0),                    // tpc references counter
  90   fNITSRef(0),                    // ITS references counter
  91   fNTRDRef(0),                    // TRD references counter
  92   fNTOFRef(0),                    // TOF references counter
  93   fTPCReferences(0),
  94   fITSReferences(0),
  95   fTRDReferences(0),
  96   fTOFReferences(0)
  97 {
  98   fTPCReferences  = new TClonesArray("AliTrackReference",10);
  99   fITSReferences  = new TClonesArray("AliTrackReference",10);
 100   fTRDReferences  = new TClonesArray("AliTrackReference",10);
 101   fTOFReferences  = new TClonesArray("AliTrackReference",10);
 102   fTRdecay.SetTrack(-1);
 103   fCharge = 0;
 104 }
 105
 106 AliMCInfo::AliMCInfo(const AliMCInfo& info):
 107   TObject(),
 108   fTrackRef(info.fTrackRef),
 109   fTrackRefOut(info.fTrackRefOut),
 110   fTRdecay(info.GetTRdecay()),
 111   fPrimPart(info.fPrimPart),
 112   fParticle(info.fParticle),
 113   fMass(info.GetMass()),
 114   fCharge(info.fCharge),
 115   fLabel(info.fLabel),
 116   fEventNr(info.fEventNr),
 117   fMCtracks(info.fMCtracks),
 118   fPdg(info.fPdg),
 119   fTPCdecay(info.fTPCdecay),
 120   fRowsWithDigitsInn(info.fRowsWithDigitsInn),
 121   fRowsWithDigits(info.fRowsWithDigits),
 122   fRowsTrackLength(info.fRowsTrackLength),
 123   fPrim(info.fPrim),
 124   fTPCRow(info.fTPCRow),
 125   fNTPCRef(info.fNTPCRef),                    // tpc references counter
 126   fNITSRef(info.fNITSRef),                    // ITS references counter
 127   fNTRDRef(info.fNTRDRef),                    // TRD references counter
 128   fNTOFRef(info.fNTOFRef),                    // TOF references counter
 129   fTPCReferences(0),
 130   fITSReferences(0),
 131   fTRDReferences(0),
 132   fTOFReferences(0)
 133 {
 134   fTPCReferences = (TClonesArray*)info.fTPCReferences->Clone();
 135   fITSReferences = (TClonesArray*)info.fITSReferences->Clone();
 136   fTRDReferences = (TClonesArray*)info.fTRDReferences->Clone();
 137   fTOFReferences = (TClonesArray*)info.fTOFReferences->Clone();
 138 }
 139
 140
 141 AliMCInfo::~AliMCInfo()
 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   //
 164   fMCtracks =1;
 165   if (!fTPCReferences) {
 166     fNTPCRef =0;
 167     return;
 168   }
 169   Float_t direction=1;
 170   //Float_t rlast=0;
 171   fNTPCRef = fTPCReferences->GetEntriesFast();
 172   fNITSRef = fITSReferences->GetEntriesFast();
 173   fNTRDRef = fTRDReferences->GetEntriesFast();
 174   fNTOFRef = fTOFReferences->GetEntriesFast();
 175
 176   for (Int_t iref =0;iref<fTPCReferences->GetEntriesFast();iref++){
 177     AliTrackReference * ref = (AliTrackReference *) fTPCReferences->At(iref);
 178     //Float_t r = (ref->X()*ref->X()+ref->Y()*ref->Y());
 179     Float_t newdirection = ref->X()*ref->Px()+ref->Y()*ref->Py(); //inside or outside
 180     if (iref==0) direction = newdirection;
 181     if ( newdirection*direction<0){
 182       //changed direction
 183       direction = newdirection;
 184       fMCtracks+=1;
 185     }
 186     //rlast=r;
 187   }
 188   //
 189   // decay info
 190   fTPCdecay=kFALSE;
 191   if (fTRdecay.GetTrack()>0){
 192     fDecayCoord[0] = fTRdecay.X();
 193     fDecayCoord[1] = fTRdecay.Y();
 194     fDecayCoord[2] = fTRdecay.Z();
 195     if ( (fTRdecay.R()<250)&&(fTRdecay.R()>85) && (TMath::Abs(fTRdecay.Z())<250) ){
 196       fTPCdecay=kTRUE;
 197     }
 198     else{
 199       fDecayCoord[0] = 0;
 200       fDecayCoord[1] = 0;
 201       fDecayCoord[2] = 0;
 202     }
 203   }
 204   //
 205   //
 206   //digits information update
 207   fRowsWithDigits    = fTPCRow.RowsOn();
 208   fRowsWithDigitsInn = fTPCRow.RowsOn(63); // 63 = number of inner rows
 209   fRowsTrackLength   = fTPCRow.Last() - fTPCRow.First();
 210   //
 211   //
 212   // calculate primary ionization per cm
 213   if (fParticle.GetPDG()){
 214     fMass = fParticle.GetMass();
 215     fCharge = fParticle.GetPDG()->Charge();
 216     if (fTPCReferences->GetEntriesFast()>0){
 217       fTrackRef = *((AliTrackReference*)fTPCReferences->At(0));
 218     }
 219     if (fMass>0){
 220       Float_t p = TMath::Sqrt(fTrackRef.Px()*fTrackRef.Px()+
 221                               fTrackRef.Py()*fTrackRef.Py()+
 222                               fTrackRef.Pz()*fTrackRef.Pz());
 223       if (p>0.001){
 224         Float_t betagama = p /fMass;
 225         fPrim = TPCBetheBloch(betagama);
 226       }else fPrim=0;
 227     }
 228   }else{
 229     fMass =0;
 230     fPrim =0;
 231   }
 232 }
 233
 234 /////////////////////////////////////////////////////////////////////////////////
 235 AliGenV0Info::AliGenV0Info():
 236   fMCd(),      //info about daughter particle -
 237   fMCm(),      //info about mother particle   - first particle for V0
 238   fMotherP(),   //particle info about mother particle
 239   fMCDist1(0), //info about closest distance according closest MC - linear DCA
 240   fMCDist2(0),    //info about closest distance parabolic DCA
 241   fMCRr(0),       // rec position of the vertex
 242   fMCR(0),        //exact r position of the vertex
 243   fInvMass(0),  //reconstructed invariant mass -
 244   fPointAngleFi(0), //point angle fi
 245   fPointAngleTh(0), //point angle theta
 246   fPointAngle(0)   //point angle full
 247 {
 248   for (Int_t i=0;i<3; i++){
 249    fMCPdr[i]=0;
 250    fMCX[i]=0;
 251    fMCXr[i]=0;
 252    fMCPm[i]=0;
 253    fMCAngle[i]=0;
 254    fMCPd[i]=0;
 255   }
 256   fMCPd[3]=0;
 257   for (Int_t i=0; i<2;i++){
 258     fPdg[i]=0;
 259     fLab[i]=0;
 260   }
 261 }
 262
 263 void AliGenV0Info::Update(Float_t vertex[3])
 264 {
 265   //
 266   // Update information - calculates derived variables
 267   //
 268
 269   fMCPd[0] = fMCd.GetParticle().Px();
 270   fMCPd[1] = fMCd.GetParticle().Py();
 271   fMCPd[2] = fMCd.GetParticle().Pz();
 272   fMCPd[3] = fMCd.GetParticle().P();
 273   //
 274   fMCX[0]  = fMCd.GetParticle().Vx();
 275   fMCX[1]  = fMCd.GetParticle().Vy();
 276   fMCX[2]  = fMCd.GetParticle().Vz();
 277   fMCR       = TMath::Sqrt( fMCX[0]*fMCX[0]+fMCX[1]*fMCX[1]);
 278   //
 279   fPdg[0]    = fMCd.GetParticle().GetPdgCode();
 280   fPdg[1]    = fMCm.GetParticle().GetPdgCode();
 281   //
 282   fLab[0]    = fMCd.GetParticle().GetUniqueID();
 283   fLab[1]    = fMCm.GetParticle().GetUniqueID();
 284   //
 285   //
 286   //
 287   Double_t x1[3],p1[3];
 288   TParticle& pdaughter = fMCd.GetParticle();
 289   x1[0] = pdaughter.Vx();
 290   x1[1] = pdaughter.Vy();
 291   x1[2] = pdaughter.Vz();
 292   p1[0] = pdaughter.Px();
 293   p1[1] = pdaughter.Py();
 294   p1[2] = pdaughter.Pz();
 295   Double_t sign = (pdaughter.GetPDG()->Charge()>0)? -1:1;
 296   AliHelix dhelix1(x1,p1,sign);
 297   //
 298   //
 299   Double_t x2[3],p2[3];
 300   //
 301   TParticle& pmother = fMCm.GetParticle();
 302   x2[0] = pmother.Vx();
 303   x2[1] = pmother.Vy();
 304   x2[2] = pmother.Vz();
 305   p2[0] = pmother.Px();
 306   p2[1] = pmother.Py();
 307   p2[2] = pmother.Pz();
 308   //
 309   //
 310   sign = (pmother.GetPDG()->Charge()>0) ? -1:1;
 311   AliHelix mhelix(x2,p2,sign);
 312
 313   //
 314   //
 315   //
 316   //find intersection linear
 317   //
 318   Double_t distance1, distance2;
 319   Double_t phase[2][2],radius[2];
 320   Int_t  points = dhelix1.GetRPHIintersections(mhelix, phase, radius);
 321   Double_t delta1=10000,delta2=10000;
 322   if (points>0){
 323     dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
 324     dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
 325     dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
 326   }
 327   else{
 328     fInvMass=-1;
 329     return;
 330   }
 331   if (points==2){
 332     dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
 333     dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
 334     dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
 335   }
 336   distance1 = TMath::Min(delta1,delta2);
 337   //
 338   //find intersection parabolic
 339   //
 340   points = dhelix1.GetRPHIintersections(mhelix, phase, radius);
 341   delta1=10000,delta2=10000;
 342
 343   if (points>0){
 344     dhelix1.ParabolicDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
 345   }
 346   if (points==2){
 347     dhelix1.ParabolicDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
 348   }
 349
 350   distance2 = TMath::Min(delta1,delta2);
 351   //
 352   if (delta1<delta2){
 353     //get V0 info
 354     dhelix1.Evaluate(phase[0][0],fMCXr);
 355     dhelix1.GetMomentum(phase[0][0],fMCPdr);
 356     mhelix.GetMomentum(phase[0][1],fMCPm);
 357     dhelix1.GetAngle(phase[0][0],mhelix,phase[0][1],fMCAngle);
 358     fMCRr = TMath::Sqrt(radius[0]);
 359   }
 360   else{
 361     dhelix1.Evaluate(phase[1][0],fMCXr);
 362     dhelix1.GetMomentum(phase[1][0],fMCPdr);
 363     mhelix.GetMomentum(phase[1][1],fMCPm);
 364     dhelix1.GetAngle(phase[1][0],mhelix,phase[1][1],fMCAngle);
 365     fMCRr = TMath::Sqrt(radius[1]);
 366   }
 367   //
 368   //
 369   fMCDist1 = TMath::Sqrt(distance1);
 370   fMCDist2 = TMath::Sqrt(distance2);
 371   //
 372   //
 373   //
 374   Float_t v[3] = {fMCXr[0]-vertex[0],fMCXr[1]-vertex[1],fMCXr[2]-vertex[2]};
 375   //Float_t v[3] = {fMCXr[0],fMCXr[1],fMCXr[2]};
 376   Float_t p[3] = {fMCPdr[0]+fMCPm[0], fMCPdr[1]+fMCPm[1],fMCPdr[2]+fMCPm[2]};
 377   Float_t vnorm2 = v[0]*v[0]+v[1]*v[1];
 378   Float_t vnorm3 = TMath::Sqrt(v[2]*v[2]+vnorm2);
 379   vnorm2 = TMath::Sqrt(vnorm2);
 380   Float_t pnorm2 = p[0]*p[0]+p[1]*p[1];
 381   Float_t pnorm3 = TMath::Sqrt(p[2]*p[2]+pnorm2);
 382   pnorm2 = TMath::Sqrt(pnorm2);
 383   //
 384   fPointAngleFi = (v[0]*p[0]+v[1]*p[1])/(vnorm2*pnorm2);
 385   fPointAngleTh = (v[2]*p[2]+vnorm2*pnorm2)/(vnorm3*pnorm3);
 386   fPointAngle   = (v[0]*p[0]+v[1]*p[1]+v[2]*p[2])/(vnorm3*pnorm3);
 387   Double_t mass1 = fMCd.GetMass();
 388   Double_t mass2 = fMCm.GetMass();
 389
 390
 391   Double_t e1    = TMath::Sqrt(mass1*mass1+
 392                               fMCPd[0]*fMCPd[0]+
 393                               fMCPd[1]*fMCPd[1]+
 394                               fMCPd[2]*fMCPd[2]);
 395   Double_t e2    = TMath::Sqrt(mass2*mass2+
 396                               fMCPm[0]*fMCPm[0]+
 397                               fMCPm[1]*fMCPm[1]+
 398                               fMCPm[2]*fMCPm[2]);
 399
 400   fInvMass =
 401     (fMCPm[0]+fMCPd[0])*(fMCPm[0]+fMCPd[0])+
 402     (fMCPm[1]+fMCPd[1])*(fMCPm[1]+fMCPd[1])+
 403     (fMCPm[2]+fMCPd[2])*(fMCPm[2]+fMCPd[2]);
 404
 405   //  fInvMass = TMath::Sqrt((e1+e2)*(e1+e2)-fInvMass);
 406   fInvMass = (e1+e2)*(e1+e2)-fInvMass;
 407   if (fInvMass>0) fInvMass = TMath::Sqrt(fInvMass);
 408 }
 409
 410
 411
 412 /////////////////////////////////////////////////////////////////////////////////
 413
 414 AliGenKinkInfo::AliGenKinkInfo():
 415   fMCd(),       //info about daughter particle - second particle for V0
 416   fMCm(),       //info about mother particle   - first particle for V0
 417   fMCDist1(0),  //info about closest distance according closest MC - linear DCA
 418   fMCDist2(0),  //info about closest distance parabolic DCA
 419   fMCRr(0),     // rec position of the vertex
 420   fMCR(0)       //exact r position of the vertex
 421 {
 422   //
 423   // default constructor
 424   //
 425   for (Int_t i=0;i<3;i++){
 426     fMCPdr[i]=0;
 427     fMCPd[i]=0;
 428     fMCX[i]=0;
 429     fMCPm[i]=0;
 430     fMCAngle[i]=0;
 431   }
 432   for (Int_t i=0; i<2; i++) {
 433     fPdg[i]= 0;
 434     fLab[i]=0;
 435   }
 436 }
 437
 438 void AliGenKinkInfo::Update()
 439 {
 440   //
 441   // Update information
 442   //  some redundancy - faster acces to the values in analysis code
 443   //
 444   fMCPd[0] = fMCd.GetParticle().Px();
 445   fMCPd[1] = fMCd.GetParticle().Py();
 446   fMCPd[2] = fMCd.GetParticle().Pz();
 447   fMCPd[3] = fMCd.GetParticle().P();
 448   //
 449   fMCX[0]  = fMCd.GetParticle().Vx();
 450   fMCX[1]  = fMCd.GetParticle().Vy();
 451   fMCX[2]  = fMCd.GetParticle().Vz();
 452   fMCR       = TMath::Sqrt( fMCX[0]*fMCX[0]+fMCX[1]*fMCX[1]);
 453   //
 454   fPdg[0]    = fMCd.GetParticle().GetPdgCode();
 455   fPdg[1]    = fMCm.GetParticle().GetPdgCode();
 456   //
 457   fLab[0]    = fMCd.GetParticle().GetUniqueID();
 458   fLab[1]    = fMCm.GetParticle().GetUniqueID();
 459   //
 460   //
 461   //
 462   Double_t x1[3],p1[3];
 463   TParticle& pdaughter = fMCd.GetParticle();
 464   x1[0] = pdaughter.Vx();
 465   x1[1] = pdaughter.Vy();
 466   x1[2] = pdaughter.Vz();
 467   p1[0] = pdaughter.Px();
 468   p1[1] = pdaughter.Py();
 469   p1[2] = pdaughter.Pz();
 470   Double_t sign = (pdaughter.GetPDG()->Charge()>0)? -1:1;
 471   AliHelix dhelix1(x1,p1,sign);
 472   //
 473   //
 474   Double_t x2[3],p2[3];
 475   //
 476   TParticle& pmother = fMCm.GetParticle();
 477   x2[0] = pmother.Vx();
 478   x2[1] = pmother.Vy();
 479   x2[2] = pmother.Vz();
 480   p2[0] = pmother.Px();
 481   p2[1] = pmother.Py();
 482   p2[2] = pmother.Pz();
 483   //
 484   const AliTrackReference & pdecay = fMCm.GetTRdecay();
 485   x2[0] = pdecay.X();
 486   x2[1] = pdecay.Y();
 487   x2[2] = pdecay.Z();
 488   p2[0] = pdecay.Px();
 489   p2[1] = pdecay.Py();
 490   p2[2] = pdecay.Pz();
 491   //
 492   sign = (pmother.GetPDG()->Charge()>0) ? -1:1;
 493   AliHelix mhelix(x2,p2,sign);
 494
 495   //
 496   //
 497   //
 498   //find intersection linear
 499   //
 500   Double_t distance1, distance2;
 501   Double_t phase[2][2],radius[2];
 502   Int_t  points = dhelix1.GetRPHIintersections(mhelix, phase, radius);
 503   Double_t delta1=10000,delta2=10000;
 504   if (points>0){
 505     dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
 506     dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
 507     dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
 508   }
 509   if (points==2){
 510     dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
 511     dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
 512     dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
 513   }
 514   distance1 = TMath::Min(delta1,delta2);
 515   //
 516   //find intersection parabolic
 517   //
 518   points = dhelix1.GetRPHIintersections(mhelix, phase, radius);
 519   delta1=10000,delta2=10000;
 520
 521   if (points>0){
 522     dhelix1.ParabolicDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
 523   }
 524   if (points==2){
 525     dhelix1.ParabolicDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
 526   }
 527
 528   distance2 = TMath::Min(delta1,delta2);
 529   //
 530   if (delta1<delta2){
 531     //get V0 info
 532     dhelix1.Evaluate(phase[0][0],fMCXr);
 533     dhelix1.GetMomentum(phase[0][0],fMCPdr);
 534     mhelix.GetMomentum(phase[0][1],fMCPm);
 535     dhelix1.GetAngle(phase[0][0],mhelix,phase[0][1],fMCAngle);
 536     fMCRr = TMath::Sqrt(radius[0]);
 537   }
 538   else{
 539     dhelix1.Evaluate(phase[1][0],fMCXr);
 540     dhelix1.GetMomentum(phase[1][0],fMCPdr);
 541     mhelix.GetMomentum(phase[1][1],fMCPm);
 542     dhelix1.GetAngle(phase[1][0],mhelix,phase[1][1],fMCAngle);
 543     fMCRr = TMath::Sqrt(radius[1]);
 544   }
 545   //
 546   //
 547   fMCDist1 = TMath::Sqrt(distance1);
 548   fMCDist2 = TMath::Sqrt(distance2);
 549
 550 }
 551
 552
 553 Float_t AliGenKinkInfo::GetQt(){
 554   //
 555   //
 556   Float_t momentumd = TMath::Sqrt(fMCPd[0]*fMCPd[0]+fMCPd[1]*fMCPd[1]+fMCPd[2]*fMCPd[2]);
 557   return TMath::Sin(fMCAngle[2])*momentumd;
 558 }
 559
 560
 561
 562
 563 ////////////////////////////////////////////////////////////////////////
 564 AliTPCdigitRow::AliTPCdigitRow()
 565 {
 566   Reset();
 567 }
 568 ////////////////////////////////////////////////////////////////////////
 569 AliTPCdigitRow & AliTPCdigitRow::operator=(const AliTPCdigitRow &digOld)
 570 {
 571   for (Int_t i = 0; i<kgRowBytes; i++) fDig[i] = digOld.fDig[i];
 572   return (*this);
 573 }
 574 ////////////////////////////////////////////////////////////////////////
 575 void AliTPCdigitRow::SetRow(Int_t row)
 576 {
 577   if (row >= 8*kgRowBytes) {
 578     cerr<<"AliTPCdigitRow::SetRow: index "<<row<<" out of bounds."<<endl;
 579     return;
 580   }
 581   Int_t iC = row/8;
 582   Int_t iB = row%8;
 583   SETBIT(fDig[iC],iB);
 584 }
 585
 586 ////////////////////////////////////////////////////////////////////////
 587 Bool_t AliTPCdigitRow::TestRow(Int_t row) const
 588 {
 589 //
 590 // return kTRUE if row is on
 591 //
 592   Int_t iC = row/8;
 593   Int_t iB = row%8;
 594   return TESTBIT(fDig[iC],iB);
 595 }
 596 ////////////////////////////////////////////////////////////////////////
 597 Int_t AliTPCdigitRow::RowsOn(Int_t upto) const
 598 {
 599 //
 600 // returns number of rows with a digit
 601 // count only rows less equal row number upto
 602 //
 603   Int_t total = 0;
 604   for (Int_t i = 0; i<kgRowBytes; i++) {
 605     for (Int_t j = 0; j < 8; j++) {
 606       if (i*8+j > upto) return total;
 607       if (TESTBIT(fDig[i],j))  total++;
 608     }
 609   }
 610   return total;
 611 }
 612 ////////////////////////////////////////////////////////////////////////
 613 void AliTPCdigitRow::Reset()
 614 {
 615 //
 616 // resets all rows to zero
 617 //
 618   for (Int_t i = 0; i<kgRowBytes; i++) {
 619     fDig[i] <<= 8;
 620   }
 621 }
 622 ////////////////////////////////////////////////////////////////////////
 623 Int_t AliTPCdigitRow::Last() const
 624 {
 625 //
 626 // returns the last row number with a digit
 627 // returns -1 if now digits
 628 //
 629   for (Int_t i = kgRowBytes-1; i>=0; i--) {
 630     for (Int_t j = 7; j >= 0; j--) {
 631       if TESTBIT(fDig[i],j) return i*8+j;
 632     }
 633   }
 634   return -1;
 635 }
 636 ////////////////////////////////////////////////////////////////////////
 637 Int_t AliTPCdigitRow::First() const
 638 {
 639 //
 640 // returns the first row number with a digit
 641 // returns -1 if now digits
 642 //
 643   for (Int_t i = 0; i<kgRowBytes; i++) {
 644     for (Int_t j = 0; j < 8; j++) {
 645       if (TESTBIT(fDig[i],j)) return i*8+j;
 646     }
 647   }
 648   return -1;
 649 }
 650
 651
 652 //_____________________________________________________________________________
 653 Float_t AliMCInfo::TPCBetheBloch(Float_t bg)
 654 {
 655   //
 656   // Bethe-Bloch energy loss formula
 657   //
 658   const Double_t kp1=0.76176e-1;
 659   const Double_t kp2=10.632;
 660   const Double_t kp3=0.13279e-4;
 661   const Double_t kp4=1.8631;
 662   const Double_t kp5=1.9479;
 663
 664   Double_t dbg = (Double_t) bg;
 665
 666   Double_t beta = dbg/TMath::Sqrt(1.+dbg*dbg);
 667
 668   Double_t aa = TMath::Power(beta,kp4);
 669   Double_t bb = TMath::Power(1./dbg,kp5);
 670
 671   bb=TMath::Log(kp3+bb);
 672
 673   return ((Float_t)((kp2-aa-bb)*kp1/aa));
 674 }
 675