Updated V2 stream of tracking (Yu.Belikov). The new long waited features are: 1)...
[u/mrichter/AliRoot.git] / ITS / AliITSFindClustersV2.C
1 #ifndef __CINT__
2   #include <iostream.h>
3
4   #include "AliRun.h"
5   #include "AliITS.h"
6   #include "AliITSgeom.h"
7   #include "AliITSRecPoint.h"
8   #include "AliITSclusterV2.h"
9   #include "AliITSsimulationFastPoints.h"
10
11   #include "TFile.h"
12   #include "TTree.h"
13   #include "TParticle.h"
14 #endif
15
16 Int_t AliITSFindClustersV2() {
17 /****************************************************************
18  *  Just something to start with                                *
19  ****************************************************************/
20    cerr<<"Looking for clusters...\n";
21
22    if (gAlice) {delete gAlice; gAlice=0;}
23
24    TFile *in=TFile::Open("galice.root","update");
25    if (!in->IsOpen()) {
26       cerr<<"Can't open galice.root !\n"; 
27       return 1;
28    }
29    
30    if (!(gAlice=(AliRun*)in->Get("gAlice"))) {
31       cerr<<"Can't find gAlice !\n";
32       return 2;
33    }
34
35    Int_t ev=0;
36    gAlice->GetEvent(ev);
37
38    AliITS *ITS  = (AliITS*)gAlice->GetModule("ITS");
39    if (!ITS) {
40       cerr<<"Can't find the ITS !\n";
41       return 3;
42    }
43
44    gAlice->MakeTree("R"); ITS->MakeBranch("R",0);
45 //////////////// Taken from ITSHitsToFastPoints.C ///////////////////////
46    ITS->SetSimulationModel(0,new AliITSsimulationFastPoints());
47    ITS->SetSimulationModel(1,new AliITSsimulationFastPoints());
48    ITS->SetSimulationModel(2,new AliITSsimulationFastPoints());
49    Int_t nsignal=25;
50    Int_t size=-1;
51    Int_t bgr_ev=Int_t(ev/nsignal);
52    ITS->HitsToFastRecPoints(ev,bgr_ev,size," ","All"," ");
53 //////////////////////////////////////////////////////////////////////////
54
55    delete gAlice; gAlice=0;
56    in->Close();
57
58 ///////////////// Conversion AliITSRecPoint -> AliITSclusterV2 //////////////
59    /*TFile */in=TFile::Open("galice.root");
60
61    if (gAlice) {delete gAlice; gAlice=0;}
62
63    if (!(gAlice=(AliRun*)in->Get("gAlice"))) {
64       cerr<<"Can't find gAlice !\n";
65       return 4;
66    }
67
68    gAlice->GetEvent(0);
69
70    /*AliITS */ITS  = (AliITS*)gAlice->GetModule("ITS");
71    if (!ITS) {
72       cerr<<"Can't find the ITS !\n";
73       return 5;
74    }
75    AliITSgeom *geom=ITS->GetITSgeom();
76  
77    TFile *out=TFile::Open("AliITSclustersV2.root","new");
78    if (!out->IsOpen()) {
79       cerr<<"Delete old AliITSclustersV2.root !\n"; 
80       return 6;
81    }
82    geom->Write();
83
84    TClonesArray *clusters=new TClonesArray("AliITSclusterV2",10000);
85    TTree *cTree=new TTree("TreeC_ITS_0","ITS clusters");
86    cTree->Branch("Clusters",&clusters);
87
88    TTree *pTree=gAlice->TreeR();
89    if (!pTree) {
90       cerr<<"Can't get TreeR !\n";
91       return 7;
92    }
93    TBranch *branch=pTree->GetBranch("ITSRecPoints");
94    if (!branch) { 
95       cerr<<"Can't get ITSRecPoints branch !\n";
96       return 8;
97    }
98    TClonesArray *points=new TClonesArray("AliITSRecPoint",10000);
99    branch->SetAddress(&points);
100
101    TClonesArray &cl=*clusters;
102    Int_t nclusters=0;
103    Int_t nentr=(Int_t)pTree->GetEntries();
104
105    cerr<<"Number of entries: "<<nentr<<endl;
106
107    Float_t lp[5]; Int_t lab[6]; //Why can't it be inside a loop ?
108
109    for (Int_t i=0; i<nentr; i++) {
110        points->Clear();
111        pTree->GetEvent(i);
112        Int_t ncl=points->GetEntriesFast(); if (ncl==0){cTree->Fill();continue;}
113        Int_t lay,lad,det; geom->GetModuleId(i,lay,lad,det);
114        Float_t x,y,zshift; geom->GetTrans(lay,lad,det,x,y,zshift); 
115        Double_t rot[9];    geom->GetRotMatrix(lay,lad,det,rot);
116        Double_t yshift = x*rot[0] + y*rot[1];
117        Int_t ndet=(lad-1)*geom->GetNdetectors(lay) + (det-1);
118        nclusters+=ncl;
119        for (Int_t j=0; j<ncl; j++) {
120           AliITSRecPoint *p=(AliITSRecPoint*)points->UncheckedAt(j);
121           //Float_t lp[5];
122           lp[0]=-p->GetX()-yshift; if (lay==1) lp[0]=-lp[0];
123           lp[1]=p->GetZ()+zshift;
124           lp[2]=p->GetSigmaX2();
125           lp[3]=p->GetSigmaZ2();
126           lp[4]=p->GetQ();
127           //Int_t lab[6]; 
128           lab[0]=p->GetLabel(0);lab[1]=p->GetLabel(1);lab[2]=p->GetLabel(2);
129           lab[3]=ndet;
130
131           Int_t label=lab[0];
132           TParticle *part=(TParticle*)gAlice->Particle(label);
133           label=-3;
134           while (part->P() < 0.005) {
135              Int_t m=part->GetFirstMother();
136              if (m<0) {cerr<<"Primary momentum: "<<part->P()<<endl; break;}
137              label=m;
138              part=(TParticle*)gAlice->Particle(label);
139           }
140           if      (lab[1]<0) lab[1]=label;
141           else if (lab[2]<0) lab[2]=label;
142           else cerr<<"No empty labels !\n";
143
144           new(cl[j]) AliITSclusterV2(lab,lp);
145        }
146        cTree->Fill(); clusters->Delete();
147    }
148    cTree->Write();
149
150    cerr<<"Number of clusters: "<<nclusters<<endl;
151
152    delete cTree; delete clusters; delete points;
153
154    delete gAlice; gAlice=0;
155
156    in->Close();
157    out->Close();
158
159    return 0;
160
161 }
162
163
164
165
166
167
168
169
170
171
172
173
174
175