STEER/AliESDtest.C

   1 //********************************************************************
   2 //     Example of the reconstruction that generates the ESD
   3 // Input files:
   4 //   a) file containing the ITS clusters
   5 //      (the AliITSFindClustersV2.C macro can be used to generate it)
   6 //   b) file containing the TPC clusters
   7 //      (the AliTPCFindClusters.C macro can be used to generate it)
   8 //   c) file containing the TRD clusters
   9 //      (the AliTRDdigits2cluster.C macro can be used to generate it)
  10 //   d) file containing the TOF digits
  11 //      (the AliTOFSDigits2Digits.C macro can be used to generate it)
  12 // Ouput file:
  13 //      AliESDs.root containing the ESD events
  14 //
  15 // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
  16 //********************************************************************
  17
  18 #if !defined(__CINT__) || defined(__MAKECINT__)
  19   #include <Riostream.h>
  20   #include "TFile.h"
  21   #include "TSystem.h"
  22   #include "TStopwatch.h"
  23   #include "TGeant3.h"
  24   #include "TArrayF.h"
  25
  26   #include "AliMagF.h"
  27   #include "AliRun.h"
  28   #include "AliRunLoader.h"
  29   #include "AliLoader.h"
  30   #include "AliHeader.h"
  31   #include "AliGenEventHeader.h"
  32
  33   #include "AliESD.h"
  34   #include "AliESDpid.h"
  35
  36   #include "AliITS.h"
  37   #include "AliITSgeom.h"
  38   #include "AliITStrackerV2.h"
  39   #include "AliV0vertexer.h"
  40   #include "AliCascadeVertexer.h"
  41   #include "AliITSpidESD.h"
  42   #include "AliITSLoader.h"
  43
  44   #include "AliTPCParam.h"
  45   #include "AliTPCtracker.h"
  46   #include "AliTPCpidESD.h"
  47   #include "AliTPCLoader.h"
  48
  49   #include "AliTRDtracker.h"
  50   #include "AliTRDPartID.h"
  51
  52   #include "AliTOFpidESD.h"
  53 #endif
  54
  55 extern TSystem *gSystem;
  56 extern AliRun *gAlice;
  57 extern TFile *gFile;
  58
  59 Int_t AliESDtest(Int_t nev=1,Int_t run=0) {
  60
  61 /**** Initialization of the NewIO *******/
  62
  63    if (gAlice) {
  64       delete gAlice->GetRunLoader();
  65       delete gAlice;
  66       gAlice=0;
  67    }
  68
  69    gSystem->Load("libgeant321");     // needed for the PID in TOF
  70    new TGeant3("");                  // must be re-done !
  71
  72    AliRunLoader *rl = AliRunLoader::Open("galice.root");
  73    if (rl == 0x0) {
  74       cerr<<"Can not open session"<<endl;
  75       return 1;
  76    }
  77    Int_t retval = rl->LoadgAlice();
  78    if (retval) {
  79       cerr<<"AliESDtest.C : LoadgAlice returned error"<<endl;
  80       delete rl;
  81       return 1;
  82    }
  83    retval = rl->LoadHeader();
  84    if (retval) {
  85       cerr<<"AliESDtest.C : LoadHeader returned error"<<endl;
  86       delete rl;
  87       return 2;
  88    }
  89    gAlice=rl->GetAliRun();
  90
  91
  92    AliKalmanTrack::SetConvConst(
  93       1000/0.299792458/gAlice->Field()->SolenoidField()
  94    );
  95
  96
  97 /**** The ITS corner ********************/
  98
  99    AliITSLoader* itsl = (AliITSLoader*)rl->GetLoader("ITSLoader");
 100    if (itsl == 0x0) {
 101       cerr<<"AliESDtest.C : Can not get the ITS loader"<<endl;
 102       return 3;
 103    }
 104    itsl->LoadRecPoints("read");
 105
 106    AliITS *dITS = (AliITS*)gAlice->GetDetector("ITS");
 107    if (!dITS) {
 108       cerr<<"AliESDtest.C : Can not find the ITS detector !"<<endl;
 109       return 4;
 110    }
 111    AliITSgeom *geom = dITS->GetITSgeom();
 112
 113    //An instance of the ITS tracker
 114    AliITStrackerV2 itsTracker(geom);
 115
 116    //An instance of the ITS PID maker
 117    Double_t parITS[]={34.,0.15,10.};
 118    AliITSpidESD itsPID(parITS);
 119
 120    //An instance of the V0 finder
 121    Double_t cuts[]={33,  // max. allowed chi2
 122                     0.16,// min. allowed negative daughter's impact parameter
 123                     0.05,// min. allowed positive daughter's impact parameter
 124                     0.080,// max. allowed DCA between the daughter tracks
 125                     0.998,// max. allowed cosine of V0's pointing angle
 126                     0.9,  // min. radius of the fiducial volume
 127                     2.9   // max. radius of the fiducial volume
 128                    };
 129    AliV0vertexer vtxer(cuts);
 130
 131    Double_t cts[]={33.,    // max. allowed chi2
 132                     0.05,   // min. allowed V0 impact parameter
 133                     0.008,  // window around the Lambda mass
 134                     0.035,  // min. allowed bachelor's impact parameter
 135                     0.10,   // max. allowed DCA between a V0 and a track
 136                     0.9985, //max. allowed cosine of the cascade pointing angle
 137                     0.9,    // min. radius of the fiducial volume
 138                     2.9     // max. radius of the fiducial volume
 139                    };
 140    AliCascadeVertexer cvtxer=AliCascadeVertexer(cts);
 141
 142 /**** The TPC corner ********************/
 143
 144    AliTPCLoader* tpcl = (AliTPCLoader*)rl->GetLoader("TPCLoader");
 145    if (tpcl == 0x0) {
 146       cerr<<"AliESDtest.C : can not get the TPC loader"<<endl;
 147       return 5;
 148    }
 149    tpcl->LoadRecPoints("read");
 150
 151    rl->CdGAFile();
 152    AliTPCParam *par=(AliTPCParam *)gDirectory->Get("75x40_100x60_150x60");
 153    if (!par) {
 154       cerr<<"TPC parameters have not been found !\n";
 155       return 6;
 156    }
 157
 158    //An instance of the TPC tracker
 159    AliTPCtracker tpcTracker(par);
 160
 161    //An instance of the TPC PID maker
 162    Double_t parTPC[]={47.,0.10,10.};
 163    AliTPCpidESD tpcPID(parTPC);
 164
 165
 166 /**** The TRD corner ********************/
 167
 168    AliLoader* trdl = rl->GetLoader("TRDLoader");
 169    if (trdl == 0x0) {
 170       cerr<<"AliESDtest.C : can not get the TRD loader"<<endl;
 171       return 5;
 172    }
 173    trdl->LoadRecPoints("read");
 174
 175    //An instance of the TRD tracker
 176    rl->CdGAFile();
 177    AliTRDtracker trdTracker(gFile);  //galice.root file
 178
 179 /*
 180    //An instance of the TRD PID maker
 181    TFile* pidFile = TFile::Open("pid.root");
 182    if (!pidFile->IsOpen()) {
 183      cerr << "Can't get pid.root !\n";
 184      return 7;
 185    }
 186    AliTRDPartID* trdPID = (AliTRDPartID*) pidFile->Get("AliTRDPartID");
 187    if (!trdPID) {
 188      cerr << "Can't get PID object !\n";
 189      return 8;
 190    }
 191 */
 192
 193
 194 /**** The TOF corner ********************/
 195
 196    AliLoader* tofl = rl->GetLoader("TOFLoader");
 197    if (tofl == 0x0) {
 198       cerr<<"AliESDtest.C : can not get the TOF loader"<<endl;
 199       return 5;
 200    }
 201    tofl->LoadDigits("read");
 202
 203
 204    //Instance of the TOF PID maker
 205    Double_t parTOF[]={130.,5.};
 206    AliTOFpidESD tofPID(parTOF);
 207
 208
 209    //rl->UnloadgAlice();
 210
 211
 212    TFile *ef=TFile::Open("AliESDs.root","RECREATE");
 213    if (!ef->IsOpen()) {cerr<<"Can't AliESDs.root !\n"; return 1;}
 214
 215    TStopwatch timer;
 216    Int_t rc=0;
 217    if (nev>rl->GetNumberOfEvents()) nev=rl->GetNumberOfEvents();
 218    //The loop over events
 219    for (Int_t i=0; i<nev; i++) {
 220      cerr<<"\n\nProcessing event number : "<<i<<endl;
 221      AliESD *event=new AliESD();
 222      event->SetRunNumber(run);
 223      event->SetEventNumber(i);
 224
 225      rl->GetEvent(i);
 226
 227 //***** Primary vertex reconstruction (MC vertex position, for the moment)
 228      TArrayF v(3);
 229      rl->GetHeader()->GenEventHeader()->PrimaryVertex(v);
 230      Double_t vtx[3]={v[0],v[1],v[2]};
 231      Double_t cvtx[6]={
 232        0.005,
 233        0.000, 0.005,
 234        0.000, 0.000, 0.010
 235      };
 236      event->SetVertex(vtx,cvtx);
 237      cvtx[1]=cvtx[0]; cvtx[2]=cvtx[5]; //trackers use only the diag.elements
 238
 239 //***** Initial path towards the primary vertex
 240      tpcTracker.SetVertex(vtx,cvtx);
 241      TTree *tpcTree=tpcl->TreeR();
 242      if (!tpcTree) {
 243         cerr<<"Can't get the TPC cluster tree !\n";
 244         return 4;
 245      }
 246      tpcTracker.LoadClusters(tpcTree);
 247      rc+=tpcTracker.Clusters2Tracks(event);
 248
 249      itsTracker.SetVertex(vtx,cvtx);
 250      TTree *itsTree=itsl->TreeR();
 251      if (!itsTree) {
 252         cerr<<"Can't get the ITS cluster tree !\n";
 253         return 4;
 254      }
 255      itsTracker.LoadClusters(itsTree);
 256      rc+=itsTracker.Clusters2Tracks(event);
 257
 258
 259 //***** Back propagation towards the outer barrel detectors
 260      rc+=itsTracker.PropagateBack(event);
 261      //itsPID.MakePID(event);
 262
 263      rc+=tpcTracker.PropagateBack(event);
 264      tpcPID.MakePID(event);
 265
 266      TTree *trdTree=trdl->TreeR();
 267      if (!trdTree) {
 268         cerr<<"Can't get the TRD cluster tree !\n";
 269         return 4;
 270      }
 271      trdTracker.LoadClusters(trdTree);
 272      rc+=trdTracker.PropagateBack(event);
 273 /*
 274      for (Int_t iTrack = 0; iTrack < event->GetNumberOfTracks(); iTrack++) {
 275        AliESDtrack* track = event->GetTrack(iTrack);
 276        trdPID->MakePID(track);
 277      }
 278 */
 279      TTree *tofTree=tofl->TreeD();
 280      if (!tofTree) {
 281         cerr<<"Can't get the TOF cluster tree !\n";
 282         return 4;
 283      }
 284      tofPID.LoadClusters(tofTree);
 285      tofPID.MakePID(event);
 286      tofPID.UnloadClusters();
 287
 288
 289
 290 //***** Here is the combined PID
 291      AliESDpid::MakePID(event);
 292
 293
 294
 295 //***** Now the final refit at the primary vertex...
 296      rc+=trdTracker.RefitInward(event);
 297      trdTracker.UnloadClusters();
 298
 299      rc+=tpcTracker.RefitInward(event);
 300      tpcTracker.UnloadClusters();
 301
 302      rc+=itsTracker.RefitInward(event);
 303      itsTracker.UnloadClusters();
 304
 305
 306
 307 //***** Hyperon reconstruction
 308      vtxer.SetVertex(vtx);
 309      rc+=vtxer.Tracks2V0vertices(event);            // V0 finding
 310      rc+=cvtxer.V0sTracks2CascadeVertices(event);   // cascade finding
 311
 312
 313
 314 //***** Some final manipulations with this event
 315      if (rc==0) {
 316         Char_t ename[100];
 317         sprintf(ename,"%d",i);
 318         ef->cd();
 319         if (!event->Write(ename)) rc++;
 320      }
 321      if (rc) {
 322         cerr<<"Something bad happened...\n";
 323      }
 324      delete event;
 325    }
 326    timer.Stop(); timer.Print();
 327
 328    //pidFile->Close();
 329
 330    delete par;
 331
 332    ef->Close();
 333
 334    delete rl;
 335
 336    return rc;
 337 }