Modifying the GetChainFromCollection function based on the additions of the TEntryList
[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   Int_t  collcode = 1; // pp collisions
17   Bool_t useMeanVtx = kFALSE;
18   
19   
20   if (gAlice) {
21     delete gAlice->GetRunLoader();
22     delete gAlice; 
23     gAlice=0;
24   }
25   AliRunLoader *rl = AliRunLoader::Open("galice.root");
26   if (rl == 0x0) {
27     cerr<<"Can not open session"<<endl;
28     return;
29   }
30   Int_t retval = rl->LoadgAlice();
31   if (retval) {
32     cerr<<"LoadgAlice returned error"<<endl;
33     delete rl;
34     return;
35   }
36   retval = rl->LoadHeader();
37   if (retval) {
38     cerr<<"LoadHeader returned error"<<endl;
39     delete rl;
40     return;
41   }
42   gAlice=rl->GetAliRun();
43   
44   
45   
46   // Get field from galice.root
47   AliMagF *fiel = (AliMagF*)gAlice->Field();
48   Double_t fieval=TMath::Abs((Double_t)fiel->SolenoidField()/10.);
49   // Set the conversion constant between curvature and Pt
50   AliTracker::SetFieldMap(fiel,kTRUE);
51   
52   /**** The TPC corner ********************/
53   
54   AliTPCtrackerParam tpcTrackerPar(collcode,fieval);
55   tpcTrackerPar.Init();
56   
57   //**** Switch on the PID class (mandatory!)
58   AliPID pid;
59   
60   /**** The ITS corner ********************/
61   
62   AliITSLoader* itsl = (AliITSLoader*)rl->GetLoader("ITSLoader");
63   if (itsl == 0x0) {
64     cerr<<"Cannot get the ITS loader"<<endl;
65     return;
66   }
67   itsl->LoadRecPoints("read");
68   
69   AliITS *dITS = (AliITS*)gAlice->GetDetector("ITS");
70   if (!dITS) {
71     cerr<<"Cannot find the ITS detector !"<<endl;
72     return;
73   }
74   AliITSgeom *geom = dITS->GetITSgeom();
75   
76   // Instance of the ITS tracker
77   AliITStrackerSA itsTracker(geom);
78   Int_t ITSclusters[6] = {1,1,1,1,1,1};
79   itsTracker.SetLayersNotToSkip(ITSclusters);
80   
81   
82   // Primary vertex reconstruction in pp
83   AliESDVertex *initVertex = 0;
84   TFile *invtx = new TFile("AliESDVertexMean.root");
85   if(collcode==1 && useMeanVtx) {
86     // open the mean vertex
87     initVertex = (AliESDVertex*)invtx->Get("vtxmean");
88   } else {
89     Double_t pos[3]={0.5,0.5,0.}; 
90     Double_t err[3]={3.,3.,5.};
91     initVertex = new AliESDVertex(pos,err);
92   }
93   invtx->Close();
94   delete invtx;
95   AliVertexerTracks *vertexer = new AliVertexerTracks();
96   vertexer->SetVtxStart(initVertex);
97   vertexer->SetDebug(0);
98   delete initVertex;
99   initVertex=0;
100   
101   /***** The TREE for ESD is born *****/
102   TTree *esdTree=new TTree("esdTree","Tree with ESD objects");
103   AliESD *event=0; AliESD *eventTPCin=0;
104   esdTree->Branch("ESD","AliESD",&event);
105   
106   if(firstEvent>rl->GetNumberOfEvents()) firstEvent=rl->GetNumberOfEvents()-1;
107   if(lastEvent>rl->GetNumberOfEvents())  lastEvent=rl->GetNumberOfEvents()-1;
108   cout<<" Number of events: "<<1+lastEvent-firstEvent<<endl;
109   
110   TFile *ppZ = TFile::Open("ITS.Vertex.root"); // z vertices from SPD
111   AliESDVertex *vertexSPD = new AliESDVertex();
112   Char_t zver[100];
113   Double_t vtx[3]={0,0,0};
114   Double_t sigmavtx[3]={0.07,0.07,0.1};
115
116
117   //<---------------------------------- The Loop over events begins
118   TStopwatch timer;
119   Int_t trc;
120   for(Int_t i=firstEvent; i<=lastEvent; i++) { 
121     
122     cout<<" Processing event number : "<<i<<endl;
123     AliESD *event = new AliESD(); 
124     event->SetRunNumber(gAlice->GetRunNumber());
125     event->SetEventNumber(i);
126     event->SetMagneticField(gAlice->Field()->SolenoidField());
127     rl->GetEvent(i);
128
129     //***** Primary vertex from SPD from file 
130     sprintf(zver,"Event%d/Vertex",i);
131     vertexSPD = (AliESDVertex*)ppZ->Get(zver);
132     if(!vertexSPD) {
133       esdTree->Fill(); delete event;
134       continue;
135     }      
136     event->SetVertex(vertexSPD);
137     vertexSPD->GetXYZ(vtx);
138     vertexSPD->GetSigmaXYZ(sigmavtx);
139
140     //***** TPC tracking
141     if ( (trc=tpcTrackerPar.BuildTPCtracks(event)) ) {
142       printf("exiting TPC tracker with code %d in event %d\n",trc,i);
143       esdTree->Fill(); delete event;
144       continue;
145     }
146
147     // make a copy of the ESD at this stage
148     eventTPCin = event;
149
150     //***** ITS tracking
151     itsTracker.AliTracker::SetVertex(vtx,sigmavtx);
152     TTree *itsTree=itsl->TreeR();
153     if (!itsTree) {
154       cerr<<"Can't get the ITS cluster tree !\n";
155       esdTree->Fill(); delete event;
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(); delete event;
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     delete event;
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(AliESD *event,AliESD *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 //--------------------------------------------------------------------------