#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)
+, fRecoInfo(-1)
, fStreamer(0x0)
, fInputFile(0x0)
, fTree(0x0)
, fAllClusters("AliTPCclusterMI",10000)
, fMapTrackEvent(10000)
, fMapTrackTrackInEvent(10000)
+, fHnDelta(0x0)
, fIsAC(kFALSE)
{
//
ConnectInputFile(file, nmaxEv);
if (!fTree) return;
+ Int_t maxev=fTree->GetEntries();
+ if (nmaxEv>0&&nmaxEv<maxev) maxev=nmaxEv;
+
InitStreamer(".debug");
gROOT->cd();
- static AliExternalTrackParam dummySeedT0;
- static AliExternalTrackParam dummySeed;
- static AliExternalTrackParam dummyTrack;
+ static AliExternalTrackParam resetParam;
AliExternalTrackParam t0seed;
AliExternalTrackParam seed;
AliExternalTrackParam track;
- AliExternalTrackParam trackITS;
AliExternalTrackParam tOrig;
- AliExternalTrackParam tOrigITS;
+ // at ITS
+ AliExternalTrackParam tOrigITS; // ideal track
+ AliExternalTrackParam tRealITS; // ITS track with realistic space point resolution
+ AliExternalTrackParam trackITS; // TPC refitted track
+
+ //between TPC inner wall and ITS
+ AliExternalTrackParam tOrigITS1;
+ AliExternalTrackParam tRealITS1;
+ AliExternalTrackParam trackITS1;
+
+ //at TPC inner wall
+ AliExternalTrackParam tOrigITS2;
+ AliExternalTrackParam tRealITS2;
+ AliExternalTrackParam trackITS2;
+
AliExternalTrackParam *dummy;
+
+ // 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;
+ }
+ }
- Int_t maxev=fTree->GetEntries();
- if (nmaxEv>0&&nmaxEv<maxev) maxev=nmaxEv;
+ // 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;
+
+ // residuals
+ // binning r, phi, z, delta
+ const Int_t nbins=4;
+ Int_t bins[nbins] = {16, 18*5, 50, 80};
+ Double_t xmin[nbins] = {86. , 0., -250., -2.};
+ Double_t xmax[nbins] = {250., 2*TMath::Pi(), 250., 2.};
+ fHnDelta = new THnF("hn", "hn", nbins, bins, xmin, xmax);
+
+ // fill streamer?
+ Bool_t fillStreamer=(fStreamer!=0x0);
+ if (fRecoInfo>-1 && ((fRecoInfo&kFillNoTrackInfo)==kFillNoTrackInfo)) fillStreamer=kFALSE;
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();
+
+ // set SC scaling factor
+ fTPCCorrection->SetCorrScaleFactor(fEvent->GetSCscale());
+
for (Int_t itr=0; itr<fEvent->GetNumberOfTracks(); ++itr){
// printf(" > ====== Processing Track %6d ======== \n",itr);
const AliToyMCTrack *tr=fEvent->GetTrack(itr);
tOrig = *tr;
+ // propagate original track to ITS comparison points
+ if (fRecoInfo<0 || (fRecoInfo&kFillITS) ==kFillITS ) {
+ tOrigITS = *tr;
+ AliTrackerBase::PropagateTrackTo(&tOrigITS, lastLayerITS,kMass,1,kTRUE,kMaxSnp,0,kFALSE,fUseMaterial);
+ tRealITS = resetParam;
+ }
+ if (fRecoInfo<0 || (fRecoInfo&kFillITS1)==kFillITS1 ) {
+ tOrigITS1 = *tr;
+ AliTrackerBase::PropagateTrackTo(&tOrigITS1,betweeTPCITS,kMass,1,kTRUE,kMaxSnp,0,kFALSE,fUseMaterial);
+ tRealITS1 = resetParam;
+ }
+ if (fRecoInfo<0 || (fRecoInfo&kFillITS2)==kFillITS2 ) {
+ tOrigITS2 = *tr;
+ AliTrackerBase::PropagateTrackTo(&tOrigITS2,iFCRadius, kMass,1,kTRUE,kMaxSnp,0,kFALSE,fUseMaterial);
+ tRealITS2 = resetParam;
+ }
+
+ // realistic ITS track propagated to reference points
+ dummy = GetTrackRefit(tr,kITS);
+ if (dummy){
+ // propagate realistic track to ITS comparison points
+ if (fRecoInfo<0 || (fRecoInfo&kFillITS) ==kFillITS ) {
+ tRealITS = *dummy;
+ AliTrackerBase::PropagateTrackTo(&tRealITS, lastLayerITS,kMass,1,kTRUE,kMaxSnp,0,kFALSE,fUseMaterial);
+ }
+ if (fRecoInfo<0 || (fRecoInfo&kFillITS1)==kFillITS1 ) {
+ tRealITS1 = *dummy;
+ AliTrackerBase::PropagateTrackTo(&tRealITS1,betweeTPCITS,kMass,1,kTRUE,kMaxSnp,0,kFALSE,fUseMaterial);
+ }
+ if (fRecoInfo<0 || (fRecoInfo&kFillITS2)==kFillITS2 ) {
+ tRealITS2 = *dummy;
+ AliTrackerBase::PropagateTrackTo(&tRealITS2,iFCRadius, kMass,1,kTRUE,kMaxSnp,0,kFALSE,fUseMaterial);
+ }
+ //
+ delete dummy;
+ dummy=0x0;
+ }
- // set dummy
- t0seed = dummySeedT0;
- seed = dummySeed;
- track = dummyTrack;
-
- Float_t z0=fEvent->GetZ();
- Float_t t0=fEvent->GetT0();
+
+ // resetParam
+ t0seed = resetParam;
+ seed = resetParam;
+ track = resetParam;
+
+ if (fRecoInfo<0 || (fRecoInfo&kFillITS) ==kFillITS ) trackITS = resetParam;
+ if (fRecoInfo<0 || (fRecoInfo&kFillITS1)==kFillITS1 ) trackITS1 = resetParam;
+ if (fRecoInfo<0 || (fRecoInfo&kFillITS2)==kFillITS2 ) trackITS2 = resetParam;
+
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);
if (dummy) {
seed = *dummy;
delete dummy;
+ dummy=0x0;
// create fitted track
if (fDoTrackFit){
// printf("track\n");
- dummy = GetFittedTrackFromSeed(tr, &seed);
+ dummy = GetFittedTrackFromSeed(tr, &seed, arrClustRes);
track = *dummy;
+ dummy=0x0;
delete dummy;
}
-
- // Copy original track and fitted track
- // for extrapolation to ITS last layer
- tOrigITS = *tr;
- trackITS = track;
-
+
// propagate seed to 0
- const Double_t kMaxSnp = 0.85;
- const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
AliTrackerBase::PropagateTrackTo(&seed,0,kMass,5,kTRUE,kMaxSnp,0,kFALSE,fUseMaterial);
- // propagate original track to ITS last layer
- Double_t lastLayerITS = 43.0; // same as in AliToyMCEventGenerator::MakeITSClusters (hard coded)
- AliTrackerBase::PropagateTrackTo(&tOrigITS,lastLayerITS,kMass,1,kTRUE,kMaxSnp,0,kFALSE,fUseMaterial);
+ //
+ // ITS comparison
+ //
+
+ // rotate fitted track to the frame of the original track and propagate to same reference
+ if (fRecoInfo<0 || (fRecoInfo&kFillITS) ==kFillITS ){
+ trackITS = track;
+ 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
+ if (fRecoInfo<0 || (fRecoInfo&kFillITS1)==kFillITS1 ){
+ trackITS1 = track;
+ 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
+ if (fRecoInfo<0 || (fRecoInfo&kFillITS2)==kFillITS2 ){
+ trackITS2 = track;
+ 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);
+ }
}
}
Int_t ctype(fCorrectionType);
-
- if (fStreamer) {
+
+ if (fillStreamer){
(*fStreamer) << "Tracks" <<
"iev=" << iev <<
"z0=" << z0 <<
"t0=" << t0 <<
+ "lastt0=" << lastT0 <<
"fTime0=" << fTime0 <<
"itr=" << itr <<
"clsType=" << fClusterType <<
"zLength=" << zLength <<
"t0seed.=" << &t0seed <<
"seed.=" << &seed <<
- "track.=" << &track <<
+
"tOrig.=" << &tOrig <<
- "trackITS.=" << &trackITS <<
- "tOrigITS.=" << &tOrigITS <<
+ "track.=" << &track;
+
+
+ // ITS match
+ if (fRecoInfo<0 || (fRecoInfo&kFillITS) ==kFillITS ){
+ (*fStreamer) << "Tracks" <<
+ "tOrigITS.=" << &tOrigITS <<
+ "tRealITS.=" << &tRealITS <<
+ "trackITS.=" << &trackITS;
+ }
+
+ if (fRecoInfo<0 || (fRecoInfo&kFillITS1) ==kFillITS1 ){
+ (*fStreamer) << "Tracks" <<
+ "tOrigITS1.=" << &tOrigITS1 <<
+ "tRealITS1.=" << &tRealITS1 <<
+ "trackITS1.=" << &trackITS1;
+ }
+
+ if (fRecoInfo<0 || (fRecoInfo&kFillITS) ==kFillITS ){
+ (*fStreamer) << "Tracks" <<
+ "tOrigITS2.=" << &tOrigITS2 <<
+ "tRealITS2.=" << &tRealITS2 <<
+ "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;
}
+ fStreamer->GetFile()->cd();
+ fHnDelta->Write();
+
+ delete arrClustRes;
Cleanup();
}
// DumpTrackInfo(&seedsCentral2);
//dump clusters
- (*fStreamer) << "clusters" <<
- "cl.=" << &fAllClusters << "\n";
+// (*fStreamer) << "clusters" <<
+// "cl.=" << &fAllClusters << "\n";
Cleanup();
}
+//____________________________________________________________________________________
+AliExternalTrackParam* AliToyMCReconstruction::GetSeedFromTrackIdeal(const AliToyMCTrack * const tr, EDet det )
+{
+ //
+ // crate a seed from the track points of the respective detector
+ //
+ AliTrackPoint seedPoint[3];
+
+ Int_t npoints=0;
+ switch (det) {
+ case kITS:
+ npoints=tr->GetNumberOfITSPoints();
+ break;
+ case kTPC:
+ npoints=(fClusterType == 0)?tr->GetNumberOfSpacePoints():tr->GetNumberOfDistSpacePoints();
+ break;
+ case kTRD:
+ npoints=tr->GetNumberOfTRDPoints();
+ break;
+ }
+
+ if (npoints<3) return 0x0;
+
+ Int_t pos[3]={0,npoints/2,npoints-1};
+ const AliCluster *cl=0x0;
+
+ for (Int_t ipoint=0;ipoint<3;++ipoint){
+ Int_t cluster=pos[ipoint];
+ switch (det) {
+ case kITS:
+ seedPoint[ipoint]=(*tr->GetITSPoint(cluster));
+ break;
+ case kTPC:
+ cl=tr->GetSpacePoint(cluster);
+ if (fClusterType == 1) cl=tr->GetDistortedSpacePoint(cluster);
+ AliTPCclusterMI::SetGlobalTrackPoint(*cl,seedPoint[ipoint]);
+ break;
+ case kTRD:
+ seedPoint[ipoint]=(*tr->GetTRDPoint(cluster));
+ break;
+ }
+ }
+
+ AliExternalTrackParam *seed = AliTrackerBase::MakeSeed(seedPoint[0], seedPoint[1], seedPoint[2]);
+ seed->ResetCovariance(10);
+
+ return seed;
+}
+
//____________________________________________________________________________________
AliExternalTrackParam* AliToyMCReconstruction::GetSeedFromTrack(const AliToyMCTrack * const tr, Bool_t forceSeed)
{
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::GetTrackRefit(const AliToyMCTrack * const tr, EDet det)
+{
+ //
+ // Get the ITS or TRD track refitted from the toy track
+ // type: 0=ITS; 1=TRD
+ //
+
+ AliExternalTrackParam *track=GetSeedFromTrackIdeal(tr,det);
+ if (!track) return 0x0;
+
+ const Double_t kMaxSnp = 0.85;
+ const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
+
+ Int_t npoints=0;
+ switch (det) {
+ case kITS:
+ npoints=tr->GetNumberOfITSPoints();
+ break;
+ case kTPC:
+ npoints=(fClusterType == 0)?tr->GetNumberOfSpacePoints():tr->GetNumberOfDistSpacePoints();
+ break;
+ case kTRD:
+ npoints=tr->GetNumberOfTRDPoints();
+ break;
+ }
+
+ const AliCluster *cl=0x0;
+
+ for (Int_t ipoint=0; ipoint<npoints; ++ipoint) {
+ AliTrackPoint pIn;
+
+ switch (det) {
+ case kITS:
+ pIn=(*tr->GetITSPoint(ipoint));
+ break;
+ case kTPC:
+ cl=tr->GetSpacePoint(ipoint);
+ if (fClusterType == 1) cl=tr->GetDistortedSpacePoint(ipoint);
+ AliTPCclusterMI::SetGlobalTrackPoint(*cl,pIn);
+ break;
+ case kTRD:
+ pIn=(*tr->GetTRDPoint(ipoint));
+ break;
+ }
+
+
+ const Double_t angle=pIn.GetAngle();
+ track->Rotate(angle);
+ AliTrackPoint prot = pIn.Rotate(track->GetAlpha()); // rotate to the local frame - non distoted point
+
+ if (!AliTrackerBase::PropagateTrackTo(track,prot.GetX(),kMass,5,kFALSE,kMaxSnp,0,kFALSE,fUseMaterial)) {
+ AliInfo(Form("Could not propagate track to x=%.2f (a=%.2f) for det %d",prot.GetX(),angle,det));
+ }
+ //
+
+ Double_t pointPos[2]={0,0};
+ Double_t pointCov[3]={0,0,0};
+ pointPos[0]=prot.GetY();//local y
+ pointPos[1]=prot.GetZ();//local z
+ pointCov[0]=prot.GetCov()[3];//simay^2
+ pointCov[1]=prot.GetCov()[4];//sigmayz
+ pointCov[2]=prot.GetCov()[5];//sigmaz^2
+
+ if (!track->Update(pointPos,pointCov)) {
+ AliInfo(Form("no update: det: %d",det));
+ break;
+ }
+
+ }
+
+ return track;
+}
+
//____________________________________________________________________________________
void AliToyMCReconstruction::SetTrackPointFromCluster(const AliTPCclusterMI *cl, AliTrackPoint &p )
// AliTrackPoint *tp=const_cast<AliTPCclusterMI*>(cl)->MakePoint(p);
// p=*tp;
// delete tp;
- const_cast<AliTPCclusterMI*>(cl)->MakePoint(p);
+// const_cast<AliTPCclusterMI*>(cl)->MakePoint(p);
+ AliTPCclusterMI::SetGlobalTrackPoint(*cl,p);
// cl->Print();
// p.Print();
p.SetVolumeID(cl->GetDetector());
}
//____________________________________________________________________________________
-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);
+ // track copy for propagation
+ AliExternalTrackParam trCopy(*tr);
Int_t ncls=(fClusterType == 0)?tr->GetNumberOfSpacePoints():tr->GetNumberOfDistSpacePoints();
const Double_t refX = tr->GetX();
const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
+
+ // parametrised track resolution
+ Double_t trackRes=gRandom->Gaus();
// loop over all other points and add to the track
for (Int_t ipoint=ncls-1; ipoint>=0; --ipoint){
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()) ) );
+ }
+
+ // fill cluster residuals to ideal track for calibration studies
+ // ideal cluster position
+ // require at least 2 TRD points
+ if (fRecoInfo<0 || (fRecoInfo&kFillDeltas) ==kFillDeltas ) {
+ trCopy.Rotate(track->GetAlpha());
+ AliTrackerBase::PropagateTrackTo(&trCopy,prot.GetX(),kMass,5,kFALSE,kMaxSnp,0,kFALSE,fUseMaterial);
+ // binning r, phi, z, delta (0=rphi, 1=z)
+ // resolution parametrisation
+ Float_t soneOverPt= trCopy.GetSigned1Pt();
+ Float_t oneOverPt = TMath::Abs(soneOverPt);
+ Float_t radius = trCopy.GetX();
+ Float_t trackY = trCopy.GetY();
+ Float_t trackZ = trCopy.GetZ();
+ Float_t trackPhi = trCopy.Phi();
+ Float_t alpha = trCopy.GetAlpha();
+
+ Float_t pointY = prot.GetY();
+ Float_t pointZ = prot.GetZ();
+
+ Float_t resRphi = 0.004390 + oneOverPt*(-0.136403) + oneOverPt*radius*(0.002266) + oneOverPt*radius*radius*(-0.000006);
+
+ Float_t resRphiRandom = resRphi*trackRes;
+ Float_t deviation = trackY+resRphiRandom-pointY;
+ Short_t npTRD = tr->GetNumberOfTRDPoints();
+
+ // rphi residuals
+ Double_t xx[4]={radius, trackPhi, trackZ ,deviation};
+ if (npTRD>=2){
+ fHnDelta->Fill(xx);
+ }
+
+ Short_t event=fTree->GetReadEntry();
+
+ if (fStreamer) {
+ (*fStreamer) << "delta" <<
+ "soneOverPt=" << soneOverPt <<
+ "r=" << radius <<
+ "trackPhi=" << trackPhi <<
+ "trackY=" << trackY <<
+ "trackZ=" << trackZ <<
+ "alpha=" << alpha <<
+ "resRphi=" << resRphi <<
+ "trackRes=" << trackRes <<
+ "pointY=" << pointY <<
+ "pointZ=" << pointZ <<
+ "npTRD=" << npTRD <<
+ "event=" << event <<
+ "\n";
+// "point.=" << &prot <<
+// "track.=" << track <<
+ }
+ }
+
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)
{
// Init the space charge map
//
- TString filename="$ALICE_ROOT/TPC/Calib/maps/SC_NeCO2_eps5_50kHz_precal.root";
+// TString filename="$ALICE_ROOT/TPC/Calib/maps/SC_NeCO2_eps5_50kHz_precal.root";
+ TString filename;
if (fTree) {
TList *l=fTree->GetUserInfo();
for (Int_t i=0; i<l->GetEntries(); ++i) {
- TString s(l->At(i)->GetName());
- if (s.Contains("SC_")) {
- filename=s;
- break;
+ TObject *o=l->At(i);
+ if (o->IsA() == TObjString::Class()){
+ TString s(o->GetName());
+ if (s.Contains("lookup.root")) {
+ filename=s;
+ break;
+ }
}
}
}
- printf("Initialising the space charge map using the file: '%s'\n",filename.Data());
+ if (filename.IsNull()) {
+ AliFatal("No SC map provided in the Userinfo of the simulation tree");
+ return;
+ }
+
+ AliInfo(Form("Initialising the space charge map using the file: '%s'\n",filename.Data()));
TFile f(filename.Data());
fTPCCorrection=(AliTPCCorrection*)f.Get("map");
}
}
- tFirst->GetListOfFriends()->Print();
+ if (tFirst->GetListOfFriends()) tFirst->GetListOfFriends()->Print();
return tFirst;
}
//track->Local2GlobalPosition(xd,track->GetAlpha());
// use TPCParams to get the sector number
- Float_t xyz[3] = {xd[0],xd[1],xd[2]};
+ Float_t xyz[3] = {static_cast<Float_t>(xd[0]),static_cast<Float_t>(xd[1]),static_cast<Float_t>(xd[2])};
Int_t i[3] = {0,0,0};
if(fTPCParam){
sector = fTPCParam->Transform0to1(xyz,i);
// 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,
delete fEvent;
fEvent = 0x0;
-
+
+ delete fHnDelta;
+ fHnDelta=0x0;
// delete fTree;
fTree=0x0;
AliExternalTrackParam tOrigITS(*toyTrack);
// propagate original track to ITS last layer
- Double_t lastLayerITS = 43.0; // same as in AliToyMCEventGenerator::MakeITSClusters (hard coded)
+// Double_t lastLayerITS = 43.0; // same as in AliToyMCEventGenerator::MakeITSClusters (hard coded)
+ const Double_t iFCRadius = 83.5; //radius constants found in AliTPCCorrection.cxx
+ Double_t lastLayerITS = iFCRadius; // its track propgated to inner TPC wall
AliTrackerBase::PropagateTrackTo(&tOrigITS,lastLayerITS,kMass,1,kTRUE,kMaxSnp,0,kFALSE,fUseMaterial);
AliExternalTrackParam dummyParam;
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