/*
$Log$
+Revision 1.20 2001/11/21 14:47:45 barbera
+Some unuseful print-out commented out
+
+Revision 1.19 2001/11/21 10:49:07 barbera
+Bug correction suggested by Rene done
+
+Revision 1.18 2001/11/20 15:46:17 barbera
+Point coordinated are calculated in cylindrical reference frame once and for all at the beginning of tracking V1
+
Revision 1.10.2.1 2001/10/24 07:26:04 hristov
All the changes from the head are merged with the release
// imposition respectively. The authors thank Mariana Bondila to have help
// them to resolve some problems. July-2000
+#include <iostream.h>
#include <fstream.h>
#include <TMath.h>
#include <TBranch.h>
//Default constructor
fITS = 0;
fresult = 0;
+ fPtref=0.;
+ fChi2max=0.;
+ //fepsphi=0.;
+ //fepsz=0.;
frecPoints = 0;
fvettid = 0;
+ fflagvert=0;
frl = 0;
+ Int_t ia;
+ for(ia=0; ia<6; ia++) {
+ fNlad[ia]=0;
+ fNdet[ia]=0;
+ fAvrad[ia]=0.;
+ fDetx[ia]=0.;
+ fDetz[ia]=0.;
+ } // end for ia
fzmin = 0;
fzmax = 0;
fphimin = 0;
fphimax = 0;
fphidet = 0;
+ fNRecPoints=0;
+ fRecCylR=0;
+ fRecCylPhi=0;
+ fRecCylZ=0;
+ fFieldFactor=0;
}
//______________________________________________________________________
-AliITSTrackerV1::AliITSTrackerV1(AliITS* IITTSS, Bool_t flag) {
+AliITSTrackerV1::AliITSTrackerV1(AliITS* IITTSS, Int_t evnumber, Bool_t flag) {
//Origin A. Badala' and G.S. Pappalardo:
// e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it
// Class constructor. It does some initializations.
//PH Initialisation taken from the default constructor
fITS = IITTSS;
fresult = 0;
- frecPoints = 0;
+ fPtref = 0.;
+ fChi2max =0.;
+ frecPoints = 0;
fvettid = 0;
+ fflagvert = flag;
frl = 0;
fzmin = 0;
fzmax = 0;
fphimin = 0;
fphimax = 0;
fphidet = 0;
-
- fPtref = 0.;
- fChi2max =0.;
- fflagvert =flag;
+
Int_t imax = 200,jmax = 450;
frl = new AliITSRad(imax,jmax);
////////// gets information on geometry /////////////////////////////
- AliITSgeom *g1 = ((AliITS*)gAlice->GetDetector("ITS"))->GetITSgeom();
- // Why not AliITS *g1 = fITS->GetITSgeom(); // ?? BSN
+ AliITSgeom *g1 = fITS->GetITSgeom();
Int_t ll=1, dd=1;
TVector det(9);
- //cout<<" nlad ed ndet \n";
Int_t ia;
for(ia=0; ia<6; ia++) {
fNlad[ia]=g1->GetNladders(ia+1);
fNdet[ia]=g1->GetNdetectors(ia+1);
//cout<<fNlad[i]<<" "<<fNdet[i]<<"\n";
} // end for ia
- //getchar();
- //cout<<" raggio medio = ";
+ //cout<<" mean radius = ";
Int_t ib;
for(ib=0; ib<6; ib++) {
g1->GetCenterThetaPhi(ib+1,ll,dd,det);
//for(Int_t la=0; la<6; la++) cout<<" "<<fDetx[la]<<" "<<
// fDetz[la]<<endl;
//getchar();
+
// allocate memory and define matrices fzmin, fzmax, fphimin and fphimax //
Double_t epsz=1.2;
Double_t epszdrift=0.05;
fphidet[im1] = new Double_t[im2max];
} // end for im1
- Float_t global[3],local[3];
+ //Float_t global[3],local[3];
+ Double_t global[3],local[3];
Double_t pigre=TMath::Pi();
Double_t xmin,ymin,xmax,ymax;
if(fphimax[im1][im2]<0.) fphimax[im1][im2]+=2.*pigre;
} // end for im2
} // end for im1
+//////////////////////////////////////////////////////////////////////////////////////////////////////////
+/////////////// allocate memory and define vector fNRecPoints and matrices fRecCylR, fRecCylPhi, fRecCylZ /////////////
+ gAlice->GetEvent(evnumber);
+ Int_t NumOfModules = g1->GetIndexMax();
+ //fRecCylR = new Float_t *[NumOfModules];
+ fRecCylR = new Double_t *[NumOfModules];
+ //fRecCylPhi = new Float_t *[NumOfModules];
+ fRecCylPhi = new Double_t *[NumOfModules];
+ //fRecCylZ = new Float_t *[NumOfModules];
+ fRecCylZ = new Double_t *[NumOfModules];
+ AliITSRecPoint *recp;
+ fNRecPoints = new Int_t[NumOfModules];
+
+ for(Int_t module=0; module<NumOfModules; module++) {
+ fITS->ResetRecPoints();
+ gAlice->TreeR()->GetEvent(module);
+ frecPoints=fITS->RecPoints();
+ Int_t nRecPoints=fNRecPoints[module]=frecPoints->GetEntries();
+ /*
+ fRecCylR[module] = new Float_t[nRecPoints];
+ fRecCylPhi[module] = new Float_t[nRecPoints];
+ fRecCylZ[module] = new Float_t[nRecPoints];
+ */
+ fRecCylR[module] = new Double_t[nRecPoints];
+ fRecCylPhi[module] = new Double_t[nRecPoints];
+ fRecCylZ[module] = new Double_t[nRecPoints];
+ Int_t ind;
+ for(ind=0; ind<fNRecPoints[module]; ind++) {
+ recp=(AliITSRecPoint*)frecPoints->UncheckedAt(ind);
+ // Float_t global[3], local[3];
+ Double_t global[3], local[3];
+ local[0]=recp->GetX();
+ local[1]=0.;
+ local[2]= recp->GetZ();
+ g1->LtoG(module,local,global);
+ /*
+ Float_t r = TMath::Sqrt(global[0]*global[0]+global[1]*global[1]); // r hit
+ Float_t phi = TMath::ATan2(global[1],global[0]); if(phi<0.) phi+=2.*TMath::Pi(); // phi hit
+ Float_t z = global[2]; // z hit
+ */
+
+ Double_t r = TMath::Sqrt(global[0]*global[0]+global[1]*global[1]); // r hit
+ Double_t phi = TMath::ATan2(global[1],global[0]); if(phi<0.) phi+=2.*TMath::Pi(); // phi hit
+ Double_t z = global[2]; // z hit
+
+ fRecCylR[module][ind]=r;
+ fRecCylPhi[module][ind]=phi;
+ fRecCylZ[module][ind]=z;
+ }
+ }
+ //}
+ //}
+////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
////////// gets magnetic field factor //////////////////////////////
// Origin A. Badala' and G.S. Pappalardo:
// e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it
// copy constructor
-
- *fITS = *cobj.fITS;
+
+ *fITS = *cobj.fITS;
*fresult = *cobj.fresult;
fPtref = cobj.fPtref;
fChi2max = cobj.fChi2max;
fphidet[im1][im2]=cobj.fphidet[im1][im2];
} // end for im2
} // end for im2
+
+
+ AliITSgeom *g1 = fITS->GetITSgeom();
+ Int_t NumOfModules = g1->GetIndexMax();
+ /*
+ fRecCylR = new Float_t *[NumOfModules];
+ fRecCylPhi = new Float_t *[NumOfModules];
+ fRecCylZ = new Float_t *[NumOfModules];
+ */
+ fRecCylR = new Double_t *[NumOfModules];
+ fRecCylPhi = new Double_t *[NumOfModules];
+ fRecCylZ = new Double_t *[NumOfModules];
+ fNRecPoints = new Int_t[NumOfModules];
+ for(Int_t module=0; module<NumOfModules; module++) {
+ Int_t nRecPoints=fNRecPoints[module]=cobj.fNRecPoints[module];
+ /*
+ fRecCylR[module] = new Float_t[nRecPoints];
+ fRecCylPhi[module] = new Float_t[nRecPoints];
+ fRecCylZ[module] = new Float_t[nRecPoints];
+ */
+ fRecCylR[module] = new Double_t[nRecPoints];
+ fRecCylPhi[module] = new Double_t[nRecPoints];
+ fRecCylZ[module] = new Double_t[nRecPoints];
+ Int_t ind;
+ for(ind=0; ind<nRecPoints; ind++) {
+ fRecCylR[module][ind]=cobj.fRecCylR[module][ind];
+ fRecCylPhi[module][ind]=cobj.fRecCylPhi[module][ind];
+ fRecCylZ[module][ind]=cobj.fRecCylZ[module][ind];
+ }
+ }
+
+}
+void AliITSTrackerV1::DelMatrix(Int_t NumOfModules) {
+ for(Int_t mod=0; mod<NumOfModules; mod++) {
+ delete fRecCylR[mod];
+ delete fRecCylPhi[mod];
+ delete fRecCylZ[mod];
+ }
+ delete fRecCylR;
+ delete fRecCylPhi;
+ delete fRecCylZ;
}
//______________________________________________________________________
AliITSTrackerV1::~AliITSTrackerV1(){
// Origin A. Badala' and G.S. Pappalardo:
// e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it
- // class destructor
-
- //delete fTimerKalman;
- //delete fTimerIntersection;
+ // class destructor
+ delete frl;
+ delete fNRecPoints;
+ for(Int_t i=0; i<6; i++) {
+ delete fzmin[i];
+ delete fzmax[i];
+ delete fphimin[i];
+ delete fphimax[i];
+ delete fphidet[i];
+ }
+
+ delete fzmin;
+ delete fzmax;
+ delete fphimin;
+ delete fphimax;
+ delete fphidet;
+
}
//______________________________________________________________________
AliITSTrackerV1 &AliITSTrackerV1::operator=(AliITSTrackerV1 obj) {
// e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it
// assignement operator
- *fITS = *obj.fITS;
+ *fITS = *obj.fITS;
*fresult = *obj.fresult;
fPtref = obj.fPtref;
fChi2max = obj.fChi2max;
} // end for im2
} // end for im1
+ AliITSgeom *g1 = fITS->GetITSgeom();
+ Int_t NumOfModules = g1->GetIndexMax();
+ /*
+ fRecCylR = new Float_t *[NumOfModules];
+ fRecCylPhi = new Float_t *[NumOfModules];
+ fRecCylZ = new Float_t *[NumOfModules];
+ */
+ fRecCylR = new Double_t *[NumOfModules];
+ fRecCylPhi = new Double_t *[NumOfModules];
+ fRecCylZ = new Double_t *[NumOfModules];
+ fNRecPoints = new Int_t[NumOfModules];
+ for(Int_t module=0; module<NumOfModules; module++) {
+ Int_t nRecPoints=fNRecPoints[module]=obj.fNRecPoints[module];
+ /*
+ fRecCylR[module] = new Float_t[nRecPoints];
+ fRecCylPhi[module] = new Float_t[nRecPoints];
+ fRecCylZ[module] = new Float_t[nRecPoints];
+ */
+ fRecCylR[module] = new Double_t[nRecPoints];
+ fRecCylPhi[module] = new Double_t[nRecPoints];
+ fRecCylZ[module] = new Double_t[nRecPoints];
+ Int_t ind;
+ for(ind=0; ind<nRecPoints; ind++) {
+ fRecCylR[module][ind]=obj.fRecCylR[module][ind];
+ fRecCylPhi[module][ind]=obj.fRecCylPhi[module][ind];
+ fRecCylZ[module][ind]=obj.fRecCylZ[module][ind];
+ }
+ }
+
+
return *this;
}
//______________________________________________________________________
tracker->LoadOuterSectors();
// Load tracks
- TFile *tf=TFile::Open("AliTPCtracksSorted.root"); //modificato per hbt
+ TFile *tf=TFile::Open("AliTPCtracksSorted.root");
if (!tf->IsOpen()) {
cerr<<"Can't open AliTPCtracksSorted.root !\n";
return ;
///////////////////////////////////////////////
- /*
+ /*
////// propagation to the end of TPC //////////////
Double_t xk=77.415;
track->PropagateTo(xk, 28.94, 1.204e-3,mass); //Ne
xk -=0.5;
track->PropagateTo(xk, 41.28, 0.029,mass); //Nomex
////////////////////////////////////////////////////////////////////
- */
+ */
// new propagation to the end of TPC
- //Double_t xk=80.;
- //track->PropagateTo(xk,0.,0.); //Ne if it's still there
- Double_t xk=77.415;
+ Double_t xk=80.;
+ // track->PropagateTo(xk,0.,0.); //Ne if it's still there //attenzione funziona solo se modifica in TPC
+ // Double_t xk=77.415;
track->PropagateTo(xk, 28.94, 1.204e-3);
xk-=0.005;
track->PropagateTo(xk, 44.77,1.71); //Tedlar
track->PropagateTo(xk, 44.77, 1.71); //Tedlar
xk=61.;
- //track->PropagateTo(xk,0.,0.); //C02
- track->PropagateTo(xk,36.2,1.98e-3); //C02
+ // track->PropagateTo(xk,0.,0.); //C02
+ track->PropagateTo(xk,36.2,1.98e-3); //C02 //attenzione funziona solo se modifica in TPC
xk -=0.005;
track->PropagateTo(xk, 24.01, 2.7); //Al
xk -=0.005;
track->PropagateTo(xk, 44.77, 1.71); //Tedlar
xk -=0.005;
- track->PropagateTo(xk, 24.01, 2.7); //Al
-
-
+ track->PropagateTo(xk, 24.01, 2.7); //Al
+ ////////////////////////////////////////////////////////////////////////////////////////////////////////
AliITSTrackV1 trackITS(*track);
trackITS.PutMass(mass); //new to add mass to track
if(fresult){ delete fresult; fresult=0;}
if(pcode==11) vecLabRef(k)=p->GetFirstMother();
} // end if
itot++; vecTotLabRef(itot)=vecLabRef(k);
- if(vecLabRef(k)==0. && clustZ == 0.) vecTotLabRef(itot) =-3.;
+ if(vecLabRef(k)==0. && clustZ == -1.) vecTotLabRef(itot) =-3.;
} // end for k
} // end for lay
Long_t labref;
if(numClustNow<numClustRef && chi2Now>fresult->GetChi2()) continue;
//cout<<" chi2Now = "<<chi2Now<<"\n";
- // commentato il 30-7-2001
+
chi2Now/=numClustNow;
if(fPtref > 1.0 && chi2Now > 30.) continue;
if((fPtref >= 0.6 && fPtref<=1.0) && chi2Now > 40.) continue;
Int_t ladp, ladm, detp,detm,ladinters,detinters;
Int_t layerfin=layerInit-1;
// cout<<"Prima di intersection \n";
- //if(!fTimerIntersection) fTimerIntersection = new TStopwatch();
- // timer
- //fTimerIntersection->Continue();
- // timer
Int_t outinters=Intersection(*trackITS,layerfin,ladinters,detinters);
- //fTimerIntersection->Stop();
- // timer
// cout<<" outinters = "<<outinters<<"\n";
// cout<<" Layer ladder detector intersection ="
// <<layerfin<<" "<<ladinters<<" "<<detinters<<"\n";
*/
Float_t epsphi=5.0, epsz=5.0;
if(fPtref<0.2) {epsphi=3.; epsz=3.;}
- // nuova definizione idetot e toucLad e toucDet to be
+ // new definition of idetot e toucLad e toucDet to be
// transformed in a method
// these values could be modified
Float_t pigre=TMath::Pi();
Float_t distz = 0.0;
Float_t phicm, phicp, distphim, distphip;
phicm=phinters;
- if(phinters>fphimax[layerfin-1][ladm]) phicm=phinters-2*pigre;
- distphim=TMath::Abs(phicm-fphimax[layerfin-1][ladm]);
+ if(phinters>fphimax[layerfin-1][ladm-1]) phicm=phinters-2*pigre; //corretto il 20-11-2001
+ distphim=TMath::Abs(phicm-fphimax[layerfin-1][ladm-1]); //corretto il 20-11-2001
phicp=phinters;
- if(phinters>fphimin[layerfin-1][ladp]) phicp=phinters-2.*pigre;
- distphip=TMath::Abs(phicp-fphimin[layerfin-1][ladp]);
+ //cout<<" fNlad[layerfin-1] e ladp = "<<fNlad[layerfin-1]<<" "<<ladp<<endl;
+ if(phinters>fphimin[layerfin-1][ladp-1]) phicp=phinters-2.*pigre; //corretto il 20-11-2001
+ distphip=TMath::Abs(phicp-fphimin[layerfin-1][ladp-1]); //corretto il 20-11-2001
Int_t flagzmin=0;
Int_t flagzmax=0;
idetot=1;
///////////////////////////////////////////////////////
/*** Rec points sorted by module *****/
/**************************************/
- Int_t indexmod; //mod ott
- // AliITSRecPoint *recp;
+ Int_t indexmod;
indexmod = g1->GetModuleIndex(lycur,(Int_t)toucLad(iriv),
- (Int_t)toucDet(iriv));
- //mod ott
+ (Int_t)toucDet(iriv));
fITS->ResetRecPoints();
gAlice->TreeR()->GetEvent(indexmod);
Int_t npoints=frecPoints->GetEntries();
- /* mod ott
- Int_t *indlist=new Int_t[npoints+1];
- Int_t counter=0;
- Int_t ind;
- for (ind=0; ind<=npoints; ind++) {
- indlist[ind]=-1;
- if (*(fvettid[index]+ind)==0) {
- indlist[counter]=ind;
- counter++;
- } // end if
- } // end for ind
- ind=-1;
- for(;;) {
- ind++;
- if(indlist[ind] < 0) recp=0;
- else
- recp = (AliITSRecPoint*)frecPoints->
- UncheckedAt(indlist[ind]);
- if((!recp) ) break;
- } // end for ;;
- */
- ///////////////////////// new //////////////////////////
+
Int_t indnew;
for(indnew=0; indnew<npoints; indnew++){
if (*(fvettid[indexmod]+indnew)==0)
recp =(AliITSRecPoint*)frecPoints->UncheckedAt(indnew);
else
continue;
- ////////////////////////////////////////////////////
TVector cluster(3),vecclust(9);
- vecclust(6)=vecclust(7)=vecclust(8)=-1.;
- Double_t sigma[2];
- // set veclust in global
- Float_t global[3], local[3];
- local[0]=recp->GetX();
- local[1]=0.;
- local[2]= recp->GetZ();
- Int_t play = lycur;
- Int_t plad = TMath::Nint(toucLad(iriv));
- Int_t pdet = TMath::Nint(toucDet(iriv));
- g1->LtoG(play,plad,pdet,local,global);
- vecclust(0)=global[0];
- vecclust(1)=global[1];
- vecclust(2)=global[2];
-
+ //vecclust(6)=vecclust(7)=vecclust(8)=-1.;
+ Double_t sigma[2];
+ // now vecclust is with cylindrical cohordinates
+ vecclust(0)=(Float_t)fRecCylR[indexmod][indnew];
+ vecclust(1)=(Float_t)fRecCylPhi[indexmod][indnew];
+ vecclust(2)=(Float_t)fRecCylZ[indexmod][indnew];
vecclust(3) = (Float_t)recp->fTracks[0];
- //vecclust(4) = (Float_t)indlist[ind];
vecclust(4) = (Float_t)indnew;
- vecclust(5) = (Float_t)indexmod; //mod ott
+ vecclust(5) = (Float_t)indexmod;
vecclust(6) = (Float_t)recp->fTracks[0];
vecclust(7) = (Float_t)recp->fTracks[1];
vecclust(8) = (Float_t)recp->fTracks[2];
sigma[0] = (Double_t) recp->GetSigmaX2();
sigma[1] = (Double_t) recp->GetSigmaZ2();
- //now we are in r,phi,z in global
- cluster(0) = TMath::Sqrt(vecclust(0)*vecclust(0)+
- vecclust(1)*vecclust(1));//r hit
- cluster(1) = TMath::ATan2(vecclust(1),vecclust(0));
- if(cluster(1)<0.) cluster(1)+=2.*TMath::Pi();
- cluster(2) = vecclust(2); // z hit
+
+ cluster(0)=fRecCylR[indexmod][indnew];
+ cluster(1)=fRecCylPhi[indexmod][indnew];
+ cluster(2)=fRecCylZ[indexmod][indnew];
+
// cout<<" layer = "<<play<<"\n";
// cout<<" cluster prima = "<<vecclust(0)<<" "
// <<vecclust(1)<<" "
//matrix to the covariance matrix
(*newTrack).AddEL(frl,1.,0);
- //if(!fTimerKalman) fTimerKalman = new TStopwatch();//timer
- //fTimerKalman->Continue(); // timer
if(fflagvert){
KalmanFilterVert(newTrack,cluster,sigmanew);
//KalmanFilterVert(newTrack,cluster,sigmanew,chi2pred);
}else{
KalmanFilter(newTrack,cluster,sigmanew);
} // end if
- //fTimerKalman->Stop(); // timer
(*newTrack).PutCluster(layernew, vecclust);
newTrack->AddClustInTrack();
listoftrack.AddLast(newTrack);
} // end for indnew
- // delete [] indlist; //mod ott
} // end of for on detectors (iriv)
}//end if(outinters==0)
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff