Vertexing with tracks included (Andrea)
[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 AliITStrackerV2
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 - LNL
14   //
15
16    Int_t  collcode = 1; // pp collisions
17    Bool_t useMeanVtx = kTRUE;
18
19    if (gAlice) {
20       delete gAlice->GetRunLoader();
21       delete gAlice; 
22       gAlice=0;
23    }
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    TDatabasePDG *DataBase = TDatabasePDG::Instance();
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    //An instance of the ITS tracker
78    AliITStrackerV2 itsTracker(geom);
79    Int_t ITSclusters[6] = {1,1,1,1,1,1};
80    itsTracker.SetLayersNotToSkip(ITSclusters);
81
82    // Primary vertex reconstruction in pp
83    AliESDVertex *initVertex = 0;
84    if(collcode==1 && useMeanVtx) {
85      // open the mean vertex
86      TFile *invtx = new TFile("AliESDVertexMean.root");
87      initVertex = (AliESDVertex*)invtx->Get("vtxmean");
88      invtx->Close();
89      delete invtx;
90    } else {
91      Double_t pos[3]={0.5,0.5,0.}; 
92      Double_t err[3]={3.,3.,5.};
93      initVertex = new AliESDVertex(pos,err);
94    }
95    AliVertexerTracks *vertexer = new AliVertexerTracks;
96    vertexer->SetVtxStart(initVertex);
97
98
99    /***** The TREE for ESD is born *****/
100    TTree *esdTree=new TTree("esdTree","Tree with ESD objects");
101    AliESD *event=0;
102    esdTree->Branch("ESD","AliESD",&event);
103    
104    if(firstEvent>rl->GetNumberOfEvents()) firstEvent=rl->GetNumberOfEvents()-1;
105    if(lastEvent>rl->GetNumberOfEvents())  lastEvent=rl->GetNumberOfEvents()-1;
106    cout<<" Number of events: "<<1+lastEvent-firstEvent<<endl;
107    
108    TFile *ppZ = TFile::Open("ITS.Vertex.root"); // z vertices from SPD
109    AliESDVertex *vertexSPD = new AliESDVertex();
110    Char_t zver[100];
111    Double_t vtx[3];
112    Double_t cvtx[6];
113
114
115    //<---------------------------------- The Loop over events begins
116    TStopwatch timer;
117    Int_t trc;
118    for(Int_t i=firstEvent; i<=lastEvent; i++) { 
119      
120      cerr<<" Processing event number : "<<i<<endl;
121      AliESD *event = new AliESD(); 
122      event->SetRunNumber(gAlice->GetRunNumber());
123      event->SetEventNumber(i);
124      event->SetMagneticField(gAlice->Field()->SolenoidField());
125      rl->GetEvent(i);
126
127      //***** Primary vertex from SPD from file 
128      sprintf(zver,"Event%d/Vertex",i);
129      vertexSPD = (AliESDVertex*)ppZ->Get(zver);
130      if(!vertexSPD) {
131        esdTree->Fill(); delete event;
132        continue;
133      }      
134      event->SetVertex(vertexSPD);
135      vertexSPD->GetXYZ(vtx);
136      vertexSPD->GetCovMatrix(cvtx);
137
138      //***** TPC tracking
139      if ( (trc=tpcTrackerPar.BuildTPCtracks(event)) ) {
140        printf("exiting TPC tracker with code %d in event %d\n",trc,i);
141        esdTree->Fill(); delete event;
142        continue;
143      }
144
145      //***** ITS tracking
146      itsTracker.SetVertex(vtx,cvtx);
147      TTree *itsTree=itsl->TreeR();
148      if (!itsTree) {
149         cerr<<"Can't get the ITS cluster tree !\n";
150         esdTree->Fill(); delete event;
151         return;
152      }     
153      itsTracker.UnloadClusters();
154      itsTracker.LoadClusters(itsTree);
155      if ( (trc=itsTracker.Clusters2Tracks(event)) ) {
156        printf("exiting ITS tracker with code %d in event %d\n",trc,i);
157        esdTree->Fill(); delete event;
158        continue;
159      }
160
161      //***** Propagate kITSin tracks to local x = 0 (through beam pipe)
162      ToLocalX0(event);
163
164
165      //***** Vertex from ESD tracks
166      if(collcode==1) { // pp
167        AliESDVertex *vertexTrks = 
168          (AliESDVertex*)vertexer->FindPrimaryVertex(event);
169        event->SetPrimaryVertex(vertexTrks); 
170      }
171
172      esdTree->Fill();
173      delete event;
174
175    }//<-----------------------------------The Loop over events ends here
176    timer.Stop(); timer.Print();
177
178    // The AliESDs.root is born
179    TFile *ef = TFile::Open("AliESDs.root","RECREATE"); 
180    if (!ef || !ef->IsOpen()) {cerr<<"Can't open AliESDs.root !\n"; return;}
181
182    //Write the TREE and close everything
183    esdTree->Write();
184    delete esdTree;
185    ef->Close();
186
187    delete vertexer;
188    delete initVertex;
189    delete rl;
190
191    return;
192 }
193 //--------------------------------------------------------------------------
194 Int_t ToLocalX0(AliESD *esd) {
195
196   Int_t ntracks = esd->GetNumberOfTracks();
197   AliESDtrack *esdTrack = 0;
198
199   // loop on tracks
200   for(Int_t tr = 0; tr<ntracks; tr++) {
201     esdTrack = (AliESDtrack*)esd->GetTrack(tr);
202     // ask for kITSin
203     if(!(esdTrack->GetStatus()&AliESDtrack::kITSin)) continue;
204     AliITStrackV2 *itsTrack = 0;
205     try {
206       itsTrack = new AliITStrackV2(*esdTrack);
207       esdTrack = 0;
208     }
209     catch (const Char_t *msg) {
210         Warning("FindPrimaryVertexForCurrentEvent",msg);
211         continue;
212     }
213     // propagate track to beam pipe (0.8 mm of Be) 
214     if(itsTrack->GetX()>3.) itsTrack->PropagateTo(3.,0.0023,65.19);
215     // propagate track to (0,0)
216     itsTrack->PropagateTo(0.,0.,0.);
217     itsTrack->UpdateESDtrack(AliESDtrack::kITSin);
218     esdTrack = new AliESDtrack(*(itsTrack->GetESDtrack()));
219     delete itsTrack;
220   } // end loop on tracks
221
222   return 0;
223 }