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