Compatible with the new tracking V1
[u/mrichter/AliRoot.git] / ITS / AliITSTrackerV1.cxx
CommitLineData
13888096 1#include <iostream.h>
2#include <TROOT.h>
3#include <TMath.h>
4#include <TRandom.h>
5#include <TBranch.h>
6#include <TVector.h>
7#include <TClonesArray.h>
8#include <TFile.h>
9#include <TTree.h>
10
11
12#include "TParticle.h"
13#include "AliRun.h"
14#include "AliITS.h"
15#include "AliITSgeomSPD.h"
16#include "AliITSgeomSDD.h"
17#include "AliITSgeomSSD.h"
18#include "AliITSgeom.h"
19#include "AliITSmodule.h"
20#include "AliITSRecPoint.h"
21#include "AliMC.h"
22#include "AliKalmanTrack.h"
23#include "AliMagF.h"
24
25
26#include "AliITSRad.h"
27#include "AliITStrack.h"
28#include "AliITSiotrack.h"
29#include "AliITStracking.h"
30#include "../TPC/AliTPC.h"
31#include "../TPC/AliTPCParam.h"
32#include "../TPC/AliTPCtracker.h"
33
34#include "AliITSgeoinfo.h"
35#include "AliITSTrackerV1.h"
36
37
38
39ClassImp(AliITSTrackerV1)
40
41
42//________________________________________________________________
43
44AliITSTrackerV1::AliITSTrackerV1(AliITS* IITTSS) {
45//Origin A. Badala' and G.S. Pappalardo: e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it
46// default constructor
47 ITS = IITTSS;
48
49}
50
51//________________________________________________________________
52
53AliITStrack AliITSTrackerV1::Tracking(AliITStrack &track, AliITStrack *reference,TObjArray *fastpoints,
54Int_t **vettid, Bool_t flagvert, AliITSRad *rl, AliITSgeoinfo *geoinfo) {
55
56//Origin A. Badala' and G.S. Pappalardo: e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it
57
58 TList *list= new TList();
59
60 AliITStrack tr(track);
61
62 list->AddLast(&tr);
63
64 Double_t Pt=(tr).GetPt();
65 //cout << "\n Pt = " << Pt <<"\n"; //stampa
66
67 AliITStracking obj(list, reference, ITS, fastpoints,TMath::Abs(Pt),vettid, flagvert, rl, geoinfo); // nuova ITS
68 list->Delete();
69 delete list;
70
71 Int_t itot=-1;
72 TVector VecTotLabref(18);
73 Int_t lay, k;
74 for(lay=5; lay>=0; lay--) {
75 TVector VecLabref(3);
76 VecLabref=(*reference).GetLabTrack(lay);
77 Float_t ClustZ=(*reference).GetZclusterTrack( lay);
78 for(k=0; k<3; k++){
79 Int_t lpp=(Int_t)VecLabref(k);
80 if(lpp>=0) {
81 TParticle *p=(TParticle*) gAlice->Particle(lpp);
82 Int_t pcode=p->GetPdgCode();
83 if(pcode==11) VecLabref(k)=p->GetFirstMother();
84 }
85 itot++; VecTotLabref(itot)=VecLabref(k);
86 if(VecLabref(k)==0. && ClustZ == 0.) VecTotLabref(itot) =-3.; }
87 }
88 Long_t labref;
89 Int_t freq;
90 (*reference).Search(VecTotLabref, labref, freq);
91
92 //if(freq < 6) labref=-labref; // cinque - sei
93 if(freq < 5) labref=-labref; // cinque - sei
94 (*reference).SetLabel(labref);
95
96 return *reference;
97
98}
99
100
101
102//________________________________________________________________
103
104
105
106void AliITSTrackerV1::DoTracking(Int_t evNumber, Int_t min_t, Int_t max_t, TFile *file, Bool_t flagvert) {
107//Origin A. Badala' and G.S. Pappalardo: e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it
108// ex macro for tracking ITS
109
583731c7 110 //Loop variables
111
112 Int_t i;
113
13888096 114 printf("begin DoTracking - file %p\n",file);
115
116/////////////////////////////////////// gets information on geometry ///////////////////////////////////
117 AliITSgeoinfo *geoinfo = new AliITSgeoinfo;
118
119 AliITSgeom *g1 = ((AliITS*)gAlice->GetDetector("ITS"))->GetITSgeom();
120
121 Int_t ll=1, dd=1;
122 TVector det(9);
123
124 //cout<<" nlad ed ndet \n";
583731c7 125 for(i=0; i<6; i++) {
13888096 126 geoinfo->Nlad[i]=g1->GetNladders(i+1);
127 geoinfo->Ndet[i]=g1->GetNdetectors(i+1);
128 //cout<<geoinfo->Nlad[i]<<" "<<geoinfo->Ndet[i]<<"\n";
129 }
130 //getchar();
131
132 //cout<<" raggio medio = ";
583731c7 133 for(i=0; i<6; i++) {
13888096 134 g1->GetCenterThetaPhi(i+1,ll,dd,det);
135 geoinfo->Avrad[i]=TMath::Sqrt(det(0)*det(0)+det(1)*det(1));
136 //cout<<geoinfo->Avrad[i]<<" ";
137 }
138 //cout<<"\n"; getchar();
139
140 geoinfo->Detx[0] = ((AliITSgeomSPD*)(g1->GetShape(1, ll, dd)))->GetDx();
141 geoinfo->Detz[0] = ((AliITSgeomSPD*)(g1->GetShape(1, ll, dd)))->GetDz();
142
143 geoinfo->Detx[1] = ((AliITSgeomSPD*)(g1->GetShape(2, ll, dd)))->GetDx();
144 geoinfo->Detz[1] = ((AliITSgeomSPD*)(g1->GetShape(2, ll, dd)))->GetDz();
145
146 geoinfo->Detx[2] = ((AliITSgeomSDD*)(g1->GetShape(3, ll, dd)))->GetDx();
147 geoinfo->Detz[2] = ((AliITSgeomSDD*)(g1->GetShape(3, ll, dd)))->GetDz();
148
149 geoinfo->Detx[3] = ((AliITSgeomSDD*)(g1->GetShape(4, ll, dd)))->GetDx();
150 geoinfo->Detz[3] = ((AliITSgeomSDD*)(g1->GetShape(4, ll, dd)))->GetDz();
151
152 geoinfo->Detx[4] = ((AliITSgeomSSD*)(g1->GetShape(5, ll, dd)))->GetDx();
153 geoinfo->Detz[4] = ((AliITSgeomSSD*)(g1->GetShape(5, ll, dd)))->GetDz();
154
155 geoinfo->Detx[5] = ((AliITSgeomSSD*)(g1->GetShape(6, ll, dd)))->GetDx();
156 geoinfo->Detz[5] = ((AliITSgeomSSD*)(g1->GetShape(6, ll, dd)))->GetDz();
157
158 //cout<<" Detx Detz\n";
159 //for(Int_t la=0; la<6; la++) cout<<" "<<geoinfo->Detx[la]<<" "<<geoinfo->Detz[la]<<"\n";
160 //getchar();
161//////////////////////////////////////////////////////////////////////////////////////////////////////////
162
163 //const char *pname="75x40_100x60";
164
165 Int_t imax=200,jmax=450;
166 AliITSRad *rl = new AliITSRad(imax,jmax);
167 //cout<<" dopo costruttore AliITSRad\n"; getchar();
168
169 struct GoodTrack {
170 Int_t lab,code;
171 Float_t px,py,pz,x,y,z,pxg,pyg,pzg,ptg;
172 Bool_t flag;
173 };
174
175
176 gAlice->GetEvent(0);
177
178 AliKalmanTrack *kkprov;
179 kkprov->SetConvConst(100/0.299792458/0.2/gAlice->Field()->Factor());
180
181 TFile *cf=TFile::Open("AliTPCclusters.root");
182 AliTPCParam *digp= (AliTPCParam*)cf->Get("75x40_100x60");
183 if (!digp) { cerr<<"TPC parameters have not been found !\n"; getchar();}
184
185 AliTPCtracker *tracker = new AliTPCtracker(digp);
186
187 // Load clusters
188 tracker->LoadInnerSectors();
189 tracker->LoadOuterSectors();
190
191
192 GoodTrack gt[15000];
193 Int_t ngood=0;
194 ifstream in("itsgood_tracks");
195
196 cerr<<"Reading itsgood tracks...\n";
197 while (in>>gt[ngood].lab>>gt[ngood].code
198 >>gt[ngood].px >>gt[ngood].py>>gt[ngood].pz
199 >>gt[ngood].x >>gt[ngood].y >>gt[ngood].z
200 >>gt[ngood].pxg >>gt[ngood].pyg >>gt[ngood].pzg
201 >>gt[ngood].ptg >>gt[ngood].flag) {
202 ngood++;
203 cerr<<ngood<<'\r';
204 if (ngood==15000) {
205 cerr<<"Too many good tracks !\n";
206 break;
207 }
208 }
209 if (!in.eof()) cerr<<"Read error (itsgood_tracks) !\n";
210
211
212// Load tracks
213 TFile *tf=TFile::Open("AliTPCtracks.root");
214 if (!tf->IsOpen()) {cerr<<"Can't open AliTPCtracks.root !\n"; return ;}
215 TObjArray tracks(200000);
216 TTree *tracktree=(TTree*)tf->Get("TPCf");
217 if (!tracktree) {cerr<<"Can't get a tree with TPC tracks !\n";}
218 TBranch *tbranch=tracktree->GetBranch("tracks");
219 Int_t nentr=(Int_t)tracktree->GetEntries();
220 Int_t kk;
221 AliTPCtrack *iotracktpc=0;
222 for (kk=0; kk<nentr; kk++) {
223 iotracktpc=new AliTPCtrack;
224 tbranch->SetAddress(&iotracktpc);
225 tracktree->GetEvent(kk);
226 tracker->CookLabel(iotracktpc,0.1);
227 tracks.AddLast(iotracktpc);
228 }
229 delete tracker;
230 tf->Close();
231
232
233 Int_t nt = tracks.GetEntriesFast();
234 cerr<<"Number of found tracks "<<nt<<endl;
235
236 TVector DataOut(9);
237 Int_t kkk=0;
238
239 Double_t ptg=0.,pxg=0.,pyg=0.,pzg=0.;
240
241 ////////////////////////////// good tracks definition in TPC ////////////////////////////////
242
243 ofstream out1 ("AliITSTrag.out");
13888096 244 for (i=0; i<ngood; i++) out1 << gt[i].ptg << "\n";
245 out1.close();
246
247
248 TVector vec(5);
249 TTree *TR=gAlice->TreeR();
250 Int_t nent=(Int_t)TR->GetEntries();
251 TClonesArray *recPoints = ITS->RecPoints();
252
253 Int_t numbpoints;
254 Int_t totalpoints=0;
255 Int_t *np = new Int_t[nent];
256 Int_t **vettid = new Int_t* [nent];
257 Int_t mod;
258
259 for (mod=0; mod<nent; mod++) {
260 vettid[mod]=0;
261 ITS->ResetRecPoints();
262 //gAlice->TreeR()->GetEvent(mod+1); //first entry in TreeR is empty
263 gAlice->TreeR()->GetEvent(mod); //first entry in TreeR is empty
264 numbpoints = recPoints->GetEntries();
265 totalpoints+=numbpoints;
266 np[mod] = numbpoints;
267 //cout<<" mod = "<<mod<<" numbpoints = "<<numbpoints<<"\n"; getchar();
268 vettid[mod] = new Int_t[numbpoints];
269 Int_t ii;
270 for (ii=0;ii<numbpoints; ii++) *(vettid[mod]+ii)=0;
271 }
272
273 AliTPCtrack *track=0;
274
275
276 if(min_t < 0) {min_t = 0; max_t = nt-1;}
277
278/*
279 ///////////////////////////////// Definition of vertex end its error ////////////////////////////
280 ////////////////////////// In the future it will be given by a method ///////////////////////////
281 Double_t Vx=0.;
282 Double_t Vy=0.;
283 Double_t Vz=0.;
284
285 Float_t sigmavx=0.0050; // 50 microns
286 Float_t sigmavy=0.0050; // 50 microns
287 Float_t sigmavz=0.010; // 100 microns
288
289 //Vx+=gRandom->Gaus(0,sigmavx); Vy+=gRandom->Gaus(0,sigmavy); Vz+=gRandom->Gaus(0,sigmavz);
290 TVector vertex(3), ervertex(3)
291 vertex(0)=Vx; vertex(1)=Vy; vertex(2)=Vz;
292 ervertex(0)=sigmavx; ervertex(1)=sigmavy; ervertex(2)=sigmavz;
293 /////////////////////////////////////////////////////////////////////////////////////////////////
294*/
295
296
297 TTree tracktree1("TreeT","Tree with ITS tracks");
298 AliITSiotrack *iotrack=0;
299 tracktree1.Branch("ITStracks","AliITSiotrack",&iotrack,32000,0);
300
301 ofstream out ("AliITSTra.out");
302
303 Int_t j;
304 for (j=min_t; j<=max_t; j++) {
305 track=(AliTPCtrack*)tracks.UncheckedAt(j);
306 Int_t flaglab=0;
307 if (!track) continue;
308 ////// elimination of not good tracks ////////////
309 Int_t ilab=TMath::Abs(track->GetLabel());
310 Int_t iii;
311 for (iii=0;iii<ngood;iii++) {
312 //cout<<" ilab, gt[iii].lab = "<<ilab<<" "<<gt[iii].lab<<"\n"; getchar();
313 if (ilab==gt[iii].lab) {
314 flaglab=1;
315 ptg=gt[iii].ptg;
316 pxg=gt[iii].pxg;
317 pyg=gt[iii].pyg;
318 pzg=gt[iii].pzg;
319 break;
320 }
321 }
322 //cout<<" j flaglab = " <<j<<" "<<flaglab<<"\n"; getchar();
323 if (!flaglab) continue;
324 //cout<<" j = " <<j<<"\n"; getchar();
325
326
327 ////// new propagation to the end of TPC //////////////
328 Double_t xk=77.415;
329 track->PropagateTo(xk, 28.94, 1.204e-3); //Ne
330 xk -=0.01;
331 track->PropagateTo(xk, 44.77, 1.71); //Tedlar
332 xk -=0.04;
333 track->PropagateTo(xk, 44.86, 1.45); //Kevlar
334 xk -=2.0;
335 track->PropagateTo(xk, 41.28, 0.029); //Nomex
336 xk-=16;
337 track->PropagateTo(xk,36.2,1.98e-3); //C02
338 xk -=0.01;
339 track->PropagateTo(xk, 24.01, 2.7); //Al
340 xk -=0.01;
341 track->PropagateTo(xk, 44.77, 1.71); //Tedlar
342 xk -=0.04;
343 track->PropagateTo(xk, 44.86, 1.45); //Kevlar
344 xk -=0.5;
345 track->PropagateTo(xk, 41.28, 0.029); //Nomex
346
347 ///////////////////////////////////////////////////////////////
348
349 ///////////////////////////////////////////////////////////////
350 AliITStrack trackITS(*track);
351 AliITStrack result(*track);
352 AliITStrack primarytrack(*track);
353
354///////////////////////////////////////////////////////////////////////////////////////////////
355 TVector Vgeant(3);
356 Vgeant=result.GetVertex();
357
358 // Definition of Dv and Zv for vertex constraint
359 Double_t sigmaDv=0.0050; Double_t sigmaZv=0.010;
360 //Double_t sigmaDv=0.0015; Double_t sigmaZv=0.0015;
361 Double_t uniform= gRandom->Uniform();
362 Double_t signdv;
363 if(uniform<=0.5) signdv=-1.;
364 else
365 signdv=1.;
366
367 Double_t Vr=TMath::Sqrt(Vgeant(0)*Vgeant(0)+ Vgeant(1)*Vgeant(1));
368 Double_t Dv=gRandom->Gaus(signdv*Vr,(Float_t)sigmaDv);
369 Double_t Zv=gRandom->Gaus(Vgeant(2),(Float_t)sigmaZv);
370
371 //cout<<" Dv e Zv = "<<Dv<<" "<<Zv<<"\n";
372 trackITS.SetDv(Dv); trackITS.SetZv(Zv);
373 trackITS.SetsigmaDv(sigmaDv); trackITS.SetsigmaZv(sigmaZv);
374 result.SetDv(Dv); result.SetZv(Zv);
375 result.SetsigmaDv(sigmaDv); result.SetsigmaZv(sigmaZv);
376 primarytrack.SetDv(Dv); primarytrack.SetZv(Zv);
377 primarytrack.SetsigmaDv(sigmaDv); primarytrack.SetsigmaZv(sigmaZv);
378
379/////////////////////////////////////////////////////////////////////////////////////////////////
380
381 primarytrack.PrimaryTrack(rl);
382 TVector d2=primarytrack.Getd2();
383 TVector tgl2=primarytrack.Gettgl2();
384 TVector dtgl=primarytrack.Getdtgl();
385 trackITS.Setd2(d2); trackITS.Settgl2(tgl2); trackITS.Setdtgl(dtgl);
386 result.Setd2(d2); result.Settgl2(tgl2); result.Setdtgl(dtgl);
387 /*
388 trackITS.SetVertex(vertex); trackITS.SetErrorVertex(ervertex);
389 result.SetVertex(vertex); result.SetErrorVertex(ervertex);
390 */
391
392 Tracking(trackITS,&result,recPoints,vettid, flagvert,rl,geoinfo);
393
394 // cout<<" progressive track number = "<<j<<"\r";
395 // cout<<j<<"\r";
396 Int_t NumofCluster=result.GetNumClust();
397 // cout<<" progressive track number = "<<j<<"\n"; // stampa
398 Long_t labITS=result.GetLabel();
399 // cout << " ITS track label = " << labITS << "\n"; // stampa
400 int lab=track->GetLabel();
401 //cout << " TPC track label = " << lab <<"\n"; // stampa
402
403
404//propagation to vertex
405
406 Double_t rbeam=3.;
407
408 result.Propagation(rbeam);
409
410 Double_t C00,C10,C11,C20,C21,C22,C30,C31,C32,C33,C40,C41,C42,C43,C44;
411 result.GetCElements(C00,C10,C11,C20,C21,C22,C30,C31,C32,C33,C40,C41,C42,C43,C44);
412
413 Double_t pt=TMath::Abs(result.GetPt());
414 Double_t Dr=result.GetD();
415 Double_t Z=result.GetZ();
416 Double_t tgl=result.GetTgl();
417 Double_t C=result.GetC();
418 Double_t Cy=C/2.;
419 Double_t Dz=Z-(tgl/Cy)*TMath::ASin(result.arga(rbeam));
420 Dz-=Vgeant(2);
421
422 // cout<<" Dr e dz alla fine = "<<Dr<<" "<<Dz<<"\n"; getchar();
423 Double_t phi=result.Getphi();
424 Double_t phivertex = phi - TMath::ASin(result.argA(rbeam));
425 Double_t duepi=2.*TMath::Pi();
426 if(phivertex>duepi) phivertex-=duepi;
427 if(phivertex<0.) phivertex+=duepi;
428 Double_t Dtot=TMath::Sqrt(Dr*Dr+Dz*Dz);
429
430//////////////////////////////////////////////////////////////////////////////////////////
431
432 Int_t idmodule,idpoint;
433 if(NumofCluster >=5) { // cinque - sei
434 //if(NumofCluster ==6) { // cinque - sei
435
436
437 AliITSiotrack outtrack;
438
439 iotrack=&outtrack;
440
441 iotrack->SetStatePhi(phi);
442 iotrack->SetStateZ(Z);
443 iotrack->SetStateD(Dr);
444 iotrack->SetStateTgl(tgl);
445 iotrack->SetStateC(C);
446 Double_t radius=result.Getrtrack();
447 iotrack->SetRadius(radius);
448 Int_t charge;
449 if(C>0.) charge=-1; else charge=1;
450 iotrack->SetCharge(charge);
451
452
453
454 iotrack->SetCovMatrix(C00,C10,C11,C20,C21,C22,C30,C31,C32,C33,C40,C41,C42,C43,C44);
455
456 Double_t px=pt*TMath::Cos(phivertex);
457 Double_t py=pt*TMath::Sin(phivertex);
458 Double_t pz=pt*tgl;
459
460 Double_t xtrack=Dr*TMath::Sin(phivertex);
461 Double_t ytrack=Dr*TMath::Cos(phivertex);
462 Double_t ztrack=Dz+Vgeant(2);
463
464
465 iotrack->SetPx(px);
466 iotrack->SetPy(py);
467 iotrack->SetPz(pz);
468 iotrack->SetX(xtrack);
469 iotrack->SetY(ytrack);
470 iotrack->SetZ(ztrack);
471 iotrack->SetLabel(labITS);
472
473 Int_t il;
474 for(il=0;il<6; il++){
475 iotrack->SetIdPoint(il,result.GetIdPoint(il));
476 iotrack->SetIdModule(il,result.GetIdModule(il));
477 }
478 tracktree1.Fill();
479
480 //cout<<" labITS = "<<labITS<<"\n";
481 //cout<<" phi z Dr tgl C = "<<phi<<" "<<Z<<" "<<Dr<<" "<<tgl<<" "<<C<<"\n"; getchar();
482
483 DataOut(kkk) = ptg; kkk++; DataOut(kkk)=labITS; kkk++; DataOut(kkk)=lab; kkk++;
484
485 for (il=0;il<6;il++) {
486 idpoint=result.GetIdPoint(il);
487 idmodule=result.GetIdModule(il);
488 *(vettid[idmodule]+idpoint)=1;
489 iotrack->SetIdPoint(il,idpoint);
490 iotrack->SetIdModule(il,idmodule);
491 }
492
493 Double_t difpt= (pt-ptg)/ptg*100.;
494 DataOut(kkk)=difpt; kkk++;
495 Double_t lambdag=TMath::ATan(pzg/ptg);
496 Double_t lam=TMath::ATan(tgl);
497 Double_t diflam = (lam - lambdag)*1000.;
498 DataOut(kkk) = diflam; kkk++;
499 Double_t phig=TMath::ATan2(pyg,pxg); if(phig<0) phig=2.*TMath::Pi()+phig;
500 Double_t phi=phivertex;
501
502 Double_t difphi = (phi - phig)*1000.;
503 DataOut(kkk)=difphi; kkk++;
504 DataOut(kkk)=Dtot*1.e4; kkk++;
505 DataOut(kkk)=Dr*1.e4; kkk++;
506 DataOut(kkk)=Dz*1.e4; kkk++;
507 Int_t r;
508 for (r=0; r<9; r++) { out<<DataOut(r)<<" ";}
509 out<<"\n";
510 kkk=0;
511
512
513 } // end if on NumofCluster
514 //gObjectTable->Print(); // stampa memoria
515 } // end for (int j=min_t; j<=max_t; j++)
516
517 out.close();
518
519
520 static Bool_t first=kTRUE;
521 static TFile *tfile;
522
523 if(first) {
524 tfile=new TFile("itstracks.root","RECREATE");
525 //cout<<"I have opened itstracks.root file "<<endl;
526 }
527 first=kFALSE;
528 tfile->cd();
529 tfile->ls();
530
531 char hname[30];
532 sprintf(hname,"TreeT%d",evNumber);
533
534 tracktree1.Write(hname);
535
536
537
538 TTree *fAli=gAlice->TreeK();
539 TFile *fileAli=0;
540
541 if (fAli) fileAli =fAli->GetCurrentFile();
542 fileAli->cd();
543
544 ////////////////////////////////////////////////////////////////////////////////////////////////
545
546 printf("delete vectors\n");
547 if(np) delete [] np;
548 if(vettid) delete [] vettid;
549
550 if(geoinfo) delete geoinfo;
551
552}