Test for Coverity
[u/mrichter/AliRoot.git] / TPC / AliBarrelRec_TPCparam.C
1 void AliBarrelRec_TPCparam(Int_t firstEvent=0,Int_t lastEvent=0) {
2   //
3   // Macro to create AliESDs.root using parametrized TPC tracking
4   // and AliITStrackerSA (MI + SA)
5   //
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)
12   //
13   // A. Dainese - INFN Legnaro
14   //
15
16
17   Int_t  collcode = 1; // pp collisions
18   Bool_t useMeanVtx = kFALSE;
19   
20   AliGeomManager::LoadGeometry("geometry.root");
21   
22   if (gAlice) {
23     delete AliRunLoader::Instance();
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   
71   AliITSRecoParam * itsRecoParam = AliITSRecoParam::GetLowFluxParam();
72   AliITSReconstructor::SetRecoParam(itsRecoParam);
73   
74   // Instance of the ITS tracker
75   AliITStrackerSA itsTracker(0);
76   Int_t ITSclusters[6] = {1,1,1,1,1,1};
77   itsTracker.SetLayersNotToSkip(ITSclusters);  
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;
92   AliVertexerTracks *vertexer = new AliVertexerTracks(AliTracker::GetBz());
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");
100   AliESDEvent *event=0; AliESDEvent *eventTPCin=0;
101   event = new AliESDEvent();
102   event->CreateStdContent();
103   event->WriteToTree(esdTree);
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     
121     cout<<" Processing event number : "<<i<<endl;
122     //AliESDEvent *event = new AliESDEvent(); 
123     event->SetRunNumber(gAlice->GetRunNumber());
124     event->SetEventNumberInFile(i);
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) {
132       esdTree->Fill(); event->Reset();
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);
142       esdTree->Fill(); event->Reset();
143       continue;
144     }
145
146     // make a copy of the ESD at this stage
147     eventTPCin = event;
148
149     //***** ITS tracking
150     itsTracker.AliTracker::SetVertex(vtx,sigmavtx);
151     //    itsl->LoadRecPoints("read");
152     TTree *itsTree=itsl->TreeR();
153     if (!itsTree) {
154       cerr<<"Can't get the ITS cluster tree !\n";
155       esdTree->Fill(); event->Reset();
156       return;
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);
162       esdTree->Fill(); event->Reset();
163       continue;
164     }
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);
173     
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();
183     event->Reset();
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
192   //Write the tree and close everything
193   esdTree->Write();
194   delete esdTree;
195   ef->Close();
196
197   delete vertexer;
198   delete rl;
199
200   return;
201 }
202 //--------------------------------------------------------------------------
203 void BackToTPCInnerWall(AliESDEvent *event,AliESDEvent *eventTPC) {
204
205   Int_t ntracks = eventTPC->GetNumberOfTracks();
206   AliESDtrack *esdTrackTPC = 0;
207
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();
216   AliESDtrack *esdTrack = 0;
217   esdTrackTPC = 0;
218   Int_t indexTPC;
219
220   // loop on tracks
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
231     if(!(esdTrack->GetStatus()&AliESDtrack::kITSin)) continue;
232
233     indexTPC = labelsTPC[TMath::Abs(esdTrack->GetLabel())];
234     esdTrackTPC = (AliESDtrack*)eventTPC->GetTrack(indexTPC);
235
236     AliITStrackMI *itsTrack = 0;
237     try {
238       itsTrack = new AliITStrackMI(*esdTrackTPC);
239       esdTrack = 0;
240     }
241     catch (const Char_t *msg) {
242         Warning("ToTPCInnerWall",msg);
243         continue;
244     }
245     itsTrack->UpdateESDtrack(AliESDtrack::kITSout);
246     esdTrack = new AliESDtrack(*(itsTrack->GetESDtrack()));
247
248     delete itsTrack;
249   } // end loop on tracks
250
251   return;
252 }
253 //--------------------------------------------------------------------------