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