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