-
-
-//________________________________________________________________
-
-
-
-void AliITS::DoTracking(Int_t evNumber, Int_t min_t, Int_t max_t, TFile *file, Bool_t flagvert) {
-
-// ex macro for tracking ITS
-
- printf("begin DoTracking - file %p\n",file);
-
- const char *pname="75x40_100x60";
-
- struct GoodTrack {
- Int_t lab,code;
- Float_t px,py,pz,x,y,z,pxg,pyg,pzg,ptg;
- Bool_t flag;
- };
-
-
- gAlice->GetEvent(0);
-
- AliTPC *TPC=(AliTPC*)gAlice->GetDetector("TPC");
- AliTPCParam *digp = (AliTPCParam*)file->Get(pname);
- if (digp!=0) TPC->SetParam(digp);
-
- GoodTrack gt[7000];
- Int_t ngood=0;
- ifstream in("good_tracks");
-
- cerr<<"Reading good tracks...\n";
- while (in>>gt[ngood].lab>>gt[ngood].code
- >>gt[ngood].px >>gt[ngood].py>>gt[ngood].pz
- >>gt[ngood].x >>gt[ngood].y >>gt[ngood].z
- >>gt[ngood].pxg >>gt[ngood].pyg >>gt[ngood].pzg
- >>gt[ngood].ptg >>gt[ngood].flag) {
- ngood++;
- cerr<<ngood<<'\r';
- if (ngood==7000) {
- cerr<<"Too many good tracks !\n";
- break;
- }
- }
- if (!in.eof()) cerr<<"Read error (good_tracks) !\n";
-
-
-// Load tracks
- TFile *tf=TFile::Open("tpctracks.root");
- if (!tf->IsOpen()) {cerr<<"Can't open tpctracks.root !\n"; return ;}
- TObjArray tracks(200000);
- TTree *tracktree=(TTree*)tf->Get("TreeT");
- TBranch *tbranch=tracktree->GetBranch("tracks");
- Int_t nentr=(Int_t)tracktree->GetEntries();
- Int_t kk;
- for (kk=0; kk<nentr; kk++) {
- AliTPCtrack *iotrack=new AliTPCtrack;
- tbranch->SetAddress(&iotrack);
- tracktree->GetEvent(kk);
- tracks.AddLast(iotrack);
- }
- tf->Close();
-
-
- Int_t nt = tracks.GetEntriesFast();
- cerr<<"Number of found tracks "<<nt<<endl;
-
- TVector DataOut(9);
- Int_t kkk=0;
-
- Double_t ptg=0.,pxg=0.,pyg=0.,pzg=0.;
-
- ////////////////////////////// good tracks definition in TPC ////////////////////////////////
-
- ofstream out1 ("AliITSTrag.out");
- Int_t i;
- for (i=0; i<ngood; i++) out1 << gt[i].ptg << "\n";
- out1.close();
-
-
- TVector vec(5);
- TTree *TR=gAlice->TreeR();
- Int_t nent=(Int_t)TR->GetEntries();
- TClonesArray *recPoints = RecPoints();
- Int_t numbpoints;
- Int_t totalpoints=0;
- Int_t *np = new Int_t[nent];
- Int_t **vettid = new Int_t* [nent];
- Int_t mod;
- for (mod=0; mod<nent; mod++) {
- vettid[mod]=0;
- this->ResetRecPoints();
- gAlice->TreeR()->GetEvent(mod+1); //first entry in TreeR is empty
- numbpoints = recPoints->GetEntries();
- totalpoints+=numbpoints;
- np[mod] = numbpoints;
- //cout<<" mod = "<<mod<<" numbpoints = "<<numbpoints<<"\n"; getchar();
- vettid[mod] = new Int_t[numbpoints];
- Int_t ii;
- for (ii=0;ii<numbpoints; ii++) *(vettid[mod]+ii)=0;
- }
-
- AliTPCtrack *track;
-
-
- if(min_t < 0) {min_t = 0; max_t = nt-1;}
-
-/*
- ///////////////////////////////// Definition of vertex end its error ////////////////////////////
- ////////////////////////// In the future it will be given by a method ///////////////////////////
- Double_t Vx=0.;
- Double_t Vy=0.;
- Double_t Vz=0.;
-
- Float_t sigmavx=0.0050; // 50 microns
- Float_t sigmavy=0.0050; // 50 microns
- Float_t sigmavz=0.010; // 100 microns
-
- //Vx+=gRandom->Gaus(0,sigmavx); Vy+=gRandom->Gaus(0,sigmavy); Vz+=gRandom->Gaus(0,sigmavz);
- TVector vertex(3), ervertex(3)
- vertex(0)=Vx; vertex(1)=Vy; vertex(2)=Vz;
- ervertex(0)=sigmavx; ervertex(1)=sigmavy; ervertex(2)=sigmavz;
- /////////////////////////////////////////////////////////////////////////////////////////////////
-*/
-
- //TDirectory *savedir=gDirectory;
-
- TTree tracktree1("TreeT","Tree with ITS tracks");
- AliITSiotrack *iotrack=0;
- tracktree1.Branch("ITStracks","AliITSiotrack",&iotrack,32000,0);
-
- ofstream out ("AliITSTra.out");
-
- Int_t j;
- for (j=min_t; j<=max_t; j++) {
- track=(AliTPCtrack*)tracks.UncheckedAt(j);
- Int_t flaglab=0;
- if (!track) continue;
- ////// elimination of not good tracks ////////////
- Int_t ilab=TMath::Abs(track->GetLabel());
- Int_t iii;
- for (iii=0;iii<ngood;iii++) {
- //cout<<" ilab, gt[iii].lab = "<<ilab<<" "<<gt[iii].lab<<"\n"; getchar();
- if (ilab==gt[iii].lab) {
- flaglab=1;
- ptg=gt[iii].ptg;
- pxg=gt[iii].pxg;
- pyg=gt[iii].pyg;
- pzg=gt[iii].pzg;
- break;
- }
- }
- //cout<<" j flaglab = " <<j<<" "<<flaglab<<"\n"; getchar();
- if (!flaglab) continue;
- //cout<<" j = " <<j<<"\n"; getchar();
- /*
- ////// old propagation to the end of TPC //////////////
- Double_t xk=76.;
- track->PropagateTo(xk);
- xk-=0.11;
- track->PropagateTo(xk,42.7,2.27); //C
- xk-=2.6;
- track->PropagateTo(xk,36.2,1.98e-3); //C02
- xk-=0.051;
- track->PropagateTo(xk,42.7,2.27); //C
- ///////////////////////////////////////////////////
- */
-
- ////// new propagation to the end of TPC //////////////
- Double_t xk=77.415;
- track->PropagateTo(xk, 28.94, 1.204e-3); //Ne
- xk -=0.01;
- track->PropagateTo(xk, 44.77, 1.71); //Tedlar
- xk -=0.04;
- track->PropagateTo(xk, 44.86, 1.45); //Kevlar
- xk -=2.0;
- track->PropagateTo(xk, 41.28, 0.029); //Nomex
- xk-=16;
- track->PropagateTo(xk,36.2,1.98e-3); //C02
- xk -=0.01;
- track->PropagateTo(xk, 24.01, 2.7); //Al
- xk -=0.01;
- track->PropagateTo(xk, 44.77, 1.71); //Tedlar
- xk -=0.04;
- track->PropagateTo(xk, 44.86, 1.45); //Kevlar
- xk -=0.5;
- track->PropagateTo(xk, 41.28, 0.029); //Nomex
-
- ///////////////////////////////////////////////////////////////
-
- ///////////////////////////////////////////////////////////////
- AliITStrack trackITS(*track);
- AliITStrack result(*track);
- AliITStrack primarytrack(*track);
-
-///////////////////////////////////////////////////////////////////////////////////////////////
- TVector Vgeant(3);
- Vgeant=result.GetVertex();
-
- // Definition of Dv and Zv for vertex constraint
- Double_t sigmaDv=0.0050; Double_t sigmaZv=0.010;
- //Double_t sigmaDv=0.0015; Double_t sigmaZv=0.0015;
- Double_t uniform= gRandom->Uniform();
- Double_t signdv;
- if(uniform<=0.5) signdv=-1.;
- else
- signdv=1.;
-
- Double_t Vr=TMath::Sqrt(Vgeant(0)*Vgeant(0)+ Vgeant(1)*Vgeant(1));
- Double_t Dv=gRandom->Gaus(signdv*Vr,(Float_t)sigmaDv);
- Double_t Zv=gRandom->Gaus(Vgeant(2),(Float_t)sigmaZv);
-
- //cout<<" Dv e Zv = "<<Dv<<" "<<Zv<<"\n";
- trackITS.SetDv(Dv); trackITS.SetZv(Zv);
- trackITS.SetsigmaDv(sigmaDv); trackITS.SetsigmaZv(sigmaZv);
- result.SetDv(Dv); result.SetZv(Zv);
- result.SetsigmaDv(sigmaDv); result.SetsigmaZv(sigmaZv);
- primarytrack.SetDv(Dv); primarytrack.SetZv(Zv);
- primarytrack.SetsigmaDv(sigmaDv); primarytrack.SetsigmaZv(sigmaZv);
-
-/////////////////////////////////////////////////////////////////////////////////////////////////
-
- primarytrack.PrimaryTrack();
- TVector d2=primarytrack.Getd2();
- TVector tgl2=primarytrack.Gettgl2();
- TVector dtgl=primarytrack.Getdtgl();
- trackITS.Setd2(d2); trackITS.Settgl2(tgl2); trackITS.Setdtgl(dtgl);
- result.Setd2(d2); result.Settgl2(tgl2); result.Setdtgl(dtgl);
- /*
- trackITS.SetVertex(vertex); trackITS.SetErrorVertex(ervertex);
- result.SetVertex(vertex); result.SetErrorVertex(ervertex);
- */
- Tracking(trackITS,&result,recPoints,vettid, flagvert);
-
- cout<<" progressive track number = "<<j<<"\r";
- // cout<<j<<"\r";
- //cout<<" progressive track number = "<<j<<"\n";
- Long_t labITS=result.GetLabel();
- //cout << " ITS track label = " << labITS << "\n";
- int lab=track->GetLabel();
- // cout << " TPC track label = " << lab <<"\n";
- //result.GetClusters(); getchar(); //to print the cluster matrix
-
-//propagation to vertex
-
- Double_t rbeam=3.;
-
- result.Propagation(rbeam);
-
- TMatrix *cov;
- cov=&result.GetCMatrix();
- Double_t pt=TMath::Abs(result.GetPt());
- Double_t Dr=result.GetD();
- Double_t Z=result.GetZ();
- Double_t tgl=result.GetTgl();
- Double_t C=result.GetC();
- Double_t Cy=C/2.;
- Double_t Dz=Z-(tgl/Cy)*TMath::ASin(result.arga(rbeam));
- Dz-=Vgeant(2);
-
- // cout<<" Dr e dz alla fine = "<<Dr<<" "<<Dz<<"\n"; getchar();
- Double_t phi=result.Getphi();
- Double_t phivertex = phi - TMath::ASin(result.argA(rbeam));
- Double_t duepi=2.*TMath::Pi();
- if(phivertex>duepi) phivertex-=duepi;
- if(phivertex<0.) phivertex+=duepi;
- Double_t Dtot=TMath::Sqrt(Dr*Dr+Dz*Dz);
-
-//////////////////////////////////////////////////////////////////////////////////////////
- Int_t NumofCluster, idmodule,idpoint;
- NumofCluster=result.GetNumClust();
- if(NumofCluster >=5) {
-
-
- AliITSiotrack outtrack;
-
- iotrack=&outtrack;
-
- iotrack->SetStatePhi(phi);
- iotrack->SetStateZ(Z);
- iotrack->SetStateD(Dr);
- iotrack->SetStateTgl(tgl);
- iotrack->SetStateC(C);
- Double_t radius=result.Getrtrack();
- iotrack->SetRadius(radius);
- Int_t charge;
- if(C>0.) charge=-1; else charge=1;
- iotrack->SetCharge(charge);
-
-
-
- iotrack->SetCovMatrix(cov);
-
- Double_t px=pt*TMath::Cos(phi);
- Double_t py=pt*TMath::Sin(phi);
- Double_t pz=pt*tgl;
-
- Double_t xtrack=Dr*TMath::Sin(phi);
- Double_t ytrack=Dr*TMath::Cos(phi);
- Double_t ztrack=Dz+Vgeant(2);
-
-
- iotrack->SetPx(px);
- iotrack->SetPy(py);
- iotrack->SetPz(pz);
- iotrack->SetX(xtrack);
- iotrack->SetY(ytrack);
- iotrack->SetZ(ztrack);
- iotrack->SetLabel(labITS);
-
- Int_t il;
- for(il=0;il<6; il++){
- iotrack->SetIdPoint(il,result.GetIdPoint(il));
- iotrack->SetIdModule(il,result.GetIdModule(il));
- }
- tracktree1.Fill();
-
- //cout<<" labITS = "<<labITS<<"\n";
- //cout<<" phi z Dr tgl C = "<<phi<<" "<<Z<<" "<<Dr<<" "<<tgl<<" "<<C<<"\n"; getchar();
-
- DataOut(kkk) = ptg; kkk++; DataOut(kkk)=labITS; kkk++; DataOut(kkk)=lab; kkk++;
-
- for (il=0;il<6;il++) {
- idpoint=result.GetIdPoint(il);
- idmodule=result.GetIdModule(il);
- *(vettid[idmodule]+idpoint)=1;
- iotrack->SetIdPoint(il,idpoint);
- iotrack->SetIdModule(il,idmodule);
- }
-
- // cout<<" +++++++++++++ pt e ptg = "<<pt<<" "<<ptg<<" ++++++++++\n";
- Double_t difpt= (pt-ptg)/ptg*100.;
- DataOut(kkk)=difpt; kkk++;
- Double_t lambdag=TMath::ATan(pzg/ptg);
- Double_t lam=TMath::ATan(tgl);
- Double_t diflam = (lam - lambdag)*1000.;
- DataOut(kkk) = diflam; kkk++;
- Double_t phig=TMath::ATan(pyg/pxg);
- Double_t phi=phivertex;
- Double_t difphi = (phi - phig)*1000.;
- DataOut(kkk)=difphi; kkk++;
- DataOut(kkk)=Dtot*1.e4; kkk++;
- DataOut(kkk)=Dr*1.e4; kkk++;
- DataOut(kkk)=Dz*1.e4; kkk++;
- Int_t r;
- for (r=0; r<9; r++) { out<<DataOut(r)<<" ";}
- out<<"\n";
- kkk=0;
-
-
- } // end if on NumofCluster
- //gObjectTable->Print(); // stampa memoria
- } // end for (int j=min_t; j<=max_t; j++)
-
- out.close();
-
-
- static Bool_t first=kTRUE;
- static TFile *tfile;
-
- if(first) {
- tfile=new TFile("itstracks.root","RECREATE");
- //cout<<"I have opened itstracks.root file "<<endl;
- }
- first=kFALSE;
- tfile->cd();
- tfile->ls();
-
- char hname[30];
- sprintf(hname,"TreeT%d",evNumber);
-
- tracktree1.Write(hname);
-
-
-
- TTree *fAli=gAlice->TreeK();
- TFile *fileAli=0;
-
- if (fAli) fileAli =fAli->GetCurrentFile();
- fileAli->cd();
-
- ////////////////////////////////////////////////////////////////////////////////////////////////
-
- printf("delete vectors\n");
- if(np) delete [] np;
- if(vettid) delete [] vettid;
-