Adding AliTPCTracklet to the repository (M.Mager)
[u/mrichter/AliRoot.git] / TPC / AliBarrelRec_TPCparam.C
CommitLineData
a654e257 1void AliBarrelRec_TPCparam(Int_t firstEvent=0,Int_t lastEvent=0) {
056891c0 2 //
a654e257 3 // Macro to create AliESDs.root using parametrized TPC tracking
970d08c3 4 // and AliITStrackerSA (MI + SA)
056891c0 5 //
a654e257 6 // Input files:
7 // - galice.root
8 // - Kinematics.root
9 // - TrackReferences.root
10 // - ITS.RecPoints.root (use AliRecontruction class)
11 // - ITS.Vertex.root (use $ALICE_ROOT/ITS/AliITSVertexerZTest.C)
e0cf7c8a 12 //
970d08c3 13 // A. Dainese - INFN Legnaro
e0cf7c8a 14 //
056891c0 15
64cf164e 16
970d08c3 17 Int_t collcode = 1; // pp collisions
225a9242 18 Bool_t useMeanVtx = kFALSE;
970d08c3 19
64cf164e 20 AliGeomManager::LoadGeometry("geometry.root");
970d08c3 21
22 if (gAlice) {
23 delete gAlice->GetRunLoader();
24 delete gAlice;
25 gAlice=0;
26 }
27 AliRunLoader *rl = AliRunLoader::Open("galice.root");
28 if (rl == 0x0) {
29 cerr<<"Can not open session"<<endl;
30 return;
31 }
32 Int_t retval = rl->LoadgAlice();
33 if (retval) {
34 cerr<<"LoadgAlice returned error"<<endl;
35 delete rl;
36 return;
37 }
38 retval = rl->LoadHeader();
39 if (retval) {
40 cerr<<"LoadHeader returned error"<<endl;
41 delete rl;
42 return;
43 }
44 gAlice=rl->GetAliRun();
45
46
47
48 // Get field from galice.root
49 AliMagF *fiel = (AliMagF*)gAlice->Field();
50 Double_t fieval=TMath::Abs((Double_t)fiel->SolenoidField()/10.);
51 // Set the conversion constant between curvature and Pt
52 AliTracker::SetFieldMap(fiel,kTRUE);
53
54 /**** The TPC corner ********************/
55
56 AliTPCtrackerParam tpcTrackerPar(collcode,fieval);
57 tpcTrackerPar.Init();
58
59 //**** Switch on the PID class (mandatory!)
60 AliPID pid;
61
62 /**** The ITS corner ********************/
63
64 AliITSLoader* itsl = (AliITSLoader*)rl->GetLoader("ITSLoader");
65 if (itsl == 0x0) {
66 cerr<<"Cannot get the ITS loader"<<endl;
67 return;
68 }
69 itsl->LoadRecPoints("read");
70
64cf164e 71 AliITSRecoParam * itsRecoParam = AliITSRecoParam::GetLowFluxParam();
72 AliITSReconstructor::SetRecoParam(itsRecoParam);
970d08c3 73
74 // Instance of the ITS tracker
64cf164e 75 AliITStrackerSA itsTracker(0);
970d08c3 76 Int_t ITSclusters[6] = {1,1,1,1,1,1};
64cf164e 77 itsTracker.SetLayersNotToSkip(ITSclusters);
970d08c3 78
79 // Primary vertex reconstruction in pp
80 AliESDVertex *initVertex = 0;
81 TFile *invtx = new TFile("AliESDVertexMean.root");
82 if(collcode==1 && useMeanVtx) {
83 // open the mean vertex
84 initVertex = (AliESDVertex*)invtx->Get("vtxmean");
85 } else {
86 Double_t pos[3]={0.5,0.5,0.};
87 Double_t err[3]={3.,3.,5.};
88 initVertex = new AliESDVertex(pos,err);
89 }
90 invtx->Close();
91 delete invtx;
64cf164e 92 AliVertexerTracks *vertexer = new AliVertexerTracks(AliTracker::GetBz());
970d08c3 93 vertexer->SetVtxStart(initVertex);
94 vertexer->SetDebug(0);
95 delete initVertex;
96 initVertex=0;
97
98 /***** The TREE for ESD is born *****/
99 TTree *esdTree=new TTree("esdTree","Tree with ESD objects");
64cf164e 100 AliESDEvent *event=0; AliESDEvent *eventTPCin=0;
101 event = new AliESDEvent();
102 event->CreateStdContent();
103 event->WriteToTree(esdTree);
970d08c3 104
105 if(firstEvent>rl->GetNumberOfEvents()) firstEvent=rl->GetNumberOfEvents()-1;
106 if(lastEvent>rl->GetNumberOfEvents()) lastEvent=rl->GetNumberOfEvents()-1;
107 cout<<" Number of events: "<<1+lastEvent-firstEvent<<endl;
108
109 TFile *ppZ = TFile::Open("ITS.Vertex.root"); // z vertices from SPD
110 AliESDVertex *vertexSPD = new AliESDVertex();
111 Char_t zver[100];
112 Double_t vtx[3]={0,0,0};
113 Double_t sigmavtx[3]={0.07,0.07,0.1};
114
115
116 //<---------------------------------- The Loop over events begins
117 TStopwatch timer;
118 Int_t trc;
119 for(Int_t i=firstEvent; i<=lastEvent; i++) {
120
225a9242 121 cout<<" Processing event number : "<<i<<endl;
64cf164e 122 //AliESDEvent *event = new AliESDEvent();
970d08c3 123 event->SetRunNumber(gAlice->GetRunNumber());
64cf164e 124 event->SetEventNumberInFile(i);
970d08c3 125 event->SetMagneticField(gAlice->Field()->SolenoidField());
126 rl->GetEvent(i);
127
128 //***** Primary vertex from SPD from file
129 sprintf(zver,"Event%d/Vertex",i);
130 vertexSPD = (AliESDVertex*)ppZ->Get(zver);
131 if(!vertexSPD) {
64cf164e 132 esdTree->Fill(); event->Reset();
970d08c3 133 continue;
134 }
135 event->SetVertex(vertexSPD);
136 vertexSPD->GetXYZ(vtx);
137 vertexSPD->GetSigmaXYZ(sigmavtx);
138
139 //***** TPC tracking
140 if ( (trc=tpcTrackerPar.BuildTPCtracks(event)) ) {
141 printf("exiting TPC tracker with code %d in event %d\n",trc,i);
64cf164e 142 esdTree->Fill(); event->Reset();
970d08c3 143 continue;
144 }
a654e257 145
225a9242 146 // make a copy of the ESD at this stage
147 eventTPCin = event;
148
970d08c3 149 //***** ITS tracking
150 itsTracker.AliTracker::SetVertex(vtx,sigmavtx);
64cf164e 151 // itsl->LoadRecPoints("read");
970d08c3 152 TTree *itsTree=itsl->TreeR();
153 if (!itsTree) {
154 cerr<<"Can't get the ITS cluster tree !\n";
64cf164e 155 esdTree->Fill(); event->Reset();
a654e257 156 return;
970d08c3 157 }
158 itsTracker.UnloadClusters();
159 itsTracker.LoadClusters(itsTree);
160 if ( (trc=itsTracker.Clusters2Tracks(event)) ) {
161 printf("exiting ITS tracker with code %d in event %d\n",trc,i);
64cf164e 162 esdTree->Fill(); event->Reset();
970d08c3 163 continue;
164 }
225a9242 165
166 // Bring kTPCin-tracks back to the TPC inner wall
167 BackToTPCInnerWall(event,eventTPCin);
168
169 // refit inward in ITS:
170 // - refit without vertex constraint
171 // - propagate through beam pipe to local x = 0
172 itsTracker.RefitInward(event);
970d08c3 173
970d08c3 174
175 //***** Vertex from ESD tracks
176 if(collcode==1) { // pp
177 AliESDVertex *vertexTrks =
178 (AliESDVertex*)vertexer->FindPrimaryVertex(event);
179 event->SetPrimaryVertex(vertexTrks);
180 }
181
182 esdTree->Fill();
64cf164e 183 event->Reset();
970d08c3 184
185 }//<-----------------------------------The Loop over events ends here
186 timer.Stop(); timer.Print();
187
188 // The AliESDs.root is born
189 TFile *ef = TFile::Open("AliESDs.root","RECREATE");
190 if (!ef || !ef->IsOpen()) {cerr<<"Can't open AliESDs.root !\n"; return;}
191
225a9242 192 //Write the tree and close everything
970d08c3 193 esdTree->Write();
194 delete esdTree;
195 ef->Close();
196
197 delete vertexer;
198 delete rl;
199
200 return;
056891c0 201}
4201c420 202//--------------------------------------------------------------------------
64cf164e 203void BackToTPCInnerWall(AliESDEvent *event,AliESDEvent *eventTPC) {
225a9242 204
205 Int_t ntracks = eventTPC->GetNumberOfTracks();
206 AliESDtrack *esdTrackTPC = 0;
4201c420 207
225a9242 208 // create relation between event and eventTPC
209 Int_t labelsTPC[100000000];
210 for(Int_t tr = 0; tr<ntracks; tr++) {
211 esdTrackTPC = (AliESDtrack*)event->GetTrack(tr);
212 labelsTPC[TMath::Abs(esdTrackTPC->GetLabel())] = tr;
213 }
214
215 ntracks = event->GetNumberOfTracks();
4201c420 216 AliESDtrack *esdTrack = 0;
225a9242 217 esdTrackTPC = 0;
218 Int_t indexTPC;
4201c420 219
220 // loop on tracks
225a9242 221 for(tr = 0; tr<ntracks; tr++) {
222 esdTrack = (AliESDtrack*)event->GetTrack(tr);
223 // set to kITSout the tracks that don't have kTPCin
224 // (they've been found by AliITStrackerSA)
225 if(!(esdTrack->GetStatus()&AliESDtrack::kTPCin)) {
226 esdTrack->SetStatus(AliESDtrack::kITSout);
227 continue;
228 }
229
230 // skip tracks that don't have kITSin
4201c420 231 if(!(esdTrack->GetStatus()&AliESDtrack::kITSin)) continue;
225a9242 232
233 indexTPC = labelsTPC[TMath::Abs(esdTrack->GetLabel())];
234 esdTrackTPC = (AliESDtrack*)eventTPC->GetTrack(indexTPC);
235
236 AliITStrackMI *itsTrack = 0;
4201c420 237 try {
225a9242 238 itsTrack = new AliITStrackMI(*esdTrackTPC);
4201c420 239 esdTrack = 0;
240 }
241 catch (const Char_t *msg) {
225a9242 242 Warning("ToTPCInnerWall",msg);
4201c420 243 continue;
244 }
225a9242 245 itsTrack->UpdateESDtrack(AliESDtrack::kITSout);
4201c420 246 esdTrack = new AliESDtrack(*(itsTrack->GetESDtrack()));
225a9242 247
4201c420 248 delete itsTrack;
249 } // end loop on tracks
250
970d08c3 251 return;
4201c420 252}
225a9242 253//--------------------------------------------------------------------------