]> git.uio.no Git - u/mrichter/AliRoot.git/blob - TPC/AliBarrelRec_TPCparam.C
AliTPCCalibSignal renamed to AliTPCCalibPulser
[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   TGeoManager::Import("geometry.root");
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");
104   AliESD *event=0; AliESD *eventTPCin=0;
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     
123     cout<<" Processing event number : "<<i<<endl;
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     }
147
148     // make a copy of the ESD at this stage
149     eventTPCin = event;
150
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;
157       return;
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     }
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);
174     
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
193   //Write the tree and close everything
194   esdTree->Write();
195   delete esdTree;
196   ef->Close();
197
198   delete vertexer;
199   delete rl;
200
201   return;
202 }
203 //--------------------------------------------------------------------------
204 void BackToTPCInnerWall(AliESD *event,AliESD *eventTPC) {
205
206   Int_t ntracks = eventTPC->GetNumberOfTracks();
207   AliESDtrack *esdTrackTPC = 0;
208
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();
217   AliESDtrack *esdTrack = 0;
218   esdTrackTPC = 0;
219   Int_t indexTPC;
220
221   // loop on tracks
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
232     if(!(esdTrack->GetStatus()&AliESDtrack::kITSin)) continue;
233
234     indexTPC = labelsTPC[TMath::Abs(esdTrack->GetLabel())];
235     esdTrackTPC = (AliESDtrack*)eventTPC->GetTrack(indexTPC);
236
237     AliITStrackMI *itsTrack = 0;
238     try {
239       itsTrack = new AliITStrackMI(*esdTrackTPC);
240       esdTrack = 0;
241     }
242     catch (const Char_t *msg) {
243         Warning("ToTPCInnerWall",msg);
244         continue;
245     }
246     itsTrack->UpdateESDtrack(AliESDtrack::kITSout);
247     esdTrack = new AliESDtrack(*(itsTrack->GetESDtrack()));
248
249     delete itsTrack;
250   } // end loop on tracks
251
252   return;
253 }
254 //--------------------------------------------------------------------------