HBTAN/AliHBTReaderTPC.cxx

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