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 ****************************************************************************/
17 #include "AliTPCtracker.h"
20 #include "AliITSgeom.h"
21 #include "AliITSRecPoint.h"
22 #include "AliITSclusterV2.h"
23 #include "AliITSsimulationFastPoints.h"
24 #include "AliITStrackerV2.h"
27 Int_t TPCFindClusters(const Char_t *inname, const Char_t *outname, Int_t n);
28 Int_t TPCFindTracks(const Char_t *inname, const Char_t *outname, Int_t n);
29 Int_t TPCSortTracks(const Char_t *inname, const Char_t *inname2, const Char_t *outname, Int_t n);
30 Int_t TPCPropagateBack(const Char_t *inname, const Char_t *outname);
32 Int_t ITSFindClusters(const Char_t *inname, const Char_t *outname, Int_t n);
33 Int_t ITSFindTracks(const Char_t *inname, const Char_t *inname2, const Char_t *outname, Int_t n);
34 Int_t ITSPropagateBack(const Char_t *inname, const Char_t *outname);
37 Int_t AliBarrelReconstruction4ITS(Int_t n=1) {
38 const Char_t *TPCdigName="rfio:galice.root";
39 const Char_t *TPCclsName="AliTPCclusters.root";
40 const Char_t *TPCtrkName="AliTPCtracks.root";
41 const Char_t *TPCtrkNameS="AliTPCtracksSorted.root";
44 const Char_t *ITSdigName="rfio:galice.root";
45 const Char_t *ITSclsName="AliITSclustersV2.root";
46 const Char_t *ITStrkName="AliITStracksV2.root";
48 AliKalmanTrack::SetConvConst(100/0.299792458/0.2/gAlice->Field()->Factor());
51 // ********** Find TPC clusters *********** //
52 if (TPCFindClusters(TPCdigName,TPCclsName, n)) {
53 cerr<<"Failed to get TPC clusters !\n";
57 // ********** Find TPC tracks *********** //
58 if (TPCFindTracks(TPCclsName,TPCtrkName,n)) {
59 cerr<<"Failed to get TPC tracks !\n";
63 // ********** Sort and label TPC tracks *********** //
64 if (TPCSortTracks(TPCclsName,TPCtrkName,TPCtrkNameS,n)) {
65 cerr<<"Failed to sort TPC tracks !\n";
69 // ********** Find ITS clusters *********** //
70 if (ITSFindClusters(ITSdigName,ITSclsName,n)) {
71 cerr<<"Failed to get ITS clusters !\n";
75 // ********** Find ITS tracks *********** //
77 //TFile *clsFile=TFile::Open(ITSclsName);
78 if (ITSFindTracks(TPCtrkNameS,ITSclsName,ITStrkName,n)) {
79 cerr<<"Failed to get ITS tracks !\n";
85 // ********** Back propagation of the ITS tracks *********** //
86 {TFile *clsFile=TFile::Open(ITSclsName);
87 if (ITSPropagateBack(ITStrkName,TPCtrkNameS)) {
88 cerr<<"Failed to propagate back the ITS tracks !\n";
94 // ********** Back propagation of the TPC tracks *********** //
95 {TFile *clsFile=TFile::Open(TPCclsName);
96 if (TPCPropagateBack(TPCtrkName,TPCtrkName)) {
97 cerr<<"Failed to propagate back the TPC tracks !\n";
106 Int_t TPCFindClusters(const Char_t *inname, const Char_t *outname, Int_t n) {
108 const Char_t *name="TPCFindClusters";
109 cerr<<'\n'<<name<<"...\n";
110 gBenchmark->Start(name);
112 TFile *out=TFile::Open(outname,"recreate");
113 TFile *in =TFile::Open(inname);
115 AliTPCParam *param=(AliTPCParam *)in->Get("75x40_100x60");
116 if (!param) {cerr<<"TPC parameters have not been found !\n"; return 1;}
120 //tpc.Digits2Clusters(out); //MI change
121 for (Int_t i=0;i<n;i++){
122 printf("Processing event %d\n",i);
123 tpc.Digits2Clusters(out,i);
124 // AliTPCclusterer::Digits2Clusters(dig, out, i);
128 gBenchmark->Stop(name);
129 gBenchmark->Show(name);
134 Int_t TPCFindTracks(const Char_t *inname, const Char_t *outname, Int_t n) {
136 const Char_t *name="TPCFindTracks";
137 cerr<<'\n'<<name<<"...\n";
138 gBenchmark->Start(name);
139 TFile *out=TFile::Open(outname,"recreate");
140 TFile *in =TFile::Open(inname);
141 AliTPCParam *param=(AliTPCParam *)in->Get("75x40_100x60");
142 if (!param) {cerr<<"TPC parameters have not been found !\n"; return 1;}
144 //AliTPCtracker *tracker=new AliTPCtracker(param);
145 //rc=tracker->Clusters2Tracks(0,out);
148 for (Int_t i=0;i<n;i++){
149 printf("Processing event %d\n",i);
150 AliTPCtracker *tracker = new AliTPCtracker(param,i);
152 tracker->Clusters2Tracks(0,out);
158 gBenchmark->Stop(name);
159 gBenchmark->Show(name);
163 Int_t TPCSortTracks(const Char_t *inname, const Char_t *inname2, const Char_t * outname, Int_t eventn){
165 const Char_t *name="TPCSortTracks";
166 cerr<<'\n'<<name<<"...\n";
167 gBenchmark->Start(name);
169 TFile *out =TFile::Open(outname,"recreate");
170 TFile *in1 =TFile::Open(inname);
171 TFile *in2 =TFile::Open(inname2);
174 AliTPCParam *param=(AliTPCParam *)in1->Get("75x40_100x60");
175 if (!param) {cerr<<"TPC parameters have not been found !\n"; return 1;}
180 for (Int_t event=0;event<eventn; event++){
181 TObjArray tarray(10000);
182 AliTPCtrack *iotrack=0;
188 sprintf(tname,"TreeT_TPC_%d",event);
190 TTree *tracktree=(TTree*)in2->Get(tname);
191 TBranch *tbranch=tracktree->GetBranch("tracks");
192 Int_t nentr=(Int_t)tracktree->GetEntries();
193 for (i=0; i<nentr; i++) {
194 iotrack=new AliTPCtrack;
195 tbranch->SetAddress(&iotrack);
196 tracktree->GetEvent(i);
197 tarray.AddLast(iotrack);
202 //assign thacks GEANT labels
204 AliTPCtracker *tracker = new AliTPCtracker(param,event);
205 tracker->LoadInnerSectors();
206 tracker->LoadOuterSectors();
207 for (i=0; i<nentr; i++) {
208 iotrack=(AliTPCtrack*)tarray.UncheckedAt(i);
209 tracker->CookLabel(iotrack,0.1);
213 //end of GEANT label assignment
215 tracktree=new TTree(tname,"Tree with TPC tracks");
216 tracktree->Branch("tracks","AliTPCtrack",&iotrack,32000,0);
217 for (i=0; i<nentr; i++) {
218 iotrack=(AliTPCtrack*)tarray.UncheckedAt(i);
229 gBenchmark->Stop(name);
230 gBenchmark->Show(name);
235 Int_t ITSFindClusters(const Char_t *inname, const Char_t *outname, Int_t n) {
237 const Char_t *name="ITSFindClusters";
238 cerr<<'\n'<<name<<"...\n";
239 gBenchmark->Start(name);
240 TFile *out=TFile::Open(outname,"recreate");
241 TFile *in =TFile::Open(inname,"update");
243 if (!(gAlice=(AliRun*)in->Get("gAlice"))) {
244 cerr<<"Can't get gAlice !\n";
248 AliITS *ITS = (AliITS*)gAlice->GetModule("ITS");
249 if (!ITS) { cerr<<"Can't get the ITS !\n"; return 1;}
250 AliITSgeom *geom=ITS->GetITSgeom();
255 for (ev = 0; ev<n; ev++){
256 in->cd(); // !!!! MI directory must point to galice. - othervise problem with Tree -connection
257 gAlice->GetEvent(ev);
258 //gAlice->TreeR()->Reset(); //reset reconstructed tree
261 TTree *pTree=gAlice->TreeR();
263 gAlice->MakeTree("R");
264 pTree = gAlice->TreeR();
266 TBranch *branch=pTree->GetBranch("ITSRecPoints");
268 //if not reconstructed ITS branch do reconstruction
269 ITS->MakeBranch("R",0);
270 //////////////// Taken from ITSHitsToFastPoints.C ///////////////////////
271 AliITSsimulationFastPoints *sim = new AliITSsimulationFastPoints();
272 for (Int_t i=0;i<3;i++) { ITS->SetSimulationModel(i,sim); }
275 Int_t bgr_ev=Int_t(ev/nsignal);
276 ITS->HitsToFastRecPoints(ev,bgr_ev,size," ","All"," ");
277 //////////////////////////////////////////////////////////////////////////
278 gAlice->GetEvent(ev); //MI comment - in HitsToFast... they reset treeR to 0
279 //they overwrite full reconstructed event ???? ... so lets connect TreeR one more
286 TClonesArray *clusters=new TClonesArray("AliITSclusterV2",10000);
288 sprintf(cname,"TreeC_ITS_%d",ev);
290 TTree *cTree=new TTree(cname,"ITS clusters");
291 cTree->Branch("Clusters",&clusters);
293 pTree=gAlice->TreeR();
294 if (!pTree) { cerr<<"Can't get TreeR !\n"; return 1; }
295 branch=pTree->GetBranch("ITSRecPoints");
296 if (!branch) { cerr<<"Can't get ITSRecPoints branch !\n"; return 1;}
297 TClonesArray *points=new TClonesArray("AliITSRecPoint",10000);
298 branch->SetAddress(&points);
300 TClonesArray &cl=*clusters;
302 Int_t nentr=(Int_t)pTree->GetEntries();
303 AliITSgeom *geom=ITS->GetITSgeom();
305 for (Int_t i=0; i<nentr; i++) {
306 if (!pTree->GetEvent(i)) {cTree->Fill(); continue;}
307 Int_t lay,lad,det; geom->GetModuleId(i,lay,lad,det);
308 Float_t x,y,zshift; geom->GetTrans(lay,lad,det,x,y,zshift);
309 Double_t rot[9]; geom->GetRotMatrix(lay,lad,det,rot);
310 Double_t yshift = x*rot[0] + y*rot[1];
311 Int_t ndet=(lad-1)*geom->GetNdetectors(lay) + (det-1);
312 Int_t ncl=points->GetEntriesFast();
314 for (Int_t j=0; j<ncl; j++) {
315 AliITSRecPoint *p=(AliITSRecPoint*)points->UncheckedAt(j);
317 lp[0]=-p->GetX()-yshift; if (lay==1) lp[0]=-lp[0];
318 lp[1]=p->GetZ()+zshift;
319 lp[2]=p->GetSigmaX2();
320 lp[3]=p->GetSigmaZ2();
323 lab[0]=p->GetLabel(0);lab[1]=p->GetLabel(1);lab[2]=p->GetLabel(2);
327 TParticle *part=(TParticle*)gAlice->Particle(label);
329 while (part->P() < 0.005) {
330 Int_t m=part->GetFirstMother();
331 if (m<0) {cerr<<"Primary momentum: "<<part->P()<<endl; break;}
333 part=(TParticle*)gAlice->Particle(label);
335 if (lab[1]<0) lab[1]=label;
336 else if (lab[2]<0) lab[2]=label;
337 else cerr<<"No empty labels !\n";
339 new(cl[j]) AliITSclusterV2(lab,lp);
341 cTree->Fill(); clusters->Delete();
345 cerr<<"Number of clusters: "<<nclusters<<endl;
346 delete cTree; delete clusters; delete points;
350 delete gAlice; gAlice=0;
353 gBenchmark->Stop(name);
354 gBenchmark->Show(name);
359 Int_t ITSFindTracks(const Char_t * inname, const Char_t *inname2, const Char_t *outname, Int_t n) {
361 const Char_t *name="ITSFindTracks";
362 cerr<<'\n'<<name<<"...\n";
363 gBenchmark->Start(name);
366 TFile *out=TFile::Open(outname,"recreate");
367 TFile *in =TFile::Open(inname);
368 TFile *in2 =TFile::Open(inname2);
370 AliITSgeom *geom=(AliITSgeom*)gFile->Get("AliITSgeom");
371 if (!geom) { cerr<<"can't get ITS geometry !\n"; return 1;}
374 for (Int_t i=0;i<n;i++){
375 AliITStrackerV2 tracker(geom,i);
376 rc=tracker.Clusters2Tracks(in,out);
383 gBenchmark->Stop(name);
384 gBenchmark->Show(name);
389 Int_t ITSPropagateBack(const Char_t *inname, const Char_t *outname) {
393 const Char_t *name="ITSPropagateBack";
394 cerr<<'\n'<<name<<"...\n";
395 gBenchmark->Start(name);
397 AliITSgeom *geom=(AliITSgeom*)gFile->Get("AliITSgeom");
398 if (!geom) { cerr<<"can't get ITS geometry !\n"; return 1;}
399 AliITStrackerV2 tracker(geom);
401 TFile *out=TFile::Open(outname,"update");
402 TFile *in =TFile::Open(inname);
403 rc=tracker.PropagateBack(in,out);
407 gBenchmark->Stop(name);
408 gBenchmark->Show(name);
413 Int_t TPCPropagateBack(const Char_t *inname, const Char_t *outname) {
415 const Char_t *name="TPCPropagateBack";
416 cerr<<'\n'<<name<<"...\n";
417 gBenchmark->Start(name);
419 AliTPCParam *param=(AliTPCParam *)gFile->Get("75x40_100x60");
420 if (!param) {cerr<<"TPC parameters have not been found !\n"; return 1;}
421 AliTPCtracker *tracker=new AliTPCtracker(param);
423 TFile *out=TFile::Open(outname,"update");
424 TFile *in =TFile::Open(inname);
425 rc=tracker->PropagateBack(in,out);
430 gBenchmark->Stop(name);
431 gBenchmark->Show(name);