358e2fc3cc06d80a60fd8fde2c7fddc4cc8d113c
[u/mrichter/AliRoot.git] / TPC / AliBarrelReconstruction.C
1 /****************************************************************************
2  * This macro is supposed to do reconstruction in the barrel ALICE trackers *
3  * (Make sure you have TPC digits and ITS hits before using this macro !!!) *
4  * It does the following steps (April 12, 2001):                            *
5  *                   1) TPC cluster finding                                 *
6  *                   2) TPC track finding                                   *
7  *                   3) ITS cluster finding V2 (via fast points !)          *
8  *                   4) ITS track finding V2                                *
9  *                   5) ITS back track propagation V2                       *
10  *                   6) TPC back track propagation                          *
11  *                (Origin: I.Belikov, CERN, Jouri.Belikov@cern.ch)          *
12  ****************************************************************************/
13
14 #ifndef __CINT__
15   #include "alles.h"
16   #include "iostream.h"
17   #include "AliRun.h"
18   #include "AliRunLoader.h"
19   #include "AliLoader.h"
20   #include "AliTPCLoader.h"
21   #include "AliITSLoader.h"
22   #include "AliMagF.h"
23   #include "AliTPCtracker.h"
24   #include "AliITS.h"
25   #include "AliITSgeom.h"
26   #include "AliITSRecPoint.h"
27   #include "AliITSclusterV2.h"
28   #include "AliITSsimulationFastPoints.h"
29   #include "AliITStrackerV2.h"
30 #endif
31
32
33
34 Int_t TPCFindClusters( Int_t n);
35 Int_t TPCFindTracks(Int_t n);
36 Int_t TPCSortTracks(const Char_t *outname,  Int_t n);
37 Int_t TPCPropagateBack();
38
39 Int_t ITSFindClusters(Int_t n);
40 Int_t ITSFindTracks(const Char_t *inname2, Int_t n);
41 Int_t ITSPropagateBack();
42
43 const char* TPCtrkNameS= "TPC.TracksSorted.root";
44
45 class AliRunLoader;
46 class AliTPCLoader;
47 class AliTPCParam;
48 AliRunLoader *rl = 0x0;
49 AliTPCLoader *tpcl = 0x0;   
50 AliTPCParam  *param = 0x0;
51 Bool_t debug = kFALSE;
52 Int_t AliBarrelReconstruction(Int_t n=1) 
53  {
54  
55    if (gAlice)
56     {
57       delete gAlice->GetRunLoader();
58       delete gAlice;//if everything was OK here it is already NULL
59       gAlice = 0x0;
60     }
61     
62    rl = AliRunLoader::Open("galice.root");
63    if (rl == 0x0)
64     {
65       cerr<<"Can not open session"<<endl;
66       return 1;
67     }
68    
69    if (rl->LoadgAlice())
70     {
71       cerr<<"Error occured while l"<<endl;
72       return 1;
73     }
74    AliKalmanTrack::SetConvConst(1000/0.299792458/rl->GetAliRun()->Field()->SolenoidField());
75    
76    tpcl = (AliTPCLoader*)rl->GetLoader("TPCLoader");
77    if (tpcl == 0x0)
78     {
79       cerr<<"Can not get TPC Loader"<<endl;
80       return 1;
81     }
82
83 // ********** Find TPC clusters *********** //
84    if ( TPCFindClusters(n) ) 
85     {
86       cerr<<"Failed to get TPC clusters !\n";
87       return 1;
88     }
89    
90 // ********** Find TPC tracks *********** //
91     if (TPCFindTracks(n)) 
92      {
93       cerr<<"Failed to get TPC tracks !\n";
94       return 1;
95      }
96
97 //   cout<<"Stopping tracking on TPC\n";
98 // ********** Sort and label TPC tracks *********** //
99    if (TPCSortTracks(TPCtrkNameS,n)) {
100       cerr<<"Failed to sort TPC tracks !\n";
101       return 1;
102     } 
103
104 //   cout<<"Stopping tracking on TPC sorted T\n";
105
106    
107 // ********** Find ITS clusters *********** //
108    if (ITSFindClusters(n)) 
109     {
110       cerr<<"Failed to get ITS clusters !\n";
111       return 1;
112     }
113    
114 // ********** Find ITS tracks *********** //
115
116    if (ITSFindTracks(TPCtrkNameS,n)) 
117     {
118       cerr<<"Failed to get ITS tracks !\n";
119       return 1;
120     }
121    //clsFile->Close();
122    cout<<"Stopping on ITSFindTracks\n";
123    delete rl;
124    rl = 0x0;
125    return 0;
126    return 1;
127
128 // ********** Back propagation of the ITS tracks *********** //
129    if ( ITSPropagateBack() ) {
130       cerr<<"Failed to propagate back the ITS tracks !\n";
131       return 1;
132    }
133
134
135 // ********** Back propagation of the TPC tracks *********** //
136    
137    if (TPCPropagateBack()) {
138       cerr<<"Failed to propagate back the TPC tracks !\n";
139       return 1;
140    }
141
142    return 0;
143 }
144
145
146 Int_t TPCFindClusters(Int_t n) 
147  {
148    Int_t rc=0;
149    const Char_t *name="TPCFindClusters";
150    cerr<<'\n'<<name<<"...\n";
151    rl->CdGAFile();
152    param=(AliTPCParam *)gDirectory->Get("75x40_100x60");
153    if (!param) 
154     {
155      param=(AliTPCParam *)gDirectory->Get("75x40_100x60_150x60");
156      if (!param) 
157       {
158         cerr<<"TPC parameters have not been found !\n";
159         return 1;
160       }
161     }
162 //   param->Dump();
163    param->Update();
164
165    gBenchmark->Start(name);
166
167    AliTPCv2 tpc;
168    tpc.SetParam(param);
169    tpc.SetLoader(tpcl);
170  
171    tpcl->LoadDigits("read");
172    tpcl->LoadRecPoints("recreate");
173   
174
175    //tpc.Digits2Clusters(out); //MI change
176    for (Int_t i=0;i<n;i++)
177     {
178       printf("Processing event %d\n",i);
179       tpc.Digits2Clusters(i);
180      //     AliTPCclusterer::Digits2Clusters(dig, out, i);
181     }
182    
183    tpcl->UnloadDigits();
184    tpcl->UnloadRecPoints();
185    gBenchmark->Stop(name);
186    gBenchmark->Show(name);
187
188    return rc;
189 }
190
191 Int_t TPCFindTracks(Int_t n) {
192
193    Int_t rc=0;
194    AliTPCtracker *tracker = 0x0;
195    const Char_t *name="TPCFindTracks";
196    rl->GetEvent(0);
197    rl->CdGAFile();
198    param=(AliTPCParam *)gDirectory->Get("75x40_100x60");
199    
200    if (!param) 
201     {
202      param=(AliTPCParam *)gDirectory->Get("75x40_100x60_150x60");
203      if (!param) 
204       {
205         cerr<<"TPC parameters have not been found !\n";
206         return 1;
207       }
208     }
209 //   param->Dump();
210    param->Update();
211
212    gBenchmark->Start(name);
213
214    tpcl->LoadRecPoints("read");
215    tpcl->LoadTracks("recreate");
216
217    //AliTPCtracker *tracker=new AliTPCtracker(param);
218    //rc=tracker->Clusters2Tracks(0,out);
219    //delete tracker;
220    
221    for (Int_t i=0;i<n;i++){
222      rl->GetEvent(i);
223      printf("Processing event %d\n",i);
224      
225      tracker = new AliTPCtracker(param, i, AliConfig::fgkDefaultEventFolderName.Data());
226      //Int_t rc=
227      tracker->Clusters2Tracks();
228      delete tracker;
229    }
230    tpcl->UnloadRecPoints();
231    tpcl->UnloadTracks();
232
233    gBenchmark->Stop(name);
234    gBenchmark->Show(name);
235
236    return rc;
237 }
238
239
240 Int_t TPCSortTracks(const Char_t * outname,  Int_t eventn){
241    Int_t rc=0;
242    const Char_t *name="TPCSortTracks";
243
244    cerr<<'\n'<<name<<"...\n";
245    rl->CdGAFile();
246    param=(AliTPCParam *)gDirectory->Get("75x40_100x60");
247    if (!param) 
248     {
249      param=(AliTPCParam *)gDirectory->Get("75x40_100x60_150x60");
250      if (!param) 
251       {
252         cerr<<"TPC parameters have not been found !\n";
253         return 1;
254       }
255     }
256    param->Update();
257
258    gBenchmark->Start(name);
259    
260  
261    AliRunLoader* rl2 = AliRunLoader::Open("galice.root","tmp");
262    
263    AliLoader* tpcl2 = (AliTPCLoader*)rl2->GetLoader("TPCLoader");
264    cout<<"tpcl2->SetTracksFileName("<<outname<<");\n";
265    tpcl2->SetTracksFileName(TString(outname));
266    tpcl2->LoadTracks("recreate");
267    
268    // loop over events 
269    for (Int_t event=0;event<eventn; event++)
270     {
271      
272      rl->GetEvent(event);
273      rl2->GetEvent(event);
274
275      TObjArray tarray(10000);
276      AliTPCtrack *iotrack=0;
277      Int_t i;
278      
279      if (tpcl->TreeT() == 0x0)   tpcl->LoadTracks("read");
280      TTree *tracktree=tpcl->TreeT();
281      if (tracktree == 0x0)
282       {
283         cerr<<"Can not get TreeT for event "<<event<<endl;
284         continue;
285       }
286      TBranch *tbranch=tracktree->GetBranch("tracks");
287      Int_t nentr=(Int_t)tracktree->GetEntries();
288      cout<<"IN Tracks nentr = "<<nentr<<endl;
289      for (i=0; i<nentr; i++) 
290       {
291        iotrack=new AliTPCtrack;
292        tbranch->SetAddress(&iotrack);
293        tracktree->GetEvent(i);
294        tarray.AddLast(iotrack);
295       }   
296      tarray.Sort();
297      
298      //assign thacks GEANT labels
299      cout<<"Running Tracker\n";
300      AliTPCtracker *tracker = new AliTPCtracker(param, event);
301
302      cout<<"Load Sectors\n";
303
304      tracker->LoadInnerSectors();
305      tracker->LoadOuterSectors();
306
307      cout<<"Cooking Labels\n";
308
309      for (i=0; i<nentr; i++) 
310       {
311        iotrack=(AliTPCtrack*)tarray.UncheckedAt(i);
312        tracker->CookLabel(iotrack,0.1);
313       }   
314      cout<<"deleting tracker\n";
315      delete tracker;
316      //in->Close();
317      //end of GEANT label assignment
318      
319
320      if (tpcl2->TreeT() == 0x0)   tpcl2->MakeTree("T");
321      tracktree= tpcl2->TreeT();
322      if (tracktree == 0x0)
323       {
324         cerr<<"Can not get TreeT for Sorted tracks\n";
325         return 1;
326       }
327      tracktree->Branch("tracks","AliTPCtrack",&iotrack,32000,0);
328      for (i=0; i<nentr; i++) {
329        iotrack=(AliTPCtrack*)tarray.UncheckedAt(i);
330        tracktree->Fill();
331      }
332
333      cout<<" Tracks File Name is "<<tpcl2->GetTracksFileName()<<endl;
334      tpcl2->WriteTracks("OVERWRITE");
335    }
336   
337
338    delete rl2;
339
340    tpcl->UnloadTracks();
341    tpcl->UnloadRecPoints(); 
342    
343    gBenchmark->Stop(name);
344    gBenchmark->Show(name);
345
346    return rc;
347 }
348
349 Int_t ITSFindClusters(Int_t n) 
350  {
351    Int_t rc=0;
352    const Char_t *name="ITSFindClusters";
353    cerr<<'\n'<<name<<"...\n";
354    gBenchmark->Start(name);
355
356    if (rl->GetEvent(0))
357     {
358       cerr<<"Problems\n";
359       return 1;
360     }
361     
362    AliITSLoader* itsl = (AliITSLoader*)rl->GetLoader("ITSLoader");
363    if (itsl == 0x0)
364     {
365       cerr<<"Can not get ITS Loader from Run Loader\n";
366       return 1;
367     } 
368   
369    gAlice=rl->GetAliRun();
370    if (!gAlice) {
371       cerr<<"Can't get gAlice !\n";
372       return 1;
373    }
374    if (rl->TreeK() == 0x0) rl->LoadKinematics();
375    if (rl->TreeE() == 0x0) rl->LoadHeader();
376    
377    itsl->LoadRawClusters("recreate");
378    
379    AliITS *ITS  = (AliITS*)gAlice->GetModule("ITS");
380    if (!ITS) { cerr<<"Can't get the ITS !\n"; return 1;}
381    AliITSgeom *geom=ITS->GetITSgeom();
382
383    Int_t ev=0;
384    for (ev = 0; ev<n; ev++)
385    {
386
387      rl->GetEvent(ev);
388      TBranch *branch = 0x0;
389      TTree *pTree = 0x0;
390      
391      pTree=itsl->TreeR();
392      if (pTree == 0x0) 
393       {
394         itsl->LoadRecPoints("read");
395         pTree=itsl->TreeR();
396         if (pTree == 0x0) 
397          {
398            ::Error("AliBarrelReonstruction.C::ITSFindClusters",
399                    "Can not get TreeR for event %d",ev);
400            return 1;
401          }
402       }
403      
404      pTree=itsl->TreeR();
405      branch=(pTree)?pTree->GetBranch("ITSRecPoints"):0x0;
406
407      if (branch== 0x0) {
408        //if not reconstructed ITS branch do reconstruction 
409        if (debug) ::Info("AliBarrelReconstruction.C","Did not get ITSRecPoints from TreeR.");
410        if (debug) ::Info("AliBarrelReconstruction.C","Making branch and Producing RecPoints");
411        itsl->SetRecPointsFileOption("recreate");
412        if (itsl->TreeR()==0x0) itsl->MakeTree("R");
413        
414        ITS->MakeBranch("RF");
415        ITS->SetTreeAddress();
416        itsl->LoadHits();
417
418        //////////////// Taken from ITSHitsToFastPoints.C ///////////////////////
419        for (Int_t i=0;i<3;i++) { 
420          ITS->SetSimulationModel(i,new AliITSsimulationFastPoints()); 
421        }
422        Int_t nsignal=25;
423        Int_t size=-1;
424        Int_t bgr_ev=Int_t(ev/nsignal);
425        
426        ITS->HitsToFastRecPoints(ev,bgr_ev,size," ","All"," ");
427        //////////////////////////////////////////////////////////////////////////
428        //MI comment  - in HitsToFast... they reset treeR to 0 
429        //they overwrite full reconstructed event ???? ... so lets connect TreeR one more
430        //time
431        itsl->UnloadHits();
432        ITS->SetTreeAddress();
433      }
434
435      TClonesArray *clusters=new TClonesArray("AliITSclusterV2",10000);
436
437      if (itsl->TreeC() == 0x0) itsl->MakeTree("C");
438      TTree *cTree = itsl->TreeC();
439      TBranch* b = cTree->GetBranch("Clusters");
440      if (b == 0x0) 
441        cTree->Branch("Clusters",&clusters);
442      else b->SetAddress(&clusters);
443      
444       
445      if (itsl->TreeR() == 0x0) itsl->LoadRecPoints();
446
447      if (branch == 0x0)
448       {//it means that we produced FastRecPoints above
449        pTree=itsl->TreeR();
450        if (!pTree) { cerr<<"Can't get TreeR !\n"; return 1; }
451        branch=pTree->GetBranch("ITSRecPointsF");
452        if (!branch) { cerr<<"Can't get Fast RecPoints branch named ITSRecPointsF  !\n"; return 1;}
453       } 
454      TClonesArray *points=new TClonesArray("AliITSRecPoint",10000);
455      branch->SetAddress(&points);
456      
457      TClonesArray &cl=*clusters;
458      Int_t nclusters=0;
459      Int_t nentr=(Int_t)branch->GetEntries();
460      AliITSgeom *geom=ITS->GetITSgeom();
461      
462      cout<<"\n\n Number of entries in TreeR = "<<nentr<<endl;
463      for (Int_t i=0; i<nentr; i++) 
464       {
465        if (!branch->GetEvent(i)) 
466         {
467           cTree->Fill(); 
468           continue;
469         }
470        
471        Int_t lay,lad,det; geom->GetModuleId(i,lay,lad,det);
472        Float_t x,y,zshift; geom->GetTrans(lay,lad,det,x,y,zshift); 
473        Double_t rot[9];    geom->GetRotMatrix(lay,lad,det,rot);
474        Double_t yshift = x*rot[0] + y*rot[1];
475        Int_t ndet=(lad-1)*geom->GetNdetectors(lay) + (det-1);
476        Int_t ncl=points->GetEntriesFast();
477        nclusters+=ncl;
478
479        if (debug) ::Info("AliBarrelReconstruction.C",
480               "i=%d lay=%d lad=%d det=%d NRP=%d",
481                i,   lay,   lad,   det,   ncl);
482
483        for (Int_t j=0; j<ncl; j++) 
484         { 
485           Float_t lp[5];
486           Int_t lab[4]; 
487
488           AliITSRecPoint *p=(AliITSRecPoint*)points->UncheckedAt(j);
489           lp[0]=-p->GetX()-yshift; if (lay==1) lp[0]=-lp[0];
490           lp[1]=p->GetZ()+zshift;
491           lp[2]=p->GetSigmaX2();
492           lp[3]=p->GetSigmaZ2();
493           lp[4]=p->GetQ();
494           lab[0]=p->GetLabel(0);
495           lab[1]=p->GetLabel(1);
496           lab[2]=p->GetLabel(2);
497           lab[3]=ndet;
498      
499           Int_t label=lab[0];
500           if (label<0) continue;
501           TParticle *part=(TParticle*)rl->Stack()->Particle(label);
502           label=-3;
503           while (part->P() < 0.005) {
504             Int_t m=part->GetFirstMother();
505             if (m<0) {cerr<<"Primary momentum: "<<part->P()<<endl; break;}
506             label=m;
507             part=(TParticle*)gAlice->Particle(label);
508           }
509           if (lab[1]<0) lab[1]=label;
510           else if (lab[2]<0) lab[2]=label;
511                else cerr<<"No empty labels !\n";
512      
513           new(cl[j]) AliITSclusterV2(lab,lp);
514          }
515        cTree->Fill(); clusters->Delete();
516        points->Delete();
517      }
518      itsl->WriteRawClusters("OVERWRITE");
519      
520      cerr<<"Number of clusters: "<<nclusters<<endl;
521      delete clusters; delete points;
522
523    }
524    itsl->UnloadRecPoints();
525    itsl->UnloadRawClusters();
526    rl->UnloadKinematics();
527
528    gBenchmark->Stop(name);
529    gBenchmark->Show(name);
530
531    return rc;
532 }
533
534 Int_t ITSFindTracks(const Char_t *inname2, Int_t n) 
535 {
536
537    Int_t rc=0;
538    const Char_t *name="ITSFindTracks";
539    cerr<<'\n'<<name<<"...\n";
540    gBenchmark->Start(name);
541
542    rl->GetEvent(0);
543
544    AliITSLoader* itsl = (AliITSLoader*)rl->GetLoader("ITSLoader");
545    if ( (itsl == 0x0) || (tpcl == 0x0))
546     {
547       cerr<<"Can not get ITS Loader from Run Loader\n";
548       return 1;
549     } 
550    
551    tpcl->SetTracksFileName(TString(inname2));
552    if(rl->GetAliRun() == 0x0) rl->LoadgAlice();
553    
554    AliITS* ITS = rl->GetAliRun()->GetModule("ITS");
555    
556    rl->CdGAFile();
557    AliITSgeom *geom=ITS->GetITSgeom();
558    if (!geom) { cerr<<"can't get ITS geometry !\n"; return 1;}
559
560    itsl->LoadTracks("recreate");
561    itsl->LoadRawClusters("read");
562    for (Int_t i=0;i<n;i++)
563     {
564       rl->GetEvent(i);
565       AliITStrackerV2 tracker(geom,i,AliConfig::fgkDefaultEventFolderName);
566       rc=tracker.Clusters2Tracks();
567     }
568
569
570    gBenchmark->Stop(name);
571    gBenchmark->Show(name);
572
573    return rc;
574 }
575
576 Int_t ITSPropagateBack() 
577 {
578    
579   Int_t rc=0;
580   /*
581    const Char_t *name="ITSPropagateBack";
582    cerr<<'\n'<<name<<"...\n";
583    gBenchmark->Start(name);
584
585    AliITSgeom *geom=(AliITSgeom*)gFile->Get("AliITSgeom");
586    if (!geom) { cerr<<"can't get ITS geometry !\n"; return 1;}
587    AliITStrackerV2 tracker(geom);
588
589    TFile *out=TFile::Open(outname,"update");
590    TFile *in =TFile::Open(inname);
591    rc=tracker.PropagateBack(in,out);
592    in->Close();
593    out->Close();
594
595    gBenchmark->Stop(name);
596    gBenchmark->Show(name);
597   */
598    return rc;
599 }
600
601 Int_t TPCPropagateBack() 
602 {
603    Int_t rc=0;
604    const Char_t *name="TPCPropagateBack";
605    cerr<<'\n'<<name<<"...\n";
606    gBenchmark->Start(name);
607
608    param=(AliTPCParam *)gDirectory->Get("75x40_100x60");
609    if (!param) 
610     {
611      param=(AliTPCParam *)gDirectory->Get("75x40_100x60_150x60");
612      if (!param) 
613       {
614         cerr<<"TPC parameters have not been found !\n";
615         return 1;
616       }
617     }
618 //   param->Dump();
619    param->Update();
620
621    AliTPCtracker *tracker = new AliTPCtracker(param, AliConfig::fgkDefaultEventFolderName);
622
623 //   TFile *out=TFile::Open(outname,"update");
624 //   TFile *in =TFile::Open(inname);
625
626    rc=tracker->PropagateBack();
627    delete tracker;
628
629    gBenchmark->Stop(name);
630    gBenchmark->Show(name);
631
632    return rc;
633 }