1 ////////////////////////////////////////////////////////////////////////
3 // AliITSandTPCrecoV2.C
6 ////////////////////////////////////////////////////////////////////////
9 /****************************************************************************
10 * This macro is supposed to do reconstruction in the barrel ALICE trackers *
11 * (Make sure you have TPC digits and ITS hits before using this macro !!!) *
12 * It does the following steps *
13 * 1) TPC cluster finding *
14 * 2) TPC track finding *
15 * 3) ITS cluster finding - fast recpoints - *
17 * TPC part taken from AliTPCTracking.C - Jiri Chudoba *
18 * The TPC part will be removed *
19 * M.Masera 30/09/2002 09:40 *
20 * Last revision: 4/3/2003 *
23 * The main function is AliITSandTPCrecoV2 All the input arguments have a *
25 * nEvents: number of events to be processed. *
26 * If nEvents<0, all the available events are used starting from *
28 * firstEvent: is the first event to be analyzed. *
29 * It is set to 0 if nEvents<0 (default) *
30 * fileNameHits: input root file with Kine and hits *
31 * fileNameDigits: input root file with ITS and TPC digits and with ITS *
32 * recpoints. Only fast points are used at the moment. *
33 * An option to use slow points should be added *
34 * fileNameClusters: output root file with TPC clusters *
35 * fileNameTracks: output root file with TPC tracks *
36 * fileNameITSClusters: output root file with ITS clusters *
37 * fileNameITSTracks: output root file with ITS tracks *
39 * AliITSandTPCrecoV2(Int_t nEvents=-1, Int_t firstEvent=0, *
40 * TString fileNameHits="galice.root", *
41 * TString fileNameDigits="galiceSDR.root", *
42 * TString fileNameClusters="tpc.clusters.root", *
43 * TString fileNameTracks="tpc.tracks.root", *
44 * TString fileNameITSClusters="its.clusters.root", *
45 * TString fileNameITSTracks="its.tracks.root"); *
49 * Int_t gDEBUG = 0; // NO verbose printouts if ==0 *
50 * TArrayD *gzvtx = 0; // Z coordinate of primary vertices *
51 * Int_t gCurrEve = 0; // current event number *
52 * Bool_t gPPMode = kTRUE; // it must be set to kTRUE for pp interactions *
53 * Bool_t gFastReco = kTRUE; // it must set to kTRUE for fast reconstr. *
54 * const Double_t kgVertexNotFound = -100; // the Z coord. of the *
55 * primary vertex is set to this value *
56 * when the vertexing fails *
57 * Bool_t gUseNewClustersV2 = kFALSE; // if kTRUE the AliITSclustereV2 class
58 * is used to produce V2 clusters. *
59 * ITS recpoints are used otherwise *
60 ****************************************************************************/
62 #if !defined(__CINT__) || defined(__MAKECINT__)
67 #include <TParticle.h>
71 #include "AliMagFCM.h"
72 #include "AliTPCtracker.h"
73 #include "AliTPCclusterer.h"
74 #include "Riostream.h"
80 #include "TParticle.h"
84 #include "TObjArray.h"
88 #include "AliITSclustererV2.h"
89 #include "AliITSgeom.h"
90 #include "AliITSRecPoint.h"
91 #include "AliITSclusterV2.h"
92 #include "AliITStrackerV2.h"
93 #include "AliITSVertexerPPZ.h"
94 #include "AliITSVertexerIons.h"
96 #include "AliHeader.h"
97 #include "AliGenEventHeader.h"
99 #include "AliGenPythiaEventHeader.h"
107 Bool_t gPPMode = kTRUE;
108 Bool_t gFastReco = kTRUE;
109 Bool_t gUseNewClustersV2 = kFALSE;
110 const Double_t kgVertexNotFound = -100.;
113 Int_t TPCFindClusters(Int_t nEvents=1, Int_t firstEvent=0,
114 TString fileNameDigits="galiceSDR.root",
115 TString fileNameClusters="tpc.clusters.root");
116 Int_t TPCFindClusters(Int_t nEvents, Int_t firstEvent,
117 TFile* fileDigits, TFile* fileClusters,
118 AliTPCParam* paramTPC=0);
119 Int_t TPCFindTracks(Int_t nEvents=1, Int_t firstEvent=0,
120 TString fileNameClusters="tpc.clusters.root",
121 TString fileNameTracks="tpc.tracks.root");
122 Int_t TPCFindTracks(Int_t nEvents, Int_t firstEvent,
123 TFile* fileClusters, TFile* fileTracks,
124 AliTPCParam* paramTPC=0);
126 Int_t AliTPCTracking(Int_t nEvents=1, Int_t firstEvent=0,
127 TString fileNameHits="galice.root",
128 TString fileNameDigits="galiceSDR.root",
129 TString fileNameClusters="tpc.clusters.root",
130 TString fileNameTracks="tpc.tracks.root");
131 Int_t AliITSandTPCrecoV2(Int_t nEvents=-1, Int_t firstEvent=0,
132 TString fileNameHits="galice.root",
133 TString fileNameDigits="galiceSDR.root",
134 TString fileNameClusters="tpc.clusters.root",
135 TString fileNameTracks="tpc.tracks.root",
136 TString fileNameITSClusters="its.clusters.root",
137 TString fileNameITSTracks="its.tracks.root");
138 Int_t ITSFindClusters(TString inname, TString recpointFileName, TString fileNameITSClusters, Int_t n, Int_t first);
139 Int_t ITSFindTracks(TString fileNameTracks, TString fileNameITSClusters, TString fileNameITSTracks, Int_t n, Int_t first, Int_t vc, Int_t vc2);
140 void FindZV(Int_t nEvents, Int_t first, TString FileNameHits, TString FileWithRecPoints);
141 Int_t DirectClusterFinder(Int_t eventn,Int_t first,TFile *in,TFile *out);
142 ////////////////////////////////////////////////////////////////////////
143 Int_t AliITSandTPCrecoV2(Int_t nEvents, Int_t firstEvent,
144 TString fileNameHits,
145 TString fileNameDigits,
146 TString fileNameClusters,
147 TString fileNameTracks,
148 TString fileNameITSClusters,
149 TString fileNameITSTracks) {
150 const Char_t *name="AliITSandTPCrecoV2";
152 gBenchmark->Start(name);
154 // Set the conversion constant for the Kalman filter
155 // and set the gAlice pointer
156 Char_t *finame = (Char_t *)fileNameHits.Data();
157 AliTracker::SetFieldFactor(finame,kFALSE);
160 nEvents = (Int_t)gAlice->TreeE()->GetEntries();
162 cout << "Processing events from " << firstEvent << " up to " << nEvents-1 <<endl;
165 if(fileNameDigits!=fileNameHits)gAlice->SetTreeDFileName(fileNameDigits.Data());
166 gzvtx = new TArrayD(nEvents);
168 // measure Z vertex coordinate for seeding
170 FindZV(nEvents,firstEvent,fileNameHits,fileNameDigits);
172 // TPC tracking - abort macro execution if tracking fails
173 if(AliTPCTracking(nEvents,firstEvent,fileNameHits,fileNameDigits,
174 fileNameClusters,fileNameTracks)) {
175 cerr << "TPC tracking failed \n";
180 if(ITSFindClusters(fileNameHits,fileNameDigits,fileNameITSClusters,nEvents,firstEvent)) {
181 cerr << "ITS clustering failed \n";
185 Int_t firstpass = 1; // 1=vert. constraint; 0=no vert. constr.; -1=skip pass
186 Int_t secondpass = 0; // as above
187 if(ITSFindTracks(fileNameTracks,fileNameITSClusters,fileNameITSTracks,nEvents,firstEvent,firstpass,secondpass)) {
188 cerr << "ITS tracking failed \n";
193 gBenchmark->Stop(name);
194 gBenchmark->Show(name);
197 ////////////////////////////////////////////////////////////////////////
198 Int_t AliTPCTracking( Int_t nEvents, Int_t firstEvent,
199 TString fileNameHits,
200 TString fileNameDigits,
201 TString fileNameClusters,
202 TString fileNameTracks) {
205 // ********** Find TPC clusters *********** //
206 if (TPCFindClusters(nEvents,firstEvent,fileNameDigits,fileNameClusters)) {
207 cerr<<"Failed to get TPC clusters: !\n";
211 // ********** Find TPC tracks *********** //
212 // if (TPCFindTracks(TPCclsName,TPCtrkName,n)) {
213 if (TPCFindTracks(nEvents,firstEvent,fileNameClusters,fileNameTracks)) {
214 cerr<<"Failed to get TPC tracks !\n";
221 ////////////////////////////////////////////////////////////////////////
222 Int_t TPCFindClusters(Int_t nEvents, Int_t firstEvent,
223 TString fileNameDigits, TString fileNameClusters) {
226 const Char_t *name="TPCFindClusters";
227 if (gDEBUG>1) cout<<name<<" starts...\n";
228 if (gDEBUG>1) gBenchmark->Start(name);
229 TFile *fileClusters = TFile::Open(fileNameClusters,"recreate");
230 TFile *fileDigits = TFile::Open(fileNameDigits);
231 if (!fileDigits->IsOpen()) {
232 cerr<<"Cannnot open "<<fileNameDigits<<" !\n";
235 if (!fileClusters->IsOpen()) {
236 cerr<<"Cannnot open "<<fileNameClusters<<" !\n";
240 rc = TPCFindClusters(nEvents,firstEvent,fileDigits,fileClusters);
243 fileClusters->Close();
246 if (gDEBUG>1) gBenchmark->Stop(name);
247 if (gDEBUG>1) gBenchmark->Show(name);
251 ////////////////////////////////////////////////////////////////////////
252 Int_t TPCFindClusters(Int_t nEvents, Int_t firstEvent,
253 TFile* fileDigits, TFile* fileClusters,
254 AliTPCParam* paramTPC) {
257 if (!paramTPC) paramTPC = AliTPC::LoadTPCParam(fileDigits);
258 if (!paramTPC) return 1;
260 for (Int_t iEvent = firstEvent; iEvent < firstEvent+nEvents; iEvent++){
261 if (gDEBUG > 2) cout<<"TPCFindClusters: event "<<iEvent<<endl;
262 AliTPCclusterer::Digits2Clusters(paramTPC, fileClusters, iEvent);
266 ////////////////////////////////////////////////////////////////////////
267 Int_t TPCFindTracks(Int_t nEvents, Int_t firstEvent,
268 TString fileNameClusters, TString fileNameTracks) {
271 const Char_t *name="TPCFindTracks";
272 if (gDEBUG>1) cout<<name<<" starts"<<endl;
273 if (gDEBUG>1) gBenchmark->Start(name);
274 TFile *fileTracks = TFile::Open(fileNameTracks,"recreate");
275 TFile *fileClusters =TFile::Open(fileNameClusters);
277 rc = TPCFindTracks(nEvents, firstEvent, fileClusters, fileTracks);
279 fileClusters->Close();
283 if (gDEBUG>1) gBenchmark->Stop(name);
284 if (gDEBUG>1) gBenchmark->Show(name);
288 ////////////////////////////////////////////////////////////////////////
289 Int_t TPCFindTracks(Int_t nEvents, Int_t firstEvent,
290 TFile *fileClusters, TFile * fileTracks,
291 AliTPCParam* paramTPC) {
294 if (!paramTPC) paramTPC = AliTPC::LoadTPCParam(fileClusters);
295 if (!paramTPC) return 1;
297 for(Int_t j=0; j<3; j++){vertex[j]=0.;}
298 AliTPCtracker *tracker = new AliTPCtracker(paramTPC);
299 for (Int_t iEvent = firstEvent; iEvent < firstEvent+nEvents; iEvent++){
300 if (gDEBUG > 2) cout<<"TPCFindTracks: event "<<iEvent<<endl;
301 tracker->SetEventNumber(iEvent);
302 if((*gzvtx)[iEvent] > -100){
303 vertex[2] = (*gzvtx)[iEvent];
308 tracker->SetVertex(vertex);
310 rc = tracker->Clusters2Tracks(0,fileTracks);
317 ////////////////////////////////////////////////////////////////////////
318 Int_t ITSFindClusters(TString inname, TString recpointFileName, TString fileNameITSClusters, Int_t n, Int_t first) {
320 const Char_t *name="ITSFindClusters";
321 cerr<<'\n'<<name<<"...\n";
322 gBenchmark->Start(name);
323 TFile *out=TFile::Open(fileNameITSClusters,"recreate");
324 TFile *in = (TFile*)gROOT->GetListOfFiles()->FindObject(inname.Data());
325 if(gUseNewClustersV2){
326 cout<<"Direct method\n";
327 return DirectClusterFinder(n,first,in,out);
330 cout<<"Clusters through rec points \n";
333 AliITS *ITS = (AliITS*)gAlice->GetModule("ITS");
334 if (!ITS) { cerr<<"Can't get the ITS !\n"; return 1;}
335 AliITSgeom *geom=ITS->GetITSgeom();
336 AliITSclustererV2 *clusterer = new AliITSclustererV2(geom);
341 if(gFastReco && gDEBUG>0)cout <<"Using fast points\n";
342 else cout<<"Using slow points\n";
343 for (ev = first; ev<n; ev++){
344 cout<<"ITSFindClusters: processing event "<<ev<<endl;
346 gAlice->GetEvent(ev);
348 TTree *pTree=gAlice->TreeR();
350 cerr << "ITSFindClusters: can't get TreeR\n";
355 branch = pTree->GetBranch("ITSRecPointsF");
358 branch = pTree->GetBranch("ITSRecPoints");
361 cerr << "ITSFindClusters: can't get ITSRecPoints branch \n";
366 TClonesArray *clusters=new TClonesArray("AliITSclusterV2",10000);
368 sprintf(cname,"TreeC_ITS_%d",ev);
370 TTree *cTree=new TTree(cname,"ITS clusters");
371 cTree->Branch("Clusters",&clusters);
373 TClonesArray *points=new TClonesArray("AliITSRecPoint",10000);
374 branch->SetAddress(&points);
376 TClonesArray &cl=*clusters;
378 Int_t nentr=(Int_t)pTree->GetEntries();
379 AliITSgeom *geom=ITS->GetITSgeom();
383 for (Int_t i=0; i<nentr; i++) {
385 clusterer->RecPoints2Clusters(points,i,clusters);
386 cTree->Fill(); clusters->Delete();
390 delete cTree; delete clusters; points->Delete(); delete points;
398 ////////////////////////////////////////////////////////////////////////
399 Int_t ITSFindTracks(TString inname, TString inname2, TString outname, Int_t n, Int_t first, Int_t vc, Int_t vc2) {
402 AliITStrackV2::SetSigmaYV(0.015);
403 AliITStrackV2::SetSigmaZV(0.017);
405 TFile *out=TFile::Open(outname.Data(),"recreate");
406 if (!out->IsOpen()) {
407 cerr<<"Can't open file "<<outname.Data()<<endl;
410 TFile *in =TFile::Open(inname.Data());
412 cerr<<"Can't open file "<<inname.Data()<<endl;
415 TFile *in2 =TFile::Open(inname2.Data());
416 if (!in2->IsOpen()) {
417 cerr<<"Can't open file "<<inname2.Data()<<endl;
421 AliITSgeom *geom=(AliITSgeom*)gFile->Get("AliITSgeom");
422 if (!geom) { cerr<<"can't get ITS geometry !\n"; return 1;}
424 AliITStrackerV2 tracker(geom);
425 for (Int_t i=first;i<n;i++){
426 tracker.SetEventNumber(i);
427 Double_t vtx=(*gzvtx)[i];
428 cout<<"++++++++++++++++++++++++++++++++++++++++++++++++++++\n";
429 cout<<"ITSFindTracks -- event "<<i<<endl;
430 cout<<"Computed vertex position: "<<vtx<<endl;
435 if(vtx != kgVertexNotFound){
436 tracker.SetVertex(vrtc); // set vertex only if it was computed
439 vc = 0; // use vertex contraint only if vertex info is available
443 tracker.SetupFirstPass(&vc);
445 tracker.SetupSecondPass(&vc2);
447 rc=tracker.Clusters2Tracks(in,out);
457 //////////////////////////////////////////////////////////////////
458 void FindZV(Int_t nEvents, Int_t first, TString FileNameHits, TString FileWithRecPoints){
459 TFile *in = (TFile*)gROOT->GetListOfFiles()->FindObject(FileNameHits.Data());
460 TFile *fo = new TFile("AliITSVertices.root","recreate");
461 if(FileNameHits!=FileWithRecPoints)gAlice->SetTreeRFileName(FileWithRecPoints);
462 AliITSVertexer *findVert;
464 findVert = new AliITSVertexerPPZ(in,fo);
467 findVert = new AliITSVertexerIons(in,fo);
469 Int_t last = first + nEvents - 1;
470 AliITSVertex *vert = 0;
471 for(Int_t i=first; i<=last; i++){
473 // The true Z coord. is fetched for comparison
474 AliHeader *header = gAlice->GetHeader();
475 AliGenEventHeader* genEventHeader = header->GenEventHeader();
476 TArrayF primaryVertex(3);
477 genEventHeader->PrimaryVertex(primaryVertex);
478 vert = findVert->FindVertexForCurrentEvent(i);
481 for(Int_t kk=0;kk<3;kk++)pos[kk]=(Double_t)primaryVertex[kk];
482 vert->SetTruePos(pos);
483 Double_t found = vert->GetZv();
484 gzvtx->AddAt(found,i);
485 findVert->WriteCurrentVertex();
488 gzvtx->AddAt(kgVertexNotFound,i);
498 Int_t DirectClusterFinder(Int_t eventn,Int_t first,TFile *in,TFile* out){
501 AliITS *ITS = (AliITS*)gAlice->GetModule("ITS");
502 if (!ITS) { cerr<<"Can't find the ITS !\n"; return 3; }
504 AliITSgeom *geom=ITS->GetITSgeom();
509 AliITSclustererV2 clusterer(geom);
510 for (Int_t i=first; i<eventn; i++) {
511 cerr<<"Processing event number: "<<i<<endl;
513 //ITS->MakeTreeC(); //To make the V1 cluster finders happy
514 clusterer.SetEvent(i);
515 if (gFastReco) clusterer.Digits2Clusters(in,out);
516 else clusterer.Hits2Clusters(in,out);
518 timer.Stop(); timer.Print();