8ea894e6cd23125b2bc8b59be04ba89a7f10fc8b
[u/mrichter/AliRoot.git] / 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 }