Update master to aliroot
[u/mrichter/AliRoot.git] / HLT / ITS / RunHLTITS.C
CommitLineData
9582ea1a 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
df7af99f 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"
df7af99f 33 #include "AliGenEventHeader.h"
34
4aa41877 35 #include "AliHLTITStrack.h"
36 #include "AliHLTITStracker.h"
37 #include "AliHLTITSVertexerZ.h"
df7af99f 38#endif
39
40//extern TSystem *gSystem;
41
9582ea1a 42Int_t RunHLTITS(Int_t nev=1,Int_t run=0) {
43
4070f709 44 // gSystem->Load("libAliHLTITS");
df7af99f 45
9582ea1a 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
df7af99f 74 AliTracker::SetFieldMap(gAlice->Field());
9582ea1a 75
9582ea1a 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 }
df7af99f 88 // AliITSgeom *geom = dITS->GetITSgeom();
89 AliITSgeom *geom = new AliITSgeom();
90 geom->ReadNewFile("$ALICE_ROOT/ITS/ITSgeometry_vPPRasymmFMD.det");
9582ea1a 91
92 //An instance of the HLT ITS tracker
4aa41877 93 AliHLTITStracker itsTracker(geom);
9582ea1a 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);
9582ea1a 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};
df7af99f 120 cout<<"MC vertex position: "<<v[2]<<endl;
121
4aa41877 122 AliHLTITSVertexerZ vertexer("null");
df7af99f 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);
9582ea1a 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}