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