X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;ds=sidebyside;f=TRD%2FAliTRDtracker.cxx;h=758d31ef2e98e2ce7383566e02668a5946d99bdd;hb=a4cf6278ec27c294e6ccd953b999c54f04585048;hp=03e0db3d4a375721eef4ead28e25df1f8c33f50e;hpb=eb37e0af138ccdf19c7527be8e99dc1dfc379590;p=u%2Fmrichter%2FAliRoot.git diff --git a/TRD/AliTRDtracker.cxx b/TRD/AliTRDtracker.cxx index 03e0db3d4a3..758d31ef2e9 100644 --- a/TRD/AliTRDtracker.cxx +++ b/TRD/AliTRDtracker.cxx @@ -26,20 +26,17 @@ // // /////////////////////////////////////////////////////////////////////////////// - -#include #include #include #include -#include -#include #include #include #include #include #include -#include "AliESD.h" +#include "AliESDEvent.h" +#include "AliESDtrack.h" #include "AliAlignObj.h" #include "AliRieman.h" #include "AliTrackPointArray.h" @@ -61,12 +58,22 @@ ClassImp(AliTRDtracker) const Float_t AliTRDtracker::fgkMinClustersInTrack = 0.5; // const Float_t AliTRDtracker::fgkLabelFraction = 0.8; // const Double_t AliTRDtracker::fgkMaxChi2 = 12.0; // -const Double_t AliTRDtracker::fgkMaxSnp = 0.95; // Corresponds to tan = 3 +const Double_t AliTRDtracker::fgkMaxSnp = 0.95; // Maximum local sine of the azimuthal angle const Double_t AliTRDtracker::fgkMaxStep = 2.0; // Maximal step size in propagation //_____________________________________________________________________________ AliTRDtracker::AliTRDtracker() :AliTracker() + ,fGeom(0) + ,fNclusters(0) + ,fClusters(0) + ,fNseeds(0) + ,fSeeds(0) + ,fNtracks(0) + ,fTracks(0) + ,fTimeBinsPerPlane(0) + ,fAddTRDseeds(kFALSE) + ,fNoTilt(kFALSE) ,fHBackfit(0x0) ,fHClSearch(0x0) ,fHRefit(0x0) @@ -79,16 +86,6 @@ AliTRDtracker::AliTRDtracker() ,fHMinD(0x0) ,fHDeltaX(0x0) ,fHXCl(0x0) - ,fGeom(0) - ,fNclusters(0) - ,fClusters(0) - ,fNseeds(0) - ,fSeeds(0) - ,fNtracks(0) - ,fTracks(0) - ,fTimeBinsPerPlane(0) - ,fAddTRDseeds(kFALSE) - ,fNoTilt(kFALSE) ,fDebugStreamer(0) { // @@ -98,11 +95,6 @@ AliTRDtracker::AliTRDtracker() for (Int_t i = 0; i < kTrackingSectors; i++) { fTrSec[i] = 0; } - for (Int_t j = 0; j < 5; j++) { - for (Int_t k = 0; k < 18; k++) { - fHoles[j][k] = kFALSE; - } - } InitLogHists(); @@ -111,6 +103,16 @@ AliTRDtracker::AliTRDtracker() //_____________________________________________________________________________ AliTRDtracker::AliTRDtracker(const AliTRDtracker &t) :AliTracker(t) + ,fGeom(0) + ,fNclusters(0) + ,fClusters(0) + ,fNseeds(0) + ,fSeeds(0) + ,fNtracks(0) + ,fTracks(0) + ,fTimeBinsPerPlane(0) + ,fAddTRDseeds(kFALSE) + ,fNoTilt(kFALSE) ,fHBackfit(0x0) ,fHClSearch(0x0) ,fHRefit(0x0) @@ -123,16 +125,6 @@ AliTRDtracker::AliTRDtracker(const AliTRDtracker &t) ,fHMinD(0x0) ,fHDeltaX(0x0) ,fHXCl(0x0) - ,fGeom(0) - ,fNclusters(0) - ,fClusters(0) - ,fNseeds(0) - ,fSeeds(0) - ,fNtracks(0) - ,fTracks(0) - ,fTimeBinsPerPlane(0) - ,fAddTRDseeds(kFALSE) - ,fNoTilt(kFALSE) ,fDebugStreamer(0) { // @@ -142,8 +134,18 @@ AliTRDtracker::AliTRDtracker(const AliTRDtracker &t) } //_____________________________________________________________________________ -AliTRDtracker::AliTRDtracker(const TFile *geomfile) +AliTRDtracker::AliTRDtracker(const TFile */*geomfile*/) :AliTracker() + ,fGeom(0) + ,fNclusters(0) + ,fClusters(new TObjArray(2000)) + ,fNseeds(0) + ,fSeeds(new TObjArray(2000)) + ,fNtracks(0) + ,fTracks(new TObjArray(1000)) + ,fTimeBinsPerPlane(0) + ,fAddTRDseeds(kFALSE) + ,fNoTilt(kFALSE) ,fHBackfit(0x0) ,fHClSearch(0x0) ,fHRefit(0x0) @@ -156,16 +158,6 @@ AliTRDtracker::AliTRDtracker(const TFile *geomfile) ,fHMinD(0x0) ,fHDeltaX(0x0) ,fHXCl(0x0) - ,fGeom(0) - ,fNclusters(0) - ,fClusters(new TObjArray(2000)) - ,fNseeds(0) - ,fSeeds(new TObjArray(2000)) - ,fNtracks(0) - ,fTracks(new TObjArray(1000)) - ,fTimeBinsPerPlane(0) - ,fAddTRDseeds(kFALSE) - ,fNoTilt(kFALSE) ,fDebugStreamer(0) { // @@ -173,34 +165,15 @@ AliTRDtracker::AliTRDtracker(const TFile *geomfile) // TDirectory *savedir = gDirectory; - TFile *in = (TFile *) geomfile; - - if (!in->IsOpen()) { - AliWarning("geometry file is not open!\n"); - AliWarning("FULL TRD geometry and DEFAULT TRD parameter will be used\n"); - } - else { - in->cd(); - fGeom = (AliTRDgeometry *) in->Get("TRDgeometry"); - } - - if (!fGeom) { - AliWarning("Cannot find TRD geometry!\n"); - fGeom = new AliTRDgeometry(); - } - fGeom->ReadGeoMatrices(); - savedir->cd(); + fGeom = new AliTRDgeometry(); for (Int_t geomS = 0; geomS < kTrackingSectors; geomS++) { Int_t trS = geomS; fTrSec[trS] = new AliTRDtrackingSector(fGeom,geomS); - for (Int_t icham = 0; icham < AliTRDgeometry::kNcham; icham++) { - fHoles[icham][trS] = fGeom->IsHole(0,icham,geomS); - } } - AliTRDpadPlane *padPlane = AliTRDCommonParam::Instance()->GetPadPlane(0,0); + AliTRDpadPlane *padPlane = fGeom->GetPadPlane(0,0); Float_t tiltAngle = TMath::Abs(padPlane->GetTiltingAngle()); if (tiltAngle < 0.1) { fNoTilt = kTRUE; @@ -241,7 +214,9 @@ AliTRDtracker::~AliTRDtracker() delete fSeeds; } - delete fGeom; + if (fGeom) { + delete fGeom; + } for (Int_t geomS = 0; geomS < kTrackingSectors; geomS++) { delete fTrSec[geomS]; @@ -340,157 +315,6 @@ Int_t AliTRDtracker::GlobalToLocalID(Int_t gid) } -//_____________________________________________________________________________ -Bool_t AliTRDtracker::Transform(AliTRDcluster *cluster) -{ - // - // Transform from cluster system to tracking system - // - - // Magic constants for geo manager transformation - const Double_t kX0shift = 2.52; - - // - // Apply alignment and calibration to transform cluster - // - Int_t detector = cluster->GetDetector(); - Int_t plane = fGeom->GetPlane(cluster->GetDetector()); - Int_t chamber = fGeom->GetChamber(cluster->GetDetector()); - Int_t sector = fGeom->GetSector(cluster->GetDetector()); - - Double_t dxAmp = (Double_t) fGeom->CamHght(); // Amplification region - Double_t driftX = TMath::Max(cluster->GetX()-dxAmp*0.5,0.0); // Drift distance - - // - // ExB correction - // - Double_t vdrift = AliTRDcalibDB::Instance()->GetVdrift(cluster->GetDetector(),0,0); - Double_t exB = AliTRDcalibDB::Instance()->GetOmegaTau(vdrift,-AliTracker::GetBz()*0.1); - - AliTRDCommonParam *commonParam = AliTRDCommonParam::Instance(); - AliTRDpadPlane *padPlane = commonParam->GetPadPlane(plane,chamber); - Double_t zshiftIdeal = 0.5*(padPlane->GetRow0()+padPlane->GetRowEnd()); - Double_t localPos[3]; - Double_t localPosTracker[3]; - localPos[0] = -cluster->GetX(); - localPos[1] = cluster->GetY() - driftX * exB; - localPos[2] = cluster->GetZ() - zshiftIdeal; - - cluster->SetY(cluster->GetY() - driftX*exB); - Double_t xplane = (Double_t) AliTRDgeometry::GetTime0(plane); - cluster->SetX(xplane - cluster->GetX()); - - TGeoHMatrix *matrix = fGeom->GetCorrectionMatrix(cluster->GetDetector()); - if (!matrix) { - // No matrix found - if somebody used geometry with holes - AliError("Invalid Geometry - Default Geometry used\n"); - return kTRUE; - } - matrix->LocalToMaster(localPos,localPosTracker); - - if (AliTRDReconstructor::StreamLevel() > 1) { - (* fDebugStreamer) << "Transform" - << "Cl.=" << cluster - << "matrix.=" << matrix - << "Detector=" << detector - << "Sector=" << sector - << "Plane=" << plane - << "Chamber=" << chamber - << "lx0=" << localPosTracker[0] - << "ly0=" << localPosTracker[1] - << "lz0=" << localPosTracker[2] - << "\n"; - } - - cluster->SetX(localPosTracker[0]+kX0shift); - cluster->SetY(localPosTracker[1]); - cluster->SetZ(localPosTracker[2]); - - return kTRUE; - -} - -//_____________________________________________________________________________ -// Bool_t AliTRDtracker::Transform(AliTRDcluster *cluster) -//{ -// // -// // Is this still needed ???? -// // -// const Double_t kDriftCorrection = 1.01; // drift coeficient correction -// const Double_t kTime0Cor = 0.32; // time0 correction -// // -// const Double_t kX0shift = 2.52; -// const Double_t kX0shift5 = 3.05; - -// // -// // apply alignment and calibration to transform cluster -// // -// // -// Int_t detector = cluster->GetDetector(); -// Int_t plane = fGeom->GetPlane(cluster->GetDetector()); -// Int_t chamber = fGeom->GetChamber(cluster->GetDetector()); -// Int_t sector = fGeom->GetSector(cluster->GetDetector()); - -// Double_t dxAmp = (Double_t) fGeom->CamHght(); // Amplification region -// Double_t driftX = TMath::Max(cluster->GetX()-dxAmp*0.5,0.); // drift distance -// // -// // ExB correction -// // -// Double_t vdrift = AliTRDcalibDB::Instance()->GetVdrift(cluster->GetDetector(),0,0); -// Double_t exB = AliTRDcalibDB::Instance()->GetOmegaTau(vdrift,-AliTracker::GetBz()*0.1); -// // - -// AliTRDCommonParam* commonParam = AliTRDCommonParam::Instance(); -// AliTRDpadPlane * padPlane = commonParam->GetPadPlane(plane,chamber); -// Double_t zshiftIdeal = 0.5*(padPlane->GetRow0()+padPlane->GetRowEnd()); -// Double_t localPos[3], globalPos[3], localPosTracker[3], localPosTracker2[3]; -// localPos[2] = -cluster->GetX(); -// localPos[0] = cluster->GetY() - driftX*exB; -// localPos[1] = cluster->GetZ() -zshiftIdeal; -// TGeoHMatrix * matrix = fGeom->GetGeoMatrix(cluster->GetDetector()); -// matrix->LocalToMaster(localPos, globalPos); - -// Double_t sectorAngle = 20.*(sector%18)+10; -// TGeoHMatrix rotSector; -// rotSector.RotateZ(sectorAngle); -// rotSector.LocalToMaster(globalPos, localPosTracker); -// // -// // -// TGeoHMatrix matrix2(*matrix); -// matrix2.MultiplyLeft(&rotSector); -// matrix2.LocalToMaster(localPos,localPosTracker2); -// // -// // -// // -// cluster->SetY(cluster->GetY() - driftX*exB); -// Double_t xplane = (Double_t) AliTRDgeometry::GetTime0(plane); -// cluster->SetX(xplane- kDriftCorrection*(cluster->GetX()-kTime0Cor)); -// (*fDebugStreamer)<<"Transform"<< -// "Cl.="<SetX(localPosTracker[0]+kX0shift5); -// else -// cluster->SetX(localPosTracker[0]+kX0shift); - -// cluster->SetY(localPosTracker[1]); -// cluster->SetZ(localPosTracker[2]); -// return kTRUE; -// } - //_____________________________________________________________________________ Bool_t AliTRDtracker::AdjustSector(AliTRDtrack *track) { @@ -502,18 +326,12 @@ Bool_t AliTRDtracker::AdjustSector(AliTRDtrack *track) Double_t y = track->GetY(); Double_t ymax = track->GetX()*TMath::Tan(0.5*alpha); - // Is this still needed ???? - //Int_t ns = AliTRDgeometry::kNsect; - //Int_t s=Int_t(track->GetAlpha()/alpha)%ns; - if (y > ymax) { - //s = (s+1) % ns; if (!track->Rotate( alpha)) { return kFALSE; } } else if (y < -ymax) { - //s = (s-1+ns) % ns; if (!track->Rotate(-alpha)) { return kFALSE; } @@ -583,7 +401,7 @@ Int_t AliTRDtracker::GetLastPlane(AliTRDtrack *track) } //_____________________________________________________________________________ -Int_t AliTRDtracker::Clusters2Tracks(AliESD *event) +Int_t AliTRDtracker::Clusters2Tracks(AliESDEvent *event) { // // Finds tracks within the TRD. The ESD event is expected to contain seeds @@ -642,7 +460,7 @@ Int_t AliTRDtracker::Clusters2Tracks(AliESD *event) } //_____________________________________________________________________________ -Int_t AliTRDtracker::PropagateBack(AliESD *event) +Int_t AliTRDtracker::PropagateBack(AliESDEvent *event) { // // Gets seeds from ESD event. The seeds are AliTPCtrack's found and @@ -655,70 +473,51 @@ Int_t AliTRDtracker::PropagateBack(AliESD *event) Int_t found = 0; // number of tracks found Float_t foundMin = 20.0; - Int_t n = event->GetNumberOfTracks(); - // Sort tracks - Float_t *quality = new Float_t[n]; - Int_t *index = new Int_t[n]; - for (Int_t i = 0; i < n; i++) { - AliESDtrack *seed = event->GetTrack(i); + // Sort tracks according to covariance of local Y and Z + Int_t nSeed = event->GetNumberOfTracks(); + Float_t *quality = new Float_t[nSeed]; + Int_t *index = new Int_t[nSeed]; + for (Int_t iSeed = 0; iSeed < nSeed; iSeed++) { + AliESDtrack *seed = event->GetTrack(iSeed); Double_t covariance[15]; seed->GetExternalCovariance(covariance); - quality[i] = covariance[0]+covariance[2]; - //quality[i] = covariance[0]; + quality[iSeed] = covariance[0] + covariance[2]; } - TMath::Sort(n,quality,index,kFALSE); + TMath::Sort(nSeed,quality,index,kFALSE); - for (Int_t i = 0; i < n; i++) { + // Backpropagate all seeds + for (Int_t iSeed = 0; iSeed < nSeed; iSeed++) { - //AliESDtrack *seed = event->GetTrack(i); - AliESDtrack *seed = event->GetTrack(index[i]); - fHBackfit->Fill(0); + // Get the seeds in sorted sequence + AliESDtrack *seed = event->GetTrack(index[iSeed]); + fHBackfit->Fill(0); // All seeds + // Check the seed status ULong_t status = seed->GetStatus(); if ((status & AliESDtrack::kTPCout) == 0) { - fHBackfit->Fill(1); + fHBackfit->Fill(1); // TPC outer edge reached continue; } - if ((status & AliESDtrack::kTRDout) != 0) { - fHBackfit->Fill(2); + fHBackfit->Fill(2); // TRD outer edge reached (does this happen ?) continue; } - Int_t lbl = seed->GetLabel(); - AliTRDtrack *track = new AliTRDtrack(*seed); + // Do the back prolongation + Int_t lbl = seed->GetLabel(); + AliTRDtrack *track = new AliTRDtrack(*seed); track->SetSeedLabel(lbl); seed->UpdateTrackParams(track,AliESDtrack::kTRDbackup); // Make backup fNseeds++; - Float_t p4 = track->GetC(); - Int_t expectedClr = FollowBackProlongation(*track); - - fHBackfit->Fill(3); + Float_t p4 = track->GetC(); + Int_t expectedClr = FollowBackProlongation(*track); + fHBackfit->Fill(3); // Back prolongation done fHX->Fill(track->GetX()); + // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - // store the last measurement - /* - fHNClTrack->Fill(track->GetNumberOfClusters()); - if (track->GetNumberOfClusters() >= foundMin) { - - fHBackfit->Fill(4); - track->CookdEdx(); - CookdEdxTimBin(*track); - CookLabel(track,1 - fgkLabelFraction); - if (track->GetBackupTrack()) { - //fHBackfit->Fill(5); - UseClusters(track->GetBackupTrack()); - seed->UpdateTrackParams(track->GetBackupTrack(),AliESDtrack::kTRDbackup); - } - } - */ - - /**/ - // inter-tracks competition ??? if ((TMath::Abs(track->GetC() - p4) / TMath::Abs(p4) < 0.2) || - (TMath::Abs(track->GetPt()) > 0.8)) { + (track->Pt() > 0.8)) { fHBackfit->Fill(4); @@ -729,7 +528,7 @@ Int_t AliTRDtracker::PropagateBack(AliESD *event) Int_t foundClr = track->GetNumberOfClusters(); if (foundClr >= foundMin) { track->CookdEdx(); - CookdEdxTimBin(*track); + track->CookdEdxTimBin(seed->GetID()); // A.Bercuci 25.07.07 CookLabel(track,1 - fgkLabelFraction); if (track->GetBackupTrack()) { UseClusters(track->GetBackupTrack()); @@ -738,7 +537,7 @@ Int_t AliTRDtracker::PropagateBack(AliESD *event) // Sign only gold tracks if (track->GetChi2() / track->GetNumberOfClusters() < 4) { if ((seed->GetKinkIndex(0) == 0) && - (TMath::Abs(track->GetPt()) < 1.5)) { + (track->Pt() < 1.5)) { UseClusters(track); } } @@ -924,10 +723,11 @@ Int_t AliTRDtracker::PropagateBack(AliESD *event) SaveLogHists(); return 0; + } //_____________________________________________________________________________ -Int_t AliTRDtracker::RefitInward(AliESD *event) +Int_t AliTRDtracker::RefitInward(AliESDEvent *event) { // // Refits tracks within the TRD. The ESD event is expected to contain seeds @@ -937,10 +737,11 @@ Int_t AliTRDtracker::RefitInward(AliESD *event) // Origin: Thomas KUHR (Thomas.Kuhr@cern.ch) // - Int_t timeBins = fTrSec[0]->GetNumberOfTimeBins(); - Float_t foundMin = fgkMinClustersInTrack * timeBins; + //Int_t timeBins = fTrSec[0]->GetNumberOfTimeBins(); + //Float_t foundMin = fgkMinClustersInTrack * timeBins; Int_t nseed = 0; Int_t found = 0; + Int_t pidQ = 0; //Int_t innerTB = fTrSec[0]->GetInnerTimeBin(); AliTRDtrack seed2; @@ -989,19 +790,15 @@ Int_t AliTRDtracker::RefitInward(AliESD *event) indexes3[i] = indexes2[i]; } - //AliTRDtrack *pt = seed2; - AliTRDtrack &t = *pt; - FollowProlongation(t); - if (t.GetNumberOfClusters() >= foundMin) { - //UseClusters(&t); - //CookLabel(pt, 1-fgkLabelFraction); - t.CookdEdx(); - CookdEdxTimBin(t); - } + FollowProlongation(*pt); + pt->CookdEdx(); + pt->CookdEdxTimBin(seed->GetID()); + pt->SetPIDMethod(AliTRDtrack::kLQ); //switch between TRD PID methods + pt->CookPID(pidQ); found++; Double_t xTPC = 250.0; - if (PropagateToX(t,xTPC,fgkMaxStep)) { + if (PropagateToX(*pt,xTPC,fgkMaxStep)) { seed->UpdateTrackParams(pt,AliESDtrack::kTRDrefit); fHRefit->Fill(5); @@ -1012,17 +809,21 @@ Int_t AliTRDtracker::RefitInward(AliESD *event) } seed->SetTRDTimBin(pt->GetPIDTimBin(i),i); } + } else { + // If not prolongation to TPC - propagate without update fHRefit->Fill(5); AliTRDtrack *seed2 = new AliTRDtrack(*seed); seed2->ResetCovariance(5.0); AliTRDtrack *pt2 = new AliTRDtrack(*seed2,seed2->GetAlpha()); delete seed2; + if (PropagateToX(*pt2,xTPC,fgkMaxStep)) { - pt2->CookdEdx( ); - CookdEdxTimBin(*pt2); + + pt2->CookdEdx(); + pt2->CookdEdxTimBin(seed->GetID()); seed->UpdateTrackParams(pt2,AliESDtrack::kTRDrefit); fHRefit->Fill(6); @@ -1032,10 +833,22 @@ Int_t AliTRDtracker::RefitInward(AliESD *event) } seed->SetTRDTimBin(pt2->GetPIDTimBin(i),i); } + + } + + // Add TRD track to ESDfriendTrack - maybe this tracks are + // not useful for post-processing - TODO make decision + if (AliTRDReconstructor::StreamLevel() > 0) { + seed->AddCalibObject(new AliTRDtrack(*pt2/*, kTRUE*/)); } delete pt2; + } + // Add TRD track to ESDfriendTrack + if (AliTRDReconstructor::StreamLevel() > 0) { + seed->AddCalibObject(new AliTRDtrack(*pt/*, kTRUE*/)); + } delete pt; } @@ -1044,7 +857,9 @@ Int_t AliTRDtracker::RefitInward(AliESD *event) AliInfo(Form("Number of found tracks from loaded seeds: %d",found)); SaveLogHists(); + return 0; + } //_____________________________________________________________________________ @@ -1060,8 +875,8 @@ Int_t AliTRDtracker::FollowProlongation(AliTRDtrack &t) Int_t sector; Int_t lastplane = GetLastPlane(&t); - Double_t radLength = 0.0; - Double_t rho = 0.0; + Double_t xx0 = 0.0; + Double_t xrho= 0.0; Int_t expectedNumberOfClusters = 0; for (Int_t iplane = lastplane; iplane >= 0; iplane--) { @@ -1069,14 +884,12 @@ Int_t AliTRDtracker::FollowProlongation(AliTRDtrack &t) Int_t row0 = GetGlobalTimeBin(0,iplane,GetTimeBinsPerPlane()-1); Int_t rowlast = GetGlobalTimeBin(0,iplane,0); - // // Propagate track close to the plane if neccessary - // Double_t currentx = fTrSec[0]->GetLayer(rowlast)->GetX(); if (currentx < (-fgkMaxStep + t.GetX())) { - // Propagate closer to chamber - safety space fgkMaxStep + // Propagate closer to chamber - safety space fgkMaxStep if (!PropagateToX(t,currentx+fgkMaxStep,fgkMaxStep)) { - break; + break; } } @@ -1084,9 +897,7 @@ Int_t AliTRDtracker::FollowProlongation(AliTRDtrack &t) break; } - // // Get material budget - // Double_t xyz0[3]; Double_t xyz1[3]; Double_t param[7]; @@ -1095,31 +906,41 @@ Int_t AliTRDtracker::FollowProlongation(AliTRDtrack &t) Double_t z; // Starting global position - t.GetXYZ(xyz0); + t.GetXYZ(xyz0); // End global position x = fTrSec[0]->GetLayer(row0)->GetX(); if (!t.GetProlongation(x,y,z)) { break; } - xyz1[0] = x * TMath::Cos(t.GetAlpha()) - y * TMath::Sin(t.GetAlpha()); + xyz1[0] = x * TMath::Cos(t.GetAlpha()) - y * TMath::Sin(t.GetAlpha()); xyz1[1] = +x * TMath::Sin(t.GetAlpha()) + y * TMath::Cos(t.GetAlpha()); xyz1[2] = z; - AliKalmanTrack::MeanMaterialBudget(xyz0,xyz1,param); - rho = param[0]; - radLength = param[1]; // Get mean propagation parameters + AliTracker::MeanMaterialBudget(xyz0,xyz1,param); + xrho= param[0]*param[4]; + xx0 = param[1]; // Get mean propagation parameters + + // Flags for marking the track momentum and direction. The track is + // marked only if it has at least 1 cluster picked up in the current + // chamber. + Bool_t kUPDATED = kFALSE; + Bool_t kMARKED = kFALSE; - // // Propagate and update - // sector = t.GetSector(); //for (Int_t itime=GetTimeBinsPerPlane()-1;itime>=0;itime--) { for (Int_t itime = 0 ; itime < GetTimeBinsPerPlane(); itime++) { + // Mark track kinematics + if (itime > 10 && kUPDATED && !kMARKED) { + t.SetTrackSegmentDirMom(iplane); + kMARKED = kTRUE; + } + Int_t ilayer = GetGlobalTimeBin(0,iplane,itime); - expectedNumberOfClusters++; + expectedNumberOfClusters++; t.SetNExpected(t.GetNExpected() + 1); if (t.GetX() > 345.0) { - t.SetNExpectedLast(t.GetNExpectedLast() + 1); + t.SetNExpectedLast(t.GetNExpectedLast() + 1); } AliTRDpropagationLayer &timeBin = *(fTrSec[sector]->GetLayer(ilayer)); AliTRDcluster *cl = 0; @@ -1131,60 +952,64 @@ Int_t AliTRDtracker::FollowProlongation(AliTRDtrack &t) AliTRDcluster *cl0 = timeBin[0]; if (!cl0) { - // No clusters in given time bin - continue; + // No clusters in given time bin + continue; } - Int_t plane = fGeom->GetPlane(cl0->GetDetector()); + Int_t plane = fGeom->GetPlane(cl0->GetDetector()); if (plane > lastplane) { - continue; + continue; } Int_t timebin = cl0->GetLocalTimeBin(); AliTRDcluster *cl2 = GetCluster(&t,plane,timebin,index); if (cl2) { - - cl = cl2; + cl = cl2; //Double_t h01 = GetTiltFactor(cl); //I.B's fix //maxChi2=t.GetPredictedChi2(cl,h01); + } - } - if (cl) { + if (cl) { - //if (cl->GetNPads()<5) + //if (cl->GetNPads()<5) Double_t dxsample = timeBin.GetdX(); - t.SetSampledEdx(TMath::Abs(cl->GetQ()/dxsample)); - Double_t h01 = GetTiltFactor(cl); - Int_t det = cl->GetDetector(); + t.SetSampledEdx(TMath::Abs(cl->GetQ()/dxsample)); + Double_t h01 = GetTiltFactor(cl); + Int_t det = cl->GetDetector(); Int_t plane = fGeom->GetPlane(det); + if (t.GetX() > 345.0) { t.SetNLast(t.GetNLast() + 1); t.SetChi2Last(t.GetChi2Last() + maxChi2); } Double_t xcluster = cl->GetX(); - t.PropagateTo(xcluster,radLength,rho); - - if (!AdjustSector(&t)) { - break; //I.B's fix + t.PropagateTo(xcluster,xx0,xrho); + if (!AdjustSector(&t)) { + break; //I.B's fix } - maxChi2 = t.GetPredictedChi2(cl,h01); - - if (maxChi2<1e+10) + + maxChi2 = t.GetPredictedChi2(cl,h01); + if (maxChi2 < 1e+10) { if (!t.UpdateMI(cl,maxChi2,index,h01,plane)) { // ???? + } + else { + //SetCluster(cl, GetNumberOfClusters()-1); + kUPDATED = kTRUE; } + } - } + } } - } + } } - return expectedNumberOfClusters; + return expectedNumberOfClusters; } @@ -1193,23 +1018,31 @@ Int_t AliTRDtracker::FollowBackProlongation(AliTRDtrack &t) { // // Starting from current radial position of track this function - // extrapolates the track up to outer timebin and in the sensitive + // extrapolates the track up to the outer timebin and in the sensitive // layers confirms prolongation if a close cluster is found. // Returns the number of clusters expected to be found in sensitive layers - // Use GEO manager for material Description + // Uses the geomanager for material description // // return number of assigned clusters ? // - Int_t sector; - Int_t clusters[1000]; - Double_t radLength = 0.0; - Double_t rho = 0.0; + + Double_t xx0 = 0.0; + Double_t xrho = 0.0; + + Float_t ratio0 = 0.0; + Int_t expectedNumberOfClusters = 0; - Float_t ratio0 = 0.0; + AliTRDtracklet tracklet; + const Int_t kNclusters = 1000; + Int_t clusters[kNclusters]; + for (Int_t i = 0; i < kNclusters; i++) { + clusters[i] = -1; + } + // Calibration fill 2D AliTRDCalibraFillHisto *calibra = AliTRDCalibraFillHisto::Instance(); if (!calibra) { @@ -1219,76 +1052,82 @@ Int_t AliTRDtracker::FollowBackProlongation(AliTRDtrack &t) calibra->ResetTrack(); } - for (Int_t i = 0; i < 1000; i++) { - clusters[i] = -1; - } - + // Loop through the TRD planes for (Int_t iplane = 0; iplane < AliESDtrack::kNPlane; iplane++) { - int hb = iplane * 10; + Int_t hb = iplane * 10; fHClSearch->Fill(hb); + // Get the global time bin numbers for the first an last + // local time bin of the given plane Int_t row0 = GetGlobalTimeBin(0,iplane,GetTimeBinsPerPlane()-1); Int_t rowlast = GetGlobalTimeBin(0,iplane,0); + + // Get the X coordinates of the propagation layer for the first time bin Double_t currentx = fTrSec[0]->GetLayer(row0)->GetX(); if (currentx < t.GetX()) { fHClSearch->Fill(hb+1); continue; } - // - // Propagate closer to chamber if neccessary - // + // Propagate closer to the current chamber if neccessary if (currentx > (fgkMaxStep + t.GetX())) { if (!PropagateToX(t,currentx-fgkMaxStep,fgkMaxStep)) { fHClSearch->Fill(hb+2); break; } } + + // Rotate track to adjacent sector if neccessary if (!AdjustSector(&t)) { fHClSearch->Fill(hb+3); break; } + + // Check whether azimuthal angle is getting too large if (TMath::Abs(t.GetSnp()) > fgkMaxSnp) { fHClSearch->Fill(hb+4); break; } - // - // Get material budget inside of chamber - // Double_t xyz0[3]; Double_t xyz1[3]; Double_t param[7]; Double_t x; Double_t y; Double_t z; - // Starting global position + // Global start position (beginning of chamber) t.GetXYZ(xyz0); - // End global position + // X-position of the end of the chamber x = fTrSec[0]->GetLayer(rowlast)->GetX(); + // Get local Y and Z at the X-position of the end of the chamber if (!t.GetProlongation(x,y,z)) { fHClSearch->Fill(hb+5); break; } + // Global end position (end of chamber) xyz1[0] = x * TMath::Cos(t.GetAlpha()) - y * TMath::Sin(t.GetAlpha()); xyz1[1] = +x * TMath::Sin(t.GetAlpha()) + y * TMath::Cos(t.GetAlpha()); xyz1[2] = z; - AliKalmanTrack::MeanMaterialBudget(xyz0,xyz1,param); - rho = param[0]; - radLength = param[1]; // Get mean propagation parameters - // - // Find clusters - // - sector = t.GetSector(); + // Calculate the mean material budget along the path inside the chamber + AliTracker::MeanMaterialBudget(xyz0,xyz1,param); + // The mean propagation parameters (density*length and radiation length) + xrho = param[0]*param[4]; + xx0 = param[1]; + + // Find the clusters and tracklet along the path inside the chamber + sector = t.GetSector(); Float_t ncl = FindClusters(sector,row0,rowlast,&t,clusters,tracklet); fHNCl->Fill(tracklet.GetN()); + // Discard if in less than 1/3 of the available timebins + // clusters are found if (tracklet.GetN() < GetTimeBinsPerPlane()/3) { fHClSearch->Fill(hb+6); continue; } + //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! // // Propagate and update track @@ -1329,7 +1168,7 @@ Int_t AliTRDtracker::FollowBackProlongation(AliTRDtrack &t) t.SetChi2Last(t.GetChi2Last() + maxChi2); } Double_t xcluster = cl->GetX(); - t.PropagateTo(xcluster,radLength,rho); + t.PropagateTo(xcluster,xx0,xrho); maxChi2 = t.GetPredictedChi2(cl,h01); if (maxChi2<1e+10) @@ -1337,7 +1176,8 @@ Int_t AliTRDtracker::FollowBackProlongation(AliTRDtrack &t) if (!t.Update(cl,maxChi2,index,h01)) { // ???? } - } + } // else SetCluster(cl, GetNumberOfClusters()-1); // A.Bercuci 25.07.07 + if (calibra->GetMITracking()) { calibra->UpdateHistograms(cl,&t); @@ -1372,25 +1212,27 @@ Int_t AliTRDtracker::FollowBackProlongation(AliTRDtrack &t) } return expectedNumberOfClusters; + } //_____________________________________________________________________________ Int_t AliTRDtracker::PropagateToX(AliTRDtrack &t, Double_t xToGo, Double_t maxStep) { // - // Starting from current radial position of track this function + // Starting from current X-position of track this function // extrapolates the track up to radial position . // Returns 1 if track reaches the plane, and 0 otherwise // const Double_t kEpsilon = 0.00001; - //Double_t tanmax = TMath::Tan(0.5*AliTRDgeometry::GetAlpha()); + + // Current track X-position Double_t xpos = t.GetX(); - Double_t dir = (xpos kEpsilon) { + // Direction: inward or outward + Double_t dir = (xpos < xToGo) ? 1.0 : -1.0; - Double_t step = dir * TMath::Min(TMath::Abs(xToGo-xpos),maxStep); + while (((xToGo - xpos) * dir) > kEpsilon) { Double_t xyz0[3]; Double_t xyz1[3]; @@ -1398,23 +1240,39 @@ Int_t AliTRDtracker::PropagateToX(AliTRDtrack &t, Double_t xToGo, Double_t maxSt Double_t x; Double_t y; Double_t z; - // Starting global position - t.GetXYZ(xyz0); + + // The next step size + Double_t step = dir * TMath::Min(TMath::Abs(xToGo-xpos),maxStep); + + // Get the global position of the starting point + t.GetXYZ(xyz0); + + // X-position after next step x = xpos + step; + // Get local Y and Z at the X-position of the next step if (!t.GetProlongation(x,y,z)) { - return 0; // No prolongation + return 0; // No prolongation possible } + // The global position of the end point of this prolongation step xyz1[0] = x * TMath::Cos(t.GetAlpha()) - y * TMath::Sin(t.GetAlpha()); xyz1[1] = +x * TMath::Sin(t.GetAlpha()) + y * TMath::Cos(t.GetAlpha()); xyz1[2] = z; - AliKalmanTrack::MeanMaterialBudget(xyz0,xyz1,param); - if (!t.PropagateTo(x,param[1],param[0])) { + // Calculate the mean material budget between start and + // end point of this prolongation step + AliTracker::MeanMaterialBudget(xyz0,xyz1,param); + + // Propagate the track to the X-position after the next step + if (!t.PropagateTo(x,param[1],param[0]*param[4])) { return 0; } + + // Rotate the track if necessary AdjustSector(&t); + + // New track X-position xpos = t.GetX(); } @@ -1433,20 +1291,14 @@ Int_t AliTRDtracker::LoadClusters(TTree *cTree) // if (ReadClusters(fClusters,cTree)) { - AliError("Problem with reading the clusters !"); - return 1; + AliError("Problem with reading the clusters !"); + return 1; } - Int_t ncl = fClusters->GetEntriesFast(); + Int_t ncl = fClusters->GetEntriesFast(); fNclusters = ncl; AliInfo(Form("Sorting %d clusters",ncl)); UInt_t index; - for (Int_t ichamber = 0; ichamber < 5; ichamber++) { - for (Int_t isector = 0; isector < 18; isector++) { - fHoles[ichamber][isector] = kTRUE; - } - } - while (ncl--) { AliTRDcluster *c = (AliTRDcluster *) fClusters->UncheckedAt(ncl); @@ -1456,13 +1308,11 @@ Int_t AliTRDtracker::LoadClusters(TTree *cTree) Int_t plane = fGeom->GetPlane(detector); Int_t trackingSector = sector; - //if (c->GetLabel(0) > 0) { - if (c->GetQ() > 10) { - Int_t chamber = fGeom->GetChamber(detector); - fHoles[chamber][trackingSector] = kFALSE; - } + //if (c->GetQ() > 10) { + // Int_t chamber = fGeom->GetChamber(detector); + //} - Int_t gtb = fTrSec[trackingSector]->CookTimeBinIndex(plane,localTimeBin); + Int_t gtb = fTrSec[trackingSector]->CookTimeBinIndex(plane,localTimeBin); if (gtb < 0) { continue; } @@ -1470,12 +1320,11 @@ Int_t AliTRDtracker::LoadClusters(TTree *cTree) index = ncl; - // Apply pos correction - Transform(c); fHXCl->Fill(c->GetX()); fTrSec[trackingSector]->GetLayer(layer)->SetX(c->GetX()); fTrSec[trackingSector]->GetLayer(layer)->InsertCluster(c,index); + } return 0; @@ -1518,7 +1367,7 @@ void AliTRDtracker::UnloadClusters() } //_____________________________________________________________________________ -void AliTRDtracker::MakeSeedsMI(Int_t /*inner*/, Int_t /*outer*/, AliESD *esd) +void AliTRDtracker::MakeSeedsMI(Int_t /*inner*/, Int_t /*outer*/, AliESDEvent *esd) { // // Creates seeds using clusters between position inner plane and outer plane @@ -1571,7 +1420,7 @@ void AliTRDtracker::MakeSeedsMI(Int_t /*inner*/, Int_t /*outer*/, AliESD *esd) } } - AliTRDpadPlane *padPlane = AliTRDCommonParam::Instance()->GetPadPlane(0,0); + AliTRDpadPlane *padPlane = fGeom->GetPadPlane(0,0); Double_t h01 = TMath::Tan(-TMath::Pi() / 180.0 * padPlane->GetTiltingAngle()); Double_t hL[6]; // Tilting angle Double_t xcl[6]; // X - position of reference cluster @@ -2688,8 +2537,8 @@ void AliTRDtracker::MakeSeedsMI(Int_t /*inner*/, Int_t /*outer*/, AliESD *esd) Int_t AliTRDtracker::ReadClusters(TObjArray *array, TTree *clusterTree) const { // - // Reads AliTRDclusters (option >= 0) or AliTRDrecPoints (option < 0) - // from the file. The names of the cluster tree and branches + // Reads AliTRDclusters from the file. + // The names of the cluster tree and branches // should match the ones used in AliTRDclusterizer::WriteClusters() // @@ -3035,31 +2884,24 @@ AliTRDtracker::AliTRDtrackingSector Double_t *zc = new Double_t[kNchambers]; Double_t *zmax = new Double_t[kNchambers]; Double_t *zmaxsensitive = new Double_t[kNchambers]; - - AliTRDCommonParam *commonParam = AliTRDCommonParam::Instance(); - if (!commonParam) { - AliErrorGeneral("AliTRDtrackingSector::Ctor" - ,"Could not get common parameters\n"); - return; - } for (Int_t plane = 0; plane < AliTRDgeometry::Nplan(); plane++) { ymax = fGeom->GetChamberWidth(plane) / 2.0; - padPlane = commonParam->GetPadPlane(plane,0); + padPlane = fGeom->GetPadPlane(plane,0); ymaxsensitive = (padPlane->GetColSize(1) * padPlane->GetNcols() - 4.0) / 2.0; for (Int_t ch = 0; ch < kNchambers; ch++) { zmax[ch] = fGeom->GetChamberLength(plane,ch) / 2.0; Float_t pad = padPlane->GetRowSize(1); - Float_t row0 = commonParam->GetRow0(plane,ch,0); - Int_t nPads = commonParam->GetRowMax(plane,ch,0); + Float_t row0 = fGeom->GetRow0(plane,ch,0); + Int_t nPads = fGeom->GetRowMax(plane,ch,0); zmaxsensitive[ch] = Float_t(nPads) * pad / 2.0; zc[ch] = -(pad * nPads) / 2.0 + row0; } dx = AliTRDcalibDB::Instance()->GetVdrift(0,0,0) - / commonParam->GetSamplingFrequency(); + / AliTRDCommonParam::Instance()->GetSamplingFrequency(); rho = 0.00295 * 0.85; //???? radLength = 11.0; @@ -3106,14 +2948,28 @@ AliTRDtracker::AliTRDtrackingSector } +//_____________________________________________________________________________ +AliTRDtracker::AliTRDtrackingSector + ::~AliTRDtrackingSector() +{ + // + // Destructor + // + + for (Int_t i = 0; i < fN; i++) { + delete fLayers[i]; + } + +} + //_____________________________________________________________________________ Int_t AliTRDtracker::AliTRDtrackingSector ::CookTimeBinIndex(Int_t plane, Int_t localTB) const { // - // depending on the digitization parameters calculates "global" - // time bin index for timebin in plane - // + // Depending on the digitization parameters calculates global + // (i.e. for a whole TRD stack of 6 planes) time bin index for + // timebin in plane // Int_t tbPerPlane = AliTRDcalibDB::Instance()->GetNumberOfTimeBins(); @@ -3446,7 +3302,7 @@ Double_t AliTRDtracker::GetTiltFactor(const AliTRDcluster *c) Int_t det = c->GetDetector(); Int_t plane = fGeom->GetPlane(det); - AliTRDpadPlane *padPlane = AliTRDCommonParam::Instance()->GetPadPlane(plane,0); + AliTRDpadPlane *padPlane = fGeom->GetPadPlane(plane,0); Double_t h01 = TMath::Tan(-TMath::Pi() / 180.0 * padPlane->GetTiltingAngle()); if (fNoTilt) { @@ -3457,134 +3313,56 @@ Double_t AliTRDtracker::GetTiltFactor(const AliTRDcluster *c) } -//_____________________________________________________________________________ -void AliTRDtracker::CookdEdxTimBin(AliTRDtrack &TRDtrack) -{ - // - // This is setting fdEdxPlane and fTimBinPlane - // Sums up the charge in each plane for track TRDtrack and also get the - // Time bin for Max. Cluster - // Prashant Shukla (shukla@physi.uni-heidelberg.de) - // - - Double_t clscharge[AliESDtrack::kNPlane][AliESDtrack::kNSlice]; - Double_t maxclscharge[AliESDtrack::kNPlane]; - Int_t nCluster[AliESDtrack::kNPlane][AliESDtrack::kNSlice]; - Int_t timebin[AliESDtrack::kNPlane]; - - // Initialization of cluster charge per plane. - for (Int_t iPlane = 0; iPlane < AliESDtrack::kNPlane; iPlane++) { - for (Int_t iSlice = 0; iSlice < AliESDtrack::kNSlice; iSlice++) { - clscharge[iPlane][iSlice] = 0.0; - nCluster[iPlane][iSlice] = 0; - } - } - - // Initialization of cluster charge per plane. - for (Int_t iPlane = 0; iPlane < AliESDtrack::kNPlane; iPlane++) { - timebin[iPlane] = -1; - maxclscharge[iPlane] = 0.0; - } - - // Loop through all clusters associated to track TRDtrack - Int_t nClus = TRDtrack.GetNumberOfClusters(); // from Kalmantrack - for (Int_t iClus = 0; iClus < nClus; iClus++) { - Double_t charge = TRDtrack.GetClusterdQdl(iClus); - Int_t index = TRDtrack.GetClusterIndex(iClus); - AliTRDcluster *pTRDcluster = (AliTRDcluster *) GetCluster(index); - if (!pTRDcluster) { - continue; - } - Int_t tb = pTRDcluster->GetLocalTimeBin(); - if (!tb) { - continue; - } - Int_t detector = pTRDcluster->GetDetector(); - Int_t iPlane = fGeom->GetPlane(detector); - if (iPlane >= AliESDtrack::kNPlane) { - AliError(Form("Wrong plane %d",iPlane)); - continue; - } - Int_t iSlice = tb * AliESDtrack::kNSlice - / AliTRDcalibDB::Instance()->GetNumberOfTimeBins(); - if (iSlice >= AliESDtrack::kNSlice) { - AliError(Form("Wrong slice %d",iSlice)); - continue; - } - clscharge[iPlane][iSlice] = clscharge[iPlane][iSlice] + charge; - if (charge > maxclscharge[iPlane]) { - maxclscharge[iPlane] = charge; - timebin[iPlane] = tb; - } - nCluster[iPlane][iSlice]++; - } // End of loop over cluster - - // Setting the fdEdxPlane and fTimBinPlane variabales - Double_t totalCharge = 0.0; - - for (Int_t iPlane = 0; iPlane < AliESDtrack::kNPlane; iPlane++) { - for (Int_t iSlice = 0; iSlice < AliESDtrack::kNSlice; iSlice++) { - if (nCluster[iPlane][iSlice]) { - clscharge[iPlane][iSlice] /= nCluster[iPlane][iSlice]; - } - TRDtrack.SetPIDsignals(clscharge[iPlane][iSlice],iPlane,iSlice); - totalCharge = totalCharge+clscharge[iPlane][iSlice]; - } - TRDtrack.SetPIDTimBin(timebin[iPlane],iPlane); - } - - // Still needed ???? - // Int_t i; - // Int_t nc=TRDtrack.GetNumberOfClusters(); - // Float_t dedx=0; - // for (i=0; i and . + // Also the corresponding tracklet is calculated // Correction coeficients - depend on TRD parameters - to be changed accordingly // - Double_t x[100]; - Double_t yt[100]; - Double_t zt[100]; - Double_t xmean = 0.0; // Reference x - Double_t dz[10][100]; - Double_t dy[10][100]; - Float_t zmean[100]; - Float_t nmean[100]; + const Int_t kN1 = 100; + const Int_t kN2 = 10; + + Double_t x[kN1]; + Double_t yt[kN1]; + Double_t zt[kN1]; + Float_t zmean[kN1]; + Float_t nmean[kN1]; + + Double_t dz[kN2][kN1]; + Double_t dy[kN2][kN1]; + Int_t indexes[kN2][kN1]; // Indexes of the clusters in the road + Int_t best[kN2][kN1]; // Index of best matching cluster + AliTRDcluster *cl[kN2][kN1]; // Pointers to the clusters in the road + + Double_t xmean = 0.0; // Reference x Int_t clfound = 0; - Int_t indexes[10][100]; // Indexes of the clusters in the road - Int_t best[10][100]; // Index of best matching cluster - AliTRDcluster *cl[10][100]; // Pointers to the clusters in the road - - for (Int_t it = 0; it < 100; it++) { - x[it] = 0.0; - yt[it] = 0.0; - zt[it] = 0.0; + + // Initialize the arrays + for (Int_t it = 0; it < kN1; it++) { + + x[it] = 0.0; + yt[it] = 0.0; + zt[it] = 0.0; clusters[it] = -2; - zmean[it] = 0.0; - nmean[it] = 0.0; - for (Int_t ih = 0; ih < 10;ih++) { - indexes[ih][it] = -2; // Reset indexes1 + zmean[it] = 0.0; + nmean[it] = 0.0; + + for (Int_t ih = 0; ih < kN2; ih++) { + indexes[ih][it] = -2; cl[ih][it] = 0; dz[ih][it] = -100.0; dy[ih][it] = -100.0; best[ih][it] = 0; } + } Double_t x0 = track->GetX(); @@ -3595,6 +3373,7 @@ Int_t AliTRDtracker::FindClusters(Int_t sector, Int_t t0, Int_t t1 Int_t plane = -1; Int_t detector = -1; Float_t padlength = 0.0; + AliTRDtrack track2(* track); Float_t snpy = track->GetSnp(); Float_t tany = TMath::Sqrt(snpy*snpy / (1.0 - snpy*snpy)); @@ -3602,13 +3381,12 @@ Int_t AliTRDtracker::FindClusters(Int_t sector, Int_t t0, Int_t t1 tany *= -1.0; } - Double_t sy2 = ExpectedSigmaY2(x0,track->GetTgl(),track->GetPt()); + Double_t sy2 = ExpectedSigmaY2(x0,track->GetTgl(),track->GetSignedPt()); Double_t sz2 = ExpectedSigmaZ2(x0,track->GetTgl()); Double_t road = 15.0 * TMath::Sqrt(track->GetSigmaY2() + sy2); if (road > 6.0) { road = 6.0; } - //road = 20.0; for (Int_t it = 0; it < t1-t0; it++) { @@ -3632,12 +3410,11 @@ Int_t AliTRDtracker::FindClusters(Int_t sector, Int_t t0, Int_t t1 // // Find 2 nearest cluster at given time bin // - int checkPoint[4] = {0,0,0,0}; - double minY = 123456789; - double minD[2] = {1,1}; + Int_t checkPoint[4] = { 0, 0, 0, 0 }; + Double_t minY = 123456789.0; + Double_t minD[2] = { 1.0, 1.0 }; for (Int_t i = timeBin.Find(y - road); i < maxn; i++) { - //for (Int_t i = 0; i < maxn; i++) { AliTRDcluster *c = (AliTRDcluster *) (timeBin[i]); h01 = GetTiltFactor(c); @@ -3647,18 +3424,16 @@ Int_t AliTRDtracker::FindClusters(Int_t sector, Int_t t0, Int_t t1 padlength = TMath::Sqrt(c->GetSigmaZ2() * 12.0); } - //if (c->GetLocalTimeBin()==0) continue; if (c->GetY() > (y + road)) { break; } fHDeltaX->Fill(c->GetX() - x[it]); - //printf("%f\t%f\t%f \n", c->GetX(), x[it], c->GetX()-x[it]); if (TMath::Abs(c->GetY()-y) < TMath::Abs(minY)) { - minY = c->GetY()-y; - minD[0] = c->GetY()-y; - minD[1] = c->GetZ()-z; + minY = c->GetY() - y; + minD[0] = c->GetY() - y; + minD[1] = c->GetZ() - z; } checkPoint[0]++; @@ -3670,13 +3445,13 @@ Int_t AliTRDtracker::FindClusters(Int_t sector, Int_t t0, Int_t t1 checkPoint[1]++; Double_t dist = TMath::Abs(c->GetZ() - z); - if (dist > (0.5 * padlength + 6.0 * sigmaz)) { // 0.5 + if (dist > (0.5 * padlength + 6.0 * sigmaz)) { continue; // 6 sigma boundary cut } checkPoint[2]++; - Double_t cost = 0.0; // Sigma boundary cost function + Double_t cost = 0.0; if (dist> (0.5 * padlength - sigmaz)){ cost = (dist - 0.5*padlength) / (2.0 * sigmaz); if (cost > -1) { @@ -3686,8 +3461,6 @@ Int_t AliTRDtracker::FindClusters(Int_t sector, Int_t t0, Int_t t1 cost = 0.0; } } - //Int_t label = TMath::Abs(track->GetLabel()); - //if (c->GetLabel(0)!=label && c->GetLabel(1)!=label&&c->GetLabel(2)!=label) continue; chi2 = track2.GetPredictedChi2(c,h01) + cost; clfound++; @@ -3696,8 +3469,8 @@ Int_t AliTRDtracker::FindClusters(Int_t sector, Int_t t0, Int_t t1 } checkPoint[3]++; - detector = c->GetDetector(); // Store the clusters in the road + detector = c->GetDetector(); for (Int_t ih = 2; ih < 9; ih++) { if (cl[ih][it] == 0) { cl[ih][it] = c; @@ -3721,15 +3494,19 @@ Int_t AliTRDtracker::FindClusters(Int_t sector, Int_t t0, Int_t t1 } - for(int iCheckPoint = 0; iCheckPoint<4; iCheckPoint++) + for(int iCheckPoint = 0; iCheckPoint<4; iCheckPoint++) { fHFindCl[iCheckPoint]->Fill(checkPoint[iCheckPoint]); + } if (checkPoint[3]) { - if (track->GetPt() > 0) fHMinYPos->Fill(minY); - else fHMinYNeg->Fill(minY); - - fHMinD->Fill(minD[0], minD[1]); - } + if (track->GetSignedPt() > 0) { + fHMinYPos->Fill(minY); + } + else { + fHMinYNeg->Fill(minY); + } + fHMinD->Fill(minD[0],minD[1]); + } if (cl[0][it]) { nfound++; @@ -3747,7 +3524,6 @@ Int_t AliTRDtracker::FindClusters(Int_t sector, Int_t t0, Int_t t1 // // Choose one of the variants // - Int_t changes[10]; Float_t sumz = 0.0; Float_t sum = 0.0; Double_t sumdy = 0.0; @@ -3757,36 +3533,38 @@ Int_t AliTRDtracker::FindClusters(Int_t sector, Int_t t0, Int_t t1 Double_t sumx2 = 0.0; Double_t mpads = 0.0; - Int_t ngood[10]; - Int_t nbad[10]; + Int_t changes[kN2]; + + Int_t ngood[kN2]; + Int_t nbad[kN2]; - Double_t meanz[10]; - Double_t moffset[10]; // Mean offset - Double_t mean[10]; // Mean value - Double_t angle[10]; // Angle + Double_t meanz[kN2]; + Double_t moffset[kN2]; // Mean offset + Double_t mean[kN2]; // Mean value + Double_t angle[kN2]; // Angle - Double_t smoffset[10]; // Sigma of mean offset - Double_t smean[10]; // Sigma of mean value - Double_t sangle[10]; // Sigma of angle - Double_t smeanangle[10]; // Correlation + Double_t smoffset[kN2]; // Sigma of mean offset + Double_t smean[kN2]; // Sigma of mean value + Double_t sangle[kN2]; // Sigma of angle + Double_t smeanangle[kN2]; // Correlation - Double_t sigmas[10]; - Double_t tchi2s[10]; // Chi2s for tracklet + Double_t sigmas[kN2]; + Double_t tchi2s[kN2]; // Chi2s for tracklet - for (Int_t it = 0; it < 10; it++) { + for (Int_t it = 0; it < kN2; it++) { ngood[it] = 0; nbad[it] = 0; meanz[it] = 0.0; - moffset[it] = 0.0; // Mean offset - mean[it] = 0.0; // Mean value - angle[it] = 0.0; // Angle + moffset[it] = 0.0; // Mean offset + mean[it] = 0.0; // Mean value + angle[it] = 0.0; // Angle smoffset[it] = 1.0e5; // Sigma of mean offset smean[it] = 1.0e5; // Sigma of mean value sangle[it] = 1.0e5; // Sigma of angle - smeanangle[it] = 0.0; // Correlation + smeanangle[it] = 0.0; // Correlation sigmas[it] = 1.0e5; tchi2s[it] = 1.0e5; // Chi2s for tracklet @@ -3894,7 +3672,8 @@ Int_t AliTRDtracker::FindClusters(Int_t sector, Int_t t0, Int_t t1 changes[iter]++; } - Double_t dx = x[it]-xmean; // Distance to reference x + // Distance to reference x + Double_t dx = x[it]-xmean; sumz += cl[best[iter][it]][it]->GetZ(); sum++; sumdy += dy[best[iter][it]][it]; @@ -3915,12 +3694,12 @@ Int_t AliTRDtracker::FindClusters(Int_t sector, Int_t t0, Int_t t1 } // - // calculates line parameters + // Calculates line parameters // Double_t det = sum*sumx2 - sumx*sumx; angle[iter] = (sum*sumxy - sumx*sumdy) / det; mean[iter] = (sumx2*sumdy - sumx*sumxy) / det; - meanz[iter] = sumz / sum; + meanz[iter] = sumz / sum; moffset[iter] = sumdy / sum; mpads /= sum; // Mean number of pads @@ -3964,7 +3743,7 @@ Int_t AliTRDtracker::FindClusters(Int_t sector, Int_t t0, Int_t t1 Double_t mindist = 100000.0; Int_t ihbest = 0; - for (Int_t ih = 0; ih < 10; ih++) { + for (Int_t ih = 0; ih < kN2; ih++) { if (!cl[ih][it]) { break; } @@ -3975,7 +3754,9 @@ Int_t AliTRDtracker::FindClusters(Int_t sector, Int_t t0, Int_t t1 ihbest = ih; } } + best[iter+1][it] = ihbest; + } // @@ -3984,8 +3765,6 @@ Int_t AliTRDtracker::FindClusters(Int_t sector, Int_t t0, Int_t t1 Double_t sy2 = smean[iter] + track->GetSigmaY2(); Double_t sa2 = sangle[iter] + track->GetSigmaSnp2(); // track->fCee; Double_t say = track->GetSigmaSnpY(); // track->fCey; - //Double_t chi20 = mean[bestiter]*mean[bestiter ] / sy2+angle[bestiter]*angle[bestiter]/sa2; - //Double_t chi21 = mean[iter]*mean[iter] / sy2+angle[iter]*angle[iter]/sa2; Double_t detchi = sy2*sa2 - say*say; Double_t invers[3] = {sa2/detchi,sy2/detchi,-say/detchi}; // Inverse value of covariance matrix @@ -4016,9 +3795,6 @@ Int_t AliTRDtracker::FindClusters(Int_t sector, Int_t t0, Int_t t1 Short_t maxpos5 = -1; Float_t maxcharge5 = 0.0; - //if (tchi2s[bestiter]>25.) sigma2*=tchi2s[bestiter]/25.; - //if (tchi2s[bestiter]>25.) sigma2=1000.; // dont'accept - Double_t exB = AliTRDcalibDB::Instance()->GetOmegaTau(AliTRDcalibDB::Instance()->GetVdrift(0,0,0) ,-AliTracker::GetBz()*0.1); Double_t expectederr = sigma2*sigma2 + 0.01*0.01; @@ -4029,8 +3805,6 @@ Int_t AliTRDtracker::FindClusters(Int_t sector, Int_t t0, Int_t t1 expectederr += changes[bestiter] * 0.01; } expectederr += (0.03 * (tany-exB)*(tany-exB)) * 15.0; - //if (tchi2s[bestiter]>18.) expectederr*= tchi2s[bestiter]/18.; - //expectederr+=10000; for (Int_t it = 0; it < t1 - t0; it++) { @@ -4038,10 +3812,10 @@ Int_t AliTRDtracker::FindClusters(Int_t sector, Int_t t0, Int_t t1 continue; } - cl[best[bestiter][it]][it]->SetSigmaY2(expectederr); // Set cluster error + // Set cluster error + cl[best[bestiter][it]][it]->SetSigmaY2(expectederr); if (!cl[best[bestiter][it]][it]->IsUsed()) { cl[best[bestiter][it]][it]->SetY(cl[best[bestiter][it]][it]->GetY()); - //cl[best[bestiter][it]][it]->Use(); } // Time bins with maximal charge @@ -4086,16 +3860,9 @@ Int_t AliTRDtracker::FindClusters(Int_t sector, Int_t t0, Int_t t1 clusters[it+t0] = indexes[best[bestiter][it]][it]; - // Still needed ???? - //if (cl[best[bestiter][it]][it]->GetLocalTimeBin()>4 && - //cl[best[bestiter][it]][it]->GetLocalTimeBin()<18) clusters[it+t0] - // = indexes[best[bestiter][it]][it]; //Test - } - // // Set tracklet parameters - // Double_t trackleterr2 = smoffset[bestiter] + 0.01*0.01; if (mpads > 3.5) { trackleterr2 += (mpads - 3.5) * 0.04; @@ -4104,7 +3871,6 @@ Int_t AliTRDtracker::FindClusters(Int_t sector, Int_t t0, Int_t t1 trackleterr2 *= TMath::Max(14.0 - nfound,1.0); trackleterr2 += 0.2 * (tany-exB)*(tany-exB); - // Set tracklet parameters tracklet.Set(xmean ,track2.GetY() + moffset[bestiter] ,meanz[bestiter] @@ -4150,6 +3916,7 @@ Int_t AliTRDtracker::FindClusters(Int_t sector, Int_t t0, Int_t t1 new(array1[it]) AliTRDcluster(dummy); } } + TGraph graph0(t1-t0,x,dy0); TGraph graph1(t1-t0,x,dyb); TGraph graphy(t1-t0,x,yt); @@ -4221,6 +3988,10 @@ Int_t AliTRDtracker::Freq(Int_t n, const Int_t *inlist // The size of output array has is 2*n // + if (n <= 0) { + return 0; + } + Int_t *sindexS = new Int_t[n]; // Temporary array for sorting Int_t *sindexF = new Int_t[2*n]; for (Int_t i = 0; i < n; i++) { @@ -4311,7 +4082,8 @@ AliTRDtrack *AliTRDtracker::RegisterSeed(AliTRDseed *seeds, Double_t *params) ,c ,params[0] ,params[6]*alpha+shift); - track->PropagateTo(params[0]-5.0); + // SetCluster(cl, 0); // A. Bercuci + track->PropagateTo(params[0]-5.0); track->ResetCovariance(1); Int_t rc = FollowBackProlongation(*track); @@ -4321,77 +4093,103 @@ AliTRDtrack *AliTRDtracker::RegisterSeed(AliTRDseed *seeds, Double_t *params) } else { track->CookdEdx(); - CookdEdxTimBin(*track); + track->CookdEdxTimBin(-1); CookLabel(track,0.9); } return track; -} -////////////////////////////////////////////////////////////////////////////////////////// +} -void AliTRDtracker::InitLogHists() { +//_____________________________________________________________________________ +void AliTRDtracker::InitLogHists() +{ + // + // Create the log histograms + // - fHBackfit = new TH1D("logTRD_backfit", "", 40, -0.5, 39.5); - fHRefit = new TH1D("logTRD_refit", "", 40, -0.5, 39.5); - fHClSearch = new TH1D("logTRD_clSearch", "", 60, -0.5, 59.5); + fHBackfit = new TH1D("logTRD_backfit" ,"" + , 40,-0.5, 39.5); + fHRefit = new TH1D("logTRD_refit" ,"" + , 40,-0.5, 39.5); + fHClSearch = new TH1D("logTRD_clSearch","" + , 60,-0.5, 59.5); - fHX = new TH1D("logTRD_X", ";x (cm)", 200, 50, 400); - fHNCl = new TH1D("logTRD_ncl", "", 40, -0.5, 39.5); - fHNClTrack = new TH1D("logTRD_nclTrack", "", 180, -0.5, 179.5); + fHX = new TH1D("logTRD_X" ,";x (cm)" + , 200, 50, 400); + fHNCl = new TH1D("logTRD_ncl" ,"" + , 40,-0.5, 39.5); + fHNClTrack = new TH1D("logTRD_nclTrack","" + , 180,-0.5,179.5); - fHMinYPos = new TH1D("logTRD_minYPos", ";#delta Y (cm)", 400, -6, 6); - fHMinYNeg = new TH1D("logTRD_minYNeg", ";#delta Y (cm)", 400, -6, 6); - fHMinZ = new TH1D("logTRD_minZ", ";#delta Z (cm)", 400, -20, 20); - fHMinD = new TH2D("logTRD_minD", ";#delta Y (cm);#delta Z (cm)", 100, -6, 6, 100, -50, 50); + fHMinYPos = new TH1D("logTRD_minYPos" ,";#delta Y (cm)" + , 400, -6, 6); + fHMinYNeg = new TH1D("logTRD_minYNeg" ,";#delta Y (cm)" + , 400, -6, 6); + fHMinZ = new TH1D("logTRD_minZ" ,";#delta Z (cm)" + , 400, -20, 20); + fHMinD = new TH2D("logTRD_minD" ,";#delta Y (cm);#delta Z (cm)" + , 100, -6, 6 + , 100, -50, 50); - fHDeltaX = new TH1D("logTRD_deltaX", ";#delta X (cm)", 100, -5, 5); - fHXCl = new TH1D("logTRD_xCl", ";cluster x position (cm)", 1000, 280, 380); + fHDeltaX = new TH1D("logTRD_deltaX" ,";#delta X (cm)" + , 100, -5, 5); + fHXCl = new TH1D("logTRD_xCl" ,";cluster x position (cm)" + ,1000, 280, 380); - const char *nameFindCl[4] = {"logTRD_clY", "logTRD_clZ", "logTRD_clB", "logTRD_clG"}; + const Char_t *nameFindCl[4] = { "logTRD_clY" + , "logTRD_clZ" + , "logTRD_clB" + , "logTRD_clG" }; - for(int i=0; i<4; i++) { - fHFindCl[i] = new TH1D(nameFindCl[i], "", 30, -0.5, 29.5); + for (Int_t i = 0; i < 4; i++) { + fHFindCl[i] = new TH1D(nameFindCl[i],"",30,-0.5,29.5); } -} -////////////////////////////////////////////////////////////////////////////////////////// +} -void AliTRDtracker::SaveLogHists() { +//_____________________________________________________________________________ +void AliTRDtracker::SaveLogHists() +{ + // + // Save the log histograms in AliESDs.root + // - TDirectory *sav = gDirectory; - TFile *logFile = 0; + TDirectory *sav = gDirectory; + TFile *logFile = 0; TSeqCollection *col = gROOT->GetListOfFiles(); - int N = col->GetEntries(); - for(int i=0; iAt(i); - if (strstr(logFile->GetName(), "AliESDs.root")) break; + Int_t nn = col->GetEntries(); + for (Int_t i = 0; i < nn; i++) { + logFile = (TFile *) col->At(i); + if (strstr(logFile->GetName(),"AliESDs.root")) { + break; + } } logFile->cd(); - fHBackfit->Write(fHBackfit->GetName(), TObject::kOverwrite); - fHRefit->Write(fHRefit->GetName(), TObject::kOverwrite); - fHClSearch->Write(fHClSearch->GetName(), TObject::kOverwrite); - fHX->Write(fHX->GetName(), TObject::kOverwrite); - fHNCl->Write(fHNCl->GetName(), TObject::kOverwrite); - fHNClTrack->Write(fHNClTrack->GetName(), TObject::kOverwrite); - - fHMinYPos->Write(fHMinYPos->GetName(), TObject::kOverwrite); - fHMinYNeg->Write(fHMinYNeg->GetName(), TObject::kOverwrite); - fHMinD->Write(fHMinD->GetName(), TObject::kOverwrite); - fHMinZ->Write(fHMinZ->GetName(), TObject::kOverwrite); - - fHDeltaX->Write(fHDeltaX->GetName(), TObject::kOverwrite); - fHXCl->Write(fHXCl->GetName(), TObject::kOverwrite); + fHBackfit->Write(fHBackfit->GetName(),TObject::kOverwrite); + fHRefit->Write(fHRefit->GetName(),TObject::kOverwrite); + fHClSearch->Write(fHClSearch->GetName(),TObject::kOverwrite); + fHX->Write(fHX->GetName(),TObject::kOverwrite); + fHNCl->Write(fHNCl->GetName(),TObject::kOverwrite); + fHNClTrack->Write(fHNClTrack->GetName(),TObject::kOverwrite); + + fHMinYPos->Write(fHMinYPos->GetName(),TObject::kOverwrite); + fHMinYNeg->Write(fHMinYNeg->GetName(),TObject::kOverwrite); + fHMinD->Write(fHMinD->GetName(),TObject::kOverwrite); + fHMinZ->Write(fHMinZ->GetName(),TObject::kOverwrite); + + fHDeltaX->Write(fHDeltaX->GetName(),TObject::kOverwrite); + fHXCl->Write(fHXCl->GetName(),TObject::kOverwrite); - for(int i=0; i<4; i++) - fHFindCl[i]->Write(fHFindCl[i]->GetName(), TObject::kOverwrite); + for (Int_t i = 0; i < 4; i++) { + fHFindCl[i]->Write(fHFindCl[i]->GetName(),TObject::kOverwrite); + } logFile->Flush(); sav->cd(); -} -////////////////////////////////////////////////////////////////////////////////////////// +}