#include <TROOT.h>
#include <TFile.h>
#include <TPRegexp.h>
+#include <TVectorF.h>
#include <AliExternalTrackParam.h>
#include <AliTPCcalibDB.h>
, fNmaxEvents(-1)
, fTime0(-1)
, fCreateT0seed(kFALSE)
+, fLongT0seed(kTRUE)
+, fFillClusterRes(kFALSE)
+, fUseT0list(kFALSE)
+, fUseZ0list(kFALSE)
+, fForceAlpha(kFALSE)
, fStreamer(0x0)
, fInputFile(0x0)
, fTree(0x0)
ConnectInputFile(file, nmaxEv);
if (!fTree) return;
+ Int_t maxev=fTree->GetEntries();
+ if (nmaxEv>0&&nmaxEv<maxev) maxev=nmaxEv;
+
InitStreamer(".debug");
gROOT->cd();
AliExternalTrackParam trackITS2;
AliExternalTrackParam *dummy;
-
- Int_t maxev=fTree->GetEntries();
- if (nmaxEv>0&&nmaxEv<maxev) maxev=nmaxEv;
+ // prepare list of T0s
+ TVectorF t0list(maxev);
+ TVectorF z0list(maxev);
+ if (fUseT0list || fUseZ0list) {
+ for (Int_t iev=0; iev<maxev; ++iev){
+ fTree->GetEvent(iev);
+ const Float_t t0=fEvent->GetT0();
+ const Float_t z0=fEvent->GetZ();
+ t0list[iev]=t0;
+ z0list[iev]=z0;
+ }
+ }
+
+ // array with cluster residuals
+ TClonesArray *arrClustRes=0x0;
+ if (fFillClusterRes){
+ arrClustRes=new TClonesArray("AliTPCclusterMI",160);
+ }
+
const Double_t lastLayerITS = 43.0; // same as in AliToyMCEventGenerator::MakeITSClusters (hard coded)
const Double_t iFCRadius = 83.5; //radius constants found in AliTPCCorrection.cxx
const Double_t betweeTPCITS = (lastLayerITS+iFCRadius)/2.; // its track propgated to inner TPC wall
const Double_t kMaxSnp = 0.85;
const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
+ Double_t lastT0=0;
for (Int_t iev=0; iev<maxev; ++iev){
printf("============== Processing Event %6d =================\n",iev);
fTree->GetEvent(iev);
+
+ Float_t z0=fEvent->GetZ();
+ Float_t t0=fEvent->GetT0();
+
for (Int_t itr=0; itr<fEvent->GetNumberOfTracks(); ++itr){
// printf(" > ====== Processing Track %6d ======== \n",itr);
const AliToyMCTrack *tr=fEvent->GetTrack(itr);
trackITS1 = resetParam;
trackITS2 = resetParam;
- Float_t z0=fEvent->GetZ();
- Float_t t0=fEvent->GetT0();
-
Float_t vDrift=GetVDrift();
Float_t zLength=GetZLength(0);
t0seed = *dummy;
delete dummy;
- // crate real seed using the time 0 from the first seed
- // set fCreateT0seed now to false to get the seed in z coordinates
+ // Long seed
+ if (fLongT0seed){
+ dummy = GetFittedTrackFromSeed(tr,&t0seed);
+ t0seed = *dummy;
+ delete dummy;
+ }
+
+ // set the T0 from the seed
+ // in case the match with the real T0 infor is requested, find the
+ // closes T0 from the list of T0s
fTime0 = t0seed.GetZ()-zLength/vDrift;
+ if (fUseT0list || fUseZ0list) {
+ fTime0 = FindClosestT0(t0list, z0list, t0seed);
+ }
+ // create real seed using the time 0 from the first seed
+ // set fCreateT0seed now to false to get the seed in z coordinates
fCreateT0seed = kFALSE;
// printf("seed (%.2g)\n",fTime0);
dummy = GetSeedFromTrack(tr);
// create fitted track
if (fDoTrackFit){
// printf("track\n");
- dummy = GetFittedTrackFromSeed(tr, &seed);
+ dummy = GetFittedTrackFromSeed(tr, &seed, arrClustRes);
track = *dummy;
delete dummy;
}
- // Copy original track and fitted track
- // for extrapolation to ITS last layer
- trackITS = track;
+ // Copy original track and fitted track
+ // for extrapolation to ITS last layer
+ trackITS = track;
trackITS1 = track;
trackITS2 = track;
// ITS comparison
//
- // rotate fitted track to the frame of the original track and propagate to same reference
- AliTrackerBase::PropagateTrackTo(&trackITS,lastLayerITS,kMass,5,kTRUE,kMaxSnp,0,kFALSE,fUseMaterial);
- trackITS.Rotate(tOrigITS.GetAlpha());
- AliTrackerBase::PropagateTrackTo(&trackITS,lastLayerITS,kMass,1,kFALSE,kMaxSnp,0,kFALSE,fUseMaterial);
-
// rotate fitted track to the frame of the original track and propagate to same reference
AliTrackerBase::PropagateTrackTo(&trackITS,lastLayerITS,kMass,5,kTRUE,kMaxSnp,0,kFALSE,fUseMaterial);
trackITS.Rotate(tOrigITS.GetAlpha());
AliTrackerBase::PropagateTrackTo(&trackITS,lastLayerITS,kMass,1,kFALSE,kMaxSnp,0,kFALSE,fUseMaterial);
// rotate fitted track to the frame of the original track and propagate to same reference
- AliTrackerBase::PropagateTrackTo(&trackITS,lastLayerITS,kMass,5,kTRUE,kMaxSnp,0,kFALSE,fUseMaterial);
- trackITS.Rotate(tOrigITS.GetAlpha());
- AliTrackerBase::PropagateTrackTo(&trackITS,lastLayerITS,kMass,1,kFALSE,kMaxSnp,0,kFALSE,fUseMaterial);
-
-
+ AliTrackerBase::PropagateTrackTo(&trackITS1,betweeTPCITS,kMass,5,kTRUE,kMaxSnp,0,kFALSE,fUseMaterial);
+ trackITS1.Rotate(tOrigITS1.GetAlpha());
+ AliTrackerBase::PropagateTrackTo(&trackITS1,betweeTPCITS,kMass,1,kFALSE,kMaxSnp,0,kFALSE,fUseMaterial);
+
+ // rotate fitted track to the frame of the original track and propagate to same reference
+ AliTrackerBase::PropagateTrackTo(&trackITS2,iFCRadius,kMass,5,kTRUE,kMaxSnp,0,kFALSE,fUseMaterial);
+ trackITS2.Rotate(tOrigITS2.GetAlpha());
+ AliTrackerBase::PropagateTrackTo(&trackITS2,iFCRadius,kMass,1,kFALSE,kMaxSnp,0,kFALSE,fUseMaterial);
}
}
"iev=" << iev <<
"z0=" << z0 <<
"t0=" << t0 <<
+ "lastt0=" << lastT0 <<
"fTime0=" << fTime0 <<
"itr=" << itr <<
"clsType=" << fClusterType <<
"track.=" << &track <<
// ITS match
"tOrigITS.=" << &tOrigITS <<
- "tOrigITS1.=" << &tOrigITS <<
- "tOrigITS2.=" << &tOrigITS <<
+ "tOrigITS1.=" << &tOrigITS1 <<
+ "tOrigITS2.=" << &tOrigITS2 <<
+
"tRealITS.=" << &tRealITS <<
- "tRealITS1.=" << &tRealITS <<
- "tRealITS2.=" << &tRealITS <<
+ "tRealITS1.=" << &tRealITS1 <<
+ "tRealITS2.=" << &tRealITS2 <<
+
"trackITS.=" << &trackITS <<
- "trackITS1.=" << &trackITS <<
- "trackITS2.=" << &trackITS <<
+ "trackITS1.=" << &trackITS1 <<
+ "trackITS2.=" << &trackITS2;
+
+ if (arrClustRes) {
+ const Int_t nCl=arrClustRes->GetEntriesFast();
+ // fracktion of outliers from track extrapolation
+ // for 3, 3.5, 4, 4.5 and 5 sigma of the cluster resolution (~1mm)
+ Float_t fracY[5]={0.};
+ Float_t fracZ[5]={0.};
+
+ for (Int_t icl=0; icl<nCl; ++icl) {
+ AliTPCclusterMI *cl=static_cast<AliTPCclusterMI*>(arrClustRes->At(icl));
+// const Float_t sigmaY=TMath::Sqrt(cl->GetSigmaY2());
+// const Float_t sigmaZ=TMath::Sqrt(cl->GetSigmaZ2());
+ for (Int_t inSig=0; inSig<5; ++inSig) {
+ fracY[inSig] += cl->GetY()>(3+inSig*.5)/**sigmaY*/;
+ fracZ[inSig] += cl->GetZ()>(3+inSig*.5)/**sigmaZ*/;
+ }
+ }
+
+ if (nCl>0) {
+ for (Int_t inSig=0; inSig<5; ++inSig) {
+ fracY[inSig]/=nCl;
+ fracZ[inSig]/=nCl;
+ }
+ }
+
+ (*fStreamer) << "Tracks" <<
+ "clustRes.=" << arrClustRes;
+ for (Int_t inSig=0; inSig<5; ++inSig) {
+ const char* fracYname=Form("clFracY%02d=", 30+inSig*5);
+ const char* fracZname=Form("clFracZ%02d=", 30+inSig*5);
+ (*fStreamer) << "Tracks" <<
+ fracYname << fracY[inSig] <<
+ fracZname << fracZ[inSig];
+ }
+ }
+
+ (*fStreamer) << "Tracks" <<
"\n";
}
}
+ lastT0=t0;
}
+ delete arrClustRes;
Cleanup();
}
AliError(Form("Seeding failed for parameters %d, %d\n",fSeedingDist,fSeedingRow));
return 0x0;
}
+
+ // determine preliminary theta
+ Float_t xyz1[3]={0,0,0};
+ Float_t xyz2[3]={0,0,0};
+ seedPoint[0].GetXYZ(xyz1);
+ seedPoint[2].GetXYZ(xyz2);
+ Float_t prelDeltaR = TMath::Sqrt(xyz2[0]*xyz2[0]+xyz2[1]*xyz2[1]) - TMath::Sqrt(xyz1[0]*xyz1[0]+xyz1[1]*xyz1[1]) ;
+ Float_t prelDeltaZ = ( seedCluster[0]->GetTimeBin() - seedCluster[2]->GetTimeBin() ) * GetVDrift();
+ Float_t prelTheta = TMath::ATan(prelDeltaR/prelDeltaZ);
+ if(prelTheta > TMath::Pi()/2) prelTheta = TMath::Pi() - prelTheta;
// do cluster correction for fCorrectionType:
// 0 - no correction
// 1 - TPC center
// 2 - average eta
// 3 - ideal
+ // 4 - preliminary eta (needs fixing!!! Not yet in full code!!!)
// assign the cluster abs time as z component to all seeds
for (Int_t iseed=0; iseed<3; ++iseed) {
Float_t xyz[3]={0,0,0};
const Int_t sector=seedCluster[iseed]->GetDetector();
const Int_t sign=1-2*((sector/18)%2);
+
+ Float_t zBeforeCorr = xyz[2];
if ( (fClusterType == 1) && (fCorrectionType != kNoCorrection) ) {
-// printf("correction type: %d\n",(Int_t)fCorrectionType);
-
// the settings below are for the T0 seed
// for known T0 the z position is already calculated in SetTrackPointFromCluster
if ( fCreateT0seed ){
if ( fCorrectionType == kTPCCenter ) xyz[2] = 125.*sign;
//!!! TODO: is this the correct association?
if ( fCorrectionType == kAverageEta ) xyz[2] = TMath::Tan(45./2.*TMath::DegToRad())*r*sign;
+ if ( fCorrectionType == kPreliminaryEta ) xyz[2] = r/TMath::Tan(prelTheta)*sign;//(needs fixing!!! Not yet in full code!!!)
}
if ( fCorrectionType == kIdeal ) xyz[2] = seedCluster[iseed]->GetZ();
-
+
+ // Store xyz only here!!! To get the Delta z from the correction...
+ zBeforeCorr = xyz[2];
+
//!!! TODO: to be replaced with the proper correction
fTPCCorrection->CorrectPoint(xyz, seedCluster[iseed]->GetDetector());
}
// after the correction set the time bin as z-Position in case of a T0 seed
if ( fCreateT0seed )
- xyz[2]=seedCluster[iseed]->GetTimeBin();
-
+ xyz[2]=seedCluster[iseed]->GetTimeBin() + sign*( zBeforeCorr - xyz[2] )/GetVDrift();
+
seedPoint[iseed].SetXYZ(xyz);
}
const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
AliExternalTrackParam *seed = AliTrackerBase::MakeSeed(seedPoint[0], seedPoint[1], seedPoint[2]);
- seed->ResetCovariance(10);
- if (fCreateT0seed){
+ if (fCreateT0seed&&!fLongT0seed){
+ // only propagate to vertex if we don't create a long seed
// if fTime0 < 0 we assume that we create a seed for the T0 estimate
- AliTrackerBase::PropagateTrackTo(seed,0,kMass,5,kTRUE,kMaxSnp,0,kFALSE,kFALSE);
+ Int_t ret=AliTrackerBase::PropagateTrackTo2(seed,0,kMass,5,kTRUE,kMaxSnp,0,kFALSE,kFALSE);
if (TMath::Abs(seed->GetX())>3) {
-// printf("Could not propagate track to 0, %.2f, %.2f, %.2f\n",seed->GetX(),seed->GetAlpha(),seed->Pt());
+// printf("Could not propagate track to 0, x:%.2f, a:%.2f (%.2f), snp:%.2f (%.2f), pt:%.2f (%.2f), %d\n",seed->GetX(),seed->GetAlpha(),tr->GetAlpha(), seed->GetSnp(), tr->GetSnp(), seed->Pt(), tr->Pt(), ret);
+ }
+ if (fForceAlpha) {
+ seed->Rotate(tr->GetAlpha());
+ AliTrackerBase::PropagateTrackTo2(seed,0,kMass,1.,kFALSE,kMaxSnp,0,kFALSE,kFALSE);
}
}
-
+
+ seed->ResetCovariance(10);
return seed;
}
}
//____________________________________________________________________________________
-AliExternalTrackParam* AliToyMCReconstruction::GetFittedTrackFromSeed(const AliToyMCTrack *tr, const AliExternalTrackParam *seed)
+AliExternalTrackParam* AliToyMCReconstruction::GetFittedTrackFromSeed(const AliToyMCTrack *tr, const AliExternalTrackParam *seed, TClonesArray *arrClustRes)
{
//
//
//
+ if (arrClustRes) {
+ arrClustRes->Clear();
+ }
+
// create track
AliExternalTrackParam *track = new AliExternalTrackParam(*seed);
AliTrackPoint pIn;
const AliTPCclusterMI *cl=tr->GetSpacePoint(ipoint);
if (fClusterType == 1) cl=tr->GetDistortedSpacePoint(ipoint);
+ const Int_t globalRow = cl->GetRow() +(cl->GetDetector() >35)*63;
+ if ( fCreateT0seed ){
+ if ( globalRow<fSeedingRow || globalRow>fSeedingRow+2*fSeedingDist ) continue;
+ }
+
SetTrackPointFromCluster(cl, pIn);
+
+ Float_t xyz[3]={0,0,0};
+ pIn.GetXYZ(xyz);
+ Float_t zBeforeCorr = xyz[2];
+
+ const Int_t sector=cl->GetDetector();
+ const Int_t sign=1-2*((sector/18)%2);
+
if (fCorrectionType != kNoCorrection){
- Float_t xyz[3]={0,0,0};
- pIn.GetXYZ(xyz);
-// if ( fCorrectionType == kIdeal ) xyz[2] = cl->GetZ();
+
+ const Float_t r=TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
+
+ if ( fCreateT0seed ){
+ if ( fCorrectionType == kTPCCenter ) xyz[2] = 125.*sign;
+ //!!! TODO: is this the correct association?
+ if ( fCorrectionType == kAverageEta ) xyz[2] = TMath::Tan(45./2.*TMath::DegToRad())*r*sign;
+ if ( fCorrectionType == kIdeal ) xyz[2] = cl->GetZ();
+ }
+
+ // Store xyz only here!!! To get the Delta z from the correction...
+ zBeforeCorr = xyz[2];
+
fTPCCorrection->CorrectPoint(xyz, cl->GetDetector());
- pIn.SetXYZ(xyz);
}
+
+ if ( fCreateT0seed )
+ xyz[2]=cl->GetTimeBin() + sign*( zBeforeCorr - xyz[2] )/GetVDrift();
+ // xyz[2]=cl->GetTimeBin();
+ pIn.SetXYZ(xyz);
+
// rotate the cluster to the local detector frame
track->Rotate(((cl->GetDetector()%18)*20+10)*TMath::DegToRad());
AliTrackPoint prot = pIn.Rotate(track->GetAlpha()); // rotate to the local frame - non distoted point
if (TMath::Abs(track->GetX())>kMaxR) break;
// if (TMath::Abs(track->GetZ())<kZcut)continue;
//
+
+ // add residuals
+ if (arrClustRes) {
+ TClonesArray &arrDummy=*arrClustRes;
+ AliTPCclusterMI *clRes = new(arrDummy[arrDummy.GetEntriesFast()]) AliTPCclusterMI(*cl);
+ clRes->SetX(prot.GetX());
+ // residuals in terms of sigma cl and track
+ clRes->SetY((track->GetY()-prot.GetY())/( sqrt ( prot.GetCov()[3] + track->GetSigmaY2()) ) );
+ clRes->SetZ((track->GetZ()-prot.GetZ())/( sqrt ( prot.GetCov()[5] + track->GetSigmaZ2()) ) );
+ }
+
Double_t pointPos[2]={0,0};
Double_t pointCov[3]={0,0,0};
pointPos[0]=prot.GetY();//local y
AliTrackerBase::PropagateTrackTo2(track,refX,kMass,5.,kTRUE,kMaxSnp,0,kFALSE,fUseMaterial);
- // rotate fittet track to the frame of the original track and propagate to same reference
- track->Rotate(tr->GetAlpha());
+ if (!fCreateT0seed || fForceAlpha){
+ // rotate fittet track to the frame of the original track and propagate to same reference
+ track->Rotate(tr->GetAlpha());
- AliTrackerBase::PropagateTrackTo2(track,refX,kMass,1.,kFALSE,kMaxSnp,0,kFALSE,fUseMaterial);
+ AliTrackerBase::PropagateTrackTo2(track,refX,kMass,1.,kFALSE,kMaxSnp,0,kFALSE,fUseMaterial);
+ }
return track;
}
-
//____________________________________________________________________________________
AliExternalTrackParam* AliToyMCReconstruction::GetFittedTrackFromSeedAllClusters(const AliToyMCTrack *tr, const AliExternalTrackParam *seed, Int_t &nClus)
{
}
}
- tFirst->GetListOfFriends()->Print();
+ if (tFirst->GetListOfFriends()) tFirst->GetListOfFriends()->Print();
return tFirst;
}
// a 'valid' position in z is needed for the seeding procedure
Double_t sign=1;
if (((sec/18)%2)==1) sign=-1;
- cl->SetZ(cl->GetTimeBin()*GetVDrift()*sign);
+ cl->SetZ(cl->GetTimeBin()*GetVDrift());
//mark cluster to be time*vDrift by setting the type to 1
cl->SetType(1);
// cl->SetZ(cl->GetTimeBin());
if (!fTree) return;
- TString debugName=fInputFile->GetName();
+ TString debugName=gSystem->BaseName(fInputFile->GetName());
debugName.ReplaceAll(".root","");
debugName.Append(Form(".%1d.%1d_%1d_%1d_%03d_%02d",
fUseMaterial,fIdealTracking,fClusterType,
for (Int_t i=0;i<from.GetN();++i) to.AddPoint(from.GetX()[i],from.GetY()[i],from.GetZ()[i],from.GetSy()[i],from.GetSz()[i]);
}
+//____________________________________________________________________________________
+Float_t AliToyMCReconstruction::FindClosestT0(const TVectorF &t0list, const TVectorF &z0list, AliExternalTrackParam &t0seed)
+{
+ //
+ // find closes T0 in a list of T0s
+ //
+
+ Long64_t size=t0list.GetNrows();
+ const Float_t *array=t0list.GetMatrixArray();
+
+ Float_t vDrift=GetVDrift();
+ Float_t zLength=GetZLength(0);
+
+ Float_t sign=(1-2*(t0seed.GetTgl()>0));
+
+ Float_t vtxCorr=0.;
+ Float_t t0=t0seed.GetZ()-zLength/vDrift;
+
+ Int_t index=0;
+ Float_t minDist=1e20;
+ for (Int_t it0=0; it0<size; ++it0) {
+ if (fUseZ0list) vtxCorr=sign*z0list[it0]/vDrift;
+// printf("vtxcorr %d: %.2g, %.2g\n",it0, vtxCorr, z0list[it0]);
+ const Float_t dist=TMath::Abs(t0list[it0]-t0-vtxCorr);
+ if (dist<minDist) {
+ index=it0;
+ minDist=dist;
+ }
+ }
+
+ return array[index];
+}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff