HBTAN/AliHBTReaderITSv2.cxx

   1 #include "AliHBTReaderITSv2.h"
   2
   3 #include <Varargs.h>
   4
   5 #include <TString.h>
   6 #include <TObjString.h>
   7 #include <TTree.h>
   8 #include <TParticle.h>
   9
  10 #include <AliRun.h>
  11 #include <AliRunLoader.h>
  12 #include <AliLoader.h>
  13 #include <AliStack.h>
  14 #include <AliMagF.h>
  15 #include <AliITStrackV2.h>
  16
  17 #include "AliHBTRun.h"
  18 #include "AliHBTEvent.h"
  19 #include "AliHBTParticle.h"
  20 #include "AliHBTParticleCut.h"
  21
  22
  23 ClassImp(AliHBTReaderITSv2)
  24
  25 AliHBTReaderITSv2::AliHBTReaderITSv2():
  26  fFileName("galice.root"),
  27  fRunLoader(0x0),
  28  fITSLoader(0x0),
  29  fMagneticField(0.0),
  30  fUseMagFFromRun(kTRUE)
  31 {
  32   //constructor,
  33   //Defaults:
  34   //  galicefilename = "galice.root"
  35 }
  36 /********************************************************************/
  37
  38 AliHBTReaderITSv2::AliHBTReaderITSv2(const Char_t* galicefilename):
  39  fFileName(galicefilename),
  40  fRunLoader(0x0),
  41  fITSLoader(0x0),
  42  fMagneticField(0.0),
  43  fUseMagFFromRun(kTRUE)
  44 {
  45   //constructor,
  46   //Defaults:
  47   //  galicefilename = "galice.root"
  48 }
  49 /********************************************************************/
  50
  51 AliHBTReaderITSv2::AliHBTReaderITSv2(TObjArray* dirs, const Char_t* galicefilename):
  52  AliHBTReader(dirs),
  53  fFileName(galicefilename),
  54  fRunLoader(0x0),
  55  fITSLoader(0x0),
  56  fMagneticField(0.0),
  57  fUseMagFFromRun(kTRUE)
  58 {
  59   //constructor,
  60   //Defaults:
  61   //  galicefilename = "galice.root"
  62
  63 }
  64 /********************************************************************/
  65
  66 void AliHBTReaderITSv2::Rewind()
  67 {
  68   //rewinds reading
  69   delete fRunLoader;
  70   fRunLoader = 0x0;
  71   fCurrentDir = 0;
  72   fNEventsRead= 0;
  73 }
  74 /********************************************************************/
  75
  76 AliHBTReaderITSv2::~AliHBTReaderITSv2()
  77 {
  78   //dtor
  79   delete fRunLoader;
  80 }
  81 /********************************************************************/
  82
  83 Int_t AliHBTReaderITSv2::ReadNext()
  84 {
  85 //reads data from next event
  86
  87  register Int_t i = 0; //iterator
  88
  89 // AliITStrackerV2 *tracker; // ITS tracker - used for cooking labels
  90  TTree *tracktree = 0x0; // tree for tracks
  91  AliITStrackV2 *iotrack = 0x0;
  92
  93  Double_t xk;
  94  Double_t par[5]; //Kalman track parameters
  95  Float_t phi, lam, pt;//angles and transverse momentum
  96  Int_t label; //label of the current track
  97
  98
  99  if (fParticlesEvent == 0x0)  fParticlesEvent = new AliHBTEvent();
 100  if (fTracksEvent == 0x0)  fTracksEvent = new AliHBTEvent();
 101
 102  fParticlesEvent->Reset();
 103  fTracksEvent->Reset();
 104  do //do while is good even if Ndirs==0 (than read from current directory)
 105    {
 106     if (fRunLoader == 0x0)
 107       if (OpenNextFile()) continue;//directory counter is increased inside in case of error
 108
 109     if (fCurrentEvent == fRunLoader->GetNumberOfEvents())
 110      {
 111        //read next directory
 112        delete fRunLoader;//close current session
 113        fRunLoader = 0x0;//assure pointer is null
 114        fCurrentDir++;//go to next dir
 115        continue;//directory counter is increased inside in case of error
 116      }
 117
 118     Info("ReadNext","Reading Event %d",fCurrentEvent);
 119
 120     fRunLoader->GetEvent(fCurrentEvent);
 121
 122     tracktree=fITSLoader->TreeT();
 123     if (!tracktree)
 124      {
 125        Error("ReadNext","Can't get a tree with ITS tracks");
 126        fCurrentEvent++;
 127        continue;
 128      }
 129
 130     TBranch *tbranch=tracktree->GetBranch("tracks");
 131     if (!tbranch)
 132      {
 133        Error("ReadNext","Can't get a branch with ITS tracks");
 134        fCurrentEvent++;
 135        continue;
 136      }
 137
 138     AliStack* stack = fRunLoader->Stack();
 139     if (stack == 0x0)
 140      {
 141        Error("ReadNext","Can not get stack for current event",fCurrentEvent);
 142        fCurrentEvent++;
 143        continue;
 144      }
 145
 146     //must be here because on the beginning conv. const. is not set yet
 147     if (iotrack == 0x0) iotrack = new AliITStrackV2(); //buffer track for reading data from tree
 148
 149     Int_t ntr = (Int_t)tracktree->GetEntries();
 150
 151     for (i=0; i < ntr; i++) //loop over all tpc tracks
 152      {
 153        tbranch->SetAddress(&iotrack);
 154        tracktree->GetEvent(i);
 155
 156        label=iotrack->GetLabel();
 157        if (label < 0)
 158         {
 159           continue;
 160         }
 161
 162        TParticle *p = stack->Particle(label);
 163        if(p == 0x0) continue; //if returned pointer is NULL
 164        if(p->GetPDG() == 0x0) continue; //if particle has crezy PDG code (not known to our database)
 165
 166        if(Pass(p->GetPdgCode())) continue; //check if we are intersted with particles of this type
 167                                            //if not take next partilce
 168
 169        AliHBTParticle* part = new AliHBTParticle(*p,i);
 170        if(Pass(part)) { delete part; continue;}//check if meets all criteria of any of our cuts
 171                                                //if it does not delete it and take next good track
 172
 173        iotrack->PropagateTo(3.,0.0028,65.19);
 174        iotrack->PropagateToVertex();
 175
 176        iotrack->GetExternalParameters(xk,par);     //get properties of the track
 177        phi=TMath::ASin(par[2]) + iotrack->GetAlpha();
 178        if (phi<-TMath::Pi()) phi+=2*TMath::Pi();
 179        if (phi>=TMath::Pi()) phi-=2*TMath::Pi();
 180        lam=par[3];
 181        pt=1.0/TMath::Abs(par[4]);
 182
 183        Double_t tpx = pt * TMath::Cos(phi); //track x coordinate of momentum
 184        Double_t tpy = pt * TMath::Sin(phi); //track y coordinate of momentum
 185        Double_t tpz = pt * lam; //track z coordinate of momentum
 186
 187        Double_t mass = p->GetMass();
 188        Double_t tEtot = TMath::Sqrt( tpx*tpx + tpy*tpy + tpz*tpz + mass*mass);//total energy of the track
 189
 190        AliHBTParticle* track = new AliHBTParticle(p->GetPdgCode(), i, tpx, tpy , tpz, tEtot, 0., 0., 0., 0.);
 191        if(Pass(track))//check if meets all criteria of any of our cuts
 192                       //if it does not delete it and take next good track
 193         {
 194          delete track;
 195          delete part;
 196          continue;
 197         }
 198
 199        fParticlesEvent->AddParticle(part);
 200        fTracksEvent->AddParticle(track);
 201      }//end of loop over tracks in the event
 202
 203     Info("ReadNext","Read %d tracks and %d particles from event %d (event %d in dir %d).",
 204             fParticlesEvent->GetNumberOfParticles(), fTracksEvent->GetNumberOfParticles(),
 205             fNEventsRead,fCurrentEvent,fCurrentDir);
 206
 207     fCurrentEvent++;
 208     fNEventsRead++;
 209     delete iotrack;
 210     return 0;
 211    }while(fCurrentDir < GetNumberOfDirs());//end of loop over directories specified in fDirs Obj Array
 212
 213  delete iotrack;
 214  return 1;
 215 }
 216
 217 /********************************************************************/
 218 Int_t AliHBTReaderITSv2::OpenNextFile()
 219 {
 220   //opens next file
 221   TString filename = GetDirName(fCurrentDir);
 222   if (filename.IsNull())
 223    {
 224      DoOpenError("Can not get directory name");
 225      return 1;
 226    }
 227   filename = filename +"/"+ fFileName;
 228   fRunLoader = AliRunLoader::Open(filename,AliConfig::GetDefaultEventFolderName());
 229   if( fRunLoader == 0x0)
 230    {
 231      DoOpenError("Can not open session.");
 232      return 1;
 233    }
 234
 235   if (fRunLoader->GetNumberOfEvents() <= 0)
 236    {
 237      DoOpenError("There is no events in this directory.");
 238      return 1;
 239    }
 240
 241   if (fRunLoader->LoadKinematics())
 242    {
 243      DoOpenError("Error occured while loading kinematics.");
 244      return 1;
 245    }
 246   fITSLoader = fRunLoader->GetLoader("ITSLoader");
 247   if ( fITSLoader == 0x0)
 248    {
 249      DoOpenError("Exiting due to problems with opening files.");
 250      return 1;
 251    }
 252
 253   Info("OpenNextSession","________________________________________________________");
 254   Info("OpenNextSession","Found %d event(s) in directory %s",
 255         fRunLoader->GetNumberOfEvents(),GetDirName(fCurrentDir).Data());
 256   Float_t mf;
 257   if (fUseMagFFromRun)
 258    {
 259      if (fRunLoader->LoadgAlice())
 260       {
 261         DoOpenError("Error occured while loading AliRun.");
 262         return 1;
 263       }
 264      mf = fRunLoader->GetAliRun()->Field()->SolenoidField();
 265      Info("OpenNextSession","Setting Magnetic Field from run: B=%fT",mf/10.);
 266      fRunLoader->UnloadgAlice();
 267    }
 268   else
 269    {
 270      Info("OpenNextSession","Setting Own Magnetic Field: B=%fT",fMagneticField);
 271      if (fMagneticField == 0x0)
 272       {
 273         Fatal("OpenNextSession","Magnetic field can not be 0.");
 274         return 1;//pro forma
 275       }
 276      mf = fMagneticField*10.;
 277    }
 278   AliKalmanTrack::SetConvConst(1000/0.299792458/mf);
 279
 280   if (fITSLoader->LoadTracks())
 281    {
 282      DoOpenError("Error occured while loading TPC tracks.");
 283      return 1;
 284    }
 285
 286   fCurrentEvent = 0;
 287   return 0;
 288 }
 289 /********************************************************************/
 290
 291 void AliHBTReaderITSv2::DoOpenError( const char *va_(fmt), ...)
 292 {
 293   // Does error display and clean-up in case error caught on Open Next Session
 294
 295    va_list ap;
 296    va_start(ap,va_(fmt));
 297    Error("OpenNextFile", va_(fmt), ap);
 298    va_end(ap);
 299
 300    delete fRunLoader;
 301    fRunLoader = 0x0;
 302    fITSLoader = 0x0;
 303    fCurrentDir++;
 304 }
 305
 306 /********************************************************************/