[u/mrichter/AliRoot.git] / ITS / AliITSTrackerV1.cxx

#include <iostream.h>
#include <TROOT.h>
#include <TMath.h>
#include <TRandom.h>
#include <TBranch.h>
#include <TVector.h>
#include <TClonesArray.h>
#include <TFile.h>
#include <TTree.h>


#include "TParticle.h"
#include "AliRun.h"
#include "AliITS.h"
#include "AliITSgeomSPD.h"
#include "AliITSgeomSDD.h"
#include "AliITSgeomSSD.h"
#include "AliITSgeom.h"
#include "AliITSmodule.h"
#include "AliITSRecPoint.h"
#include "AliMC.h"
#include "AliKalmanTrack.h" 
#include "AliMagF.h"


#include "AliITSRad.h" 
#include "AliITStrack.h"
#include "AliITSiotrack.h"
#include "AliITStracking.h"   
#include "../TPC/AliTPC.h"
#include "../TPC/AliTPCParam.h"
#include "../TPC/AliTPCtracker.h"

#include "AliITSgeoinfo.h"
#include "AliITSTrackerV1.h"


ClassImp(AliITSTrackerV1)


//________________________________________________________________

AliITSTrackerV1::AliITSTrackerV1(AliITS* IITTSS) {
//Origin  A. Badala' and G.S. Pappalardo:  e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it 
// default constructor   
  ITS = IITTSS;

}

//________________________________________________________________

AliITStrack  AliITSTrackerV1::Tracking(AliITStrack &track, AliITStrack *reference,TObjArray *fastpoints, 
Int_t **vettid, Bool_t flagvert,  AliITSRad *rl, AliITSgeoinfo *geoinfo) { 
										
//Origin  A. Badala' and G.S. Pappalardo:  e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it 
  										    
  TList *list= new TList();   

  AliITStrack tr(track);
  
  list->AddLast(&tr);
  
  Double_t Pt=(tr).GetPt();
  //cout << "\n Pt = " << Pt <<"\n";  //stampa

  AliITStracking obj(list, reference, ITS, fastpoints,TMath::Abs(Pt),vettid, flagvert, rl, geoinfo);  // nuova ITS 
  list->Delete();
  delete list;

  Int_t itot=-1;
  TVector VecTotLabref(18);
  Int_t lay, k;
  for(lay=5; lay>=0; lay--) {
    TVector VecLabref(3); 
    VecLabref=(*reference).GetLabTrack(lay);
    Float_t ClustZ=(*reference).GetZclusterTrack( lay);   
    for(k=0; k<3; k++){ 
		Int_t lpp=(Int_t)VecLabref(k);
		if(lpp>=0) {
		  TParticle *p=(TParticle*) gAlice->Particle(lpp);
		  Int_t pcode=p->GetPdgCode();
		  if(pcode==11) VecLabref(k)=p->GetFirstMother();
		}    
    itot++; VecTotLabref(itot)=VecLabref(k);
    if(VecLabref(k)==0. && ClustZ == 0.) VecTotLabref(itot) =-3.; }  
  }
  Long_t labref;
  Int_t freq;  
  (*reference).Search(VecTotLabref, labref, freq);
    
  //if(freq < 6) labref=-labref;        // cinque - sei
  if(freq < 5) labref=-labref;        // cinque - sei	
  (*reference).SetLabel(labref);

  return *reference; 

}


//________________________________________________________________


void AliITSTrackerV1::DoTracking(Int_t evNumber, Int_t min_t, Int_t max_t, TFile *file, Bool_t flagvert) {
//Origin  A. Badala' and G.S. Pappalardo:  e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it 
//   ex macro for tracking ITS

  printf("begin DoTracking - file %p\n",file);
  
///////////////////////////////////////  gets information on geometry ///////////////////////////////////  
  AliITSgeoinfo *geoinfo = new AliITSgeoinfo;

  AliITSgeom *g1 = ((AliITS*)gAlice->GetDetector("ITS"))->GetITSgeom(); 
  
  Int_t ll=1, dd=1;
  TVector det(9);
  
  //cout<<" nlad ed ndet \n";
  for(Int_t i=0; i<6; i++) {
    geoinfo->Nlad[i]=g1->GetNladders(i+1);
    geoinfo->Ndet[i]=g1->GetNdetectors(i+1);
	 //cout<<geoinfo->Nlad[i]<<" "<<geoinfo->Ndet[i]<<"\n"; 
  }
  //getchar();

  //cout<<" raggio medio = ";
  for(Int_t i=0; i<6; i++) {  
    g1->GetCenterThetaPhi(i+1,ll,dd,det);
    geoinfo->Avrad[i]=TMath::Sqrt(det(0)*det(0)+det(1)*det(1));
	 //cout<<geoinfo->Avrad[i]<<" ";
  }
  //cout<<"\n"; getchar();
  
  geoinfo->Detx[0] = ((AliITSgeomSPD*)(g1->GetShape(1, ll, dd)))->GetDx();
  geoinfo->Detz[0] = ((AliITSgeomSPD*)(g1->GetShape(1, ll, dd)))->GetDz();
  
  geoinfo->Detx[1] = ((AliITSgeomSPD*)(g1->GetShape(2, ll, dd)))->GetDx();
  geoinfo->Detz[1] = ((AliITSgeomSPD*)(g1->GetShape(2, ll, dd)))->GetDz();
  
  geoinfo->Detx[2] = ((AliITSgeomSDD*)(g1->GetShape(3, ll, dd)))->GetDx();
  geoinfo->Detz[2] = ((AliITSgeomSDD*)(g1->GetShape(3, ll, dd)))->GetDz();
  
  geoinfo->Detx[3] = ((AliITSgeomSDD*)(g1->GetShape(4, ll, dd)))->GetDx();
  geoinfo->Detz[3] = ((AliITSgeomSDD*)(g1->GetShape(4, ll, dd)))->GetDz();
  
  geoinfo->Detx[4] = ((AliITSgeomSSD*)(g1->GetShape(5, ll, dd)))->GetDx();
  geoinfo->Detz[4] = ((AliITSgeomSSD*)(g1->GetShape(5, ll, dd)))->GetDz();
  
  geoinfo->Detx[5] = ((AliITSgeomSSD*)(g1->GetShape(6, ll, dd)))->GetDx();
  geoinfo->Detz[5] = ((AliITSgeomSSD*)(g1->GetShape(6, ll, dd)))->GetDz();
  
  //cout<<"    Detx     Detz\n";
  //for(Int_t la=0; la<6; la++) cout<<"    "<<geoinfo->Detx[la]<<"     "<<geoinfo->Detz[la]<<"\n";
  //getchar();
//////////////////////////////////////////////////////////////////////////////////////////////////////////  
  
  //const char *pname="75x40_100x60";
  
 Int_t imax=200,jmax=450;
  AliITSRad *rl = new AliITSRad(imax,jmax);
  //cout<<" dopo costruttore AliITSRad\n"; getchar();
    
  struct GoodTrack {
    Int_t lab,code;
    Float_t px,py,pz,x,y,z,pxg,pyg,pzg,ptg;
    Bool_t flag;
  };
  

  gAlice->GetEvent(0);
 
    AliKalmanTrack *kkprov;
    kkprov->SetConvConst(100/0.299792458/0.2/gAlice->Field()->Factor());  

  TFile *cf=TFile::Open("AliTPCclusters.root");  
  AliTPCParam *digp= (AliTPCParam*)cf->Get("75x40_100x60");
  if (!digp) { cerr<<"TPC parameters have not been found !\n"; getchar();}
  
   AliTPCtracker *tracker = new AliTPCtracker(digp);  
 
 // Load clusters
   tracker->LoadInnerSectors();
   tracker->LoadOuterSectors();
       
  
  GoodTrack gt[15000];
  Int_t ngood=0;
  ifstream in("itsgood_tracks");

  cerr<<"Reading itsgood 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==15000) {
      cerr<<"Too many good tracks !\n";
      break;
    }
  }
  if (!in.eof()) cerr<<"Read error (itsgood_tracks) !\n";
  
  
// Load tracks
  TFile *tf=TFile::Open("AliTPCtracks.root"); 
  if (!tf->IsOpen()) {cerr<<"Can't open AliTPCtracks.root !\n"; return ;}
  TObjArray tracks(200000);
   TTree *tracktree=(TTree*)tf->Get("TPCf"); 
   if (!tracktree) {cerr<<"Can't get a tree with TPC tracks !\n";}   
  TBranch *tbranch=tracktree->GetBranch("tracks");
  Int_t nentr=(Int_t)tracktree->GetEntries();
  Int_t kk;
   AliTPCtrack *iotracktpc=0;    
  for (kk=0; kk<nentr; kk++) {
    iotracktpc=new AliTPCtrack; 
    tbranch->SetAddress(&iotracktpc);
    tracktree->GetEvent(kk);    
    tracker->CookLabel(iotracktpc,0.1);       
    tracks.AddLast(iotracktpc);         
  }  
   delete tracker;      
  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 = ITS->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;
    ITS->ResetRecPoints();  
    //gAlice->TreeR()->GetEvent(mod+1); //first entry in TreeR is empty
    gAlice->TreeR()->GetEvent(mod); //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=0;

     
  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;
  /////////////////////////////////////////////////////////////////////////////////////////////////
*/      
 

  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();
 
 	 	 
	 ////// 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(rl);
    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,rl,geoinfo);  
	     
    // cout<<" progressive track number = "<<j<<"\r";
   // cout<<j<<"\r";
    Int_t NumofCluster=result.GetNumClust();  
   // cout<<" progressive track number = "<<j<<"\n";    // stampa
    Long_t labITS=result.GetLabel();
   // cout << " ITS track label = " << labITS << "\n"; 	// stampa	    
    int lab=track->GetLabel();		    
    //cout << " TPC track label = " << lab <<"\n";      // stampa
	 
	     
//propagation to vertex
	
    Double_t rbeam=3.;
     
    result.Propagation(rbeam);
       	 
	 Double_t C00,C10,C11,C20,C21,C22,C30,C31,C32,C33,C40,C41,C42,C43,C44;
	 result.GetCElements(C00,C10,C11,C20,C21,C22,C30,C31,C32,C33,C40,C41,C42,C43,C44);
	 	 
    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 idmodule,idpoint;
	 if(NumofCluster >=5)  {            // cinque - sei
	 //if(NumofCluster ==6)  {            // cinque - sei 


      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(C00,C10,C11,C20,C21,C22,C30,C31,C32,C33,C40,C41,C42,C43,C44);  

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