[u/mrichter/AliRoot.git] / PWG1 / AliGenInfo.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/


///////////////////////////////////////////////////////////////////////////
/*

Origin: marian.ivanov@cern.ch
Container classes with MC infomation

*/

#if !defined(__CINT__) || defined(__MAKECINT__)
#include <stdio.h>
#include <string.h>
//ROOT includes
#include "TROOT.h"
#include "Rtypes.h"
#include "TFile.h"
#include "TTree.h"
#include "TChain.h"
#include "TCut.h"
#include "TString.h"
#include "TStopwatch.h"
#include "TParticle.h"
#include "TSystem.h"
#include "TCanvas.h"
#include "TGeometry.h"
#include "TPolyLine3D.h"

//ALIROOT includes
#include "AliRun.h"
#include "AliStack.h"
#include "AliSimDigits.h"
#include "AliTPCParam.h"
#include "AliTPC.h"
#include "AliTPCLoader.h"
#include "AliDetector.h"
#include "AliTrackReference.h"
#include "AliTPCParamSR.h"
#include "AliTracker.h"
#include "AliMagF.h"
#include "AliHelix.h"
#include "AliTrackPointArray.h"

#endif
#include "AliGenInfo.h" 
//
// 

ClassImp(AliTPCdigitRow);
ClassImp(AliMCInfo);
ClassImp(AliGenV0Info)
ClassImp(AliGenKinkInfo)


////////////////////////////////////////////////////////////////////////
AliMCInfo::AliMCInfo():
  fTrackRef(),
  fTrackRefOut(),
  fTRdecay(),
  fPrimPart(0),
  fParticle(),
  fMass(0),
  fCharge(0),
  fLabel(0),
  fEventNr(),
  fMCtracks(),
  fPdg(0),
  fTPCdecay(0),
  fRowsWithDigitsInn(0),
  fRowsWithDigits(0),
  fRowsTrackLength(0),
  fPrim(0),
  fTPCRow(), 
  fNTPCRef(0),                    // tpc references counter
  fNITSRef(0),                    // ITS references counter
  fNTRDRef(0),                    // TRD references counter
  fNTOFRef(0),                    // TOF references counter
  fTPCReferences(0),
  fITSReferences(0),
  fTRDReferences(0),
  fTOFReferences(0)
{
  fTPCReferences  = new TClonesArray("AliTrackReference",10);
  fITSReferences  = new TClonesArray("AliTrackReference",10);
  fTRDReferences  = new TClonesArray("AliTrackReference",10);
  fTOFReferences  = new TClonesArray("AliTrackReference",10);
  fTRdecay.SetTrack(-1);
  fCharge = 0;
}

AliMCInfo::AliMCInfo(const AliMCInfo& info):
  TObject(),
  fTrackRef(info.fTrackRef),
  fTrackRefOut(info.fTrackRefOut),
  fTRdecay(info.GetTRdecay()),
  fPrimPart(info.fPrimPart),
  fParticle(info.fParticle),
  fMass(info.GetMass()),
  fCharge(info.fCharge),
  fLabel(info.fLabel),
  fEventNr(info.fEventNr),
  fMCtracks(info.fMCtracks),
  fPdg(info.fPdg),
  fTPCdecay(info.fTPCdecay),
  fRowsWithDigitsInn(info.fRowsWithDigitsInn),
  fRowsWithDigits(info.fRowsWithDigits),
  fRowsTrackLength(info.fRowsTrackLength),
  fPrim(info.fPrim),
  fTPCRow(info.fTPCRow), 
  fNTPCRef(info.fNTPCRef),                    // tpc references counter
  fNITSRef(info.fNITSRef),                    // ITS references counter
  fNTRDRef(info.fNTRDRef),                    // TRD references counter
  fNTOFRef(info.fNTOFRef),                    // TOF references counter
  fTPCReferences(0),
  fITSReferences(0),
  fTRDReferences(0),
  fTOFReferences(0)
{
  fTPCReferences = (TClonesArray*)info.fTPCReferences->Clone();
  fITSReferences = (TClonesArray*)info.fITSReferences->Clone();
  fTRDReferences = (TClonesArray*)info.fTRDReferences->Clone();
  fTOFReferences = (TClonesArray*)info.fTOFReferences->Clone();
}


AliMCInfo::~AliMCInfo()
{
  if (fTPCReferences) {
    delete fTPCReferences;
  }
  if (fITSReferences){
    delete fITSReferences;
  }
  if (fTRDReferences){    
    delete fTRDReferences;  
  }
  if (fTOFReferences){    
    delete fTOFReferences;  
  }

}


void AliMCInfo::Update()
{
  //
  //
  fMCtracks =1;
  if (!fTPCReferences) {
    fNTPCRef =0;
    return;
  }
  Float_t direction=1;
  //Float_t rlast=0;
  fNTPCRef = fTPCReferences->GetEntriesFast();
  fNITSRef = fITSReferences->GetEntriesFast();
  fNTRDRef = fTRDReferences->GetEntriesFast();
  fNTOFRef = fTOFReferences->GetEntriesFast();
  
  for (Int_t iref =0;iref<fTPCReferences->GetEntriesFast();iref++){
    AliTrackReference * ref = (AliTrackReference *) fTPCReferences->At(iref);
    //Float_t r = (ref->X()*ref->X()+ref->Y()*ref->Y());
    Float_t newdirection = ref->X()*ref->Px()+ref->Y()*ref->Py(); //inside or outside
    if (iref==0) direction = newdirection;
    if ( newdirection*direction<0){
      //changed direction
      direction = newdirection;
      fMCtracks+=1;
    }
    //rlast=r;			    
  }
  //
  // decay info
  fTPCdecay=kFALSE;
  if (fTRdecay.GetTrack()>0){
    fDecayCoord[0] = fTRdecay.X();
    fDecayCoord[1] = fTRdecay.Y();
    fDecayCoord[2] = fTRdecay.Z();
    if ( (fTRdecay.R()<250)&&(fTRdecay.R()>85) && (TMath::Abs(fTRdecay.Z())<250) ){
      fTPCdecay=kTRUE;     
    }
    else{
      fDecayCoord[0] = 0;
      fDecayCoord[1] = 0;
      fDecayCoord[2] = 0;
    }
  }
  //
  //
  //digits information update
  fRowsWithDigits    = fTPCRow.RowsOn();    
  fRowsWithDigitsInn = fTPCRow.RowsOn(63); // 63 = number of inner rows
  fRowsTrackLength   = fTPCRow.Last() - fTPCRow.First();
  //
  //
  // calculate primary ionization per cm  
  if (fParticle.GetPDG()){
    fMass = fParticle.GetMass();  
    fCharge = fParticle.GetPDG()->Charge();
    if (fTPCReferences->GetEntriesFast()>0){
      fTrackRef = *((AliTrackReference*)fTPCReferences->At(0));
    }
    if (fMass>0){
      Float_t p = TMath::Sqrt(fTrackRef.Px()*fTrackRef.Px()+
			      fTrackRef.Py()*fTrackRef.Py()+
			      fTrackRef.Pz()*fTrackRef.Pz());    
      if (p>0.001){
	Float_t betagama = p /fMass;
	fPrim = TPCBetheBloch(betagama);
      }else fPrim=0;
    }
  }else{
    fMass =0;
    fPrim =0;
  }  
}

/////////////////////////////////////////////////////////////////////////////////
AliGenV0Info::AliGenV0Info():
  fMCd(),      //info about daughter particle - 
  fMCm(),      //info about mother particle   - first particle for V0
  fMotherP(),   //particle info about mother particle
  fMCDist1(0), //info about closest distance according closest MC - linear DCA
  fMCDist2(0),    //info about closest distance parabolic DCA
  fMCRr(0),       // rec position of the vertex 
  fMCR(0),        //exact r position of the vertex
  fInvMass(0),  //reconstructed invariant mass -
  fPointAngleFi(0), //point angle fi
  fPointAngleTh(0), //point angle theta
  fPointAngle(0)   //point angle full
{
  for (Int_t i=0;i<3; i++){
   fMCPdr[i]=0;    
   fMCX[i]=0;     
   fMCXr[i]=0;    
   fMCPm[i]=0;    
   fMCAngle[i]=0; 
   fMCPd[i]=0;    
  }
  fMCPd[3]=0;     
  for (Int_t i=0; i<2;i++){
    fPdg[i]=0;   
    fLab[i]=0;  
  }
}

void AliGenV0Info::Update(Float_t vertex[3])
{
  //
  // Update information - calculates derived variables
  //
  
  fMCPd[0] = fMCd.GetParticle().Px();
  fMCPd[1] = fMCd.GetParticle().Py();
  fMCPd[2] = fMCd.GetParticle().Pz();
  fMCPd[3] = fMCd.GetParticle().P();
  //
  fMCX[0]  = fMCd.GetParticle().Vx();
  fMCX[1]  = fMCd.GetParticle().Vy();
  fMCX[2]  = fMCd.GetParticle().Vz();
  fMCR       = TMath::Sqrt( fMCX[0]*fMCX[0]+fMCX[1]*fMCX[1]);
  //
  fPdg[0]    = fMCd.GetParticle().GetPdgCode();
  fPdg[1]    = fMCm.GetParticle().GetPdgCode();
  //
  fLab[0]    = fMCd.GetParticle().GetUniqueID();
  fLab[1]    = fMCm.GetParticle().GetUniqueID();
  //
  //
  //
  Double_t x1[3],p1[3];
  TParticle& pdaughter = fMCd.GetParticle();
  x1[0] = pdaughter.Vx();      
  x1[1] = pdaughter.Vy();
  x1[2] = pdaughter.Vz();
  p1[0] = pdaughter.Px();
  p1[1] = pdaughter.Py();
  p1[2] = pdaughter.Pz();
  Double_t sign = (pdaughter.GetPDG()->Charge()>0)? -1:1;
  AliHelix dhelix1(x1,p1,sign);
  //
  //
  Double_t x2[3],p2[3];            
  //
  TParticle& pmother = fMCm.GetParticle();
  x2[0] = pmother.Vx();      
  x2[1] = pmother.Vy();      
  x2[2] = pmother.Vz();      
  p2[0] = pmother.Px();
  p2[1] = pmother.Py();
  p2[2] = pmother.Pz();
  //
  //
  sign = (pmother.GetPDG()->Charge()>0) ? -1:1;
  AliHelix mhelix(x2,p2,sign);
  
  //
  //
  //
  //find intersection linear
  //
  Double_t distance1, distance2;
  Double_t phase[2][2],radius[2];
  Int_t  points = dhelix1.GetRPHIintersections(mhelix, phase, radius);
  Double_t delta1=10000,delta2=10000;    
  if (points>0){
    dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
    dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
    dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
  }
  else{
    fInvMass=-1;
    return;
  }
  if (points==2){    
    dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
    dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
    dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
  }
  distance1 = TMath::Min(delta1,delta2);
  //
  //find intersection parabolic
  //
  points = dhelix1.GetRPHIintersections(mhelix, phase, radius);
  delta1=10000,delta2=10000;  
  
  if (points>0){
    dhelix1.ParabolicDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
  }
  if (points==2){    
    dhelix1.ParabolicDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
  }
  
  distance2 = TMath::Min(delta1,delta2);
  //
  if (delta1<delta2){
    //get V0 info
    dhelix1.Evaluate(phase[0][0],fMCXr);
    dhelix1.GetMomentum(phase[0][0],fMCPdr);
    mhelix.GetMomentum(phase[0][1],fMCPm);
    dhelix1.GetAngle(phase[0][0],mhelix,phase[0][1],fMCAngle);
    fMCRr = TMath::Sqrt(radius[0]);
  }
  else{
    dhelix1.Evaluate(phase[1][0],fMCXr);
    dhelix1.GetMomentum(phase[1][0],fMCPdr);
    mhelix.GetMomentum(phase[1][1],fMCPm);
    dhelix1.GetAngle(phase[1][0],mhelix,phase[1][1],fMCAngle);
    fMCRr = TMath::Sqrt(radius[1]);
  }     
  //            
  //
  fMCDist1 = TMath::Sqrt(distance1);
  fMCDist2 = TMath::Sqrt(distance2);      
  //
  //
  // 
  Float_t v[3] = {fMCXr[0]-vertex[0],fMCXr[1]-vertex[1],fMCXr[2]-vertex[2]};
  //Float_t v[3] = {fMCXr[0],fMCXr[1],fMCXr[2]};
  Float_t p[3] = {fMCPdr[0]+fMCPm[0], fMCPdr[1]+fMCPm[1],fMCPdr[2]+fMCPm[2]};
  Float_t vnorm2 = v[0]*v[0]+v[1]*v[1];
  Float_t vnorm3 = TMath::Sqrt(v[2]*v[2]+vnorm2);
  vnorm2 = TMath::Sqrt(vnorm2);
  Float_t pnorm2 = p[0]*p[0]+p[1]*p[1];
  Float_t pnorm3 = TMath::Sqrt(p[2]*p[2]+pnorm2);
  pnorm2 = TMath::Sqrt(pnorm2);
  //
  fPointAngleFi = (v[0]*p[0]+v[1]*p[1])/(vnorm2*pnorm2);
  fPointAngleTh = (v[2]*p[2]+vnorm2*pnorm2)/(vnorm3*pnorm3);  
  fPointAngle   = (v[0]*p[0]+v[1]*p[1]+v[2]*p[2])/(vnorm3*pnorm3);
  Double_t mass1 = fMCd.GetMass();
  Double_t mass2 = fMCm.GetMass();

  
  Double_t e1    = TMath::Sqrt(mass1*mass1+
			      fMCPd[0]*fMCPd[0]+
			      fMCPd[1]*fMCPd[1]+
			      fMCPd[2]*fMCPd[2]);
  Double_t e2    = TMath::Sqrt(mass2*mass2+
			      fMCPm[0]*fMCPm[0]+
			      fMCPm[1]*fMCPm[1]+
			      fMCPm[2]*fMCPm[2]);
  
  fInvMass =  
    (fMCPm[0]+fMCPd[0])*(fMCPm[0]+fMCPd[0])+
    (fMCPm[1]+fMCPd[1])*(fMCPm[1]+fMCPd[1])+
    (fMCPm[2]+fMCPd[2])*(fMCPm[2]+fMCPd[2]);
  
  //  fInvMass = TMath::Sqrt((e1+e2)*(e1+e2)-fInvMass);
  fInvMass = (e1+e2)*(e1+e2)-fInvMass;
  if (fInvMass>0) fInvMass = TMath::Sqrt(fInvMass);    
}


/////////////////////////////////////////////////////////////////////////////////

AliGenKinkInfo::AliGenKinkInfo():
  fMCd(),       //info about daughter particle - second particle for V0
  fMCm(),       //info about mother particle   - first particle for V0
  fMCDist1(0),  //info about closest distance according closest MC - linear DCA
  fMCDist2(0),  //info about closest distance parabolic DCA
  fMCRr(0),     // rec position of the vertex 
  fMCR(0)       //exact r position of the vertex
{
  //
  // default constructor
  //
  for (Int_t i=0;i<3;i++){
    fMCPdr[i]=0;
    fMCPd[i]=0;
    fMCX[i]=0;
    fMCPm[i]=0;
    fMCAngle[i]=0;
  }
  for (Int_t i=0; i<2; i++) {
    fPdg[i]= 0;
    fLab[i]=0;
  }
}

void AliGenKinkInfo::Update()
{
  //
  // Update information
  //  some redundancy - faster acces to the values in analysis code
  //
  fMCPd[0] = fMCd.GetParticle().Px();
  fMCPd[1] = fMCd.GetParticle().Py();
  fMCPd[2] = fMCd.GetParticle().Pz();
  fMCPd[3] = fMCd.GetParticle().P();
  //
  fMCX[0]  = fMCd.GetParticle().Vx();
  fMCX[1]  = fMCd.GetParticle().Vy();
  fMCX[2]  = fMCd.GetParticle().Vz();
  fMCR       = TMath::Sqrt( fMCX[0]*fMCX[0]+fMCX[1]*fMCX[1]);
  //
  fPdg[0]    = fMCd.GetParticle().GetPdgCode();
  fPdg[1]    = fMCm.GetParticle().GetPdgCode();
  //
  fLab[0]    = fMCd.GetParticle().GetUniqueID();
  fLab[1]    = fMCm.GetParticle().GetUniqueID();
  //
  //
  //
  Double_t x1[3],p1[3];
  TParticle& pdaughter = fMCd.GetParticle();
  x1[0] = pdaughter.Vx();      
  x1[1] = pdaughter.Vy();
  x1[2] = pdaughter.Vz();
  p1[0] = pdaughter.Px();
  p1[1] = pdaughter.Py();
  p1[2] = pdaughter.Pz();
  Double_t sign = (pdaughter.GetPDG()->Charge()>0)? -1:1;
  AliHelix dhelix1(x1,p1,sign);
  //
  //
  Double_t x2[3],p2[3];            
  //
  TParticle& pmother = fMCm.GetParticle();
  x2[0] = pmother.Vx();      
  x2[1] = pmother.Vy();      
  x2[2] = pmother.Vz();      
  p2[0] = pmother.Px();
  p2[1] = pmother.Py();
  p2[2] = pmother.Pz();
  //
  const AliTrackReference & pdecay = fMCm.GetTRdecay();
  x2[0] = pdecay.X();      
  x2[1] = pdecay.Y();      
  x2[2] = pdecay.Z();      
  p2[0] = pdecay.Px();
  p2[1] = pdecay.Py();
  p2[2] = pdecay.Pz();
  //
  sign = (pmother.GetPDG()->Charge()>0) ? -1:1;
  AliHelix mhelix(x2,p2,sign);
  
  //
  //
  //
  //find intersection linear
  //
  Double_t distance1, distance2;
  Double_t phase[2][2],radius[2];
  Int_t  points = dhelix1.GetRPHIintersections(mhelix, phase, radius);
  Double_t delta1=10000,delta2=10000;    
  if (points>0){
    dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
    dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
    dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
  }
  if (points==2){    
    dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
    dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
    dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
  }
  distance1 = TMath::Min(delta1,delta2);
  //
  //find intersection parabolic
  //
  points = dhelix1.GetRPHIintersections(mhelix, phase, radius);
  delta1=10000,delta2=10000;  
  
  if (points>0){
    dhelix1.ParabolicDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
  }
  if (points==2){    
    dhelix1.ParabolicDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
  }
  
  distance2 = TMath::Min(delta1,delta2);
  //
  if (delta1<delta2){
    //get V0 info
    dhelix1.Evaluate(phase[0][0],fMCXr);
    dhelix1.GetMomentum(phase[0][0],fMCPdr);
    mhelix.GetMomentum(phase[0][1],fMCPm);
    dhelix1.GetAngle(phase[0][0],mhelix,phase[0][1],fMCAngle);
    fMCRr = TMath::Sqrt(radius[0]);
  }
  else{
    dhelix1.Evaluate(phase[1][0],fMCXr);
    dhelix1.GetMomentum(phase[1][0],fMCPdr);
    mhelix.GetMomentum(phase[1][1],fMCPm);
    dhelix1.GetAngle(phase[1][0],mhelix,phase[1][1],fMCAngle);
    fMCRr = TMath::Sqrt(radius[1]);
  }     
  //            
  //
  fMCDist1 = TMath::Sqrt(distance1);
  fMCDist2 = TMath::Sqrt(distance2);      
    
}


Float_t AliGenKinkInfo::GetQt(){
  //
  //
  Float_t momentumd = TMath::Sqrt(fMCPd[0]*fMCPd[0]+fMCPd[1]*fMCPd[1]+fMCPd[2]*fMCPd[2]);
  return TMath::Sin(fMCAngle[2])*momentumd;
}


////////////////////////////////////////////////////////////////////////
AliTPCdigitRow::AliTPCdigitRow()
{
  Reset();
}
////////////////////////////////////////////////////////////////////////
AliTPCdigitRow & AliTPCdigitRow::operator=(const AliTPCdigitRow &digOld)
{
  for (Int_t i = 0; i<kgRowBytes; i++) fDig[i] = digOld.fDig[i];
  return (*this);
}
////////////////////////////////////////////////////////////////////////
void AliTPCdigitRow::SetRow(Int_t row) 
{
  if (row >= 8*kgRowBytes) {
    cerr<<"AliTPCdigitRow::SetRow: index "<<row<<" out of bounds."<<endl;
    return;
  }
  Int_t iC = row/8;
  Int_t iB = row%8;
  SETBIT(fDig[iC],iB);
}

////////////////////////////////////////////////////////////////////////
Bool_t AliTPCdigitRow::TestRow(Int_t row) const
{
//
// return kTRUE if row is on
//
  Int_t iC = row/8;
  Int_t iB = row%8;
  return TESTBIT(fDig[iC],iB);
}
////////////////////////////////////////////////////////////////////////
Int_t AliTPCdigitRow::RowsOn(Int_t upto) const
{
//
// returns number of rows with a digit  
// count only rows less equal row number upto
//
  Int_t total = 0;
  for (Int_t i = 0; i<kgRowBytes; i++) {
    for (Int_t j = 0; j < 8; j++) {
      if (i*8+j > upto) return total;
      if (TESTBIT(fDig[i],j))  total++;
    }
  }
  return total;
}
////////////////////////////////////////////////////////////////////////
void AliTPCdigitRow::Reset()
{
//
// resets all rows to zero
//
  for (Int_t i = 0; i<kgRowBytes; i++) {
    fDig[i] <<= 8;
  }
}
////////////////////////////////////////////////////////////////////////
Int_t AliTPCdigitRow::Last() const
{
//
// returns the last row number with a digit
// returns -1 if now digits 
//
  for (Int_t i = kgRowBytes-1; i>=0; i--) {
    for (Int_t j = 7; j >= 0; j--) {
      if TESTBIT(fDig[i],j) return i*8+j;
    }
  }
  return -1;
}
////////////////////////////////////////////////////////////////////////
Int_t AliTPCdigitRow::First() const
{
//
// returns the first row number with a digit
// returns -1 if now digits 
//
  for (Int_t i = 0; i<kgRowBytes; i++) {
    for (Int_t j = 0; j < 8; j++) {
      if (TESTBIT(fDig[i],j)) return i*8+j;
    }
  }
  return -1;
}


//_____________________________________________________________________________
Float_t AliMCInfo::TPCBetheBloch(Float_t bg)
{
  //
  // Bethe-Bloch energy loss formula
  //
  const Double_t kp1=0.76176e-1;
  const Double_t kp2=10.632;
  const Double_t kp3=0.13279e-4;
  const Double_t kp4=1.8631;
  const Double_t kp5=1.9479;

  Double_t dbg = (Double_t) bg;

  Double_t beta = dbg/TMath::Sqrt(1.+dbg*dbg);

  Double_t aa = TMath::Power(beta,kp4);
  Double_t bb = TMath::Power(1./dbg,kp5);

  bb=TMath::Log(kp3+bb);
  
  return ((Float_t)((kp2-aa-bb)*kp1/aa));
}
Commit	Line	Data
c92725b7	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
d92975ba	21	Container classes with MC infomation
c92725b7	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"
c92725b7	36	#include "TStopwatch.h"
	37	#include "TParticle.h"
	38	#include "TSystem.h"
c92725b7	39	#include "TCanvas.h"
c92725b7	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"
c92725b7	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)
c92725b7	67
	68
	69
	70	////////////////////////////////////////////////////////////////////////
022044bf	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)
c92725b7	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
022044bf	106	AliMCInfo::AliMCInfo(const AliMCInfo& info):
	107	TObject(),
	108	fTrackRef(info.fTrackRef),
	109	fTrackRefOut(info.fTrackRefOut),
d92975ba	110	fTRdecay(info.GetTRdecay()),
022044bf	111	fPrimPart(info.fPrimPart),
022044bf	112	fParticle(info.fParticle),
d92975ba	113	fMass(info.GetMass()),
022044bf	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
c92725b7	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	/////////////////////////////////////////////////////////////////////////////////
022044bf	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
c92725b7	247	{
022044bf	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	}
c92725b7	261	}
c92725b7	262
	263	void AliGenV0Info::Update(Float_t vertex[3])
	264	{
022044bf	265	//
	266	// Update information - calculates derived variables
	267	//
	268
d92975ba	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();
c92725b7	277	fMCR = TMath::Sqrt( fMCX[0]fMCX[0]+fMCX[1]fMCX[1]);
c92725b7	278	//
d92975ba	279	fPdg[0] = fMCd.GetParticle().GetPdgCode();
d92975ba	280	fPdg[1] = fMCm.GetParticle().GetPdgCode();
c92725b7	281	//
d92975ba	282	fLab[0] = fMCd.GetParticle().GetUniqueID();
d92975ba	283	fLab[1] = fMCm.GetParticle().GetUniqueID();
c92725b7	284	//
	285	//
	286	//
	287	Double_t x1[3],p1[3];
d92975ba	288	TParticle& pdaughter = fMCd.GetParticle();
c92725b7	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	//
d92975ba	301	TParticle& pmother = fMCm.GetParticle();
c92725b7	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]+vnorm2pnorm2)/(vnorm3*pnorm3);
386	fPointAngle = (v[0]p[0]+v[1]p[1]+v[2]p[2])/(vnorm3pnorm3);
d92975ba	387	Double_t mass1 = fMCd.GetMass();
d92975ba	388	Double_t mass2 = fMCm.GetMass();
c92725b7	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;
022044bf	407	if (fInvMass>0) fInvMass = TMath::Sqrt(fInvMass);
c92725b7	408	}
	409
	410
	411
	412	/////////////////////////////////////////////////////////////////////////////////
022044bf	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
c92725b7	438	void AliGenKinkInfo::Update()
c92725b7	439	{
022044bf	440	//
	441	// Update information
	442	// some redundancy - faster acces to the values in analysis code
	443	//
d92975ba	444	fMCPd[0] = fMCd.GetParticle().Px();
	445	fMCPd[1] = fMCd.GetParticle().Py();
	446	fMCPd[2] = fMCd.GetParticle().Pz();
	447	fMCPd[3] = fMCd.GetParticle().P();
c92725b7	448	//
d92975ba	449	fMCX[0] = fMCd.GetParticle().Vx();
	450	fMCX[1] = fMCd.GetParticle().Vy();
	451	fMCX[2] = fMCd.GetParticle().Vz();
c92725b7	452	fMCR = TMath::Sqrt( fMCX[0]fMCX[0]+fMCX[1]fMCX[1]);
c92725b7	453	//
d92975ba	454	fPdg[0] = fMCd.GetParticle().GetPdgCode();
d92975ba	455	fPdg[1] = fMCm.GetParticle().GetPdgCode();
c92725b7	456	//
d92975ba	457	fLab[0] = fMCd.GetParticle().GetUniqueID();
d92975ba	458	fLab[1] = fMCm.GetParticle().GetUniqueID();
c92725b7	459	//
	460	//
	461	//
	462	Double_t x1[3],p1[3];
d92975ba	463	TParticle& pdaughter = fMCd.GetParticle();
c92725b7	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	//
d92975ba	476	TParticle& pmother = fMCm.GetParticle();
c92725b7	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	//
d92975ba	484	const AliTrackReference & pdecay = fMCm.GetTRdecay();
c92725b7	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
c92725b7	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	////////////////////////////////////////////////////////////////////////
022044bf	587	Bool_t AliTPCdigitRow::TestRow(Int_t row) const
c92725b7	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	////////////////////////////////////////////////////////////////////////
022044bf	597	Int_t AliTPCdigitRow::RowsOn(Int_t upto) const
c92725b7	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	////////////////////////////////////////////////////////////////////////
022044bf	623	Int_t AliTPCdigitRow::Last() const
c92725b7	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	////////////////////////////////////////////////////////////////////////
022044bf	637	Int_t AliTPCdigitRow::First() const
c92725b7	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
022044bf	651
d92975ba	652	//_____________________________________________________________________________
d92975ba	653	Float_t AliMCInfo::TPCBetheBloch(Float_t bg)
c92725b7	654	{
c92725b7	655	//
d92975ba	656	// Bethe-Bloch energy loss formula
c92725b7	657	//
d92975ba	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;
c92725b7	663
d92975ba	664	Double_t dbg = (Double_t) bg;
c92725b7	665
d92975ba	666	Double_t beta = dbg/TMath::Sqrt(1.+dbg*dbg);
c92725b7	667
d92975ba	668	Double_t aa = TMath::Power(beta,kp4);
d92975ba	669	Double_t bb = TMath::Power(1./dbg,kp5);
c92725b7	670
d92975ba	671	bb=TMath::Log(kp3+bb);
c92725b7	672
d92975ba	673	return ((Float_t)((kp2-aa-bb)*kp1/aa));
c92725b7	674	}
c92725b7	675