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 "AliL3ITStrack.h"
  37   #include "AliL3ITStracker.h"
  38   #include "AliL3ITSVertexerZ.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("libAliL3ITS.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    AliKalmanTrack::SetConvConst(
  78       1000/0.299792458/gAlice->Field()->SolenoidField()
  79    );
  80
  81    AliITSLoader* itsl = (AliITSLoader*)rl->GetLoader("ITSLoader");
  82    if (itsl == 0x0) {
  83       cerr<<"AliESDtest.C : Can not get the ITS loader"<<endl;
  84       return 3;
  85    }
  86    itsl->LoadRecPoints("read");
  87
  88    AliITS *dITS = (AliITS*)gAlice->GetDetector("ITS");
  89    if (!dITS) {
  90       cerr<<"AliESDtest.C : Can not find the ITS detector !"<<endl;
  91       return 4;
  92    }
  93    //   AliITSgeom *geom = dITS->GetITSgeom();
  94    AliITSgeom *geom = new AliITSgeom();
  95    geom->ReadNewFile("$ALICE_ROOT/ITS/ITSgeometry_vPPRasymmFMD.det");
  96
  97    //An instance of the HLT ITS tracker
  98    AliL3ITStracker itsTracker(geom);
  99
 100    TFile *ef=TFile::Open("AliESDs.root");
 101    if (!ef || !ef->IsOpen()) {cerr<<"Can't AliESDs.root !\n"; return 1;}
 102    AliESD* event = new AliESD;
 103    TTree* tree = (TTree*) ef->Get("esdTree");
 104    if (!tree) {cerr<<"no ESD tree found\n"; return 1;};
 105    tree->SetBranchAddress("ESD", &event);
 106
 107    TFile *itsf=TFile::Open("AliESDits.root","RECREATE");
 108    if ((!itsf)||(!itsf->IsOpen())) {
 109       cerr<<"Can't AliESDits.root !\n"; return 1;
 110    }
 111
 112    Int_t rc=0;
 113    if (nev>rl->GetNumberOfEvents()) nev=rl->GetNumberOfEvents();
 114    //The loop over events
 115    for (Int_t i=0; i<nev; i++) {
 116
 117      cerr<<"\n\nProcessing event number : "<<i<<endl;
 118      tree->GetEvent(i);
 119      rl->GetEvent(i);
 120
 121      TArrayF v(3);
 122      rl->GetHeader()->GenEventHeader()->PrimaryVertex(v);
 123      Double_t vtx[3]={v[0],v[1],v[2]};
 124      Double_t cvtx[3]={0.005,0.005,0.010};
 125      cout<<"MC vertex position: "<<v[2]<<endl;
 126
 127      AliL3ITSVertexerZ vertexer("null");
 128      AliESDVertex* vertex = NULL;
 129      TStopwatch timer2;
 130      timer2.Start();
 131      TTree* treeClusters = itsl->TreeR();
 132      //     vertex = vertexer.FindVertexForCurrentEvent(i);
 133      //     AliESDVertex *vertex = vertexer.FindVertexForCurrentEvent(geom,treeClusters);
 134      vertex = new AliESDVertex(vtx,cvtx);
 135      timer2.Stop();
 136      timer2.Print();
 137      if(!vertex){
 138        cerr<<"Vertex not found"<<endl;
 139        vertex = new AliESDVertex(vtx,cvtx);
 140      }
 141      else {
 142        vertex->SetTruePos(vtx);  // store also the vertex from MC
 143      }
 144
 145      event->SetVertex(vertex);
 146
 147      Double_t vtxPos[3];
 148      Double_t vtxErr[3];
 149      vertex->GetXYZ(vtxPos);
 150      vertex->GetSigmaXYZ(vtxErr);
 151      itsTracker.SetVertex(vtxPos,vtxErr);
 152
 153      TTree *itsTree=itsl->TreeR();
 154      if (!itsTree) {
 155        cerr<<"Can't get the ITS cluster tree !\n";
 156        return 4;
 157      }
 158      itsTracker.LoadClusters(itsTree);
 159      rc+=itsTracker.Clusters2Tracks(event);
 160      //     rc+=itsTracker.PropagateBack(event);
 161      itsTracker.UnloadClusters();
 162
 163      if (rc==0) {
 164        TTree* tree = new TTree("esdTree", "Tree with ESD objects");
 165        tree->Branch("ESD", "AliESD", &event);
 166        tree->Fill();
 167        itsf->cd();
 168        tree->Write();
 169      }
 170      if (rc) {
 171        cerr<<"Something bad happened...\n";
 172      }
 173
 174    }
 175    delete event;
 176
 177    itsf->Close();
 178    ef->Close();
 179
 180    //   delete rl;
 181
 182    timer.Stop();
 183    timer.Print();
 184
 185    return rc;
 186 }