From b48af42837326738183c8897f1d4c55c382d00c2 Mon Sep 17 00:00:00 2001 From: nilsen Date: Fri, 19 Oct 2001 21:32:35 +0000 Subject: [PATCH] Minor changes to remove compliation warning on gcc 2.92.2 compiler, and cleanded up a little bit of code. --- ITS/AliITSClusterFinderSDD.cxx | 105 +- ITS/AliITSClusterFinderSPD.cxx | 8 +- ITS/AliITSRawCluster.cxx | 407 ++--- ITS/AliITSTrackerV1.cxx | 2852 ++++++++++++++++---------------- ITS/AliITSTrackerV1.h | 119 +- ITS/AliITSgeomSDD.cxx | 17 +- ITS/AliITSgeomSPD.cxx | 22 +- ITS/AliITSsimulationSSD.cxx | 8 +- ITS/AliITSvPPRsymm.cxx | 5 +- 9 files changed, 1801 insertions(+), 1742 deletions(-) diff --git a/ITS/AliITSClusterFinderSDD.cxx b/ITS/AliITSClusterFinderSDD.cxx index fe2a41b8c50..b26aecf4e1f 100644 --- a/ITS/AliITSClusterFinderSDD.cxx +++ b/ITS/AliITSClusterFinderSDD.cxx @@ -36,11 +36,11 @@ AliITSClusterFinderSDD::AliITSClusterFinderSDD(AliITSsegmentation *seg, TClonesArray *recp){ // standard constructor - fSegmentation=seg; - fResponse=response; - fDigits=digits; - fClusters=recp; - fNclusters= fClusters->GetEntriesFast(); + fSegmentation = seg; + fResponse = response; + fDigits = digits; + fClusters = recp; + fNclusters = fClusters->GetEntriesFast(); SetCutAmplitude(); SetDAnode(); SetDTime(); @@ -49,19 +49,19 @@ AliITSClusterFinderSDD::AliITSClusterFinderSDD(AliITSsegmentation *seg, SetMaxNCells(); SetTimeCorr(); SetMinCharge(); - fMap=new AliITSMapA1(fSegmentation,fDigits,fCutAmplitude); + fMap = new AliITSMapA1(fSegmentation,fDigits,fCutAmplitude); } //______________________________________________________________________ AliITSClusterFinderSDD::AliITSClusterFinderSDD(){ // default constructor - fSegmentation=0; - fResponse=0; - fDigits=0; - fClusters=0; - fNclusters=0; - fMap=0; - fCutAmplitude=0; + fSegmentation = 0; + fResponse = 0; + fDigits = 0; + fClusters = 0; + fNclusters = 0; + fMap = 0; + fCutAmplitude = 0; SetDAnode(); SetDTime(); SetMinPeak(); @@ -83,22 +83,22 @@ void AliITSClusterFinderSDD::SetCutAmplitude(Float_t nsigma){ fResponse->GetNoiseParam(noise,baseline); noise_after_el = ((AliITSresponseSDD*)fResponse)->GetNoiseAfterElectronics(); - fCutAmplitude=(Int_t)((baseline + nsigma*noise_after_el)); + fCutAmplitude = (Int_t)((baseline + nsigma*noise_after_el)); } //______________________________________________________________________ void AliITSClusterFinderSDD::Find1DClusters(){ // find 1D clusters - static AliITS *iTS=(AliITS*)gAlice->GetModule("ITS"); + static AliITS *iTS = (AliITS*)gAlice->GetModule("ITS"); // retrieve the parameters - Int_t fNofMaps = fSegmentation->Npz(); + Int_t fNofMaps = fSegmentation->Npz(); Int_t fMaxNofSamples = fSegmentation->Npx(); - Int_t fNofAnodes = fNofMaps/2; - Int_t dummy=0; - Float_t fTimeStep = fSegmentation->Dpx(dummy); - Float_t fSddLength = fSegmentation->Dx(); - Float_t fDriftSpeed = fResponse->DriftSpeed(); - Float_t anodePitch = fSegmentation->Dpz(dummy); + Int_t fNofAnodes = fNofMaps/2; + Int_t dummy = 0; + Float_t fTimeStep = fSegmentation->Dpx(dummy); + Float_t fSddLength = fSegmentation->Dx(); + Float_t fDriftSpeed = fResponse->DriftSpeed(); + Float_t anodePitch = fSegmentation->Dpz(dummy); // map the signal fMap->SetThreshold(fCutAmplitude); @@ -112,7 +112,7 @@ void AliITSClusterFinderSDD::Find1DClusters(){ Int_t i; Float_t **dfadc = new Float_t*[fNofAnodes]; for(i=0;i=fMaxNofSamples) break; @@ -154,7 +154,7 @@ void AliITSClusterFinderSDD::Find1DClusters(){ } // end if if(dfadc[k][id] > dfadcmax) { dfadcmax = dfadc[k][id]; - imaxd = id; + imaxd = id; } // end if } // end for m it = imaxd; @@ -166,7 +166,7 @@ void AliITSClusterFinderSDD::Find1DClusters(){ if(lthrt >= lthrmint && lthra >= lthrmina) ilcl = 1; if(ilcl) { nofFoundClusters++; - Int_t tstop = tstart; + Int_t tstop = tstart; Float_t dfadcmin = 10000.; Int_t ij; for(ij=0; ij<20; ij++) { @@ -181,17 +181,17 @@ void AliITSClusterFinderSDD::Find1DClusters(){ } // end for ij Float_t clusterCharge = 0.; - Float_t clusterAnode = k+0.5; - Float_t clusterTime = 0.; - Float_t clusterMult = 0.; + Float_t clusterAnode = k+0.5; + Float_t clusterTime = 0.; + Int_t clusterMult = 0; Float_t clusterPeakAmplitude = 0.; - Int_t its,peakpos=-1; + Int_t its,peakpos = -1; Float_t n, baseline; fResponse->GetNoiseParam(n,baseline); for(its=tstart; its<=tstop; its++) { fadc=(float)fMap->GetSignal(idx,its); - if(fadc>baseline) fadc-=baseline; - else fadc=0.; + if(fadc>baseline) fadc -= baseline; + else fadc = 0.; clusterCharge += fadc; // as a matter of fact we should take the peak // pos before FFT @@ -199,32 +199,35 @@ void AliITSClusterFinderSDD::Find1DClusters(){ if(fadc > clusterPeakAmplitude) { clusterPeakAmplitude = fadc; //peakpos=fMap->GetHitIndex(idx,its); - Int_t shift=(int)(fTimeCorr/fTimeStep); + Int_t shift = (int)(fTimeCorr/fTimeStep); if(its>shift && its<(fMaxNofSamples-shift)) - peakpos=fMap->GetHitIndex(idx,its+shift); - else peakpos=fMap->GetHitIndex(idx,its); - if(peakpos<0) peakpos=fMap->GetHitIndex(idx,its); + peakpos = fMap->GetHitIndex(idx,its+shift); + else peakpos = fMap->GetHitIndex(idx,its); + if(peakpos<0) peakpos =fMap->GetHitIndex(idx,its); } // end if clusterTime += fadc*its; if(fadc > 0) clusterMult++; if(its == tstop) { clusterTime /= (clusterCharge/fTimeStep); // ns - if(clusterTime>fTimeCorr) clusterTime-=fTimeCorr; + if(clusterTime>fTimeCorr) clusterTime -=fTimeCorr; //ns } // end if } // end for its Float_t clusteranodePath = (clusterAnode - fNofAnodes/2)* - anodePitch; + anodePitch; Float_t clusterDriftPath = clusterTime*fDriftSpeed; clusterDriftPath = fSddLength-clusterDriftPath; if(clusterCharge <= 0.) break; - AliITSRawClusterSDD clust(j+1,clusterAnode,clusterTime, - clusterCharge, - clusterPeakAmplitude, - peakpos,0.,0.,clusterDriftPath, - clusteranodePath,clusterMult,0,0, - 0,0,0,0,0); + AliITSRawClusterSDD clust(j+1,//i + clusterAnode,clusterTime,//ff + clusterCharge, //f + clusterPeakAmplitude, //f + peakpos, //i + 0.,0.,clusterDriftPath,//fff + clusteranodePath, //f + clusterMult, //i + 0,0,0,0,0,0,0);//7*i iTS->AddCluster(1,&clust); it = tstop; } // ilcl diff --git a/ITS/AliITSClusterFinderSPD.cxx b/ITS/AliITSClusterFinderSPD.cxx index 1ec3b830121..b8d53a3faa3 100644 --- a/ITS/AliITSClusterFinderSPD.cxx +++ b/ITS/AliITSClusterFinderSPD.cxx @@ -382,7 +382,13 @@ void AliITSClusterFinderSPD::ClusterFinder(Int_t ndigits, zcenter[i] = zcenter[i] - fSegmentation->Dz()/2.; - AliITSRawClusterSPD *clust = new AliITSRawClusterSPD(zcenter[i],xcenter[i],ndig,ndz,ndx,0.,0.,0.,0.,0.,0.,0.); + AliITSRawClusterSPD *clust = new AliITSRawClusterSPD(zcenter[i], //f + xcenter[i], //f + ndig, //f + ndz,ndx, //ii + 0,0,0,0, //iiii + 0.,0., //ff + 0); //i iTS->AddCluster(0,clust); delete clust; }//end loop on clusters diff --git a/ITS/AliITSRawCluster.cxx b/ITS/AliITSRawCluster.cxx index 3b9e14aafb2..8698ca668c8 100644 --- a/ITS/AliITSRawCluster.cxx +++ b/ITS/AliITSRawCluster.cxx @@ -1,214 +1,247 @@ +/************************************************************************** + * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * * + * Author: The ALICE Off-line Project. * + * Contributors are mentioned in the code where appropriate. * + * * + * Permission to use, copy, modify and distribute this software and its * + * documentation strictly for non-commercial purposes is hereby granted * + * without fee, provided that the above copyright notice appears in all * + * copies and that both the copyright notice and this permission notice * + * appear in the supporting documentation. The authors make no claims * + * about the suitability of this software for any purpose. It is * + * provided "as is" without express or implied warranty. * + **************************************************************************/ + +/* +$Log$ +*/ #include #include #include "AliITSRawCluster.h" ClassImp(AliITSRawCluster) - ClassImp(AliITSRawClusterSDD) -//-------------------------------------- -AliITSRawClusterSDD::AliITSRawClusterSDD(Int_t wing, Float_t Anode,Float_t Time, Float_t Charge, - Float_t PeakAmplitude,Int_t PeakPosition, Float_t Asigma, Float_t Tsigma, Float_t DriftPath, Float_t AnodeOffset, Int_t Samples, - Int_t Tstart, Int_t Tstop, Int_t Tstartf, Int_t Tstopf, Int_t Anodes, Int_t Astart, Int_t Astop) -{ - // constructor - fWing = wing; - fAnode = Anode; - fTime = Time; - fQ = Charge; - fPeakAmplitude = PeakAmplitude; - fPeakPosition = PeakPosition; - fAsigma = Asigma; - fTsigma = Tsigma; - fNanodes = Anodes; - fTstart = Tstart; - fTstop = Tstop; - fTstartf = Tstartf; - fTstopf = Tstopf; - fAstart = Astart; - fAstop = Astop; - fMultiplicity = Samples; - fSumAmplitude = 0; - Int_t sign = 1; - for(Int_t i=0;iA()*clJ->Q())/(fQ+clJ->Q()); - fTime = (fTime*fQ + clJ->T()*clJ->Q())/(fQ+clJ->Q()); - fX = (fX*fQ + clJ->X()*clJ->Q())/(fQ+clJ->Q()); - fZ = (fZ*fQ + clJ->Z()*clJ->Q())/(fQ+clJ->Q()); - fQ += clJ->Q(); - if(fSumAmplitude == 0) fSumAmplitude += fPeakAmplitude; - /* - fAnode = (fAnode*fSumAmplitude+clJ->A()*clJ->PeakAmpl())/(fSumAmplitude+clJ->PeakAmpl()); - fTime = (fTime*fSumAmplitude +clJ->T()*clJ->PeakAmpl())/(fSumAmplitude+clJ->PeakAmpl()); - fX = (fX*fSumAmplitude +clJ->X()*clJ->PeakAmpl())/(fSumAmplitude+clJ->PeakAmpl()); - fZ = (fZ*fSumAmplitude +clJ->Z()*clJ->PeakAmpl())/(fSumAmplitude+clJ->PeakAmpl()); - */ - fSumAmplitude += clJ->PeakAmpl(); - - fTstart = clJ->Tstart(); - fTstop = clJ->Tstop(); - if(fTstartf > clJ->Tstartf()) fTstartf = clJ->Tstartf(); - if(fTstopf < clJ->Tstopf()) fTstopf = clJ->Tstopf(); - if(fAstop < clJ->Astop()) fAstop = clJ->Astop(); - - fMultiplicity += (Int_t) (clJ->Samples()); - (fNanodes)++; - if(clJ->PeakAmpl() > fPeakAmplitude) { - fPeakAmplitude = clJ->PeakAmpl(); - fPeakPosition = clJ->PeakPos(); - } - - return; -} - -//-------------------------------------- -Bool_t AliITSRawClusterSDD::Brother(AliITSRawClusterSDD* cluster,Float_t danode,Float_t dtime) { + // add + + fAnode = (fAnode*fQ + clJ->A()*clJ->Q())/(fQ+clJ->Q()); + fTime = ( fTime*fQ + clJ->T()*clJ->Q())/(fQ+clJ->Q()); + fX = ( fX*fQ + clJ->X()*clJ->Q())/(fQ+clJ->Q()); + fZ = ( fZ*fQ + clJ->Z()*clJ->Q())/(fQ+clJ->Q()); + fQ += clJ->Q(); + if(fSumAmplitude == 0) fSumAmplitude += fPeakAmplitude; + /* + fAnode = (fAnode*fSumAmplitude+clJ->A()*clJ->PeakAmpl())/ + (fSumAmplitude+clJ->PeakAmpl()); + fTime = (fTime*fSumAmplitude +clJ->T()*clJ->PeakAmpl())/ + (fSumAmplitude+clJ->PeakAmpl()); + fX = (fX*fSumAmplitude +clJ->X()*clJ->PeakAmpl())/ + (fSumAmplitude+clJ->PeakAmpl()); + fZ = (fZ*fSumAmplitude +clJ->Z()*clJ->PeakAmpl())/ + (fSumAmplitude+clJ->PeakAmpl()); + */ + fSumAmplitude += clJ->PeakAmpl(); + + fTstart = clJ->Tstart(); + fTstop = clJ->Tstop(); + if(fTstartf > clJ->Tstartf()) fTstartf = clJ->Tstartf(); + if( fTstopf < clJ->Tstopf() ) fTstopf = clJ->Tstopf(); + if( fAstop < clJ->Astop() ) fAstop = clJ->Astop(); + + fMultiplicity += (Int_t) (clJ->Samples()); + (fNanodes)++; + if(clJ->PeakAmpl() > fPeakAmplitude) { + fPeakAmplitude = clJ->PeakAmpl(); + fPeakPosition = clJ->PeakPos(); + } // end if + + return; +} +//______________________________________________________________________ +Bool_t AliITSRawClusterSDD::Brother(AliITSRawClusterSDD* cluster, + Float_t danode,Float_t dtime) { - Bool_t brother = kFALSE; + Bool_t brother = kFALSE; - Bool_t test2 = kFALSE; - Bool_t test3 = kFALSE; - Bool_t test4 = kFALSE; - Bool_t test5 = kFALSE; + Bool_t test2 = kFALSE; + Bool_t test3 = kFALSE; + Bool_t test4 = kFALSE; + Bool_t test5 = kFALSE; - if(fWing != cluster->W()) return brother; + if(fWing != cluster->W()) return brother; - if(fTstopf >= cluster->Tstart() && fTstartf <= cluster->Tstop()) test2 = kTRUE; - if(cluster->Astop() == (fAstop+1)) test3 = kTRUE; + if(fTstopf >= cluster->Tstart() && + fTstartf <= cluster->Tstop()) test2 = kTRUE; + if(cluster->Astop() == (fAstop+1)) test3 = kTRUE; - if(TMath::Abs(fTime-cluster->T()) < dtime) test4 = kTRUE; - if(TMath::Abs(fAnode-cluster->A()) < danode) test5 = kTRUE; + if(TMath::Abs(fTime-cluster->T()) < dtime) test4 = kTRUE; + if(TMath::Abs(fAnode-cluster->A()) < danode) test5 = kTRUE; - if((test2 && test3) || (test4 && test5) ) { - return brother = kTRUE; - } + if((test2 && test3) || (test4 && test5) ) { + return brother = kTRUE; + } // end if - return brother; + return brother; } - -//-------------------------------------- +//______________________________________________________________________ void AliITSRawClusterSDD::PrintInfo() { - // print - cout << ", Anode " << fAnode << ", Time: " << fTime << ", Charge: " << fQ; - cout << ", Samples: " << fMultiplicity; - cout << ", X: " << fX << ", Z: " << fZ << "tstart " << fTstart << "tstop "<< fTstop <fZStop < clJ->ZStop()) this->fZStop = clJ->ZStop(); - - this->fZ = this->fZ + clJ->Z(); - this->fX = (this->fX + clJ->X())/2.; - this->fQ = this->fQ + clJ->Q(); - - this->fXStart = clJ->XStart(); // for a comparison with the next - this->fXStop = clJ->XStop(); // z column - - if(this->fXStartf > clJ->XStartf()) this->fXStartf = clJ->XStartf(); - if(this->fXStopf < clJ->XStopf()) this->fXStopf = clJ->XStopf(); - if(this->fZend < clJ->Zend()) this->fZend = clJ->Zend(); - this->fNClX = this->fXStopf - this->fXStartf + 1; - (this->fNClZ)++; - - return; -} - -//-------------------------------------- -Bool_t AliITSRawClusterSPD::Brother(AliITSRawClusterSPD* cluster,Float_t dz,Float_t dx) { - // fXStart, fXstop and fZend information is used now instead of dz and dx - // to check an absent (or a present) of the gap between two pixels in - // both x and z directions. The increasing order of fZend is used. - - Bool_t brother = kFALSE; - Bool_t test2 = kFALSE; - Bool_t test3 = kFALSE; - - // Diagonal clusters are included: - if(fXStop >= (cluster->XStart() -1) && fXStart <= (cluster->XStop()+1)) test2 = kTRUE; - - // Diagonal clusters are excluded: - // if(fXStop >= cluster->XStart() && fXStart <= cluster->XStop()) test2 = kTRUE; - if(cluster->Zend() == (fZend + 1)) test3 = kTRUE; - if(test2 && test3) { - // cout<<"test 2,3 0k, brother = true "<fZStop < clJ->ZStop()) this->fZStop = clJ->ZStop(); + this->fZ = this->fZ + clJ->Z(); + this->fX = (this->fX + clJ->X())/2.; + this->fQ = this->fQ + clJ->Q(); + this->fXStart = clJ->XStart(); // for a comparison with the next + this->fXStop = clJ->XStop(); // z column + if(this->fXStartf > clJ->XStartf()) this->fXStartf = clJ->XStartf(); + if(this->fXStopf < clJ->XStopf()) this->fXStopf = clJ->XStopf(); + if(this->fZend < clJ->Zend()) this->fZend = clJ->Zend(); + this->fNClX = this->fXStopf - this->fXStartf + 1; + (this->fNClZ)++; + + return; } - -//-------------------------------------- -void AliITSRawClusterSPD::PrintInfo() -{ - // print - cout << ", Z: " << fZ << ", X: " << fX << ", Charge: " << fQ<= (cluster->XStart() -1) && + fXStart <= (cluster->XStop()+1)) test2 = kTRUE; + + // Diagonal clusters are excluded: + // if(fXStop >= cluster->XStart() && + // fXStart <= cluster->XStop()) test2 = kTRUE; + if(cluster->Zend() == (fZend + 1)) test3 = kTRUE; + if(test2 && test3) { + // cout<<"test 2,3 0k, brother = true "< #include @@ -40,1450 +58,1420 @@ #include "AliITSVertex.h" ClassImp(AliITSTrackerV1) - - -//________________________________________________________________ - +//______________________________________________________________________ AliITSTrackerV1::AliITSTrackerV1(AliITS* IITTSS, 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. - - fITS = IITTSS; - 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(); - - 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<GetShape(1, ll, dd)))->GetDx(); - fDetz[0] = ((AliITSgeomSPD*)(g1->GetShape(1, ll, dd)))->GetDz(); - - fDetx[1] = ((AliITSgeomSPD*)(g1->GetShape(2, ll, dd)))->GetDx(); - fDetz[1] = ((AliITSgeomSPD*)(g1->GetShape(2, ll, dd)))->GetDz(); - - fDetx[2] = ((AliITSgeomSDD*)(g1->GetShape(3, ll, dd)))->GetDx(); - fDetz[2] = ((AliITSgeomSDD*)(g1->GetShape(3, ll, dd)))->GetDz(); - - fDetx[3] = ((AliITSgeomSDD*)(g1->GetShape(4, ll, dd)))->GetDx(); - fDetz[3] = ((AliITSgeomSDD*)(g1->GetShape(4, ll, dd)))->GetDz(); - - fDetx[4] = ((AliITSgeomSSD*)(g1->GetShape(5, ll, dd)))->GetDx(); - fDetz[4] = ((AliITSgeomSSD*)(g1->GetShape(5, ll, dd)))->GetDz(); - - fDetx[5] = ((AliITSgeomSSD*)(g1->GetShape(6, ll, dd)))->GetDx(); - fDetz[5] = ((AliITSgeomSSD*)(g1->GetShape(6, ll, dd)))->GetDz(); - - //cout<<" Detx Detz\n"; - //for(Int_t la=0; la<6; la++) cout<<" "<GetCenterThetaPhi(im1+1,1,im2+1,det); - if(im2!=0) fzmin[im1][im2]=det(2)-fDetz[im1]; - else - fzmin[im1][im2]=det(2)-(fDetz[im1])*epsz; - if(im2!=(im2max-1)) fzmax[im1][im2]=det(2)+fDetz[im1]; - else + //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. + + fITS = IITTSS; + 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 + 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<GetShape(1, ll, dd)))->GetDx(); + fDetz[0] = ((AliITSgeomSPD*)(g1->GetShape(1, ll, dd)))->GetDz(); + + fDetx[1] = ((AliITSgeomSPD*)(g1->GetShape(2, ll, dd)))->GetDx(); + fDetz[1] = ((AliITSgeomSPD*)(g1->GetShape(2, ll, dd)))->GetDz(); + + fDetx[2] = ((AliITSgeomSDD*)(g1->GetShape(3, ll, dd)))->GetDx(); + fDetz[2] = ((AliITSgeomSDD*)(g1->GetShape(3, ll, dd)))->GetDz(); + + fDetx[3] = ((AliITSgeomSDD*)(g1->GetShape(4, ll, dd)))->GetDx(); + fDetz[3] = ((AliITSgeomSDD*)(g1->GetShape(4, ll, dd)))->GetDz(); + + fDetx[4] = ((AliITSgeomSSD*)(g1->GetShape(5, ll, dd)))->GetDx(); + fDetz[4] = ((AliITSgeomSSD*)(g1->GetShape(5, ll, dd)))->GetDz(); + + fDetx[5] = ((AliITSgeomSSD*)(g1->GetShape(6, ll, dd)))->GetDx(); + fDetz[5] = ((AliITSgeomSSD*)(g1->GetShape(6, ll, dd)))->GetDz(); + //cout<<" Detx Detz\n"; + //for(Int_t la=0; la<6; la++) cout<<" "<GetCenterThetaPhi(im1+1,1,im2+1,det); + if(im2!=0) fzmin[im1][im2]=det(2)-fDetz[im1]; + else + fzmin[im1][im2]=det(2)-(fDetz[im1])*epsz; + if(im2!=(im2max-1)) fzmax[im1][im2]=det(2)+fDetz[im1]; + else fzmax[im1][im2]=det(2)+fDetz[im1]*epsz; - if(im1==2 || im1==3) {fzmin[im1][im2]-=epszdrift; fzmax[im1][im2]+=epszdrift;} - } - } - - fphimin = new Double_t*[6]; fphimax = new Double_t*[6]; - for(im1=0;im1<6;im1++) { - im2max=fNlad[im1]; - fphimin[im1] = new Double_t[im2max]; fphimax[im1] = new Double_t[im2max]; - } - - fphidet = new Double_t*[6]; - for(im1=0; im1<6; im1++) { - im2max=fNlad[im1]; - fphidet[im1] = new Double_t[im2max]; - } - - Float_t global[3],local[3]; - Double_t pigre=TMath::Pi(); - Double_t xmin,ymin,xmax,ymax; - - for(im1=0; im1<6; im1++) { - im2max=fNlad[im1]; - for(im2=0; im2GetCenterThetaPhi(im1+1,im2+1,idet,det); - fphidet[im1][im2]= TMath::ATan2(Double_t(det(1)),Double_t(det(0))); - if(fphidet[im1][im2]<0.) fphidet[im1][im2]+=2.*pigre; - local[1]=local[2]=0.; - local[0]= -(fDetx[im1]); - if(im1==0) local[0]= (fDetx[im1]); //to take into account different reference system - g1->LtoG(im1+1,im2+1,idet,local,global); - xmax=global[0]; ymax=global[1]; - local[0]= (fDetx[im1]); - if(im1==0) local[0]= -(fDetx[im1]); //take into account different reference system - g1->LtoG(im1+1,im2+1,idet,local,global); - xmin=global[0]; ymin=global[1]; - fphimin[im1][im2]= TMath::ATan2(ymin,xmin); if(fphimin[im1][im2]<0.) fphimin[im1][im2]+=2.*pigre; - fphimax[im1][im2]= TMath::ATan2(ymax,xmax); if(fphimax[im1][im2]<0.) fphimax[im1][im2]+=2.*pigre; - } - } - -////////////////////////////////////////////////////////////////////////////////////////////////////////// - -//////////////////////////////////////// gets magnetic field factor //////////////////////////////// - - AliMagF * fieldPointer = gAlice->Field(); - fFieldFactor = (Double_t)fieldPointer->Factor(); - //cout<< " field factor = "<GetCenterThetaPhi(im1+1,im2+1,idet,det); + fphidet[im1][im2] = TMath::ATan2(Double_t(det(1)), + Double_t(det(0))); + if(fphidet[im1][im2]<0.) fphidet[im1][im2]+=2.*pigre; + local[1]=local[2]=0.; + local[0]= -(fDetx[im1]); + if(im1==0) local[0]= (fDetx[im1]); //to take into account + // different reference system + g1->LtoG(im1+1,im2+1,idet,local,global); + xmax=global[0]; ymax=global[1]; + local[0]= (fDetx[im1]); + if(im1==0) local[0]= -(fDetx[im1]);//take into account different + // reference system + g1->LtoG(im1+1,im2+1,idet,local,global); + xmin=global[0]; ymin=global[1]; + fphimin[im1][im2]= TMath::ATan2(ymin,xmin); + if(fphimin[im1][im2]<0.) fphimin[im1][im2]+=2.*pigre; + fphimax[im1][im2]= TMath::ATan2(ymax,xmax); + if(fphimax[im1][im2]<0.) fphimax[im1][im2]+=2.*pigre; + } // end for im2 + } // end for im1 + + ////////// gets magnetic field factor ////////////////////////////// + + AliMagF * fieldPointer = gAlice->Field(); + fFieldFactor = (Double_t)fieldPointer->Factor(); + //cout<< " field factor = "<GetEvent(evNumber); //modificato per gestire hbt - - AliKalmanTrack *kkprov; - kkprov->SetConvConst(100/0.299792458/0.2/fFieldFactor); - - - TFile *cf=TFile::Open("AliTPCclusters.root"); - AliTPCParam *digp= (AliTPCParam*)cf->Get("75x40_100x60"); - if (!digp) { cerr<<"TPC parameters have not been found !\n"; getchar();} - - cf->cd(); - AliTPCtracker *tracker = new AliTPCtracker(digp,evNumber); - - // Load clusters - tracker->LoadInnerSectors(); - tracker->LoadOuterSectors(); - - -// Load tracks - TFile *tf=TFile::Open("AliTPCtracksSorted.root"); //modificato per hbt - if (!tf->IsOpen()) {cerr<<"Can't open AliTPCtracksSorted.root !\n"; return ;} - TObjArray tracks(200000); - char tname[100]; - sprintf(tname,"TreeT_TPC_%d",evNumber); - - TTree *tracktree=(TTree*)tf->Get(tname); - 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; kkSetAddress(&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 "<TreeR(); - Int_t nent=(Int_t)tr->GetEntries(); - frecPoints = fITS->RecPoints(); - - Int_t numbpoints; - Int_t totalpoints=0; - Int_t *np = new Int_t[nent]; - fvettid = new Int_t* [nent]; - Int_t mod; - - for (mod=0; modResetRecPoints(); - gAlice->TreeR()->GetEvent(mod); - numbpoints = frecPoints->GetEntries(); - totalpoints+=numbpoints; - np[mod] = numbpoints; - //cout<<" mod = "<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 - - /////////////////////////////////////////////////////////////// +//______________________________________________________________________ +void AliITSTrackerV1::DoTracking(Int_t evNumber,Int_t minTr,Int_t maxTr, + TFile *file) { + // Origin A. Badala' and G.S. Pappalardo: + // e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it + // The method needs the event number, the minimum and maximum order + // number of TPC tracks that + // are to be tracked trough the ITS, and the file where the recpoints + // are registered. + // The method can be called by a macro. It preforms the tracking for + // all good TPC tracks + + printf("begin DoTracking - file %p\n",file); + + gAlice->GetEvent(evNumber); //modificato per gestire hbt - - AliITSTrackV1 trackITS(*track); - - if(fresult){ delete fresult; fresult=0;} - fresult = new AliITSTrackV1(trackITS); - - AliITSTrackV1 primaryTrack(trackITS); - - vgeant=(*fresult).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; + AliKalmanTrack *kkprov; + kkprov->SetConvConst(100/0.299792458/0.2/fFieldFactor); + + TFile *cf=TFile::Open("AliTPCclusters.root"); + AliTPCParam *digp= (AliTPCParam*)cf->Get("75x40_100x60"); + if (!digp) { cerr<<"TPC parameters have not been found !\n"; getchar();} + + cf->cd(); + AliTPCtracker *tracker = new AliTPCtracker(digp,evNumber); + + // Load clusters + tracker->LoadInnerSectors(); + tracker->LoadOuterSectors(); + + // Load tracks + TFile *tf=TFile::Open("AliTPCtracksSorted.root"); //modificato per hbt + if (!tf->IsOpen()) { + cerr<<"Can't open AliTPCtracksSorted.root !\n"; + return ; + } // end if + TObjArray tracks(200000); + char tname[100]; + sprintf(tname,"TreeT_TPC_%d",evNumber); + + TTree *tracktree=(TTree*)tf->Get(tname); + 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; kkSetAddress(&ioTrackTPC); + tracktree->GetEvent(kk); + tracker->CookLabel(ioTrackTPC,0.1); + tracks.AddLast(ioTrackTPC); + } // end for kk + delete tracker; + tf->Close(); + + Int_t nt = tracks.GetEntriesFast(); + cerr<<"Number of found tracks "<TreeR(); + Int_t nent=(Int_t)tr->GetEntries(); + frecPoints = fITS->RecPoints(); + + Int_t numbpoints; + Int_t totalpoints=0; + Int_t *np = new Int_t[nent]; + fvettid = new Int_t* [nent]; + Int_t mod; + + for (mod=0; modResetRecPoints(); + gAlice->TreeR()->GetEvent(mod); + numbpoints = frecPoints->GetEntries(); + totalpoints+=numbpoints; + np[mod] = numbpoints; + //cout<<" mod = "<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 + //////////////////////////////////////////////////////////////////// + + AliITSTrackV1 trackITS(*track); + if(fresult){ delete fresult; fresult=0;} + fresult = new AliITSTrackV1(trackITS); + + AliITSTrackV1 primaryTrack(trackITS); + vgeant=(*fresult).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.; - + 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 = "<AddLast(&trackITS); - - fPtref=TMath::Abs( (trackITS).GetPt() ); - - if(fPtref>1.0) fChi2max=40.; - if(fPtref<=1.0) fChi2max=20.; - if(fPtref<0.4 ) fChi2max=100.; - if(fPtref<0.2 ) fChi2max=40.; - // if(fPtref<0.4 ) fChi2max=30.; - // if(fPtref<0.2 ) fChi2max=20.; - //if(fPtref<0.2 ) fChi2max=10.; - //if(fPtref<0.1 ) fChi2max=5.; - - - cout << "\n Pt = " << fPtref <<"\n"; //stampa - RecursiveTracking(list); - 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=(*fresult).GetLabTrack(lay); - Float_t clustZ=(*fresult).GetZclusterTrack( lay); - for(k=0; k<3; k++){ + Double_t dv=gRandom->Gaus(signdv*vr,(Float_t)sigmaDv); + Double_t zv=gRandom->Gaus(vgeant(2),(Float_t)sigmaZv); + //cout<<" Dv e Zv = "<AddLast(&trackITS); + + fPtref=TMath::Abs( (trackITS).GetPt() ); + if(fPtref>1.0) fChi2max=40.; + if(fPtref<=1.0) fChi2max=20.; + if(fPtref<0.4 ) fChi2max=100.; + if(fPtref<0.2 ) fChi2max=40.; + // if(fPtref<0.4 ) fChi2max=30.; + // if(fPtref<0.2 ) fChi2max=20.; + //if(fPtref<0.2 ) fChi2max=10.; + //if(fPtref<0.1 ) fChi2max=5.; + cout << "\n Pt = " << fPtref <<"\n"; //stampa + RecursiveTracking(list); + 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=(*fresult).GetLabTrack(lay); + Float_t clustZ=(*fresult).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; - (*fresult).Search(vecTotLabRef, labref, freq); - - - //if(freq < 6) labref=-labref; // cinque - sei - if(freq < 5) labref=-labref; // cinque - sei - (*fresult).SetLabel(labref); - - // cout<<" progressive track number = "<GetLabel(); - cout << " TPC track label = " << lab <<"\n"; // stampa - - -//propagation to vertex - - Double_t rbeam=3.; - - (*fresult).Propagation(rbeam); - - Double_t c00,c10,c11,c20,c21,c22,c30,c31,c32,c33,c40,c41,c42,c43,c44; - (*fresult).GetCElements(c00,c10,c11,c20,c21,c22,c30,c31,c32,c33,c40,c41,c42,c43,c44); + TParticle *p=(TParticle*) gAlice->Particle(lpp); + Int_t pcode=p->GetPdgCode(); + if(pcode==11) vecLabRef(k)=p->GetFirstMother(); + } // end if + itot++; vecTotLabRef(itot)=vecLabRef(k); + if(vecLabRef(k)==0. && clustZ == 0.) vecTotLabRef(itot) =-3.; + } // end for k + } // end for lay + Long_t labref; + Int_t freq; + (*fresult).Search(vecTotLabRef, labref, freq); + + //if(freq < 6) labref=-labref; // cinque - sei + if(freq < 5) labref=-labref; // cinque - sei + (*fresult).SetLabel(labref); + + // cout<<" progressive track number = "<GetLabel(); + cout << " TPC track label = " << lab <<"\n"; // stampa + //propagation to vertex + + Double_t rbeam=3.; + (*fresult).Propagation(rbeam); + Double_t c00,c10,c11,c20,c21,c22,c30,c31,c32,c33,c40,c41,c42,c43,c44; + (*fresult).GetCElements(c00, + c10,c11, + c20,c21,c22, + c30,c31,c32,c33, + c40,c41,c42,c43,c44); - Double_t pt=TMath::Abs((*fresult).GetPt()); - Double_t dr=(*fresult).GetD(); - Double_t z=(*fresult).GetZ(); - Double_t tgl=(*fresult).GetTgl(); - Double_t c=(*fresult).GetC(); - Double_t cy=c/2.; - Double_t dz=z-(tgl/cy)*TMath::ASin((*fresult).Arga(rbeam)); - dz-=vgeant(2); - - // cout<<" dr e dz alla fine = "<duepi) phivertex-=duepi; - if(phivertex<0.) phivertex+=duepi; - -////////////////////////////////////////////////////////////////////////////////////////// - - 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=(*fresult).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); - ioTrack->SetTPCLabel(lab); - - Int_t il; - for(il=0;il<6; il++){ - ioTrack->SetIdPoint(il,(*fresult).GetIdPoint(il)); - ioTrack->SetIdModule(il,(*fresult).GetIdModule(il)); - } - tracktree1.Fill(); - - //cout<<" labITS = "<0.) *(fvettid[idmodule]+idpoint)=1; //modificata angela - ioTrack->SetIdPoint(il,idpoint); - ioTrack->SetIdModule(il,idmodule); - } - - - } // end if on numOfCluster - //gObjectTable->Print(); // stampa memoria - } // end for (int j=minTr; j<=maxTr; j++) - - static Bool_t first=kTRUE; - static TFile *tfile; - - if(first) { - tfile=new TFile("itstracks.root","RECREATE"); - //cout<<"I have opened itstracks.root file "<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(fvettid) delete [] fvettid; - if(fresult) {delete fresult; fresult=0;} - + Double_t pt=TMath::Abs((*fresult).GetPt()); + Double_t dr=(*fresult).GetD(); + Double_t z=(*fresult).GetZ(); + Double_t tgl=(*fresult).GetTgl(); + Double_t c=(*fresult).GetC(); + Double_t cy=c/2.; + Double_t dz=z-(tgl/cy)*TMath::ASin((*fresult).Arga(rbeam)); + dz-=vgeant(2); + // cout<<" dr e dz alla fine = "<duepi) phivertex-=duepi; + if(phivertex<0.) phivertex+=duepi; + ///////////////////////////////////////////////////////////// + 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=(*fresult).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); + ioTrack->SetTPCLabel(lab); + Int_t il; + for(il=0;il<6; il++){ + ioTrack->SetIdPoint(il,(*fresult).GetIdPoint(il)); + ioTrack->SetIdModule(il,(*fresult).GetIdModule(il)); + } // end for il + tracktree1.Fill(); + for (il=0;il<6;il++) { + idpoint=(*fresult).GetIdPoint(il); + idmodule=(*fresult).GetIdModule(il); + //*(fvettid[idmodule]+idpoint)=1; + if(idmodule>0.) *(fvettid[idmodule]+idpoint)=1;//modificata + // angela + ioTrack->SetIdPoint(il,idpoint); + ioTrack->SetIdModule(il,idmodule); + } // end for il + } // end if on numOfCluster + //gObjectTable->Print(); // stampa memoria + } // end for (int j=minTr; j<=maxTr; j++) + static Bool_t first=kTRUE; + static TFile *tfile; + if(first) { + tfile=new TFile("itstracks.root","RECREATE"); + //cout<<"I have opened itstracks.root file "<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(fvettid) delete [] fvettid; + if(fresult) {delete fresult; fresult=0;} } - - - +//______________________________________________________________________ void AliITSTrackerV1::RecursiveTracking(TList *trackITSlist) { - -/////////////////////// This function perform the recursive tracking in ITS detectors ///////////////////// -/////////////////////// reference is a pointer to the final best track ///////////////////// -//Origin A. Badala' and G.S. Pappalardo: e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it -// The authors thank Mariana Bondila to have help them to resolve some problems. July-2000 - - //Rlayer[0]=4.; Rlayer[1]=7.; Rlayer[2]=14.9; Rlayer[3]=23.8; Rlayer[4]=39.1; Rlayer[5]=43.6; //vecchio - - ////////////////////// - Float_t sigmaphil[6], sigmazl[6]; - sigmaphil[0]=1.44e-6/(fAvrad[0]*fAvrad[0]); - sigmaphil[1]=1.44e-6/(fAvrad[1]*fAvrad[1]); - sigmaphil[2]=1.444e-5/(fAvrad[2]*fAvrad[2]); - sigmaphil[3]=1.444e-5/(fAvrad[3]*fAvrad[3]); - sigmaphil[4]=4e-6/(fAvrad[4]*fAvrad[4]); - sigmaphil[5]=4e-6/(fAvrad[5]*fAvrad[5]); - - sigmazl[0]=1e-2; - sigmazl[1]=1e-2; - sigmazl[2]=7.84e-4; - sigmazl[3]=7.84e-4; - sigmazl[4]=0.6889; - sigmazl[5]=0.6889; - /////////////////////////////////////////////////////////// - - - Int_t index; - AliITSgeom *g1 = fITS->GetITSgeom(); - AliITSRecPoint *recp; - for(index =0; indexGetSize(); index++) { - AliITSTrackV1 *trackITS = (AliITSTrackV1 *) trackITSlist->At(index); - - if((*trackITS).GetLayer()==7) fresult->SetChi2(10.223e140); - // cout <<" Layer inizio = "<<(*trackITS).GetLayer()<<"\n"; - // cout<<"fvtrack =" <<"\n"; - // cout << (*trackITS)(0) << " "<<(*trackITS)(1)<<" "<<(*trackITS)(2)<<" "<<(*trackITS)(3)<<" "<<(*trackITS)(4)<<"\n"; - // cout<< " rtrack = "<<(*trackITS).Getrtrack()<<"\n"; - // cout<< " Pt = "<<(*trackITS).GetPt()<<"\n"; - // getchar(); - Double_t chi2Now, chi2Ref; - Float_t numClustRef = fresult->GetNumClust(); - if((*trackITS).GetLayer()==1 ) { - chi2Now = trackITS->GetChi2(); - Float_t numClustNow = trackITS->GetNumClust(); - if(trackITS->GetNumClust()) chi2Now /= (Double_t )trackITS->GetNumClust(); - chi2Ref = fresult->GetChi2(); - - if(fresult->GetNumClust()) chi2Ref /= (Double_t )fresult->GetNumClust(); - //cout<<" chi2Now and chi2Ref = "< numClustRef ) {*fresult = *trackITS;} - if((numClustNow == numClustRef )&& (chi2Now < chi2Ref)) {*fresult = *trackITS;} - continue; - } - - if(trackITS->Getfnoclust()>=2) continue; - Float_t numClustNow = trackITS->GetNumClust(); - if(numClustNow) { - chi2Now = trackITS->GetChi2(); - - - if(numClustNowfresult->GetChi2()) continue; - - //cout<<" chi2Now = "< 1.0 && chi2Now > 30.) continue; - if((fPtref >= 0.6 && fPtref<=1.0) && chi2Now > 40.) continue; - // if((fPtref <= 0.6 && fPtref>0.2)&& chi2Now > 40.) continue; - // if(fPtref <= 0.2 && chi2Now > 8.) continue; - if((fPtref <= 0.6 && fPtref>0.2)&& chi2Now > 30.) continue; - if(fPtref <= 0.2 && chi2Now > 7.) continue; - - - ///////////////////////////// - - } - - Int_t layerInit = (*trackITS).GetLayer(); - Int_t layernew = layerInit - 2; // -1 for new layer, -1 for matrix index - - TList listoftrack; - 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 = "< fNlad[layerfin-1]) ladp=1; - detp=detinters+1; - detm=detinters-1; - Int_t idetot=1; - /* - toucLad(0)=ladinters; toucLad(1)=ladm; toucLad(2)=ladp; - toucLad(3)=ladinters; toucLad(4)=ladm; toucLad(5)=ladp; - toucLad(6)=ladinters; toucLad(7)=ladm; toucLad(8)=ladp; - toucDet(0)=detinters; toucDet(1)=detinters; toucDet(2)=detinters; - if(detm > 0 && detp <= fNdet[layerfin-1]) { - idetot=9; - toucDet(3)=detm; toucDet(4)=detm; toucDet(5)=detm; - toucDet(6)=detp; toucDet(7)=detp; toucDet(8)=detp; - } - - if(detm > 0 && detp > fNdet[layerfin-1]) { - idetot=6; - toucDet(3)=detm; toucDet(4)=detm; toucDet(5)=detm; - } - - if(detm <= 0 && detp <= fNdet[layerfin-1]) { - idetot=6; - toucDet(3)=detp; toucDet(4)=detp; toucDet(5)=detp; - } - */ + // This function perform the recursive tracking in ITS detectors + // reference is a pointer to the final best track + // Origin A. Badala' and G.S. Pappalardo: + // e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it + // The authors thank Mariana Bondila to have help them to resolve some + // problems. July-2000 + + //Rlayer[0]=4.; Rlayer[1]=7.; Rlayer[2]=14.9; + // Rlayer[3]=23.8; Rlayer[4]=39.1; Rlayer[5]=43.6; //vecchio + + ////////////////////// + Float_t sigmaphil[6], sigmazl[6]; + sigmaphil[0]=1.44e-6/(fAvrad[0]*fAvrad[0]); + sigmaphil[1]=1.44e-6/(fAvrad[1]*fAvrad[1]); + sigmaphil[2]=1.444e-5/(fAvrad[2]*fAvrad[2]); + sigmaphil[3]=1.444e-5/(fAvrad[3]*fAvrad[3]); + sigmaphil[4]=4e-6/(fAvrad[4]*fAvrad[4]); + sigmaphil[5]=4e-6/(fAvrad[5]*fAvrad[5]); + sigmazl[0]=1e-2; + sigmazl[1]=1e-2; + sigmazl[2]=7.84e-4; + sigmazl[3]=7.84e-4; + sigmazl[4]=0.6889; + sigmazl[5]=0.6889; + /////////////////////////////////////////////////////////// + Int_t index; + AliITSgeom *g1 = fITS->GetITSgeom(); + AliITSRecPoint *recp; + for(index =0; indexGetSize(); index++) { + AliITSTrackV1 *trackITS = (AliITSTrackV1 *) trackITSlist->At(index); + if((*trackITS).GetLayer()==7) fresult->SetChi2(10.223e140); + // cout <<" Layer inizio = "<<(*trackITS).GetLayer()<<"\n"; + // cout<<"fvtrack =" <<"\n"; + // cout << (*trackITS)(0) << " "<<(*trackITS)(1)<<" " + // <<(*trackITS)(2)<<" "<<(*trackITS)(3)<<" " + // <<(*trackITS)(4)<<"\n"; + // cout<< " rtrack = "<<(*trackITS).Getrtrack()<<"\n"; + // cout<< " Pt = "<<(*trackITS).GetPt()<<"\n"; + // getchar(); + Double_t chi2Now, chi2Ref; + Float_t numClustRef = fresult->GetNumClust(); + if((*trackITS).GetLayer()==1 ) { + chi2Now = trackITS->GetChi2(); + Float_t numClustNow = trackITS->GetNumClust(); + if(trackITS->GetNumClust()) + chi2Now /= (Double_t)trackITS->GetNumClust(); + chi2Ref = fresult->GetChi2(); + if(fresult->GetNumClust()) + chi2Ref /= (Double_t)fresult->GetNumClust(); + //cout<<" chi2Now and chi2Ref = "< numClustRef ) {*fresult = *trackITS;} + if((numClustNow == numClustRef )&& + (chi2Now < chi2Ref)) { + *fresult = *trackITS; + } // end if + continue; + } // end if + + if(trackITS->Getfnoclust()>=2) continue; + Float_t numClustNow = trackITS->GetNumClust(); + if(numClustNow) { + chi2Now = trackITS->GetChi2(); + + if(numClustNowfresult->GetChi2()) continue; + //cout<<" chi2Now = "< 1.0 && chi2Now > 30.) continue; + if((fPtref >= 0.6 && fPtref<=1.0) && chi2Now > 40.) continue; + // if((fPtref <= 0.6 && fPtref>0.2)&& chi2Now > 40.) continue; + // if(fPtref <= 0.2 && chi2Now > 8.) continue; + if((fPtref <= 0.6 && fPtref>0.2)&& chi2Now > 30.) continue; + if(fPtref <= 0.2 && chi2Now > 7.) continue; + ///////////////////////////// + } // end if + + Int_t layerInit = (*trackITS).GetLayer(); + Int_t layernew = layerInit - 2;// -1 for new layer, -1 for matrix index + TList listoftrack; + 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 = "< fNlad[layerfin-1]) ladp=1; + detp=detinters+1; + detm=detinters-1; + Int_t idetot=1; + /* + toucLad(0)=ladinters; toucLad(1)=ladm; toucLad(2)=ladp; + toucLad(3)=ladinters; toucLad(4)=ladm; toucLad(5)=ladp; + toucLad(6)=ladinters; toucLad(7)=ladm; toucLad(8)=ladp; + toucDet(0)=detinters; toucDet(1)=detinters; toucDet(2)=detinters; + if(detm > 0 && detp <= fNdet[layerfin-1]) { + idetot=9; + toucDet(3)=detm; toucDet(4)=detm; toucDet(5)=detm; + toucDet(6)=detp; toucDet(7)=detp; toucDet(8)=detp; + } // end if + if(detm > 0 && detp > fNdet[layerfin-1]) { + idetot=6; + toucDet(3)=detm; toucDet(4)=detm; toucDet(5)=detm; + } // end if + if(detm <= 0 && detp <= fNdet[layerfin-1]) { + idetot=6; + toucDet(3)=detp; toucDet(4)=detp; toucDet(5)=detp; + } // end if + */ 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 transformed in a method -// these values could be modified - Float_t pigre=TMath::Pi(); - Float_t rangephi=5., rangez=5.; - if(layerfin==1 || layerfin ==2){ - rangephi=40.*epsphi*TMath::Sqrt(sigmaphil[layerfin-1]+(*trackITS).GetSigmaphi()); - rangez = 40.*epsz*TMath::Sqrt(sigmazl[layerfin-1]+(*trackITS).GetSigmaZ()); - } - if(layerfin==3 || layerfin ==4){ - //rangephi=30.*fepsphi*TMath::Sqrt(sigmaphil[layerfin-1]+(*trackITS).GetSigmaphi()); - //rangez = 40.*fepsz*TMath::Sqrt(sigmazl[layerfin-1]+(*trackITS).GetSigmaZ()); - rangephi=40.*epsphi*TMath::Sqrt(sigmaphil[layerfin-1]+(*trackITS).GetSigmaphi()); - rangez = 50.*epsz*TMath::Sqrt(sigmazl[layerfin-1]+(*trackITS).GetSigmaZ()); - } - if(layerfin==5 || layerfin ==6){ - rangephi=20.*epsphi*TMath::Sqrt(sigmaphil[layerfin-1]+(*trackITS).GetSigmaphi()); - rangez =5.*epsz*TMath::Sqrt(sigmazl[layerfin-1]+(*trackITS).GetSigmaZ()); - } - Float_t phinters, zinters; - phinters=(*trackITS).Getphi(); - zinters=(*trackITS).GetZ(); - Float_t distz; - 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]); - phicp=phinters; - if(phinters>fphimin[layerfin-1][ladp]) phicp=phinters-2.*pigre; - distphip=TMath::Abs(phicp-fphimin[layerfin-1][ladp]); - Int_t flagzmin=0; - Int_t flagzmax=0; - idetot=1; - toucLad(0)=ladinters; toucDet(0)=detinters; - if(detm>0) distz=TMath::Abs(zinters-fzmax[layerfin-1][detm-1]); - if(detm>0 && rangez>=distz){ - flagzmin=1; - idetot++; toucLad(idetot-1)=ladinters; toucDet(idetot-1)=detm; - if(rangephi>=distphim){ - idetot++; toucLad(idetot-1)=ladm; toucDet(idetot-1)=detinters; - idetot++; toucLad(idetot-1)=ladm; toucDet(idetot-1)=detm; - } - if(rangephi>=distphip){ - idetot++; toucLad(idetot-1)=ladp; toucDet(idetot-1)=detinters; - idetot++; toucLad(idetot-1)=ladp; toucDet(idetot-1)=detm; - } - } //end detm>0.... - if(detp<=fNdet[layerfin-1]) distz=TMath::Abs(zinters-fzmin[layerfin-1][detp-1]); - if(detp<=fNdet[layerfin-1] && rangez>=distz){ - flagzmax=1; - idetot++; toucLad(idetot-1)=ladinters; toucDet(idetot-1)=detp; - if(rangephi>=distphim){ - idetot++; toucLad(idetot-1)=ladm; toucDet(idetot-1)=detp; - if(flagzmin == 0) {idetot++; toucLad(idetot-1)=ladm; toucDet(idetot-1)=detinters;} - } - if(rangephi>=distphip){ - idetot++; toucLad(idetot-1)=ladp; toucDet(idetot-1)=detp; - if(flagzmin == 0) {idetot++; toucLad(idetot-1)=ladp; toucDet(idetot-1)=detinters;} - } - } //end detm=distphim){idetot++; toucLad(idetot-1)=ladm; toucDet(idetot-1)=detinters;} - if(rangephi>=distphip){idetot++; toucLad(idetot-1)=ladp; toucDet(idetot-1)=detinters;} - } - -/////////////////////////////////////////////////////////////////////////////////////////////////// - - - Int_t iriv; - for (iriv=0; irivGetModuleIndex(lycur,toucLad(iriv),toucDet(iriv)); - //mod ott - 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++; - } - } - - ind=-1; - - for(;;) { - ind++; - if(indlist[ind] < 0) recp=0; - else recp = (AliITSRecPoint*)frecPoints->UncheckedAt(indlist[ind]); - - if((!recp) ) break; - */ -/////////////////////////new - Int_t indnew; - for(indnew=0; indnewUncheckedAt(indnew); - else + if(fPtref<0.2) {epsphi=3.; epsz=3.;} + // nuova definizione idetot e toucLad e toucDet to be + // transformed in a method + // these values could be modified + Float_t pigre=TMath::Pi(); + Float_t rangephi=5., rangez=5.; + if(layerfin==1 || layerfin ==2){ + rangephi=40.*epsphi*TMath::Sqrt(sigmaphil[layerfin-1]+ + (*trackITS).GetSigmaphi()); + rangez = 40.*epsz*TMath::Sqrt(sigmazl[layerfin-1]+ + (*trackITS).GetSigmaZ()); + } // end if + if(layerfin==3 || layerfin ==4){ + //rangephi=30.*fepsphi*TMath::Sqrt(sigmaphil[layerfin-1]+ + // (*trackITS).GetSigmaphi()); + //rangez = 40.*fepsz*TMath::Sqrt(sigmazl[layerfin-1]+ + // (*trackITS).GetSigmaZ()); + rangephi=40.*epsphi*TMath::Sqrt(sigmaphil[layerfin-1]+ + (*trackITS).GetSigmaphi()); + rangez = 50.*epsz*TMath::Sqrt(sigmazl[layerfin-1]+ + (*trackITS).GetSigmaZ()); + } // end if + if(layerfin==5 || layerfin ==6){ + rangephi=20.*epsphi*TMath::Sqrt(sigmaphil[layerfin-1]+ + (*trackITS).GetSigmaphi()); + rangez =5.*epsz*TMath::Sqrt(sigmazl[layerfin-1]+ + (*trackITS).GetSigmaZ()); + } // end if + Float_t phinters, zinters; + phinters=(*trackITS).Getphi(); + zinters=(*trackITS).GetZ(); + 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]); + phicp=phinters; + if(phinters>fphimin[layerfin-1][ladp]) phicp=phinters-2.*pigre; + distphip=TMath::Abs(phicp-fphimin[layerfin-1][ladp]); + Int_t flagzmin=0; + Int_t flagzmax=0; + idetot=1; + toucLad(0)=ladinters; toucDet(0)=detinters; + if(detm>0) distz=TMath::Abs(zinters-fzmax[layerfin-1][detm-1]); + if(detm>0 && rangez>=distz){ + flagzmin=1; + idetot++; toucLad(idetot-1)=ladinters; toucDet(idetot-1)=detm; + if(rangephi>=distphim){ + idetot++; + toucLad(idetot-1)=ladm; + toucDet(idetot-1)=detinters; + idetot++; + toucLad(idetot-1)=ladm; + toucDet(idetot-1)=detm; + } // end if + if(rangephi>=distphip){ + idetot++; + toucLad(idetot-1)=ladp; + toucDet(idetot-1)=detinters; + idetot++; + toucLad(idetot-1)=ladp; + toucDet(idetot-1)=detm; + } // end if + } //end detm>0.... + if(detp<=fNdet[layerfin-1]) + distz=TMath::Abs(zinters-fzmin[layerfin-1][detp-1]); + if(detp<=fNdet[layerfin-1] && rangez>=distz){ + flagzmax=1; + idetot++; toucLad(idetot-1)=ladinters; toucDet(idetot-1)=detp; + if(rangephi>=distphim){ + idetot++; toucLad(idetot-1)=ladm; toucDet(idetot-1)=detp; + if(flagzmin == 0) { + idetot++; + toucLad(idetot-1)=ladm; + toucDet(idetot-1)=detinters; + } // end if + } // end if + if(rangephi>=distphip){ + idetot++; + toucLad(idetot-1)=ladp; + toucDet(idetot-1)=detp; + if(flagzmin == 0) { + idetot++; + toucLad(idetot-1)=ladp; + toucDet(idetot-1)=detinters; + } // end if + } // end if + } //end detm=distphim){ + idetot++; + toucLad(idetot-1)=ladm; + toucDet(idetot-1)=detinters; + } // end if + if(rangephi>=distphip){ + idetot++; + toucLad(idetot-1)=ladp; + toucDet(idetot-1)=detinters; + } // end if + } // end if + //////////////////////////////////////////////////////////// + Int_t iriv; + for (iriv=0; irivGetModuleIndex(lycur,(Int_t)toucLad(iriv), + (Int_t)toucDet(iriv)); + //mod ott + 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; indnewUncheckedAt(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(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(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 - - // cout<<" layer = "<6.) cluster(1)=cluster(1)+(2.*TMath::Pi()); - if(cluster(1)>6. && (*trackITS).Getphi()<6.) cluster(1)=cluster(1)-(2.*TMath::Pi()); - if(TMath::Abs(cluster(1)-(*trackITS).Getphi()) > sigmatotphi) continue; - // cout<<" supero sigmaphi \n"; - AliITSTrackV1 *newTrack = new AliITSTrackV1((*trackITS)); - (*newTrack).SetLayer((*trackITS).GetLayer()-1); - - if (TMath::Abs(rTrack-cluster(0))/rTrack>1e-6) - (*newTrack).Correct(Double_t(cluster(0))); - //cout<<" cluster(2) e (*newTrack).GetZ() = "< sigmatotz){ - delete newTrack; - continue;} - - Double_t sigmanew[2]; - sigmanew[0]= sigmaphi; - sigmanew[1]=sigma[1]; - - Double_t m[2]; - m[0]=cluster(1); - m[1]=cluster(2); -// Double_t chi2pred=newTrack->GetPredChi2(m,sigmanew); - // cout<<" chi2pred = "<fChi2max) continue; //aggiunto il 30-7-2001 - - if(iriv == 0) flaghit=1; - - (*newTrack).AddMS(frl); // add the multiple scattering 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); - } - //fTimerKalman->Stop(); // timer - - - (*newTrack).PutCluster(layernew, vecclust); - newTrack->AddClustInTrack(); - - listoftrack.AddLast(newTrack); - - } // end of for(;;) on rec points - - // delete [] indlist; //mod ott - - } // end of for on detectors - - }//end if(outinters==0) - - if(flaghit==0 || outinters==-2) { - AliITSTrackV1 *newTrack = new AliITSTrackV1(*trackITS); - (*newTrack).Setfnoclust(); - (*newTrack).SetLayer((*trackITS).GetLayer()-1); - (*newTrack).AddMS(frl); // add the multiple scattering matrix to the covariance matrix - (*newTrack).AddEL(frl,1.,0); - - listoftrack.AddLast(newTrack); - } - - - //gObjectTable->Print(); // stampa memoria + //////////////////////////////////////////////////// + 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(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(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 + // cout<<" layer = "<6.) + cluster(1)=cluster(1)+(2.*TMath::Pi()); + if(cluster(1)>6. && (*trackITS).Getphi()<6.) + cluster(1)=cluster(1)-(2.*TMath::Pi()); + if(TMath::Abs(cluster(1)-(*trackITS).Getphi())>sigmatotphi) + continue; + // cout<<" supero sigmaphi \n"; + AliITSTrackV1 *newTrack = new AliITSTrackV1((*trackITS)); + (*newTrack).SetLayer((*trackITS).GetLayer()-1); + if (TMath::Abs(rTrack-cluster(0))/rTrack>1e-6) + (*newTrack).Correct(Double_t(cluster(0))); + //cout<<" cluster(2) e(*newTrack).GetZ()="< sigmatotz){ + delete newTrack; + continue; + } // end if + Double_t sigmanew[2]; + sigmanew[0]= sigmaphi; + sigmanew[1]=sigma[1]; + Double_t m[2]; + m[0]=cluster(1); + m[1]=cluster(2); + // Double_t chi2pred=newTrack->GetPredChi2(m,sigmanew); + // cout<<" chi2pred = "<fChi2max) continue; //aggiunto il 30-7-2001 + if(iriv == 0) flaghit=1; + (*newTrack).AddMS(frl); // add the multiple scattering + //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) + + if(flaghit==0 || outinters==-2) { + AliITSTrackV1 *newTrack = new AliITSTrackV1(*trackITS); + (*newTrack).Setfnoclust(); + (*newTrack).SetLayer((*trackITS).GetLayer()-1); + (*newTrack).AddMS(frl); // add the multiple scattering matrix + // to the covariance matrix + (*newTrack).AddEL(frl,1.,0); + listoftrack.AddLast(newTrack); + } // end if + + //gObjectTable->Print(); // stampa memoria - RecursiveTracking(&listoftrack); - listoftrack.Delete(); - } // end of for on tracks - - //gObjectTable->Print(); // stampa memoria - -} + RecursiveTracking(&listoftrack); + listoftrack.Delete(); + } // end of for on tracks (index) -Int_t AliITSTrackerV1::Intersection(AliITSTrackV1 &track, Int_t layer, Int_t &ladder, Int_t &detector) { -//Origin A. Badala' and G.S. Pappalardo: e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it -// Found the intersection and the detector - - Double_t rk=fAvrad[layer-1]; - if(track.DoNotCross(rk)){ /*cout<< " Do not cross \n";*/ return -1;} - track.Propagation(rk); - Double_t zinters=track.GetZ(); - Double_t phinters=track.Getphi(); - //cout<<"zinters = "< fzmin[layer-1][iD-1] && zinters <= fzmax[layer-1][iD-1]) { - if(iz>1) { - cout<< " Errore su iz in Intersection \n"; getchar(); - } - else { - listDet(iz)= iD; distZCenter(iz)=TMath::Abs(zinters-det(2)); iz++; - } - } - } - - if(iz==0) {/* cout<< " No detector along Z \n";*/ return -2;} - detector=Int_t (listDet(0)); - if(iz>1 && (distZCenter(0)>distZCenter(1))) detector=Int_t (listDet(1)); - - TVector listLad(2); - TVector distPhiCenter(2); - Int_t ip=0; - Double_t pigre=TMath::Pi(); - - Int_t iLd; - for(iLd = 1; iLd<= fNlad[layer-1]; iLd++) { - Double_t phimin=fphimin[layer-1][iLd-1]; - Double_t phimax=fphimax[layer-1][iLd-1]; - Double_t phidet=fphidet[layer-1][iLd-1]; - Double_t phiconfr=phinters; - if(phimin>phimax) { - //if(phimin <5.5) {cout<<" Error in Intersection for phi \n"; getchar();} - phimin-=(2.*pigre); - if(phinters>(1.5*pigre)) phiconfr=phinters-(2.*pigre); - if(phidet>(1.5*pigre)) phidet-=(2.*pigre); - } - if(phiconfr>phimin && phiconfr<= phimax) { - if(ip>1) { - cout<< " Errore su ip in Intersection \n"; getchar(); - } - else { - listLad(ip)= iLd; distPhiCenter(ip)=TMath::Abs(phiconfr-phidet); ip++; - } - } - } - if(ip==0) { cout<< " No detector along phi \n"; getchar();} - ladder=Int_t (listLad(0)); - if(ip>1 && (distPhiCenter(0)>distPhiCenter(1))) ladder=Int_t (listLad(1)); - - return 0; + //gObjectTable->Print(); // stampa memoria } - -void AliITSTrackerV1::KalmanFilter(AliITSTrackV1 *newTrack,TVector &cluster,Double_t sigma[2]){ -//Origin A. Badala' and G.S. Pappalardo: e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it -// Kalman filter without vertex constraint - - - ////////////////////////////// Evaluation of the measurement vector ///////////////////////////////////// - - Double_t m[2]; - Double_t rk,phik,zk; - rk=cluster(0); phik=cluster(1); zk=cluster(2); - m[0]=phik; m[1]=zk; - - ///////////////////////////////////// Evaluation of the error matrix V /////////////////////////////// - - Double_t v00=sigma[0]; - Double_t v11=sigma[1]; - - /////////////////////////////////////////////////////////////////////////////////////////// - - - Double_t cin00,cin10,cin20,cin30,cin40,cin11,cin21,cin31,cin41,cin22,cin32,cin42,cin33,cin43,cin44; - - newTrack->GetCElements(cin00,cin10,cin11,cin20,cin21,cin22,cin30,cin31,cin32,cin33,cin40, - cin41,cin42,cin43,cin44); //get C matrix - - Double_t rold00=cin00+v00; - Double_t rold10=cin10; - Double_t rold11=cin11+v11; - -//////////////////////////////////// R matrix inversion /////////////////////////////////////////////// - - Double_t det=rold00*rold11-rold10*rold10; - Double_t r00=rold11/det; - Double_t r10=-rold10/det; - Double_t r11=rold00/det; - -//////////////////////////////////////////////////////////////////////////////////////////////////////// - - Double_t k00=cin00*r00+cin10*r10; - Double_t k01=cin00*r10+cin10*r11; - Double_t k10=cin10*r00+cin11*r10; - Double_t k11=cin10*r10+cin11*r11; - Double_t k20=cin20*r00+cin21*r10; - Double_t k21=cin20*r10+cin21*r11; - Double_t k30=cin30*r00+cin31*r10; - Double_t k31=cin30*r10+cin31*r11; - Double_t k40=cin40*r00+cin41*r10; - Double_t k41=cin40*r10+cin41*r11; - - Double_t x0,x1,x2,x3,x4; - newTrack->GetXElements(x0,x1,x2,x3,x4); // get the state vector - - Double_t savex0=x0, savex1=x1; - - x0+=k00*(m[0]-savex0)+k01*(m[1]-savex1); - x1+=k10*(m[0]-savex0)+k11*(m[1]-savex1); - x2+=k20*(m[0]-savex0)+k21*(m[1]-savex1); - x3+=k30*(m[0]-savex0)+k31*(m[1]-savex1); - x4+=k40*(m[0]-savex0)+k41*(m[1]-savex1); - - Double_t c00,c10,c20,c30,c40,c11,c21,c31,c41,c22,c32,c42,c33,c43,c44; - - c00=cin00-k00*cin00-k01*cin10; - c10=cin10-k00*cin10-k01*cin11; - c20=cin20-k00*cin20-k01*cin21; - c30=cin30-k00*cin30-k01*cin31; - c40=cin40-k00*cin40-k01*cin41; - - c11=cin11-k10*cin10-k11*cin11; - c21=cin21-k10*cin20-k11*cin21; - c31=cin31-k10*cin30-k11*cin31; - c41=cin41-k10*cin40-k11*cin41; - - c22=cin22-k20*cin20-k21*cin21; - c32=cin32-k20*cin30-k21*cin31; - c42=cin42-k20*cin40-k21*cin41; - - c33=cin33-k30*cin30-k31*cin31; - c43=cin43-k30*cin40-k31*cin41; - - c44=cin44-k40*cin40-k41*cin41; - - newTrack->PutXElements(x0,x1,x2,x3,x4); // put the new state vector - - newTrack->PutCElements(c00,c10,c11,c20,c21,c22,c30,c31,c32,c33,c40,c41,c42,c43,c44); // put in track the - // new cov matrix - Double_t vmcold00=v00-c00; - Double_t vmcold10=-c10; - Double_t vmcold11=v11-c11; - -///////////////////////////////////// Matrix vmc inversion //////////////////////////////////////////////// - - det=vmcold00*vmcold11-vmcold10*vmcold10; - Double_t vmc00=vmcold11/det; - Double_t vmc10=-vmcold10/det; - Double_t vmc11=vmcold00/det; - -//////////////////////////////////////////////////////////////////////////////////////////////////////////// - - Double_t chi2=(m[0]-x0)*( vmc00*(m[0]-x0) + 2.*vmc10*(m[1]-x1) ) + - (m[1]-x1)*vmc11*(m[1]-x1); - +//______________________________________________________________________ +Int_t AliITSTrackerV1::Intersection(AliITSTrackV1 &track,Int_t layer, + Int_t &ladder,Int_t &detector) { + // Origin A. Badala' and G.S. Pappalardo + // e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it + // Found the intersection and the detector + + Double_t rk=fAvrad[layer-1]; + if(track.DoNotCross(rk)){ /*cout<< " Do not cross \n";*/ return -1;} + track.Propagation(rk); + Double_t zinters=track.GetZ(); + Double_t phinters=track.Getphi(); + //cout<<"zinters = "< fzmin[layer-1][iD-1] && zinters <= fzmax[layer-1][iD-1]) { + if(iz>1) { + cout<< " Errore su iz in Intersection \n"; + getchar(); + }else { + listDet(iz)= iD; distZCenter(iz)=TMath::Abs(zinters-det(2)); + iz++; + } // end if + } // end if + } // end for iD + + if(iz==0) {/* cout<< " No detector along Z \n";*/ return -2;} + detector=Int_t (listDet(0)); + if(iz>1 && (distZCenter(0)>distZCenter(1))) detector=Int_t (listDet(1)); + + TVector listLad(2); + TVector distPhiCenter(2); + Int_t ip=0; + Double_t pigre=TMath::Pi(); + Int_t iLd; + for(iLd = 1; iLd<= fNlad[layer-1]; iLd++) { + Double_t phimin=fphimin[layer-1][iLd-1]; + Double_t phimax=fphimax[layer-1][iLd-1]; + Double_t phidet=fphidet[layer-1][iLd-1]; + Double_t phiconfr=phinters; + if(phimin>phimax) { + //if(phimin <5.5) {cout<<" Error in Intersection for phi \n"; + // getchar();} + phimin-=(2.*pigre); + if(phinters>(1.5*pigre)) phiconfr=phinters-(2.*pigre); + if(phidet>(1.5*pigre)) phidet-=(2.*pigre); + } // end if + if(phiconfr>phimin && phiconfr<= phimax) { + if(ip>1) { + cout<< " Errore su ip in Intersection \n"; getchar(); + }else { + listLad(ip)= iLd; + distPhiCenter(ip)=TMath::Abs(phiconfr-phidet); ip++; + } // end if + } // end if + } // end for iLd + if(ip==0) { cout<< " No detector along phi \n"; getchar();} + ladder=Int_t (listLad(0)); + if(ip>1 && (distPhiCenter(0)>distPhiCenter(1))) ladder=Int_t (listLad(1)); + return 0; +} +//______________________________________________________________________ +void AliITSTrackerV1::KalmanFilter(AliITSTrackV1 *newTrack,TVector &cluster, + Double_t sigma[2]){ + //Origin A. Badala' and G.S. Pappalardo: + // e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it + // Kalman filter without vertex constraint + ////// Evaluation of the measurement vector //////////////////////// + Double_t m[2]; + Double_t rk,phik,zk; + rk=cluster(0); phik=cluster(1); zk=cluster(2); + m[0]=phik; m[1]=zk; + //////////////////////// Evaluation of the error matrix V ///////// + Double_t v00=sigma[0]; + Double_t v11=sigma[1]; + //////////////////////////////////////////////////////////////////// + Double_t cin00,cin10,cin20,cin30,cin40,cin11,cin21,cin31,cin41,cin22, + cin32,cin42,cin33,cin43,cin44; + + newTrack->GetCElements(cin00, + cin10,cin11, + cin20,cin21,cin22, + cin30,cin31,cin32,cin33, + cin40,cin41,cin42,cin43,cin44); //get C matrix + Double_t rold00=cin00+v00; + Double_t rold10=cin10; + Double_t rold11=cin11+v11; + ////////////////////// R matrix inversion ///////////////////////// + Double_t det=rold00*rold11-rold10*rold10; + Double_t r00=rold11/det; + Double_t r10=-rold10/det; + Double_t r11=rold00/det; + //////////////////////////////////////////////////////////////////// + Double_t k00=cin00*r00+cin10*r10; + Double_t k01=cin00*r10+cin10*r11; + Double_t k10=cin10*r00+cin11*r10; + Double_t k11=cin10*r10+cin11*r11; + Double_t k20=cin20*r00+cin21*r10; + Double_t k21=cin20*r10+cin21*r11; + Double_t k30=cin30*r00+cin31*r10; + Double_t k31=cin30*r10+cin31*r11; + Double_t k40=cin40*r00+cin41*r10; + Double_t k41=cin40*r10+cin41*r11; + Double_t x0,x1,x2,x3,x4; + newTrack->GetXElements(x0,x1,x2,x3,x4); // get the state vector + Double_t savex0=x0, savex1=x1; + x0+=k00*(m[0]-savex0)+k01*(m[1]-savex1); + x1+=k10*(m[0]-savex0)+k11*(m[1]-savex1); + x2+=k20*(m[0]-savex0)+k21*(m[1]-savex1); + x3+=k30*(m[0]-savex0)+k31*(m[1]-savex1); + x4+=k40*(m[0]-savex0)+k41*(m[1]-savex1); + Double_t c00,c10,c20,c30,c40,c11,c21,c31,c41,c22,c32,c42,c33,c43,c44; + c00=cin00-k00*cin00-k01*cin10; + c10=cin10-k00*cin10-k01*cin11; + c20=cin20-k00*cin20-k01*cin21; + c30=cin30-k00*cin30-k01*cin31; + c40=cin40-k00*cin40-k01*cin41; + c11=cin11-k10*cin10-k11*cin11; + c21=cin21-k10*cin20-k11*cin21; + c31=cin31-k10*cin30-k11*cin31; + c41=cin41-k10*cin40-k11*cin41; + c22=cin22-k20*cin20-k21*cin21; + c32=cin32-k20*cin30-k21*cin31; + c42=cin42-k20*cin40-k21*cin41; + c33=cin33-k30*cin30-k31*cin31; + c43=cin43-k30*cin40-k31*cin41; + c44=cin44-k40*cin40-k41*cin41; + newTrack->PutXElements(x0,x1,x2,x3,x4); // put the new state vector + newTrack->PutCElements(c00, + c10,c11, + c20,c21,c22, + c30,c31,c32,c33, + c40,c41,c42,c43,c44); // put in track the + // new cov matrix + Double_t vmcold00=v00-c00; + Double_t vmcold10=-c10; + Double_t vmcold11=v11-c11; + ////////////////////// Matrix vmc inversion /////////////////////// + det=vmcold00*vmcold11-vmcold10*vmcold10; + Double_t vmc00=vmcold11/det; + Double_t vmc10=-vmcold10/det; + Double_t vmc11=vmcold00/det; + //////////////////////////////////////////////////////////////////// + Double_t chi2=(m[0]-x0)*( vmc00*(m[0]-x0) + 2.*vmc10*(m[1]-x1) ) + + (m[1]-x1)*vmc11*(m[1]-x1); newTrack->SetChi2(newTrack->GetChi2()+chi2); - -} - - -//void AliITSTrackerV1::KalmanFilterVert(AliITSTrackV1 *newTrack,TVector &cluster,Double_t sigma[2]){ -void AliITSTrackerV1::KalmanFilterVert(AliITSTrackV1 *newTrack,TVector &cluster,Double_t sigma[2]/*, Double_t chi2pred*/){ -//Origin A. Badala' and G.S. Pappalardo: e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it -// Kalman filter with vertex constraint - - ////////////////////////////// Evaluation of the measurement vector m /////////////// - - Double_t m[4]; - Double_t rk,phik,zk; - rk=cluster(0); phik=cluster(1); zk=cluster(2); - m[0]=phik; m[1]=zk; - - Double_t cc=(*newTrack).GetC(); - Double_t zv=(*newTrack).GetZv(); - Double_t dv=(*newTrack).GetDv(); - Double_t cy=cc/2.; - Double_t tgl= (zk-zv)*cy/TMath::ASin(cy*rk); - m[2]=dv; m[3]=tgl; - - ///////////////////////////////////// Evaluation of the error matrix V ////////////// - Int_t layer=newTrack->GetLayer(); - Double_t v00=sigma[0]; - Double_t v11=sigma[1]; - Double_t v31=sigma[1]/rk; - Double_t sigmaDv=newTrack->GetsigmaDv(); - Double_t v22=sigmaDv*sigmaDv + newTrack->Getd2(layer-1); - Double_t v32=newTrack->Getdtgl(layer-1); - Double_t sigmaZv=newTrack->GetsigmaZv(); - Double_t v33=(sigma[1]+sigmaZv*sigmaZv)/(rk*rk) + newTrack->Gettgl2(layer-1); - /////////////////////////////////////////////////////////////////////////////////////// - - Double_t cin00,cin10,cin11,cin20,cin21,cin22,cin30,cin31,cin32,cin33,cin40,cin41,cin42,cin43,cin44; - - newTrack->GetCElements(cin00,cin10,cin11,cin20,cin21,cin22,cin30,cin31,cin32,cin33,cin40, - cin41,cin42,cin43,cin44); //get C matrix - - Double_t r[4][4]; - r[0][0]=cin00+v00; - r[1][0]=cin10; - r[2][0]=cin20; - r[3][0]=cin30; - r[1][1]=cin11+v11; - r[2][1]=cin21; - r[3][1]=cin31+sigma[1]/rk; - r[2][2]=cin22+sigmaDv*sigmaDv+newTrack->Getd2(layer-1); - r[3][2]=cin32+newTrack->Getdtgl(layer-1); - r[3][3]=cin33+(sigma[1]+sigmaZv*sigmaZv)/(rk*rk) + newTrack->Gettgl2(layer-1); - - r[0][1]=r[1][0]; r[0][2]=r[2][0]; r[0][3]=r[3][0]; r[1][2]=r[2][1]; r[1][3]=r[3][1]; - r[2][3]=r[3][2]; - -///////////////////// Matrix R inversion //////////////////////////////////////////// - - const Int_t kn=4; - Double_t big, hold; - Double_t d=1.; - Int_t ll[kn],mm[kn]; - - Int_t i,j,k; - - for(k=0; k k) { - for(i=0; i k ) { - for(j=0; j=0; k--) { - i=ll[k]; - if(i > k) { - for (j=0; j k) { - for (i=0; iGetXElements(x0,x1,x2,x3,x4); // get the state vector - - Double_t savex0=x0, savex1=x1, savex2=x2, savex3=x3; - - x0+=k00*(m[0]-savex0)+k01*(m[1]-savex1)+k02*(m[2]-savex2)+ - k03*(m[3]-savex3); - x1+=k10*(m[0]-savex0)+k11*(m[1]-savex1)+k12*(m[2]-savex2)+ - k13*(m[3]-savex3); - x2+=k20*(m[0]-savex0)+k21*(m[1]-savex1)+k22*(m[2]-savex2)+ - k23*(m[3]-savex3); - x3+=k30*(m[0]-savex0)+k31*(m[1]-savex1)+k32*(m[2]-savex2)+ - k33*(m[3]-savex3); - x4+=k40*(m[0]-savex0)+k41*(m[1]-savex1)+k42*(m[2]-savex2)+ - k43*(m[3]-savex3); - - Double_t c00,c10,c20,c30,c40,c11,c21,c31,c41,c22,c32,c42,c33,c43,c44; - - c00=cin00-k00*cin00-k01*cin10-k02*cin20-k03*cin30; - c10=cin10-k00*cin10-k01*cin11-k02*cin21-k03*cin31; - c20=cin20-k00*cin20-k01*cin21-k02*cin22-k03*cin32; - c30=cin30-k00*cin30-k01*cin31-k02*cin32-k03*cin33; - c40=cin40-k00*cin40-k01*cin41-k02*cin42-k03*cin43; - - c11=cin11-k10*cin10-k11*cin11-k12*cin21-k13*cin31; - c21=cin21-k10*cin20-k11*cin21-k12*cin22-k13*cin32; - c31=cin31-k10*cin30-k11*cin31-k12*cin32-k13*cin33; - c41=cin41-k10*cin40-k11*cin41-k12*cin42-k13*cin43; - - c22=cin22-k20*cin20-k21*cin21-k22*cin22-k23*cin32; - c32=cin32-k20*cin30-k21*cin31-k22*cin32-k23*cin33; - c42=cin42-k20*cin40-k21*cin41-k22*cin42-k23*cin43; - - c33=cin33-k30*cin30-k31*cin31-k32*cin32-k33*cin33; - c43=cin43-k30*cin40-k31*cin41-k32*cin42-k33*cin43; - - c44=cin44-k40*cin40-k41*cin41-k42*cin42-k43*cin43; - - newTrack->PutXElements(x0,x1,x2,x3,x4); // put the new state vector - - newTrack->PutCElements(c00,c10,c11,c20,c21,c22,c30,c31,c32,c33,c40,c41,c42,c43,c44); // put in track the - // new cov matrix - - Double_t vmc[4][4]; - - vmc[0][0]=v00-c00; vmc[1][0]=-c10; vmc[2][0]=-c20; vmc[3][0]=-c30; - vmc[1][1]=v11-c11; vmc[2][1]=-c21; vmc[3][1]=v31-c31; - vmc[2][2]=v22-c22; vmc[3][2]=v32-c32; - vmc[3][3]=v33-c33; - - vmc[0][1]=vmc[1][0]; vmc[0][2]=vmc[2][0]; vmc[0][3]=vmc[3][0]; - vmc[1][2]=vmc[2][1]; vmc[1][3]=vmc[3][1]; - vmc[2][3]=vmc[3][2]; - - -/////////////////////// vmc matrix inversion /////////////////////////////////// - - d=1.; - - for(k=0; k k) { - for(i=0; i k ) { - for(j=0; j=0; k--) { - i=ll[k]; - if(i > k) { - for (j=0; j k) { - for (i=0; iSetChi2(newTrack->GetChi2()+chi2); - // newTrack->SetChi2(newTrack->GetChi2()+chi2pred); - -} - +} +//---------------------------------------------------------------------- +//void AliITSTrackerV1::KalmanFilterVert(AliITSTrackV1 *newTrack, +// TVector &cluster,Double_t sigma[2]){ +void AliITSTrackerV1::KalmanFilterVert(AliITSTrackV1 *newTrack, + TVector &cluster,Double_t sigma[2] + /*, Double_t chi2pred*/){ + //Origin A. Badala' and G.S. Pappalardo: + // e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it + // Kalman filter with vertex constraint + ///////////////////// Evaluation of the measurement vector m //////// + Double_t m[4]; + Double_t rk,phik,zk; + rk=cluster(0); phik=cluster(1); zk=cluster(2); + m[0]=phik; m[1]=zk; + Double_t cc=(*newTrack).GetC(); + Double_t zv=(*newTrack).GetZv(); + Double_t dv=(*newTrack).GetDv(); + Double_t cy=cc/2.; + Double_t tgl= (zk-zv)*cy/TMath::ASin(cy*rk); + m[2]=dv; m[3]=tgl; + /////////////////////// Evaluation of the error matrix V ////////// + Int_t layer=newTrack->GetLayer(); + Double_t v00=sigma[0]; + Double_t v11=sigma[1]; + Double_t v31=sigma[1]/rk; + Double_t sigmaDv=newTrack->GetsigmaDv(); + Double_t v22=sigmaDv*sigmaDv + newTrack->Getd2(layer-1); + Double_t v32=newTrack->Getdtgl(layer-1); + Double_t sigmaZv=newTrack->GetsigmaZv(); + Double_t v33=(sigma[1]+sigmaZv*sigmaZv)/(rk*rk)+newTrack->Gettgl2(layer-1); + ////////////////////////////////////////////////////////////////// + Double_t cin00,cin10,cin11,cin20,cin21,cin22, + cin30,cin31,cin32,cin33,cin40,cin41,cin42,cin43,cin44; + newTrack->GetCElements(cin00, + cin10,cin11, + cin20,cin21,cin22, + cin30,cin31,cin32,cin33, + cin40,cin41,cin42,cin43,cin44); //get C matrix + Double_t r[4][4]; + r[0][0]=cin00+v00; + r[1][0]=cin10; + r[2][0]=cin20; + r[3][0]=cin30; + r[1][1]=cin11+v11; + r[2][1]=cin21; + r[3][1]=cin31+sigma[1]/rk; + r[2][2]=cin22+sigmaDv*sigmaDv+newTrack->Getd2(layer-1); + r[3][2]=cin32+newTrack->Getdtgl(layer-1); + r[3][3]=cin33+(sigma[1]+sigmaZv*sigmaZv)/(rk*rk)+ + newTrack->Gettgl2(layer-1); + r[0][1]=r[1][0]; r[0][2]=r[2][0]; r[0][3]=r[3][0]; + r[1][2]=r[2][1]; r[1][3]=r[3][1]; r[2][3]=r[3][2]; + ///////////////////// Matrix R inversion ////////////////////////// + const Int_t kn=4; + Double_t big, hold; + Double_t d=1.; + Int_t ll[kn],mm[kn]; + Int_t i,j,k; + + for(k=0; k k) { + for(i=0; i k ) { + for(j=0; j=0; k--) { + i=ll[k]; + if(i > k) { + for (j=0; j k) { + for (i=0; iGetXElements(x0,x1,x2,x3,x4); // get the state vector + Double_t savex0=x0, savex1=x1, savex2=x2, savex3=x3; + x0+=k00*(m[0]-savex0)+k01*(m[1]-savex1)+k02*(m[2]-savex2)+ + k03*(m[3]-savex3); + x1+=k10*(m[0]-savex0)+k11*(m[1]-savex1)+k12*(m[2]-savex2)+ + k13*(m[3]-savex3); + x2+=k20*(m[0]-savex0)+k21*(m[1]-savex1)+k22*(m[2]-savex2)+ + k23*(m[3]-savex3); + x3+=k30*(m[0]-savex0)+k31*(m[1]-savex1)+k32*(m[2]-savex2)+ + k33*(m[3]-savex3); + x4+=k40*(m[0]-savex0)+k41*(m[1]-savex1)+k42*(m[2]-savex2)+ + k43*(m[3]-savex3); + Double_t c00,c10,c20,c30,c40,c11,c21,c31,c41,c22,c32,c42,c33,c43,c44; + c00=cin00-k00*cin00-k01*cin10-k02*cin20-k03*cin30; + c10=cin10-k00*cin10-k01*cin11-k02*cin21-k03*cin31; + c20=cin20-k00*cin20-k01*cin21-k02*cin22-k03*cin32; + c30=cin30-k00*cin30-k01*cin31-k02*cin32-k03*cin33; + c40=cin40-k00*cin40-k01*cin41-k02*cin42-k03*cin43; + c11=cin11-k10*cin10-k11*cin11-k12*cin21-k13*cin31; + c21=cin21-k10*cin20-k11*cin21-k12*cin22-k13*cin32; + c31=cin31-k10*cin30-k11*cin31-k12*cin32-k13*cin33; + c41=cin41-k10*cin40-k11*cin41-k12*cin42-k13*cin43; + c22=cin22-k20*cin20-k21*cin21-k22*cin22-k23*cin32; + c32=cin32-k20*cin30-k21*cin31-k22*cin32-k23*cin33; + c42=cin42-k20*cin40-k21*cin41-k22*cin42-k23*cin43; + c33=cin33-k30*cin30-k31*cin31-k32*cin32-k33*cin33; + c43=cin43-k30*cin40-k31*cin41-k32*cin42-k33*cin43; + c44=cin44-k40*cin40-k41*cin41-k42*cin42-k43*cin43; + + newTrack->PutXElements(x0,x1,x2,x3,x4); // put the new state vector + newTrack->PutCElements(c00, + c10,c11, + c20,c21,c22, + c30,c31,c32,c33, + c40,c41,c42,c43,c44); // put in track the + // new cov matrix + Double_t vmc[4][4]; + vmc[0][0]=v00-c00; vmc[1][0]=-c10; vmc[2][0]=-c20; vmc[3][0]=-c30; + vmc[1][1]=v11-c11; vmc[2][1]=-c21; vmc[3][1]=v31-c31; + vmc[2][2]=v22-c22; vmc[3][2]=v32-c32; + vmc[3][3]=v33-c33; + vmc[0][1]=vmc[1][0]; vmc[0][2]=vmc[2][0]; vmc[0][3]=vmc[3][0]; + vmc[1][2]=vmc[2][1]; vmc[1][3]=vmc[3][1]; + vmc[2][3]=vmc[3][2]; + /////////////////////// vmc matrix inversion /////////////////////// + d=1.; + for(k=0; k k) { + for(i=0; i k ) { + for(j=0; j=0; k--) { + i=ll[k]; + if(i > k) { + for (j=0; jk + j=mm[k]; + if(j > k) { + for (i=0; ik + } // end for k + //////////////////////////////////////////////////////////////////// + Double_t chi2=(m[0]-x0)*( vmc[0][0]*(m[0]-x0) + 2.*vmc[1][0]*(m[1]-x1) + + 2.*vmc[2][0]*(m[2]-x2)+ 2.*vmc[3][0]*(m[3]-x3) )+ + (m[1]-x1)* ( vmc[1][1]*(m[1]-x1) + 2.*vmc[2][1]*(m[2]-x2)+ + 2.*vmc[3][1]*(m[3]-x3) ) + + (m[2]-x2)* ( vmc[2][2]*(m[2]-x2)+ 2.*vmc[3][2]*(m[3]-x3) ) + + (m[3]-x3)*vmc[3][3]*(m[3]-x3); + //cout<<" chi2 kalman = "<SetChi2(newTrack->GetChi2()+chi2); + // newTrack->SetChi2(newTrack->GetChi2()+chi2pred); +} diff --git a/ITS/AliITSTrackerV1.h b/ITS/AliITSTrackerV1.h index 274471f397b..8588b4b421e 100644 --- a/ITS/AliITSTrackerV1.h +++ b/ITS/AliITSTrackerV1.h @@ -1,21 +1,23 @@ -//Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * -//See cxx source for full Copyright notice */ +//Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. +//See cxx source for full Copyright notice // //The purpose of this class is to permorm the ITS tracking. // -//The constructor has the task to inizialize some private members. -//The method DoTracking is written to be called by a macro. It gets the event number, the minimum and maximum -//order number of TPC tracks that are to be tracked trough the ITS, and the file where the recpoints are -//registered. +// The constructor has the task to inizialize some private members. +// The method DoTracking is written to be called by a macro. It gets the event +// number, the minimum and maximum order number of TPC tracks that are to be +// tracked trough the ITS, and the file where the recpoints are registered. // -//The method AliiTStracking is a recursive function that performs the tracking trough the ITS +// The method AliiTStracking is a recursive function that performs the +// tracking trough the ITS // -//The method Intersection found the layer, ladder and detector whre the intersection take place and caluclate -//the cohordinates of this intersection. It returns an integer that is 0 if the intersection has been found -//successfully. +// The method Intersection found the layer, ladder and detector where the +// intersection take place and caluclate the cohordinates of this +// intersection. It returns an integer that is 0 if the intersection has +// been found successfully. // -//The two mwthods Kalmanfilter and kalmanfiltervert operate the kalmanfilter without and with the vertex -//imposition respectively. +// The two mwthods Kalmanfilter and kalmanfiltervert operate the +// kalmanfilter without and with the verteximposition respectively. #ifndef ALIITSTRACKERV1_H #define ALIITSTRACKERV1_H @@ -33,60 +35,49 @@ class AliITSgeoinfo; class TStopwatch; class AliITSTrackerV1 : public TObject { - - public: - - AliITSTrackerV1(AliITS* IITTSS, Bool_t flag); - - AliITSTrackerV1(const AliITSTrackerV1 &cobj); - - ~AliITSTrackerV1(); - + public: + AliITSTrackerV1(AliITS* IITTSS, Bool_t flag); + AliITSTrackerV1(const AliITSTrackerV1 &cobj); + ~AliITSTrackerV1(); AliITSTrackerV1 &operator=(AliITSTrackerV1 obj); - - void DoTracking(Int_t evNumber, Int_t minTr, Int_t maxTr, TFile *file); - - void RecursiveTracking(TList *trackITSlist); - - Int_t Intersection(AliITSTrackV1 &track, Int_t layer, Int_t &ladder, Int_t &detector); - void KalmanFilter(AliITSTrackV1 *newtrack, TVector &cluster, Double_t sigma[2]); - - void KalmanFilterVert(AliITSTrackV1 *newtrack, TVector &cluster, Double_t sigma[2]); - //void KalmanFilterVert(AliITSTrackV1 *newtrack, TVector &cluster, Double_t sigma[2], Double_t chi2pred); - - private: + void DoTracking(Int_t evNumber, Int_t minTr, Int_t maxTr, TFile *file); + void RecursiveTracking(TList *trackITSlist); + Int_t Intersection(AliITSTrackV1 &track, Int_t layer,Int_t &ladder, + Int_t &detector); + void KalmanFilter(AliITSTrackV1 *newtrack, TVector &cluster, + Double_t sigma[2]); + void KalmanFilterVert(AliITSTrackV1 *newtrack, TVector &cluster, + Double_t sigma[2]); + //void KalmanFilterVert(AliITSTrackV1 *newtrack, TVector &cluster, + // Double_t sigma[2], Double_t chi2pred); + + private: + AliITS* fITS; //! pointer to AliITS + AliITSTrackV1 *fresult; // result is a pointer to the final best track + Double_t fPtref; // transvers momentum obtained from TPC tracking + Double_t fChi2max; // cluster with chi2>chi2max are cut. It is + // pt dependend. aggiunto il 31-7-2001 + Double_t fepsphi; //eps for definition window in phi aggiunto il 1-8-2001 + Double_t fepsz; //eps for definition window in z aggiunto il 1-8-2001 + TObjArray *frecPoints; // pointer to RecPoints + Int_t **fvettid; // flag vector of used clusters + Bool_t fflagvert; // a flag to impose or not the vertex constraint + AliITSRad *frl; // pointer to get the radiation lenght matrix + Int_t fNlad[6]; // Number of ladders for a given layer + Int_t fNdet[6]; // Number of detector for a given layer + Double_t fAvrad[6]; // Average radius for a given layer + Double_t fDetx[6];// Semidimension of detectors along x axis for a given layer + Double_t fDetz[6];// Semidimension of detectors along z axis for a given layer + Double_t **fzmin;// Matrix of zmin for a given layer and a given detector + Double_t **fzmax;// Matrix of zmax for a given layer and a given detector + Double_t **fphimin;// Matrix of phimin for a given layer and a given ladder + Double_t **fphimax;// Matrix of phimax for a given layer and a given ladder + Double_t **fphidet; // azimuthal angle for a given layer and a given ladder + Double_t fFieldFactor; // Magnetic filed factor + //TStopwatch *fTimerKalman; // timer for kalman filter + //TStopwatch *fTimerIntersection; // timer for Intersection - AliITS* fITS; // pointer to AliITS - AliITSTrackV1 *fresult; // result is a pointer to the final best track - Double_t fPtref; // transvers momentum obtained from TPC tracking - Double_t fChi2max; // cluster with chi2>chi2max are cut. It is pt dependend. aggiunto il 31-7-2001 - Double_t fepsphi; //eps for definition window in phi //aggiunto il 1-8-2001 - Double_t fepsz; //eps for definition window in z //aggiunto il 1-8-2001 - TObjArray *frecPoints; // pointer to RecPoints - Int_t **fvettid; // flag vector of used clusters - Bool_t fflagvert; // a flag to impose or not the vertex constraint - AliITSRad *frl; // pointer to get the radiation lenght matrix - - Int_t fNlad[6]; // Number of ladders for a given layer - Int_t fNdet[6]; // Number of detector for a given layer - Double_t fAvrad[6]; // Average radius for a given layer - Double_t fDetx[6]; // Semidimension of detectors along x axis for a given layer - Double_t fDetz[6]; // Semidimension of detectors along z axis for a given layer - - Double_t **fzmin; // Matrix of zmin for a given layer and a given detector - Double_t **fzmax; // Matrix of zmax for a given layer and a given detector - - Double_t **fphimin; // Matrix of phimin for a given layer and a given ladder - Double_t **fphimax; // Matrix of phimax for a given layer and a given ladder - - Double_t **fphidet; // azimuthal angle for a given layer and a given ladder - - Double_t fFieldFactor; // Magnetic filed factor - - //TStopwatch *fTimerKalman; // timer for kalman filter - //TStopwatch *fTimerIntersection; // timer for Intersection - - ClassDef(AliITSTrackerV1,1) + ClassDef(AliITSTrackerV1,1) //???? }; #endif diff --git a/ITS/AliITSgeomSDD.cxx b/ITS/AliITSgeomSDD.cxx index 3905dac1212..b7ebad9012c 100644 --- a/ITS/AliITSgeomSDD.cxx +++ b/ITS/AliITSgeomSDD.cxx @@ -15,6 +15,11 @@ /* $Log$ +Revision 1.13 2001/10/12 22:07:20 nilsen +A patch for C++ io manipulation functions so that they will work both +with GNU gcc 2.96 and GNU gcc 3.01 compilers. Needs to be tested with +other platforms. + Revision 1.12 2001/08/24 21:06:37 nilsen Added more documentation, fixed up some coding violations, and some forward declorations. @@ -250,7 +255,7 @@ void AliITSgeomSDD::Print(ostream *os) const { #endif #else Int_t fmt; -#endif; +#endif fmt = os->setf(ios::scientific); // set scientific floating point output *os << "TBRIK" << " "; @@ -314,6 +319,11 @@ istream &operator>>(istream &is,AliITSgeomSDD &r){ //====================================================================== /* $Log$ +Revision 1.13 2001/10/12 22:07:20 nilsen +A patch for C++ io manipulation functions so that they will work both +with GNU gcc 2.96 and GNU gcc 3.01 compilers. Needs to be tested with +other platforms. + Revision 1.12 2001/08/24 21:06:37 nilsen Added more documentation, fixed up some coding violations, and some forward declorations. @@ -779,6 +789,11 @@ istream &operator>>(istream &is,AliITSgeomSDD256 &r){ //====================================================================== /* $Log$ +Revision 1.13 2001/10/12 22:07:20 nilsen +A patch for C++ io manipulation functions so that they will work both +with GNU gcc 2.96 and GNU gcc 3.01 compilers. Needs to be tested with +other platforms. + Revision 1.12 2001/08/24 21:06:37 nilsen Added more documentation, fixed up some coding violations, and some forward declorations. diff --git a/ITS/AliITSgeomSPD.cxx b/ITS/AliITSgeomSPD.cxx index 5517066d305..ed334699a65 100644 --- a/ITS/AliITSgeomSPD.cxx +++ b/ITS/AliITSgeomSPD.cxx @@ -15,6 +15,11 @@ /* $Log$ +Revision 1.13 2001/10/12 22:07:20 nilsen +A patch for C++ io manipulation functions so that they will work both +with GNU gcc 2.96 and GNU gcc 3.01 compilers. Needs to be tested with +other platforms. + Revision 1.12 2001/08/24 21:06:37 nilsen Added more documentation, fixed up some coding violations, and some forward declorations. @@ -186,7 +191,7 @@ void AliITSgeomSPD::Print(ostream *os) const { #endif #else Int_t fmt; -#endif; +#endif fmt = os->setf(ios::scientific); // set scientific floating point output *os << "TBRIK" << " "; @@ -244,6 +249,11 @@ istream &operator>>(istream &is,AliITSgeomSPD &r){ /* $Log$ +Revision 1.13 2001/10/12 22:07:20 nilsen +A patch for C++ io manipulation functions so that they will work both +with GNU gcc 2.96 and GNU gcc 3.01 compilers. Needs to be tested with +other platforms. + Revision 1.12 2001/08/24 21:06:37 nilsen Added more documentation, fixed up some coding violations, and some forward declorations. @@ -348,6 +358,11 @@ istream &operator>>(istream &is,AliITSgeomSPD300 &r){ //===================================================================== /* $Log$ +Revision 1.13 2001/10/12 22:07:20 nilsen +A patch for C++ io manipulation functions so that they will work both +with GNU gcc 2.96 and GNU gcc 3.01 compilers. Needs to be tested with +other platforms. + Revision 1.12 2001/08/24 21:06:37 nilsen Added more documentation, fixed up some coding violations, and some forward declorations. @@ -482,6 +497,11 @@ istream &operator>>(istream &is,AliITSgeomSPD425Short &r){ /* $Log$ +Revision 1.13 2001/10/12 22:07:20 nilsen +A patch for C++ io manipulation functions so that they will work both +with GNU gcc 2.96 and GNU gcc 3.01 compilers. Needs to be tested with +other platforms. + Revision 1.12 2001/08/24 21:06:37 nilsen Added more documentation, fixed up some coding violations, and some forward declorations. diff --git a/ITS/AliITSsimulationSSD.cxx b/ITS/AliITSsimulationSSD.cxx index 487f6a94641..410d13034e9 100644 --- a/ITS/AliITSsimulationSSD.cxx +++ b/ITS/AliITSsimulationSSD.cxx @@ -370,11 +370,11 @@ void AliITSsimulationSSD::IntegrateGaussian(Int_t k,Double_t par, Double_t w, // for the time being, signal is the charge // in ChargeToSignal signal is converted in ADC channel - fMapA2->AddSignal(k,strip,dXCharge1); + fMapA2->AddSignal(k,(Int_t)strip,dXCharge1); if(((Int_t) strip) < (GetNStrips()-1)) { // strip doesn't have to be the last (remind: last=GetNStrips()-1) // otherwise part of the charge is lost - fMapA2->AddSignal(k,(strip+1),dXCharge2); + fMapA2->AddSignal(k,((Int_t)strip+1),dXCharge2); } // end if if(dXCharge1 > 1.) { @@ -408,11 +408,11 @@ void AliITSsimulationSSD::IntegrateGaussian(Int_t k,Double_t par, Double_t w, // for the time being, signal is the charge // in ChargeToSignal signal is converted in ADC channel - fMapA2->AddSignal(k,strip,dXCharge1); + fMapA2->AddSignal(k,(Int_t)strip,dXCharge1); if(((Int_t) strip) > 0) { // strip doesn't have to be the first // otherwise part of the charge is lost - fMapA2->AddSignal(k,(strip-1),dXCharge2); + fMapA2->AddSignal(k,((Int_t)strip-1),dXCharge2); } // end if if(dXCharge1 > 1.) { diff --git a/ITS/AliITSvPPRsymm.cxx b/ITS/AliITSvPPRsymm.cxx index 7c15f8fd088..f68062ca369 100644 --- a/ITS/AliITSvPPRsymm.cxx +++ b/ITS/AliITSvPPRsymm.cxx @@ -15,6 +15,9 @@ /* $Log$ +Revision 1.35 2001/10/19 10:16:28 barbera +A bug corrected in the definition of a TPCON + Revision 1.34 2001/10/18 12:25:07 barbera Detailed geometry in BuildGeometry() commented out (450 MB needed to compile the file). Six cylinders put back but improved by comparison with the ITS coarse geometry @@ -29110,7 +29113,7 @@ void AliITSvPPRsymm::SetDefaults(){ s1->GetDz()*2.*kconv, // for now. s1->GetDy()*2.*kconv); // x,z,y full width in microns. bx[0] = 1000./((s1->GetDx()*kconv/seg1->Dpx(0))/resp1->DriftSpeed()); // clock in Mhz - seg1->SetNPads(256,bx[0]);// Use AliITSgeomSDD for now + seg1->SetNPads(256,(Int_t)(bx[0]));// Use AliITSgeomSDD for now SetSegmentationModel(kSDD,seg1); const char *kData1=(iDetType->GetResponseModel())->DataType(); const char *kopt=iDetType->GetResponseModel()->ZeroSuppOption(); -- 2.39.3