[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

  //Loop variables

  Int_t i;

  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(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(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");
  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
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(DrDr+DzDz);
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	}