HBTAN/AliHBTReaderTPC.cxx

   1 #include "AliHBTReaderTPC.h"
   2
   3 #include <TTree.h>
   4 #include <TFile.h>
   5 #include <TParticle.h>
   6
   7 #include <Varargs.h>
   8
   9 #include <AliRun.h>
  10 #include <AliLoader.h>
  11 #include <AliStack.h>
  12 #include <AliMagF.h>
  13 #include <AliTPCtrack.h>
  14 #include <AliTPCParam.h>
  15 #include <AliTPCtracker.h>
  16 #include <AliTPCLoader.h>
  17
  18 #include "AliHBTRun.h"
  19 #include "AliHBTEvent.h"
  20 #include "AliHBTParticle.h"
  21 #include "AliHBTParticleCut.h"
  22
  23 ClassImp(AliHBTReaderTPC)
  24 //______________________________________________
  25 //
  26 // class AliHBTReaderTPC
  27 //
  28 //reader for TPC tracking
  29 //needs galice.root, AliTPCtracks.root, AliTPCclusters.root
  30 //
  31 //more info: http://alisoft.cern.ch/people/skowron/analyzer/index.html
  32 //Piotr.Skowronski@cern.ch
  33 AliHBTReaderTPC::AliHBTReaderTPC():
  34  fFileName("galice.root"),
  35  fRunLoader(0x0),
  36  fTPCLoader(0x0),
  37  fMagneticField(0.0),
  38  fUseMagFFromRun(kTRUE)
  39 {
  40   //constructor,
  41   //Defaults:
  42   //  galicefilename = ""  - this means: Do not open gAlice file -
  43   //                         just leave the global pointer untouched
  44
  45 }
  46
  47 AliHBTReaderTPC::AliHBTReaderTPC(const Char_t* galicefilename):
  48  fFileName(galicefilename),
  49  fRunLoader(0x0),
  50  fTPCLoader(0x0),
  51  fMagneticField(0.0),
  52  fUseMagFFromRun(kTRUE)
  53 {
  54   //constructor,
  55   //Defaults:
  56   //  galicefilename = ""  - this means: Do not open gAlice file -
  57   //                         just leave the global pointer untouched
  58
  59 }
  60 /********************************************************************/
  61 AliHBTReaderTPC::AliHBTReaderTPC(TObjArray* dirs, const Char_t* galicefilename):
  62  AliHBTReader(dirs),
  63  fFileName(galicefilename),
  64  fRunLoader(0x0),
  65  fTPCLoader(0x0),
  66  fMagneticField(0.0),
  67  fUseMagFFromRun(kTRUE)
  68 {
  69   //constructor,
  70   //Defaults:
  71   //  galicefilename = ""  - this means: Do not open gAlice file -
  72   //                         just leave the global pointer untached
  73
  74 }
  75 /********************************************************************/
  76
  77 AliHBTReaderTPC::~AliHBTReaderTPC()
  78 {
  79  //desctructor
  80    delete fRunLoader;
  81 }
  82 /********************************************************************/
  83
  84 void AliHBTReaderTPC::Rewind()
  85 {
  86   delete fRunLoader;
  87   fRunLoader = 0x0;
  88   fCurrentDir = 0;
  89   fNEventsRead= 0;
  90 }
  91 /********************************************************************/
  92
  93 Int_t AliHBTReaderTPC::ReadNext()
  94  {
  95  //reads data and puts put to the particles and tracks objects
  96  //reurns 0 if everything is OK
  97  //
  98   Info("Read","");
  99
 100
 101   TObjArray *tarray = new TObjArray(5000); //cotainer for tpc tracks
 102   tarray->SetOwner(); //set the ownership of the objects it contains
 103                       //when array is is deleted or cleared all objects
 104                       //that it contains are deleted
 105
 106   if (fParticlesEvent == 0x0)  fParticlesEvent = new AliHBTEvent();
 107   if (fTracksEvent == 0x0)  fTracksEvent = new AliHBTEvent();
 108
 109   fParticlesEvent->Reset();
 110   fTracksEvent->Reset();
 111
 112   do  //do{}while; is OK even if 0 dirs specified. In that case we try to read from "./"
 113    {
 114
 115     if (fRunLoader == 0x0)
 116       if (OpenNextSession()) continue;//directory counter is increased inside in case of error
 117
 118     if (fCurrentEvent == fRunLoader->GetNumberOfEvents())
 119      {
 120        //read next directory
 121        delete fRunLoader;//close current session
 122        fRunLoader = 0x0;//assure pointer is null
 123        fCurrentDir++;//go to next dir
 124        continue;//directory counter is increased inside in case of error
 125      }
 126
 127     Info("ReadNext","Reading Event %d",fCurrentEvent);
 128
 129     fRunLoader->GetEvent(fCurrentEvent);
 130     TTree *tracktree = fTPCLoader->TreeT();//get the tree
 131     if (!tracktree) //check if we got the tree
 132       {//if not return with error
 133          Error("ReadNext","Can't get a tree with TPC tracks !\n");
 134          continue;
 135       }
 136
 137     TBranch *trackbranch=tracktree->GetBranch("tracks");//get the branch with tracks
 138     if (!trackbranch) ////check if we got the branch
 139       {//if not return with error
 140         Error("ReadNext","Can't get a branch with TPC tracks !\n");
 141         continue;
 142       }
 143     Int_t ntpctracks=(Int_t)tracktree->GetEntries();//get number of TPC tracks
 144     Info("ReadNext","Found %d TPC tracks.",ntpctracks);
 145     //Copy tracks to array
 146
 147     AliTPCtrack *iotrack=0;
 148
 149     for (Int_t i=0; i<ntpctracks; i++) //loop over all tpc tracks
 150      {
 151        iotrack=new AliTPCtrack;   //create new tracks
 152        trackbranch->SetAddress(&iotrack); //tell the branch ehere to put track data from tree(file)
 153        tracktree->GetEvent(i); //stream track i to the iotrack
 154        tarray->AddLast(iotrack); //put the track in the array
 155      }
 156
 157     Double_t xk;
 158     Double_t par[5];
 159     Float_t phi, lam, pt;//angles and transverse momentum
 160     Int_t label; //label of the current track
 161
 162     AliStack* stack = fRunLoader->Stack();
 163     if (stack == 0x0)
 164      {
 165        Error("ReadNext","Can not get stack for current event",fCurrentEvent);
 166        fCurrentEvent++;
 167        continue;
 168      }
 169     stack->Particles();
 170
 171     for (Int_t i=0; i<ntpctracks; i++) //loop over all good tracks
 172      {
 173        iotrack = (AliTPCtrack*)tarray->At(i);
 174        label = iotrack->GetLabel();
 175
 176        if (label < 0) continue;
 177
 178        TParticle *p = (TParticle*)stack->Particle(label);
 179
 180        if(p == 0x0) continue; //if returned pointer is NULL
 181        if(p->GetPDG() == 0x0) continue; //if particle has crezy PDG code (not known to our database)
 182
 183        if(Pass(p->GetPdgCode())) continue; //check if we are intersted with particles of this type
 184                                    //if not take next partilce
 185
 186        AliHBTParticle* part = new AliHBTParticle(*p,i);
 187        if(Pass(part)) { delete part; continue;}//check if meets all criteria of any of our cuts
 188                                                //if it does not delete it and take next good track
 189
 190 //       iotrack->PropagateTo(3.,0.0028,65.19);
 191 //       iotrack->PropagateToVertex(36.66,1.2e-3);//orig
 192        iotrack->PropagateToVertex(50.,0.0353);
 193
 194        iotrack->GetExternalParameters(xk,par);     //get properties of the track
 195        phi=TMath::ASin(par[2]) + iotrack->GetAlpha();
 196        if (phi<-TMath::Pi()) phi+=2*TMath::Pi();
 197        if (phi>=TMath::Pi()) phi-=2*TMath::Pi();
 198        lam=par[3];
 199        pt=1.0/TMath::Abs(par[4]);
 200
 201        Double_t tpx = pt * TMath::Cos(phi); //track x coordinate of momentum
 202        Double_t tpy = pt * TMath::Sin(phi); //track y coordinate of momentum
 203        Double_t tpz = pt * lam; //track z coordinate of momentum
 204
 205        Double_t mass = p->GetMass();
 206        Double_t tEtot = TMath::Sqrt( tpx*tpx + tpy*tpy + tpz*tpz + mass*mass);//total energy of the track
 207
 208        AliHBTParticle* track = new AliHBTParticle(p->GetPdgCode(),i, tpx, tpy , tpz, tEtot, 0., 0., 0., 0.);
 209        if(Pass(track))//check if meets all criteria of any of our cuts
 210                     //if it does not delete it and take next good track
 211         {
 212           delete track;
 213           delete part;
 214           continue;
 215         }
 216        fParticlesEvent->AddParticle(part);
 217        fTracksEvent->AddParticle(track);
 218       }
 219
 220     Info("ReadNext","Read %d tracks and %d particles from event %d (event %d in dir %d).",
 221             fParticlesEvent->GetNumberOfParticles(), fTracksEvent->GetNumberOfParticles(),
 222             fNEventsRead,fCurrentEvent,fCurrentDir);
 223
 224     fCurrentEvent++;
 225     fNEventsRead++;
 226     delete tarray;
 227     return 0;
 228    }while(fCurrentDir < GetNumberOfDirs());
 229
 230   delete tarray;
 231   return 1;
 232  }
 233 /********************************************************************/
 234
 235 Int_t AliHBTReaderTPC::OpenNextSession()
 236 {
 237   TString filename = GetDirName(fCurrentDir);
 238   if (filename.IsNull())
 239    {
 240      DoOpenError("Can not get directory name");
 241      return 1;
 242    }
 243   filename = filename +"/"+ fFileName;
 244   fRunLoader = AliRunLoader::Open(filename,AliConfig::fgkDefaultEventFolderName);
 245   if( fRunLoader == 0x0)
 246    {
 247      DoOpenError("Can not open session.");
 248      return 1;
 249    }
 250
 251   fRunLoader->LoadHeader();
 252   if (fRunLoader->GetNumberOfEvents() <= 0)
 253    {
 254      DoOpenError("There is no events in this directory.");
 255      return 1;
 256    }
 257
 258   if (fRunLoader->LoadKinematics())
 259    {
 260      DoOpenError("Error occured while loading kinematics.");
 261      return 1;
 262    }
 263   fTPCLoader = (AliTPCLoader*)fRunLoader->GetLoader("TPCLoader");
 264   if ( fTPCLoader == 0x0)
 265    {
 266      DoOpenError("Exiting due to problems with opening files.");
 267      return 1;
 268    }
 269   Info("OpenNextSession","________________________________________________________");
 270   Info("OpenNextSession","Found %d event(s) in directory %s",
 271         fRunLoader->GetNumberOfEvents(),GetDirName(fCurrentDir).Data());
 272   Float_t mf;
 273   if (fUseMagFFromRun)
 274    {
 275      fRunLoader->LoadgAlice();
 276      mf = fRunLoader->GetAliRun()->Field()->SolenoidField();
 277      Info("OpenNextSession","Setting Magnetic Field from run: B=%fT",mf/10.);
 278      fRunLoader->UnloadgAlice();
 279    }
 280   else
 281    {
 282      Info("OpenNextSession","Setting Own Magnetic Field: B=%fT",fMagneticField);
 283      if (fMagneticField == 0x0)
 284       {
 285         Fatal("OpenNextSession","Magnetic field can not be 0.");
 286         return 1;//pro forma
 287       }
 288      mf = fMagneticField*10.;
 289    }
 290   AliKalmanTrack::SetConvConst(1000/0.299792458/mf);
 291
 292
 293   fRunLoader->CdGAFile();
 294   AliTPCParam *TPCParam= (AliTPCParam*)gDirectory->Get("75x40_100x60");
 295   if (!TPCParam)
 296    {
 297     TPCParam=(AliTPCParam *)gDirectory->Get("75x40_100x60_150x60");
 298     if (!TPCParam)
 299      {
 300        DoOpenError("TPC parameters have not been found !\n");
 301        return 1;
 302      }
 303    }
 304
 305   if (fTPCLoader->LoadTracks())
 306    {
 307      DoOpenError("Error occured while loading TPC tracks.");
 308      return 1;
 309    }
 310
 311   fCurrentEvent = 0;
 312   return 0;
 313 }
 314 /********************************************************************/
 315
 316 void AliHBTReaderTPC::DoOpenError( const char *va_(fmt), ...)
 317 {
 318   // Does error display and clean-up in case error caught on Open Next Session
 319
 320    va_list ap;
 321    va_start(ap,va_(fmt));
 322    Error("OpenNextSession", va_(fmt), ap);
 323    va_end(ap);
 324
 325    delete fRunLoader;
 326    fRunLoader = 0x0;
 327    fTPCLoader = 0x0;
 328    fCurrentDir++;
 329 }
 330
 331 /********************************************************************/
 332 /********************************************************************/
 333 /********************************************************************/
 334