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