Updated version of the PID code (from B. Batyunya)
[u/mrichter/AliRoot.git] / ITS / AliBarrelReconstructionV2.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 "AliMagF.h"
17   #include "AliTPCtracker.h"
18
19   #include "AliITS.h"
20   #include "AliITSgeom.h"
21   #include "AliITSRecPoint.h"
22   #include "AliITSclusterV2.h"
23   #include "AliITSsimulationFastPoints.h"
24   #include "AliITStrackerV2.h"
25   #include<iostream.h>
26 #endif
27
28 Int_t TPCFindClusters(const Char_t *inname, const Char_t *outname, Int_t n);
29 Int_t TPCFindTracks(const Char_t *inname, const Char_t *outname, Int_t n);
30 Int_t TPCSortTracks(const Char_t *inname, const Char_t *inname2, const Char_t *outname,  Int_t n);
31 Int_t TPCPropagateBack(const Char_t *inname, const Char_t *outname);
32
33 Int_t ITSFindClusters(const Char_t *inname,  const Char_t *outname, Int_t n);
34 Int_t ITSFindTracks(const Char_t *inname, const Char_t *inname2, const Char_t *outname, Int_t n);
35 Int_t ITSPropagateBack(const Char_t *inname, const Char_t *outname);
36
37
38 Int_t AliBarrelReconstructionV2(Int_t n=1) {
39   //const Char_t *TPCdigName="rfio:galice.root";
40    const Char_t *TPCdigName="galice.root";
41    const Char_t *TPCclsName="AliTPCclusters.root";
42    const Char_t *TPCtrkName="AliTPCtracks.root";
43    const Char_t *TPCtrkNameS="AliTPCtracksSorted.root";
44
45
46    //const Char_t *ITSdigName="rfio:galice.root";
47    const Char_t *ITSdigName="galice.root";
48    const Char_t *ITSclsName="AliITSclustersV2.root";
49    const Char_t *ITStrkName="AliITStracksV2.root";
50
51    AliKalmanTrack::SetConvConst(100/0.299792458/0.2/gAlice->Field()->Factor());
52
53     
54 // ********** Find TPC clusters *********** //
55    if (TPCFindClusters(TPCdigName,TPCclsName, n)) {
56       cerr<<"Failed to get TPC clusters !\n";
57       return 1;
58    }      
59
60 // ********** Find TPC tracks *********** //
61     if (TPCFindTracks(TPCclsName,TPCtrkName,n)) {
62       cerr<<"Failed to get TPC tracks !\n";
63       return 1;
64     }
65   
66 // ********** Sort and label TPC tracks *********** //
67    if (TPCSortTracks(TPCclsName,TPCtrkName,TPCtrkNameS,n)) {
68       cerr<<"Failed to sort TPC tracks !\n";
69       return 1;
70     }
71
72 // ********** Find ITS clusters *********** //
73    if (ITSFindClusters(ITSdigName,ITSclsName,n)) {
74       cerr<<"Failed to get ITS clusters !\n";
75       return 1;
76    }
77
78 // ********** Find ITS tracks *********** //
79    {
80      //TFile *clsFile=TFile::Open(ITSclsName);
81    if (ITSFindTracks(TPCtrkNameS,ITSclsName,ITStrkName,n)) {
82       cerr<<"Failed to get ITS tracks !\n";
83       return 1;
84    }
85    //clsFile->Close();
86    }
87    return 1;
88 // ********** Back propagation of the ITS tracks *********** //
89    {TFile *clsFile=TFile::Open(ITSclsName);
90    if (ITSPropagateBack(ITStrkName,TPCtrkNameS)) {
91       cerr<<"Failed to propagate back the ITS tracks !\n";
92       return 1;
93    }
94    clsFile->Close();}
95
96
97 // ********** Back propagation of the TPC tracks *********** //
98    {TFile *clsFile=TFile::Open(TPCclsName);
99    if (TPCPropagateBack(TPCtrkName,TPCtrkName)) {
100       cerr<<"Failed to propagate back the TPC tracks !\n";
101       return 1;
102    }
103    clsFile->Close();}
104
105    return 0;
106 }
107
108
109 Int_t TPCFindClusters(const Char_t *inname, const Char_t *outname, Int_t n) {
110    Int_t rc=0;
111    const Char_t *name="TPCFindClusters";
112    cerr<<'\n'<<name<<"...\n";
113    gBenchmark->Start(name);
114
115    TFile *out=TFile::Open(outname,"recreate");
116    TFile *in =TFile::Open(inname);
117
118    AliTPCParam *param=(AliTPCParam *)in->Get("75x40_100x60");
119    if (!param) {cerr<<"TPC parameters have not been found !\n"; return 1;}
120    AliTPCv2 tpc;
121    tpc.SetParam(param);
122
123    //tpc.Digits2Clusters(out); //MI change
124    cout<<"TPCFindClusters: nev ="<<n<<endl;
125    for (Int_t i=0;i<n;i++){
126      printf("Processing event %d\n",i);
127      tpc.Digits2Clusters(out,i);
128      //  AliTPCclusterer::Digits2Clusters(dig, out, i);
129    }
130    in->Close();
131    out->Close();
132    gBenchmark->Stop(name);
133    gBenchmark->Show(name);
134
135    return rc;
136 }
137
138 Int_t TPCFindTracks(const Char_t *inname, const Char_t *outname, Int_t n) {
139    Int_t rc=0;
140    const Char_t *name="TPCFindTracks";
141    cerr<<'\n'<<name<<"...\n";
142    gBenchmark->Start(name);
143    TFile *out=TFile::Open(outname,"recreate");
144    TFile *in =TFile::Open(inname);
145    AliTPCParam *param=(AliTPCParam *)in->Get("75x40_100x60");
146    if (!param) {cerr<<"TPC parameters have not been found !\n"; return 1;}
147
148    //AliTPCtracker *tracker=new AliTPCtracker(param);
149    //rc=tracker->Clusters2Tracks(0,out);
150    //delete tracker;
151    cout<<"TPCFindTracks: nev ="<<n<<endl;
152
153    for (Int_t i=0;i<n;i++){
154      printf("Processing event %d\n",i);
155      AliTPCtracker *tracker = new AliTPCtracker(param,i);
156      //Int_t rc=
157        tracker->Clusters2Tracks(0,out);
158      delete tracker;
159    }
160
161    in->Close();
162    out->Close();
163    gBenchmark->Stop(name);
164    gBenchmark->Show(name);
165
166    return rc;
167 }
168 Int_t TPCSortTracks(const Char_t *inname, const Char_t *inname2, const Char_t * outname,  Int_t eventn){
169    Int_t rc=0;
170    const Char_t *name="TPCSortTracks";
171    cerr<<'\n'<<name<<"...\n";
172    gBenchmark->Start(name);
173
174    TFile *out =TFile::Open(outname,"recreate");
175    TFile *in1 =TFile::Open(inname);
176    TFile *in2 =TFile::Open(inname2);
177
178
179    AliTPCParam *param=(AliTPCParam *)in1->Get("75x40_100x60");
180    if (!param) {cerr<<"TPC parameters have not been found !\n"; return 1;}
181
182    cout<<"TPCSortedtracks: nev ="<<eventn<<endl;
183
184
185    // loop over events 
186    for (Int_t event=0;event<eventn; event++){
187      TObjArray tarray(10000);
188      AliTPCtrack *iotrack=0;
189      Int_t i;
190
191
192      in2->cd();
193      char   tname[100];
194      sprintf(tname,"TreeT_TPC_%d",event);
195
196      TTree *tracktree=(TTree*)in2->Get(tname);
197      TBranch *tbranch=tracktree->GetBranch("tracks");
198      Int_t nentr=(Int_t)tracktree->GetEntries();
199      for (i=0; i<nentr; i++) {
200        iotrack=new AliTPCtrack;
201        tbranch->SetAddress(&iotrack);
202        tracktree->GetEvent(i);
203        tarray.AddLast(iotrack);
204      }   
205      tarray.Sort();
206      // in2->Close();
207      
208      //assign thacks GEANT labels
209      in1->cd();
210      AliTPCtracker *tracker = new AliTPCtracker(param,event);
211      tracker->LoadInnerSectors();
212      tracker->LoadOuterSectors();
213      for (i=0; i<nentr; i++) {
214        iotrack=(AliTPCtrack*)tarray.UncheckedAt(i);
215        tracker->CookLabel(iotrack,0.1);
216      }   
217      delete tracker;
218      //in->Close();
219      //end of GEANT label assignment
220      
221      tracktree=new TTree(tname,"Tree with TPC tracks");
222      tracktree->Branch("tracks","AliTPCtrack",&iotrack,32000,0);
223      for (i=0; i<nentr; i++) {
224        iotrack=(AliTPCtrack*)tarray.UncheckedAt(i);
225        tracktree->Fill();
226      }
227      out->cd();
228      tracktree->Write();
229    }
230
231    out->Close();
232    in2->Close();
233    in1->Close();
234
235    gBenchmark->Stop(name);
236    gBenchmark->Show(name);
237
238    return rc;
239 }
240
241 Int_t ITSFindClusters(const Char_t *inname, const Char_t *outname, Int_t n) {
242    Int_t rc=0;
243    const Char_t *name="ITSFindClusters";
244    cerr<<'\n'<<name<<"...\n";
245    gBenchmark->Start(name);
246    TFile *out=TFile::Open(outname,"recreate");
247    TFile *in =TFile::Open(inname,"update");
248
249    if (!(gAlice=(AliRun*)in->Get("gAlice"))) {
250       cerr<<"Can't get gAlice !\n";
251       return 1;
252    }
253
254    AliITS *ITS  = (AliITS*)gAlice->GetModule("ITS");
255    if (!ITS) { cerr<<"Can't get the ITS !\n"; return 1;}
256    AliITSgeom *geom=ITS->GetITSgeom();
257    out->cd();   
258    geom->Write();
259      
260    cout<<"ITSFindClusters: nev ="<<n<<endl;
261    Int_t ev=0;
262    for (ev = 0; ev<n; ev++){
263      in->cd();   // !!!! MI directory must point to galice. - othervise problem with Tree -connection
264      gAlice->GetEvent(ev);
265      //gAlice->TreeR()->Reset();   //reset reconstructed tree
266
267      
268      TTree *pTree=gAlice->TreeR();
269      if (!pTree){
270        gAlice->MakeTree("R");
271        pTree = gAlice->TreeR();
272      }
273      TBranch *branch=pTree->GetBranch("ITSRecPoints");
274      /*
275      if (!branch) {
276        //if not reconstructed ITS branch do reconstruction 
277        ITS->MakeBranch("R",0);
278
279        //////////////// Taken from ITSHitsToFastPoints.C ///////////////////////
280        AliITSsimulationFastPoints *sim = new AliITSsimulationFastPoints();
281        for (Int_t i=0;i<3;i++) { ITS->SetSimulationModel(i,sim); }
282        Int_t nsignal=25;
283        Int_t size=-1;
284        Int_t bgr_ev=Int_t(ev/nsignal);
285        ITS->HitsToFastRecPoints(ev,bgr_ev,size," ","All"," ");
286        //////////////////////////////////////////////////////////////////////////
287        gAlice->GetEvent(ev);   //MI comment  - in HitsToFast... they reset treeR to 0 
288        //they overwrite full reconstructed event ???? ... so lets connect TreeR one more
289        //time
290      }
291      */
292
293      
294      out->cd();
295      TClonesArray *clusters=new TClonesArray("AliITSclusterV2",10000);
296      char   cname[100];
297      sprintf(cname,"TreeC_ITS_%d",ev);
298   
299      TTree *cTree=new TTree(cname,"ITS clusters");
300      cTree->Branch("Clusters",&clusters);
301      
302      pTree=gAlice->TreeR();
303      if (!pTree) { cerr<<"Can't get TreeR !\n"; return 1; }
304      branch=pTree->GetBranch("ITSRecPoints");
305      if (!branch) { cerr<<"Can't get ITSRecPoints branch !\n"; return 1;}
306      TClonesArray *points=new TClonesArray("AliITSRecPoint",10000);
307      branch->SetAddress(&points);
308      
309      TClonesArray &cl=*clusters;
310      Int_t nclusters=0;
311      Int_t nentr=(Int_t)pTree->GetEntries();
312      Int_t lab[6]; 
313      Float_t lp[5];
314      AliITSgeom *geom=ITS->GetITSgeom();
315
316      cout<<"!!!! nentr ="<<nentr<<endl;
317      for (Int_t i=0; i<nentr; i++) {
318        if (!pTree->GetEvent(i)) {cTree->Fill(); continue;}
319        Int_t lay,lad,det; geom->GetModuleId(i,lay,lad,det);
320        Float_t x,y,zshift; geom->GetTrans(lay,lad,det,x,y,zshift); 
321        Double_t rot[9];    geom->GetRotMatrix(lay,lad,det,rot);
322        Double_t yshift = x*rot[0] + y*rot[1];
323        Int_t ndet=(lad-1)*geom->GetNdetectors(lay) + (det-1);
324        Int_t ncl=points->GetEntriesFast();
325        nclusters+=ncl;
326
327      Float_t kmip=1.;
328      if(lay<5&&lay>2){kmip=280.;}; // b.b.
329      if(lay<7&&lay>4){kmip=38.;};
330
331
332        for (Int_t j=0; j<ncl; j++) {
333          AliITSRecPoint *p=(AliITSRecPoint*)points->UncheckedAt(j);
334          lp[0]=-p->GetX()-yshift; if (lay==1) lp[0]=-lp[0];
335          lp[1]=p->GetZ()+zshift;
336          lp[2]=p->GetSigmaX2();
337          lp[3]=p->GetSigmaZ2();
338          lp[4]=p->GetQ(); lp[4]/=kmip;    // b.b.
339          if(ncl==11)cout<<"mod,cl,Q ="<<i<<","<<j<<","<<lp[4]<<endl;
340          lab[0]=p->GetLabel(0);lab[1]=p->GetLabel(1);lab[2]=p->GetLabel(2);
341          lab[3]=ndet;
342          
343          Int_t label=lab[0];
344         if(label>=0) { // b.b.
345          TParticle *part=(TParticle*)gAlice->Particle(label);
346          label=-3;
347          while (part->P() < 0.005) {
348            Int_t m=part->GetFirstMother();
349            if (m<0) {cerr<<"Primary momentum: "<<part->P()<<endl; break;}
350            label=m;
351            part=(TParticle*)gAlice->Particle(label);
352          }
353         }
354          if      (lab[1]<0) lab[1]=label;
355          else if (lab[2]<0) lab[2]=label;
356          else cerr<<"No empty labels !\n";
357          
358          new(cl[j]) AliITSclusterV2(lab,lp);
359        }
360        cTree->Fill(); clusters->Delete();
361        points->Delete();
362      }
363      cTree->Write();
364      cerr<<"Number of clusters: "<<nclusters<<endl;
365      delete cTree; delete clusters; delete points;
366
367    }
368
369    delete gAlice; gAlice=0;
370    in->Close();
371    out->Close();
372    gBenchmark->Stop(name);
373    gBenchmark->Show(name);
374
375    return rc;
376 }
377
378 Int_t ITSFindTracks(const Char_t * inname, const Char_t *inname2, const Char_t *outname, Int_t n) {
379    Int_t rc=0;
380    const Char_t *name="ITSFindTracks";
381    cerr<<'\n'<<name<<"...\n";
382    gBenchmark->Start(name);
383
384
385    TFile *out=TFile::Open(outname,"recreate");
386    TFile *in =TFile::Open(inname);
387    TFile *in2 =TFile::Open(inname2);
388
389    AliITSgeom *geom=(AliITSgeom*)gFile->Get("AliITSgeom");
390    if (!geom) { cerr<<"can't get ITS geometry !\n"; return 1;}
391
392    cout<<"ITSFindtracks: nev ="<<n<<endl;
393
394    for (Int_t i=0;i<n;i++){
395      AliITStrackerV2 tracker(geom,i);
396      rc=tracker.Clusters2Tracks(in,out);
397    }
398
399    in->Close();
400    in2->Close();
401    out->Close();
402
403    gBenchmark->Stop(name);
404    gBenchmark->Show(name);
405
406    return rc;
407 }
408
409 Int_t ITSPropagateBack(const Char_t *inname, const Char_t *outname) {
410    
411   Int_t rc=0;
412   /*
413    const Char_t *name="ITSPropagateBack";
414    cerr<<'\n'<<name<<"...\n";
415    gBenchmark->Start(name);
416
417    AliITSgeom *geom=(AliITSgeom*)gFile->Get("AliITSgeom");
418    if (!geom) { cerr<<"can't get ITS geometry !\n"; return 1;}
419    AliITStrackerV2 tracker(geom);
420
421    TFile *out=TFile::Open(outname,"update");
422    TFile *in =TFile::Open(inname);
423    rc=tracker.PropagateBack(in,out);
424    in->Close();
425    out->Close();
426
427    gBenchmark->Stop(name);
428    gBenchmark->Show(name);
429   */
430    return rc;
431 }
432
433 Int_t TPCPropagateBack(const Char_t *inname, const Char_t *outname) {
434    Int_t rc=0;
435    const Char_t *name="TPCPropagateBack";
436    cerr<<'\n'<<name<<"...\n";
437    gBenchmark->Start(name);
438
439    AliTPCParam *param=(AliTPCParam *)gFile->Get("75x40_100x60");
440    if (!param) {cerr<<"TPC parameters have not been found !\n"; return 1;}
441    AliTPCtracker *tracker=new AliTPCtracker(param);
442
443    TFile *out=TFile::Open(outname,"update");
444    TFile *in =TFile::Open(inname);
445    rc=tracker->PropagateBack(in,out);
446    delete tracker;
447    in->Close();
448    out->Close();
449
450    gBenchmark->Stop(name);
451    gBenchmark->Show(name);
452
453    return rc;
454 }
455
456
457
458
459