HLT/ITS/RunHLTITS.C

   1 // The following macro runs the HLT ITS tracker over the HLT
   2 // tracks stored in the ESD and stores the output tracks in a
   3 // separate ESD file AliESDits.root
   4
   5 #if !defined(__CINT__) || defined(__MAKECINT__)
   6   #include <TMath.h>
   7   #include <TError.h>
   8   #include <Riostream.h>
   9   #include <TH1F.h>
  10   #include <TH2F.h>
  11   #include <TNtuple.h>
  12   #include <TProfile2D.h>
  13   #include <TF1.h>
  14   #include <TGraphErrors.h>
  15   #include <TTree.h>
  16   #include <TParticle.h>
  17   #include <TCanvas.h>
  18   #include <TFile.h>
  19   #include <TROOT.h>
  20   #include <TStopwatch.h>
  21   #include <TSystem.h>
  22
  23   #include "AliStack.h"
  24   #include "AliHeader.h"
  25   #include "AliTrackReference.h"
  26   #include "AliRunLoader.h"
  27   #include "AliITS.h"
  28   #include "AliITSLoader.h"
  29   #include "AliITSgeom.h"
  30   #include "AliITStrackerV2.h"
  31   #include "AliRun.h"
  32   #include "AliESD.h"
  33   #include "AliGenEventHeader.h"
  34
  35   #include "AliHLTITStrack.h"
  36   #include "AliHLTITStracker.h"
  37   #include "AliHLTITSVertexerZ.h"
  38 #endif
  39
  40 //extern TSystem *gSystem;
  41
  42 Int_t RunHLTITS(Int_t nev=1,Int_t run=0) {
  43
  44   //  gSystem->Load("libAliHLTITS.so");
  45
  46   TStopwatch timer;
  47   timer.Start();
  48
  49    if (gAlice) {
  50       delete gAlice->GetRunLoader();
  51       delete gAlice;
  52       gAlice=0;
  53    }
  54
  55    AliRunLoader *rl = AliRunLoader::Open("galice.root");
  56    if (rl == 0x0) {
  57       cerr<<"Can not open session"<<endl;
  58       return 1;
  59    }
  60    Int_t retval = rl->LoadgAlice();
  61    if (retval) {
  62       cerr<<"AliESDtest.C : LoadgAlice returned error"<<endl;
  63       delete rl;
  64       return 1;
  65    }
  66    retval = rl->LoadHeader();
  67    if (retval) {
  68       cerr<<"AliESDtest.C : LoadHeader returned error"<<endl;
  69       delete rl;
  70       return 2;
  71    }
  72    gAlice=rl->GetAliRun();
  73
  74    AliTracker::SetFieldMap(gAlice->Field());
  75
  76    AliITSLoader* itsl = (AliITSLoader*)rl->GetLoader("ITSLoader");
  77    if (itsl == 0x0) {
  78       cerr<<"AliESDtest.C : Can not get the ITS loader"<<endl;
  79       return 3;
  80    }
  81    itsl->LoadRecPoints("read");
  82
  83    AliITS *dITS = (AliITS*)gAlice->GetDetector("ITS");
  84    if (!dITS) {
  85       cerr<<"AliESDtest.C : Can not find the ITS detector !"<<endl;
  86       return 4;
  87    }
  88    //   AliITSgeom *geom = dITS->GetITSgeom();
  89    AliITSgeom *geom = new AliITSgeom();
  90    geom->ReadNewFile("$ALICE_ROOT/ITS/ITSgeometry_vPPRasymmFMD.det");
  91
  92    //An instance of the HLT ITS tracker
  93    AliHLTITStracker itsTracker(geom);
  94
  95    TFile *ef=TFile::Open("AliESDs.root");
  96    if (!ef || !ef->IsOpen()) {cerr<<"Can't AliESDs.root !\n"; return 1;}
  97    AliESD* event = new AliESD;
  98    TTree* tree = (TTree*) ef->Get("esdTree");
  99    if (!tree) {cerr<<"no ESD tree found\n"; return 1;};
 100    tree->SetBranchAddress("ESD", &event);
 101
 102    TFile *itsf=TFile::Open("AliESDits.root","RECREATE");
 103    if ((!itsf)||(!itsf->IsOpen())) {
 104       cerr<<"Can't AliESDits.root !\n"; return 1;
 105    }
 106
 107    Int_t rc=0;
 108    if (nev>rl->GetNumberOfEvents()) nev=rl->GetNumberOfEvents();
 109    //The loop over events
 110    for (Int_t i=0; i<nev; i++) {
 111
 112      cerr<<"\n\nProcessing event number : "<<i<<endl;
 113      tree->GetEvent(i);
 114      rl->GetEvent(i);
 115
 116      TArrayF v(3);
 117      rl->GetHeader()->GenEventHeader()->PrimaryVertex(v);
 118      Double_t vtx[3]={v[0],v[1],v[2]};
 119      Double_t cvtx[3]={0.005,0.005,0.010};
 120      cout<<"MC vertex position: "<<v[2]<<endl;
 121
 122      AliHLTITSVertexerZ vertexer("null");
 123      AliESDVertex* vertex = NULL;
 124      TStopwatch timer2;
 125      timer2.Start();
 126      TTree* treeClusters = itsl->TreeR();
 127      //     vertex = vertexer.FindVertexForCurrentEvent(i);
 128      //     AliESDVertex *vertex = vertexer.FindVertexForCurrentEvent(geom,treeClusters);
 129      vertex = new AliESDVertex(vtx,cvtx);
 130      timer2.Stop();
 131      timer2.Print();
 132      if(!vertex){
 133        cerr<<"Vertex not found"<<endl;
 134        vertex = new AliESDVertex(vtx,cvtx);
 135      }
 136      else {
 137        vertex->SetTruePos(vtx);  // store also the vertex from MC
 138      }
 139
 140      event->SetVertex(vertex);
 141
 142      Double_t vtxPos[3];
 143      Double_t vtxErr[3];
 144      vertex->GetXYZ(vtxPos);
 145      vertex->GetSigmaXYZ(vtxErr);
 146      itsTracker.SetVertex(vtxPos,vtxErr);
 147
 148      TTree *itsTree=itsl->TreeR();
 149      if (!itsTree) {
 150        cerr<<"Can't get the ITS cluster tree !\n";
 151        return 4;
 152      }
 153      itsTracker.LoadClusters(itsTree);
 154      rc+=itsTracker.Clusters2Tracks(event);
 155      //     rc+=itsTracker.PropagateBack(event);
 156      itsTracker.UnloadClusters();
 157
 158      if (rc==0) {
 159        TTree* tree = new TTree("esdTree", "Tree with ESD objects");
 160        tree->Branch("ESD", "AliESD", &event);
 161        tree->Fill();
 162        itsf->cd();
 163        tree->Write();
 164      }
 165      if (rc) {
 166        cerr<<"Something bad happened...\n";
 167      }
 168
 169    }
 170    delete event;
 171
 172    itsf->Close();
 173    ef->Close();
 174
 175    //   delete rl;
 176
 177    timer.Stop();
 178    timer.Print();
 179
 180    return rc;
 181 }