MC label for TRD tracks + bug fix in the stand alone tracker
authorcblume <cblume@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 29 Feb 2008 16:39:12 +0000 (16:39 +0000)
committercblume <cblume@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 29 Feb 2008 16:39:12 +0000 (16:39 +0000)
TRD/AliTRDtrackV1.cxx
TRD/AliTRDtrackV1.h
TRD/AliTRDtrackerDebug.cxx
TRD/AliTRDtrackerDebug.h
TRD/AliTRDtrackerV1.cxx
TRD/AliTRDtrackerV1.h

index 771e57a..e35452c 100644 (file)
@@ -144,6 +144,58 @@ AliTRDtrackV1::AliTRDtrackV1(AliTRDseedV1 *trklts, const Double_t p[5], const Do
 }
 
 //_______________________________________________________________
+Bool_t AliTRDtrackV1::CookLabel(Float_t wrong)
+{
+       // set MC label for this tracklet
+
+  Int_t s[kMAXCLUSTERSPERTRACK][2];
+  for (Int_t i = 0; i < kMAXCLUSTERSPERTRACK; i++) {
+    s[i][0] = -1;
+    s[i][1] =  0;
+  }
+
+  Bool_t labelAdded;
+       Int_t label;
+       AliTRDcluster *c    = 0x0;
+  for (Int_t ip = 0; ip < AliESDtrack::kNPlane; ip++) {
+               if(fTrackletIndex[ip] < 0) continue;
+               for (Int_t ic = 0; ic < AliTRDseed::knTimebins; ic++) {
+                       if(!(c = fTracklet[ip].GetClusters(ic))) continue;
+                       for (Int_t k = 0; k < 3; k++) { 
+                               label      = c->GetLabel(k);
+                               labelAdded = kFALSE; 
+                               Int_t j = 0;
+                               if (label >= 0) {
+                                       while ((!labelAdded) && (j < kMAXCLUSTERSPERTRACK)) {
+                                               if ((s[j][0] == label) || 
+                                                               (s[j][1] ==     0)) {
+                                                       s[j][0] = label; 
+                                                       s[j][1]++; 
+                                                       labelAdded = kTRUE;
+                                               }
+                                               j++;
+                                       }
+                               }
+                       }
+               }
+       }
+
+  Int_t max = 0;
+  label = -123456789;
+  for (Int_t i = 0; i < kMAXCLUSTERSPERTRACK; i++) {
+    if (s[i][1] <= max) continue;
+               max   = s[i][1]; 
+               label = s[i][0];
+  }
+
+  if ((1. - Float_t(max)/GetNumberOfClusters()) > wrong) label = -label;
+
+  SetLabel(label); 
+       
+       return kTRUE;
+}
+
+//_______________________________________________________________
 Bool_t AliTRDtrackV1::CookPID()
 {
   //
index 521b93b..51748ce 100644 (file)
@@ -30,6 +30,7 @@ class AliTRDtrackV1 : public AliTRDtrack
                                                              return *this; }
 
        Bool_t         CookPID();
+       Bool_t         CookLabel(Float_t wrong);
        Float_t        GetMomentum(Int_t plane) const;
        Double_t       GetPredictedChi2(const AliTRDseedV1 *tracklet) const;
         Double_t       GetPredictedChi2(const AliCluster* /*c*/) const                   { return 0.0; }
index bc7a53f..a89a34f 100644 (file)
@@ -31,6 +31,9 @@
 #include "TTree.h"
 #include "TTreeStream.h"
 #include "TLinearFitter.h"
+#include "TGraph.h"
+#include "TCanvas.h"
+#include "TMath.h"
 
 #include "AliLog.h"
 #include "AliTRDgeometry.h"
 #include "AliTRDseedV1.h"
 #include "AliTRDseed.h"
 #include "AliTRDcluster.h"
+#include "AliTRDgeometry.h"
 
 ClassImp(AliTRDtrackerDebug)
 
+Int_t AliTRDtrackerDebug::fgEventNumber = 0;
+Int_t AliTRDtrackerDebug::fgTrackNumber = 0;
+Int_t AliTRDtrackerDebug::fgCandidateNumber = 0;
+
 //____________________________________________________
 AliTRDtrackerDebug::AliTRDtrackerDebug() : AliTRDtrackerV1()
        ,fOutputStreamer(0x0)
@@ -76,7 +84,6 @@ void AliTRDtrackerDebug::Draw(const Option_t *)
 Bool_t AliTRDtrackerDebug::Init()
 {
 // steer linking data for various debug streams        
-       
        fTrack = new AliTRDtrackV1();
        fTree->SetBranchAddress("ncl", &fNClusters);
        fTree->SetBranchAddress("track.", &fTrack);
@@ -309,3 +316,466 @@ void AliTRDtrackerDebug::ResidualsTrackletsTrack() const
        }
 }
 
+//____________________________________________________
+void AliTRDtrackerDebug::AnalyseFindable(Char_t *treename){
+//
+// Calculates the number of findable tracklets defined as the number of tracklets
+// per track candidate where the tan phi_tracklet is below 0.15 (maximum inclination
+// in y-direction.
+// 
+// Parameters: -the treename (this method can be used for all trees which store the
+//                              tracklets
+// Output:             -void
+//
+// A new TTree containing the number of findable tracklets and the number of clusters
+// attached to the full track is stored to disk
+//
+       // Link the File
+       TFile *debfile = TFile::Open("TRD.TrackerDebug.root");
+       fTree = (TTree *)(debfile->Get(treename));
+       if(!fTree){
+               AliError(Form("Tree %s not found in file TRDdebug.root. Abborting!", treename));
+               debfile->Close();
+               return;
+       }
+       
+       AliTRDseedV1 *tracklets[kNPlanes];
+       for(Int_t iPlane = 0; iPlane < AliTRDtrackerV1::kNPlanes; iPlane++)
+               tracklets[iPlane] = 0x0;
+       for(Int_t iPlane = 0; iPlane < kNPlanes; iPlane++)
+               fTree->SetBranchAddress(Form("S%d.", iPlane), &tracklets[iPlane]);
+       fTree->SetBranchAddress("EventNumber", &fgEventNumber);
+       fTree->SetBranchAddress("CandidateNumber", &fgCandidateNumber);
+       
+       Int_t findable = 0, nClusters = 0;
+       Int_t nEntries = fTree->GetEntriesFast();
+       for(Int_t iEntry = 0; iEntry < nEntries; iEntry++){
+               printf("Entry %d\n", iEntry);
+               fTree->GetEntry(iEntry);
+               findable = 0;
+               nClusters = 0;
+               // Calculate Findable
+               for(Int_t iPlane = 0; iPlane < kNPlanes; iPlane++){
+                       if (TMath::Abs(tracklets[iPlane]->GetYref(0) / tracklets[iPlane]->GetX0()) < 0.15) findable++;
+                       if (!tracklets[iPlane]->IsOK()) continue;
+                       nClusters += tracklets[iPlane]->GetN2();
+               }
+               
+               // Fill Histogramms
+               TTreeSRedirector &cstreamer = *fOutputStreamer;
+               cstreamer << "AnalyseFindable"
+                       << "EventNumber="               << fgEventNumber
+                       << "CandidateNumber="   << fgCandidateNumber
+                       << "Findable="                  << findable
+                       << "NClusters="                 << nClusters
+                       << "\n";
+       }
+}
+//____________________________________________________
+void AliTRDtrackerDebug::AnalyseTiltedRiemanFit(){
+//
+// Creating a Data Set for the method FitTiltedRieman containing usefull variables
+// Each variable can be addressed to tracks later. Data can be processed later.
+//
+// Parameters: -
+// Output:     -
+//
+// TODO: Match everything with Init and Process
+//
+       TFile *debfile = TFile::Open("TRD.TrackerDebug.root");
+       fTree = (TTree *)(debfile->Get("MakeSeeds2"));
+       if(!fTree) return;
+       Int_t nEntries = fTree->GetEntries();
+       TLinearFitter *tiltedRiemanFitter = 0x0;
+       fTree->SetBranchAddress("FitterT.", &tiltedRiemanFitter);
+       fTree->SetBranchAddress("EventNumber", &fgEventNumber);
+       fTree->SetBranchAddress("CandidateNumber", &fgCandidateNumber);
+       for(Int_t entry = 0; entry < nEntries; entry++){
+               fTree->GetEntry(entry);
+               Double_t a = tiltedRiemanFitter->GetParameter(0);
+               Double_t b = tiltedRiemanFitter->GetParameter(1);
+               Double_t c = tiltedRiemanFitter->GetParameter(2);
+               Double_t offset = tiltedRiemanFitter->GetParameter(3);
+               Double_t slope  = tiltedRiemanFitter->GetParameter(4);
+               Float_t radius = GetTrackRadius(a, b, c);
+               Float_t curvature = GetTrackCurvature(a, b, c);
+               Float_t dca = GetDCA(a, b, c);
+               TTreeSRedirector &cstreamer = *fOutputStreamer;
+               cstreamer << "AnalyseTiltedRiemanFit"
+               << "EventNumber="               << fgEventNumber
+               << "CandidateNumber=" << fgCandidateNumber
+               << "Radius="                                    << radius
+               << "Curvature="                         << curvature
+               << "DCA="                                                       << dca
+               << "Offset="                                    << offset
+               << "Slope="                                             << slope
+               << "\n";
+       }
+}
+
+//____________________________________________________
+Float_t AliTRDtrackerDebug::GetTrackRadius(Float_t a, Float_t b, Float_t c) const {
+//
+// Calculates the track radius using the parameters given by the tilted Rieman fit 
+//
+// Parameters: The three parameters from the Rieman fit
+// Output:     The track radius
+//
+       Float_t radius = 0;
+       if(1.0 + b*b - c*a > 0.0)
+               radius = TMath::Sqrt(1.0 + b*b - c*a )/a;
+       return radius;
+}
+
+//____________________________________________________
+Float_t AliTRDtrackerDebug::GetTrackCurvature(Float_t a, Float_t b, Float_t c) const {
+//
+// Calculates the track curvature using the parameters given by the linear fitter 
+//
+// Parameters: the three parameters from the tilted Rieman fitter
+// Output:             the full track curvature
+//
+       Float_t curvature =  1.0 + b*b - c*a;
+       if (curvature > 0.0) 
+               curvature  =  a / TMath::Sqrt(curvature);
+       return curvature;
+}
+
+//____________________________________________________
+Float_t AliTRDtrackerDebug::GetDCA(Float_t a, Float_t b, Float_t c) const {
+//
+// Calculates the Distance to Clostest Approach for the Vertex using the paramters
+// given by the tilted Rieman fit 
+//
+// Parameters: the three parameters from the tilted Rieman fitter
+// Output:     the Distance to Closest Approach
+//
+       Float_t dca  =  0.0;
+       if (1.0 + b*b - c*a > 0.0) {
+               dca = -c / (TMath::Sqrt(1.0 + b*b - c*a) + TMath::Sqrt(1.0 + b*b));
+       }
+       return dca;
+}
+
+//____________________________________________________
+void AliTRDtrackerDebug::AnalyseMinMax()
+{
+//
+       TFile *debfile = TFile::Open("TRD.TrackerDebug.root");
+       if(!debfile){
+               AliError("File TRD.TrackerDebug.root not found!");
+               return; 
+       }
+       fTree = (TTree *)(debfile->Get("MakeSeeds0"));
+       if(!fTree){
+               AliError("Tree MakeSeeds0 not found in File TRD.TrackerDebug.root.");
+               return;
+       }
+       AliTRDseedV1 *cseed[4] = {0x0, 0x0, 0x0, 0x0};
+       AliTRDcluster *c[4] = {0x0, 0x0, 0x0, 0x0};
+       for(Int_t il = 0; il < 4; il++){
+               fTree->SetBranchAddress(Form("Seed%d.", il),    &cseed[il]);
+               fTree->SetBranchAddress(Form("c%d.",il), &c[il]);
+       }
+       fTree->SetBranchAddress("CandidateNumber",      &fgCandidateNumber);
+       fTree->SetBranchAddress("EventNumber",  &fgEventNumber);
+       Int_t entries = fTree->GetEntries();
+       for(Int_t ientry = 0; ientry < entries; ientry++){
+               fTree->GetEntry(ientry);
+               Float_t minmax[2] = { -100.0,  100.0 };
+               for (Int_t iLayer = 0; iLayer < 4; iLayer++) {
+                       Float_t max = c[iLayer]->GetZ() + cseed[iLayer]->GetPadLength() * 0.5 + 1.0 - cseed[iLayer]->GetZref(0);
+                       if (max < minmax[1]) minmax[1] = max;
+                       Float_t min = c[iLayer]->GetZ()-cseed[iLayer]->GetPadLength() * 0.5 - 1.0 - cseed[iLayer]->GetZref(0);
+                       if (min > minmax[0]) minmax[0] = min;
+               }
+               TTreeSRedirector &cstreamer = *fOutputStreamer;
+               cstreamer << "AnalyseMinMaxLayer"
+               << "EventNumber="                               << fgEventNumber
+               << "CandidateNumber="           << fgCandidateNumber
+               << "Min="                                                               << minmax[0]
+               << "Max="                                                               << minmax[1]
+               << "\n";
+       }
+}
+
+//____________________________________________________
+TCanvas* AliTRDtrackerDebug::PlotSeedingConfiguration(Char_t *direction, Int_t event, Int_t candidate){
+//
+// Plots the four seeding clusters, the helix fit and the reference Points for
+// a given combination consisting of eventnr. and candidatenr.
+//
+// Parameters:         -direction (y or z)
+//                             -Event Nr
+//             -Candidate that has to be plotted
+//
+       const Float_t kxmin = 280;
+       const Float_t kxmax = 380;
+       const Float_t kxdelta = (kxmax - kxmin)/1000;
+       
+       if((strcmp(direction, "y") != 0) && (strcmp(direction, "z") != 0)){
+               AliError(Form("Direction %s does not exist. Abborting!", direction));
+               return 0x0;
+       }
+
+       TFile *debfile = TFile::Open("TRD.TrackerDebug.root");
+       if(!debfile){
+               AliError("File TRD.TrackerDebug.root not found!");
+               return 0x0; 
+       }
+       fTree = (TTree *)(debfile->Get("MakeSeeds0"));
+       if(!fTree){
+               AliError("Tree MakeSeeds0 not found in File TRD.TrackerDebug.root.");
+               return 0x0;
+       }
+       
+       TGraph *seedcl = new TGraph(4);
+       TGraph *seedRef = new TGraph(4);
+       TGraph *riemanFit = new TGraph(1000);
+       seedcl->SetMarkerStyle(20);
+       seedcl->SetMarkerColor(kRed);
+       seedRef->SetMarkerStyle(2);
+
+       AliTRDcluster *c[4] = {0x0, 0x0, 0x0, 0x0};
+       AliRieman *rim = 0x0;
+       Bool_t found = kFALSE;
+       for(Int_t il = 0; il < 4; il++) fTree->SetBranchAddress(Form("c%d.",il), &c[il]);
+       fTree->SetBranchAddress("EventNumber", &fgEventNumber);
+       fTree->SetBranchAddress("CandidateNumber", &fgCandidateNumber);
+       fTree->SetBranchAddress("RiemanFitter.", &rim);
+       Int_t entries = fTree->GetEntries();
+       for(Int_t entry = 0; entry < entries; entry++){
+               fTree->GetEntry(entry);
+               if(fgEventNumber < event) continue;
+               if(fgEventNumber > event) break;
+               // EventNumber matching: Do the same for the candidate number
+               if(fgCandidateNumber < candidate) continue;
+               if(fgCandidateNumber > candidate) break;
+               found = kTRUE;
+               Int_t nPoints = 0;
+               for(Int_t il = 0; il < 4; il++){
+                       Float_t cluster = 0.0, reference = 0.0;
+                       if(!strcmp(direction, "y")){
+                               cluster = c[il]->GetY();
+                               reference = rim->GetYat(c[il]->GetX());
+                       }
+                       else{
+                               cluster = c[il]->GetZ();
+                               reference = rim->GetZat(c[il]->GetX());
+                       }
+                       seedcl->SetPoint(nPoints, cluster, c[il]->GetX());
+                       seedRef->SetPoint(nPoints, reference , c[il]->GetX());
+                       nPoints++;
+               }
+               // evaluate the fitter Function numerically
+               nPoints = 0;
+               for(Int_t ipt = 0; ipt < 1000; ipt++){
+                       Float_t x = kxmin + ipt * kxdelta;
+                       Float_t point = 0.0;
+                       if(!strcmp(direction, "y"))
+                               point = rim->GetYat(x);
+                       else
+                               point = rim->GetZat(x);
+                       riemanFit->SetPoint(nPoints++, point, x);
+               }
+               // We reached the End: break
+               break;
+       }
+       if(found){
+               seedcl->SetTitle(Form("Event %d, Candidate %d\n", fgEventNumber, fgCandidateNumber));
+               seedRef->SetTitle(Form("Event %d, Candidate %d\n", fgEventNumber, fgCandidateNumber));
+               riemanFit->SetTitle(Form("Event %d, Candidate %d\n", fgEventNumber, fgCandidateNumber));
+               TCanvas *c1 = new TCanvas();
+               seedcl->Draw("ap");
+               seedRef->Draw("psame");
+               riemanFit->Draw("lpsame");
+               return c1;
+       }
+       else{
+               AliError(Form("Combination consisting of event %d and candidate %d not found", event, candidate));
+               delete seedcl;
+               delete seedRef;
+               delete riemanFit;
+               return 0x0;
+       }
+}
+
+//____________________________________________________
+TCanvas* AliTRDtrackerDebug::PlotFullTrackFit(Int_t event, Int_t candidate, Int_t iteration, Char_t *direction){
+//
+// Plots the tracklets (clusters and reference in y direction) and the fitted function for several iterations
+// in the function ImproveSeedQuality (default is before ImproveSeedQuality)
+// 
+// Parameters: -Event Number
+//             -Candidate Number
+//             -Iteration Number in ImproveSeedQuality (default: -1 = before ImproveSeedQuality)
+//                        -direction (default: y)
+// Output:     -TCanvas (containing the Picture);
+//
+       const Float_t kxmin = 280;
+       const Float_t kxmax = 380;
+       const Float_t kxdelta = (kxmax - kxmin)/1000;
+       
+       if(strcmp(direction, "y") && strcmp(direction, "z")){
+               AliError(Form("Direction %s does not exist. Abborting!", direction));
+               return 0x0;
+       }
+
+       TFile *debfile = TFile::Open("TRD.TrackerDebug.root");
+       if(!debfile){
+               AliError("File TRD.TrackerDebug.root not found.");
+               return 0x0;
+       }
+       TString *treename = 0x0;
+       if(iteration > -1)
+               treename = new TString("ImproveSeedQuality");
+       else
+               treename = new TString("MakeSeeds1");
+       fTree = (TTree *)(debfile->Get(treename->Data()));
+       if(!fTree){
+               AliError(Form("Tree %s not found in File TRD.TrackerDebug.root.", treename->Data()));
+               delete treename;
+               return 0x0;
+       }
+       delete treename;
+
+       TGraph *fitfun = new TGraph(1000);
+       // Prepare containers
+       Float_t x0[AliTRDtrackerV1::kNPlanes],
+                       refP[AliTRDtrackerV1::kNPlanes],
+                       clx[AliTRDtrackerV1::kNPlanes * AliTRDtrackerV1::kNTimeBins],
+                       clp[AliTRDtrackerV1::kNPlanes * AliTRDtrackerV1::kNTimeBins];
+       Int_t nLayers = 0, ncls = 0;
+       
+       TLinearFitter *fitter = 0x0;
+       AliTRDseedV1 *tracklet[6] = {0x0, 0x0, 0x0, 0x0, 0x0, 0x0};
+       for(Int_t iLayer = 0; iLayer < 6; iLayer++)
+               fTree->SetBranchAddress(Form("S%d.", iLayer), &tracklet[iLayer]);
+       fTree->SetBranchAddress("FitterT.", &fitter);
+       fTree->SetBranchAddress("EventNumber", &fgEventNumber);
+       fTree->SetBranchAddress("CandidateNumber", &fgCandidateNumber);
+       
+       Int_t nEntries = fTree->GetEntriesFast();
+       Bool_t found = kFALSE;
+       for(Int_t entry = 0; entry < nEntries; entry++){
+               fTree->GetEntry(entry);
+               if(fgEventNumber < event) continue;
+               if(fgEventNumber > event) break;
+               // EventNumber matching: Do the same for the candidate number
+               if(fgCandidateNumber < candidate) continue;
+               if(fgCandidateNumber > candidate) break;
+               found = kTRUE;
+               
+               for(Int_t iLayer = 0; iLayer < 6; iLayer++){
+                       if(!tracklet[iLayer]->IsOK()) continue;
+                       x0[nLayers] = tracklet[iLayer]->GetX0();
+                       if(!strcmp(direction, "y"))
+                               refP[nLayers] = tracklet[iLayer]->GetYref(0);
+                       else
+                               refP[nLayers] = tracklet[iLayer]->GetZref(0);
+                       nLayers++;
+                       for(Int_t itb = 0; itb < 30; itb++){
+                               if(!tracklet[iLayer]->IsUsable(itb)) continue;
+                               AliTRDcluster *cl = tracklet[iLayer]->GetClusters(itb);
+                               if(!strcmp(direction, "y"))
+                                       clp[ncls] = cl->GetY();
+                               else
+                                       clp[ncls] = cl->GetZ();
+                               clx[ncls] = cl->GetX();
+                               ncls++;
+                       }
+               }
+               // Add function derived by the tilted Rieman fit (Defined by the curvature)
+               Int_t nPoints = 0;
+               if(!strcmp(direction, "y")){
+                       Double_t a = fitter->GetParameter(0);
+                       Double_t b = fitter->GetParameter(1);
+                       Double_t c = fitter->GetParameter(2);
+                       Double_t curvature =  1.0 + b*b - c*a;
+                       if (curvature > 0.0) {
+                               curvature  =  a / TMath::Sqrt(curvature);
+                       }
+                       // Numerical evaluation of the function:
+                       for(Int_t ipt = 0; ipt < 1000; ipt++){
+                               Float_t x = kxmin + ipt * kxdelta;
+                               Double_t res = (x * a + b);                                                             // = (x - x0)/y0
+                               res *= res;
+                               res  = 1.0 - c * a + b * b - res;                                       // = (R^2 - (x - x0)^2)/y0^2
+                               Double_t y = 0.;
+                               if (res >= 0) {
+                                       res = TMath::Sqrt(res);
+                                       y    = (1.0 - res) / a;
+                               }
+                               fitfun->SetPoint(nPoints++, y, x);
+                       }
+               }
+               else{
+                       Double_t offset = fitter->GetParameter(3);
+                       Double_t slope  = fitter->GetParameter(4);       
+                       // calculate the reference x (defined as medium between layer 2 and layer 3)
+                       // same procedure as in the tracker code
+                       Float_t medx = 0, xref = 0;
+                       Int_t startIndex = 5, nDistances = 0;
+                       for(Int_t iLayer = 5; iLayer > 0; iLayer--){
+                               if(tracklet[iLayer]->IsOK() && tracklet[iLayer - 1]->IsOK()){
+                                       medx += tracklet[iLayer]->GetX0() - tracklet[iLayer - 1]->GetX0();
+                                       startIndex = iLayer - 1;
+                                       nDistances++;
+                               }
+                       }
+                       if(nDistances){
+                               medx /= nDistances;
+                       }
+                       else{
+                               Float_t xpos[2];        memset(xpos, 0, sizeof(Float_t) * 2);
+                               Int_t ien = 0, idiff = 0;
+                               for(Int_t iLayer = 5; iLayer > 0; iLayer--){
+                                       if(tracklet[iLayer]->IsOK()){
+                                               xpos[ien++] = tracklet[iLayer]->GetX0();
+                                               startIndex = iLayer;
+                                       }
+                                       if(ien)
+                                               idiff++;
+                                       if(ien >=2)
+                                               break;
+                               }
+                               medx = (xpos[0] - xpos[1])/idiff;
+                       }
+                       xref = tracklet[startIndex]->GetX0() + medx * (2.5 - startIndex) - 0.5 * (AliTRDgeometry::AmThick() + AliTRDgeometry::DrThick());
+
+                       for(Int_t ipt = 0; ipt < 1000; ipt++){
+                               Float_t x = kxmin + ipt * kxdelta;
+                               Float_t z = offset + slope * (x - xref);
+                               fitfun->SetPoint(nPoints++, z, x);
+                       }
+               }
+               break;
+       }
+       if(found){
+               TGraph *trGraph         = new TGraph(ncls);
+               TGraph *refPoints       = new TGraph(nLayers);
+               trGraph->SetMarkerStyle(20);
+               trGraph->SetMarkerColor(kRed);
+               refPoints->SetMarkerStyle(21);
+               refPoints->SetMarkerColor(kBlue);
+               // fill the graphs
+               for(Int_t iLayer = 0; iLayer < nLayers; iLayer++)
+                       refPoints->SetPoint(iLayer, refP[iLayer], x0[iLayer]);
+               for(Int_t icls = 0; icls < ncls; icls++)
+                       trGraph->SetPoint(icls, clp[icls], clx[icls]);
+               TCanvas *c1 = new TCanvas();
+               trGraph->SetTitle(Form("Event %d, Candidate %d\n", fgEventNumber, fgCandidateNumber));
+               refPoints->SetTitle(Form("Event %d, Candidate %d\n", fgEventNumber, fgCandidateNumber));
+               fitfun->SetTitle(Form("Event %d, Candidate %d\n", fgEventNumber, fgCandidateNumber));
+               trGraph->Draw("ap");
+               refPoints->Draw("psame");
+               fitfun->Draw("lpsame");
+               return c1;
+       }
+       else{
+               AliError(Form("Combination consisting of event %d and candidate %d not found", event, candidate));
+               delete fitfun;
+               return 0x0;
+       }
+}
+
index 8f6ac60..7d6b70a 100644 (file)
@@ -13,6 +13,7 @@
 //  Authors:                                                              //
 //                                                                        //
 //    Alex Bercuci <A.Bercuci@gsi.de>                                     //
+//    Markus Fasel <M.Fasel@gsi.de>                                       //
 //                                                                        // 
 ////////////////////////////////////////////////////////////////////////////
 
@@ -21,6 +22,7 @@
 #endif
 
 class TTree;
+class TCanvas;
 class TTreeSRedirector;
 class AliTRDtrackV1;
 class AliTRDseedV1;
@@ -41,19 +43,40 @@ public:
        void        ResidualsClustersTracklet(const AliTRDseedV1 *tracklet) const;
        void        ResidualsClustersParametrisation(const AliTRDseedV1 *tracklet) const;
        void        ResidualsTrackletsTrack() const;
-
-
+       
+       void        AnalyseTiltedRiemanFit();
+       void        AnalyseMinMax();
+       void        AnalyseFindable(Char_t *treename);
+
+       TCanvas*    PlotSeedingConfiguration(Char_t *direction, Int_t event, Int_t Candidate);
+       TCanvas*    PlotFullTrackFit(Int_t event, Int_t candidate, Int_t iteration = -1, Char_t *direction = "y");
+       
+       static Int_t GetEventNumber(){ return fgEventNumber; }
+       static Int_t GetTrackNumber(){ return fgTrackNumber; }
+       static Int_t GetCandidateNumber(){ return fgCandidateNumber; }
+       
+       static void SetEventNumber(Int_t eventNumber){ fgEventNumber = eventNumber; }
+       static void SetTrackNumber(Int_t trackNumber){ fgTrackNumber = trackNumber; }
+       static void SetCandidateNumber(Int_t candidateNumber){ fgCandidateNumber = candidateNumber; }
+                       
 private:
        AliTRDtrackerDebug(const AliTRDtrackerDebug &);
        AliTRDtrackerDebug& operator=(const AliTRDtrackerDebug &);
 
+       Float_t     GetTrackRadius(Float_t a, Float_t b, Float_t c) const;
+       Float_t     GetTrackCurvature(Float_t a, Float_t b, Float_t c) const;
+       Float_t     GetDCA(Float_t a, Float_t b, Float_t c) const;
 
-  TTreeSRedirector *fOutputStreamer;                 //!Output streamer
+       TTreeSRedirector *fOutputStreamer;                 //!Output streamer
        TTree            *fTree;       // debug tree
        AliTRDseedV1     *fTracklet;   // current tracklet
        AliTRDtrackV1    *fTrack;      // current TRD track
        Int_t            fNClusters;   // N clusters for current track
        Float_t          fAlpha;       // sector
+       
+       static Int_t fgEventNumber;                             //  Event Number in the tracking code
+       static Int_t fgTrackNumber;       //  Track Number per Event
+       static Int_t fgCandidateNumber;   //  Candidate Number per event (Set in MakeSeeds)
 
        ClassDef(AliTRDtrackerDebug, 1) // debug suite of the TRD tracker
 };
index 99bc18d..16340c3 100644 (file)
 //                                                                           //
 ///////////////////////////////////////////////////////////////////////////////
 
-#include <Riostream.h>
-#include <stdio.h>
-#include <string.h>
+// #include <Riostream.h>
+// #include <stdio.h>
+// #include <string.h>
 
 #include <TBranch.h>
-#include <TFile.h>
-#include <TGraph.h>
-#include <TH1D.h>
-#include <TH2D.h>
+#include <TDirectory.h>
 #include <TLinearFitter.h>
-#include <TROOT.h>
 #include <TTree.h>  
 #include <TClonesArray.h>
-#include <TRandom.h>
 #include <TTreeStream.h>
 
 #include "AliLog.h"
 #include "AliESDEvent.h"
-#include "AliAlignObj.h"
+#include "AliGeomManager.h"
 #include "AliRieman.h"
 #include "AliTrackPointArray.h"
 
-#include "AliTRDtrackerV1.h"
-#include "AliTRDtrackingChamber.h"
 #include "AliTRDgeometry.h"
 #include "AliTRDpadPlane.h"
-#include "AliTRDgeometry.h"
-#include "AliTRDcluster.h" 
-#include "AliTRDtrack.h"
-#include "AliTRDseed.h"
 #include "AliTRDcalibDB.h"
-#include "AliTRDCommonParam.h"
 #include "AliTRDReconstructor.h"
 #include "AliTRDCalibraFillHisto.h"
-#include "AliTRDchamberTimeBin.h"
 #include "AliTRDrecoParam.h"
+
+#include "AliTRDcluster.h" 
 #include "AliTRDseedV1.h"
 #include "AliTRDtrackV1.h"
-#include "Cal/AliTRDCalDet.h"
+#include "AliTRDtrackerV1.h"
+#include "AliTRDtrackerDebug.h"
+#include "AliTRDtrackingChamber.h"
+#include "AliTRDchamberTimeBin.h"
+
 
 
 ClassImp(AliTRDtrackerV1)
@@ -306,7 +299,7 @@ Int_t AliTRDtrackerV1::PropagateBack(AliESDEvent *event)
                                //AliInfo(Form("Making backup track ncls [%d]...", foundClr));
                                track->CookdEdx();
                                track->CookdEdxTimBin(seed->GetID()); // A.Bercuci 25.07.07
-                               CookLabel(track,1 - fgkLabelFraction);
+                               track->CookLabel(1. - fgkLabelFraction);
                                if (track->GetBackupTrack()) UseClusters(track->GetBackupTrack());
                                
                                
@@ -627,10 +620,12 @@ Int_t AliTRDtrackerV1::FollowProlongation(AliTRDtrackV1 &t)
                        t.SetTracklet(tracklet, iplane, index);
                }
 
+               Int_t eventNumber = AliTRDtrackerDebug::GetEventNumber();
                TTreeSRedirector &cstreamer = *fgDebugStreamer;
                cstreamer << "FollowProlongation"
-                       << "ncl="      << nClustersExpected
-                       << "track.="   << &t
+                       << "EventNumber="       << eventNumber
+                       << "ncl="                                       << nClustersExpected
+                       << "track.="                    << &t
                        << "\n";
        }
 
@@ -781,15 +776,50 @@ Int_t AliTRDtrackerV1::FollowBackProlongation(AliTRDtrackV1 &t)
 
        if(AliTRDReconstructor::StreamLevel() > 1){
                TTreeSRedirector &cstreamer = *fgDebugStreamer;
+               Int_t eventNumber = AliTRDtrackerDebug::GetEventNumber();
                cstreamer << "FollowBackProlongation"
-                       << "ncl="      << nClustersExpected
-                       << "track.="   << &t
+                       << "EventNumber="                       << eventNumber
+                       << "ncl="                                                       << nClustersExpected
+                       << "track.="                                    << &t
                        << "\n";
        }
        
        return nClustersExpected;
 }
 
+//_____________________________________________________________________________
+void AliTRDtrackerV1::FitLeastSquare(Int_t nPoints, Float_t *x, Float_t *y, Float_t *error, Float_t *fitparams){
+// 
+// Performing a least square fit 
+//     
+// Parameters: - Number of Points
+//                             - x-data (array)
+//                             - y-data (array)
+//                             - error assumption in y-coordinate
+//                             - fitparams (array with length 2, as output)
+// Output:             -
+//
+// The first entry in fitparams is the offset, the second entry the slope
+//
+// @TODO: Implement error paramterisation
+//
+       Float_t sumweights = 0,
+                       sumx = 0,
+                       sumy = 0,
+                       sumxy = 0,
+                       sumx2 = 0;
+       for(Int_t ipt = 0; ipt < nPoints; ipt++){
+               sumweights += error[ipt];
+               sumx += x[ipt] * error[ipt];
+               sumy += y[ipt] * error[ipt];
+               sumxy += x[ipt] * y[ipt] * error[ipt];
+               sumx2 += x[ipt] * x[ipt] * error[ipt];
+       }
+       Float_t denominator = sumweights * sumx2 - sumx * sumx;
+       fitparams[0] = (sumy * sumx2 - sumx * sumxy) / denominator;
+       fitparams[1] = (sumweights * sumxy - sumx * sumy)/ denominator;
+}
+
 //_________________________________________________________________________
 Float_t AliTRDtrackerV1::FitRieman(AliTRDseedV1 *tracklets, Double_t *chi2, Int_t *planes){
 //
@@ -825,8 +855,8 @@ Float_t AliTRDtrackerV1::FitRieman(AliTRDseedV1 *tracklets, Double_t *chi2, Int_
                
                // chi2
                if((!tracklets[ppl[il]].IsOK()) && (!planes)) continue;
-               chi2[0] += tracklets[ppl[il]].GetChi2Z();
-               chi2[1] += tracklets[ppl[il]].GetChi2Y();
+               chi2[0] += tracklets[ppl[il]].GetChi2Y();
+               chi2[1] += tracklets[ppl[il]].GetChi2Z();
        }
        return fitter->GetC();
 }
@@ -898,7 +928,7 @@ Float_t AliTRDtrackerV1::FitTiltedRiemanConstraint(AliTRDseedV1 *tracklets, Doub
        Float_t x, y, z, w, t, error, tilt;
        Double_t uvt[2];
        Int_t nPoints = 0;
-  for(Int_t ipl = 0; ipl < AliTRDgeometry::kNplan; ipl++){
+       for(Int_t ipl = 0; ipl < AliTRDgeometry::kNplan; ipl++){
                if(!tracklets[ipl].IsOK()) continue;
                for(Int_t itb = 0; itb < fgNTimeBins; itb++){
                        if(!tracklets[ipl].IsUsable(itb)) continue;
@@ -920,7 +950,7 @@ Float_t AliTRDtrackerV1::FitTiltedRiemanConstraint(AliTRDseedV1 *tracklets, Doub
 
        // Calculate curvature
        Double_t a = fitter->GetParameter(0);
-       Double_t b = fitter->GetParameter(0);
+       Double_t b = fitter->GetParameter(1);
        Double_t curvature = a/TMath::Sqrt(b*b + 1);
 
        Float_t chi2track = fitter->GetChisquare()/Double_t(nPoints);
@@ -929,16 +959,20 @@ Float_t AliTRDtrackerV1::FitTiltedRiemanConstraint(AliTRDseedV1 *tracklets, Doub
 
        if(AliTRDReconstructor::StreamLevel() >= 5){
                //Linear Model on z-direction
-         Double_t xref = (tracklets[2].GetX0() + tracklets[3].GetX0())/2;              // Relative to the middle of the stack
+               Double_t xref = CalculateReferenceX(tracklets);         // Relative to the middle of the stack
                Double_t slope = fitter->GetParameter(2);
                Double_t zref = slope * xref;
-               Float_t chi2Z = CalculateChi2Z(tracklets, zref, slope);
+               Float_t chi2Z = CalculateChi2Z(tracklets, zref, slope, xref);
+               Int_t eventNumber = AliTRDtrackerDebug::GetEventNumber();
+               Int_t candidateNumber = AliTRDtrackerDebug::GetCandidateNumber();
                TTreeSRedirector &treeStreamer = *fgDebugStreamer;
                treeStreamer << "FitTiltedRiemanConstraint"
-                       << "Curvature=" << curvature
-                       << "Chi2Track=" << chi2track
-                       << "Chi2Z="                     << chi2Z
-                       << "zref="                      << zref
+                       << "EventNumber="               << eventNumber
+                       << "CandidateNumber="   << candidateNumber
+                       << "Curvature="                         << curvature
+                       << "Chi2Track="                         << chi2track
+                       << "Chi2Z="                                             << chi2Z
+                       << "zref="                                              << zref
                        << "\n";
        }
        return chi2track;
@@ -974,37 +1008,25 @@ Float_t AliTRDtrackerV1::FitTiltedRieman(AliTRDseedV1 *tracklets, Bool_t sigErro
 //
 // Paramters:   - Array of tracklets (connected to the track candidate)
 //              - Flag selecting the error definition
-// Output:      - Chi2 value of the track
-//
-// debug level 5
+// Output:      - Chi2 values of the track (in Parameter list)
 //
        TLinearFitter *fitter = GetTiltedRiemanFitter();
-  fitter->StoreData(kTRUE);
+       fitter->StoreData(kTRUE);
        fitter->ClearPoints();
        
-       // Calculate the reference position:
-       Int_t nDistances = 0;
-       Float_t meanDistance = 0.;
-       Int_t startIndex = 5;
-       for(Int_t il =5; il > 0; il--){
-               if(tracklets[il].IsOK() && tracklets[il -1].IsOK()){
-                       meanDistance += tracklets[il].GetX0() - tracklets[il -1].GetX0();
-                       nDistances++;
-               }
-               if(tracklets[il].IsOK()) startIndex = il;
-       }
-       meanDistance /= nDistances;
-       if(tracklets[0].IsOK()) startIndex = 0;
-       Double_t xref = tracklets[startIndex].GetX0() + (2.5 - startIndex) * meanDistance - 0.5 * (AliTRDgeometry::AmThick() + AliTRDgeometry::DrThick());
-
-       Float_t x, y, z, t, tilt, xdelta, rhs, error;
-       Float_t dzMean = 0;     Int_t dzcounter = 0;    // A reference z and a reference slope is used if the fitresults in z-direction are not acceptable
+       Double_t xref = CalculateReferenceX(tracklets);
+       Double_t x, y, z, t, tilt, xdelta, rhs, error;
        Double_t uvt[4];
        Int_t nPoints = 0;
-       for(Int_t ipl = 0; ipl < AliTRDgeometry::kNplan; ipl++){
+       // Containers for Least-square fitter
+       Float_t x0[kNPlanes], zfit[kNPlanes], errors[kNPlanes];
+       Int_t nLayers = 0;
+       for(Int_t ipl = 0; ipl < kNPlanes; ipl++){
                if(!tracklets[ipl].IsOK()) continue;
-               dzMean += tracklets[ipl].GetZfitR(1);
-               dzcounter++;
+               x0[nLayers] = tracklets[ipl].GetX0();
+               zfit[nLayers] = tracklets[ipl].GetZfit(0);
+               Double_t meanError = 0;
+               Int_t ncls = 0;
                for(Int_t itb = 0; itb < fgNTimeBins; itb++){
                        if (!tracklets[ipl].IsUsable(itb)) continue;
                        x = tracklets[ipl].GetX(itb) + tracklets[ipl].GetX0();
@@ -1024,7 +1046,11 @@ Float_t AliTRDtrackerV1::FitTiltedRieman(AliTRDseedV1 *tracklets, Bool_t sigErro
                        error *= sigError ? tracklets[ipl].GetSigmaY() : 0.2;
                        fitter->AddPoint(uvt, rhs, error);
                        nPoints++;
+                       meanError += tracklets[ipl].GetClusters(itb)->GetSigmaZ2();
+                       ncls++;
                }
+               errors[nLayers] = meanError / ncls;
+               nLayers++;
        }
        
        fitter->Eval();
@@ -1036,17 +1062,20 @@ Float_t AliTRDtrackerV1::FitTiltedRieman(AliTRDseedV1 *tracklets, Bool_t sigErro
        // Do not accept non possible z and dzdx combinations
        Bool_t acceptablez = kTRUE;
        Double_t zref = 0.0;
-       for (Int_t iLayer = 0; iLayer < AliTRDgeometry::kNplan; iLayer++) {
+       for (Int_t iLayer = 0; iLayer < kNPlanes; iLayer++) {
                if(!tracklets[iLayer].IsOK()) continue;
                zref = offset + slope * (tracklets[iLayer].GetX0() - xref);
                if (TMath::Abs(tracklets[iLayer].GetZProb() - zref) > tracklets[iLayer].GetPadLength() * 0.5 + 1.0) 
                        acceptablez = kFALSE;
        }
        if (!acceptablez) {
-               dzMean /= dzcounter;
-               Double_t zmf  = tracklets[startIndex].GetZfitR(0) + dzMean * (xref - tracklets[startIndex].GetX0()); // Z-Position of the track at the middle of a stack assuming a linear dependence on x (approximation)
+               Float_t fitparams[2];
+               FitLeastSquare(nLayers, x0, zfit, errors, fitparams);
+               Double_t dzmf   = fitparams[1];
+               Double_t zmf    = fitparams[0] + dzmf * xref;
+               //printf("In FitTilted Rieman: zmf = %f, meandz = %f\n", zmf, dzmf);
                fgTiltedRieman->FixParameter(3, zmf);
-               fgTiltedRieman->FixParameter(4, dzMean);
+               fgTiltedRieman->FixParameter(4, dzmf);
                fitter->Eval();
                fitter->ReleaseParameter(3);
                fitter->ReleaseParameter(4);
@@ -1059,11 +1088,8 @@ Float_t AliTRDtrackerV1::FitTiltedRieman(AliTRDseedV1 *tracklets, Bool_t sigErro
        Double_t b     =  fitter->GetParameter(1);
        Double_t c     =  fitter->GetParameter(2);
        Double_t curvature =  1.0 + b*b - c*a;
-       Double_t dca  =  0.0;                                                                                                                   // Distance to closest approach
-       if (curvature > 0.0) {
-               dca = -c / (TMath::Sqrt(1.0 + b*b - c*a) + TMath::Sqrt(1.0 + b*b));
+       if (curvature > 0.0) 
                curvature  =  a / TMath::Sqrt(curvature);
-       }
 
        Double_t chi2track = fitter->GetChisquare()/Double_t(nPoints);
 
@@ -1111,28 +1137,24 @@ Float_t AliTRDtrackerV1::FitTiltedRieman(AliTRDseedV1 *tracklets, Bool_t sigErro
                tracklets[iLayer].SetC(curvature);
                tracklets[iLayer].SetChi2(chi2track);
   }
-
-
-       if(AliTRDReconstructor::StreamLevel() >= 5){
-               Double_t chi2Z = CalculateChi2Z(tracklets, offset, slope);
-               TTreeSRedirector &treeStreamer = *fgDebugStreamer;
-               treeStreamer << "FitTiltedRieman"
-                       << "error="         << sigError
-                       << "Curvature="     << curvature
-                       << "Chi2track="     << chi2track
-                       << "Chi2Z="         << chi2Z
-                       << "D="             << c
-                       << "DCA="           << dca
-                       << "Offset="        << offset
-                       << "Slope="         << slope
-                       << "\n";
+       
+       if(AliTRDReconstructor::StreamLevel() >=5){
+               TTreeSRedirector &cstreamer = *fgDebugStreamer;
+               Int_t eventNumber                       = AliTRDtrackerDebug::GetEventNumber();
+               Int_t candidateNumber   = AliTRDtrackerDebug::GetCandidateNumber();
+               Double_t chi2z = CalculateChi2Z(tracklets, offset, slope, xref);
+               cstreamer << "FitTiltedRieman0"
+               << "EventNumber="                       << eventNumber
+               << "CandidateNumber="   << candidateNumber
+               << "xref="                                              << xref
+               << "Chi2Z="                                             << chi2z
+               << "\n";
        }
-
        return chi2track;
 }
 
 //_________________________________________________________________________
-Float_t AliTRDtrackerV1::CalculateChi2Z(AliTRDseedV1 *tracklets, Double_t offset, Double_t slope)
+Float_t AliTRDtrackerV1::CalculateChi2Z(AliTRDseedV1 *tracklets, Double_t offset, Double_t slope, Double_t xref)
 {
 //
 // Calculates the chi2-value of the track in z-Direction including tilting pad correction.
@@ -1141,9 +1163,9 @@ Float_t AliTRDtrackerV1::CalculateChi2Z(AliTRDseedV1 *tracklets, Double_t offset
 // Parameters: - Array of tracklets (AliTRDseedV1) related to the track candidate
 //             - the offset for the reference x
 //             - the slope
+//             - the reference x position
 // Output:     - The Chi2 value of the track in z-Direction
 //
-       Double_t xref = .5 * (tracklets[2].GetX0() + tracklets[3].GetX0());
        Float_t chi2Z = 0, nLayers = 0;
        for (Int_t iLayer = 0; iLayer < AliTRDgeometry::kNplan; iLayer++) {
                if(!tracklets[iLayer].IsOK()) continue;
@@ -1155,8 +1177,6 @@ Float_t AliTRDtrackerV1::CalculateChi2Z(AliTRDseedV1 *tracklets, Double_t offset
        return chi2Z;
 }
 
-
-
 //_____________________________________________________________________________
 Int_t AliTRDtrackerV1::PropagateToX(AliTRDtrackV1 &t, Double_t xToGo, Double_t maxStep)
 {
@@ -1261,6 +1281,7 @@ Int_t AliTRDtrackerV1::ReadClusters(TClonesArray* &array, TTree *clusterTree) co
     Int_t nCluster = clusterArray->GetEntriesFast();  
     for (Int_t iCluster = 0; iCluster < nCluster; iCluster++) { 
       if(!(c = (AliTRDcluster *) clusterArray->UncheckedAt(iCluster))) continue;
+                       c->SetInChamber();
       new((*fClusters)[ncl++]) AliTRDcluster(*c);
       clusterArray->RemoveAt(iCluster); 
     }
@@ -1326,6 +1347,8 @@ void AliTRDtrackerV1::UnloadClusters()
 
   for (int i = 0; i < AliTRDgeometry::kNsect; i++) fTrSec[i].Clear();
 
+       // Increment the Event Number
+       AliTRDtrackerDebug::SetEventNumber(AliTRDtrackerDebug::GetEventNumber()  + 1);
 }
 
 //_____________________________________________________________________________
@@ -1446,6 +1469,9 @@ Int_t AliTRDtrackerV1::Clusters2TracksSM(Int_t sector, AliESDEvent *esd)
        for(int itrack=0; itrack<nTracks; itrack++) 
           esd->AddTrack((AliESDtrack*)esdTrackList[itrack]);
 
+       // Reset Track and Candidate Number
+       AliTRDtrackerDebug::SetCandidateNumber(0);
+       AliTRDtrackerDebug::SetTrackNumber(0);
        return nTracks;
 }
 
@@ -1682,34 +1708,40 @@ Int_t AliTRDtrackerV1::Clusters2TracksStack(AliTRDtrackingChamber **stack, TClon
                                        }
                                        //Int_t eventNrInFile = esd->GetEventNumberInFile();
                                        //AliInfo(Form("Number of clusters %d.", nclusters));
+                                       Int_t eventNumber = AliTRDtrackerDebug::GetEventNumber();
+                                       Int_t trackNumber = AliTRDtrackerDebug::GetTrackNumber();
+                                       Int_t candidateNumber = AliTRDtrackerDebug::GetCandidateNumber();
                                        TTreeSRedirector &cstreamer = *fgDebugStreamer;
                                        cstreamer << "Clusters2TracksStack"
-                                               << "Iter="      << fSieveSeeding
-                                               << "Like="      << fTrackQuality[trackIndex]
-                                               << "S0.="       << dseed[0]
-                                               << "S1.="       << dseed[1]
-                                               << "S2.="       << dseed[2]
-                                               << "S3.="       << dseed[3]
-                                               << "S4.="       << dseed[4]
-                                               << "S5.="       << dseed[5]
-                                               << "p0=" << trackParams[0]
-                                               << "p1=" << trackParams[1]
-                                               << "p2=" << trackParams[2]
-                                               << "p3=" << trackParams[3]
-                                               << "p4=" << trackParams[4]
-                                               << "p5=" << trackParams[5]
-                                               << "p6=" << trackParams[6]
-                                               << "Label="     << label
-                                               << "Label1="    << label1
-                                               << "Label2="    << label2
-                                               << "FakeRatio=" << fakeratio
-                                               << "Freq="      << frequency
-                                               << "Ncl="       << ncl
-                                               << "NLayers="   << nlayers
-                                               << "Findable="  << findable
-                                               << "NUsed="     << nused
+                                               << "EventNumber="               << eventNumber
+                                               << "TrackNumber="               << trackNumber
+                                               << "CandidateNumber="   << candidateNumber
+                                               << "Iter="                              << fSieveSeeding
+                                               << "Like="                              << fTrackQuality[trackIndex]
+                                               << "S0.="                               << dseed[0]
+                                               << "S1.="                               << dseed[1]
+                                               << "S2.="                               << dseed[2]
+                                               << "S3.="                               << dseed[3]
+                                               << "S4.="                               << dseed[4]
+                                               << "S5.="                               << dseed[5]
+                                               << "p0="                                << trackParams[0]
+                                               << "p1="                                << trackParams[1]
+                                               << "p2="                                << trackParams[2]
+                                               << "p3="                                << trackParams[3]
+                                               << "p4="                                << trackParams[4]
+                                               << "p5="                                << trackParams[5]
+                                               << "p6="                                << trackParams[6]
+                                               << "Label="                             << label
+                                               << "Label1="                    << label1
+                                               << "Label2="                    << label2
+                                               << "FakeRatio="                 << fakeratio
+                                               << "Freq="                              << frequency
+                                               << "Ncl="                               << ncl
+                                               << "NLayers="                   << nlayers
+                                               << "Findable="                  << findable
+                                               << "NUsed="                             << nused
                                                << "\n";
-                               }
+                                       }
                        
                                AliTRDtrackV1 *track = MakeTrack(&sseed[trackIndex*kNPlanes], trackParams);
                                if(!track){
@@ -1722,6 +1754,7 @@ Int_t AliTRDtrackerV1::Clusters2TracksStack(AliTRDtrackingChamber **stack, TClon
                                esdTrack.SetLabel(track->GetLabel());
                                new ((*esdTrackList)[ntracks0++]) AliESDtrack(esdTrack);
                                ntracks1++;
+                               AliTRDtrackerDebug::SetTrackNumber(AliTRDtrackerDebug::GetTrackNumber() + 1);
                        }
 
                        jSieve++;
@@ -1952,72 +1985,58 @@ Int_t AliTRDtrackerV1::MakeSeeds(AliTRDtrackingChamber **stack, AliTRDseedV1 *ss
                                        if (c[0]->GetLabel(0) != c[3]->GetLabel(0)) isFake = kTRUE;
                                        if (c[1]->GetLabel(0) != c[3]->GetLabel(0)) isFake = kTRUE;
                                        if (c[2]->GetLabel(0) != c[3]->GetLabel(0)) isFake = kTRUE;
-                                       
+
+                                       Double_t xpos[4];
+                                       for(Int_t l = 0; l < kNSeedPlanes; l++) xpos[l] = layer[l]->GetX();
                                        Float_t yref[4];
                                        for(int il=0; il<4; il++) yref[il] = cseed[planes[il]].GetYref(0);
                                        Int_t ll = c[3]->GetLabel(0);
+                                       Int_t eventNumber = AliTRDtrackerDebug::GetEventNumber();
+                                       Int_t candidateNumber = AliTRDtrackerDebug::GetCandidateNumber();
+                                       AliRieman *rim = GetRiemanFitter();
                                        TTreeSRedirector &cs0 = *fgDebugStreamer;
-                                                       cs0 << "MakeSeeds0"
-                                                       <<"isFake=" << isFake
-                                                       <<"label=" << ll
-                                                       <<"chi2z=" << chi2[0]
-                                                       <<"chi2y=" << chi2[1]
-                                                       <<"yref0=" << yref[0]
-                                                       <<"yref1=" << yref[1]
-                                                       <<"yref2=" << yref[2]
-                                                       <<"yref3=" << yref[3]
-                                                       <<"c0.="   << c[0]
-                                                       <<"c1.="   << c[1]
-                                                       <<"c2.="   << c[2]
-                                                       <<"c3.="   << c[3]
-                                                       <<"\n";
+                                       cs0 << "MakeSeeds0"
+                                               <<"EventNumber="                << eventNumber
+                                               <<"CandidateNumber="    << candidateNumber
+                                               <<"isFake="                             << isFake
+                                               <<"config="                             << config
+                                               <<"label="                              << ll
+                                               <<"chi2z="                              << chi2[0]
+                                               <<"chi2y="                              << chi2[1]
+                                               <<"Y2exp="                              << cond2[0]     
+                                               <<"Z2exp="                              << cond2[1]
+                                           <<"X0="                                     << xpos[0] //layer[sLayer]->GetX()
+                                           <<"X1="                                     << xpos[1] //layer[sLayer + 1]->GetX()
+                                           <<"X2="                                     << xpos[2] //layer[sLayer + 2]->GetX()
+                                           <<"X3="                                     << xpos[3] //layer[sLayer + 3]->GetX()
+                                               <<"yref0="                              << yref[0]
+                                               <<"yref1="                              << yref[1]
+                                               <<"yref2="                              << yref[2]
+                                               <<"yref3="                              << yref[3]
+                                               <<"c0.="                                << c[0]
+                                               <<"c1.="                                << c[1]
+                                               <<"c2.="                                << c[2]
+                                               <<"c3.="                                << c[3]
+                                               <<"Seed0.="                             << &cseed[planes[0]]
+                                               <<"Seed1.="                             << &cseed[planes[1]]
+                                               <<"Seed2.="                             << &cseed[planes[2]]
+                                               <<"Seed3.="                             << &cseed[planes[3]]
+                                               <<"RiemanFitter.="              << rim
+                                               <<"\n";
                                }
 
                                if(chi2[0] > AliTRDReconstructor::RecoParam()->GetChi2Z()/*7./(3. - sLayer)*//*iter*/){
                                        //AliInfo(Form("Failed chi2 filter on chi2Z [%f].", chi2[0]));
+                                       AliTRDtrackerDebug::SetCandidateNumber(AliTRDtrackerDebug::GetCandidateNumber() + 1);
                                        continue;
                                }
                                if(chi2[1] > AliTRDReconstructor::RecoParam()->GetChi2Y()/*1./(3. - sLayer)*//*iter*/){
                                        //AliInfo(Form("Failed chi2 filter on chi2Y [%f].", chi2[1]));
+                                       AliTRDtrackerDebug::SetCandidateNumber(AliTRDtrackerDebug::GetCandidateNumber() + 1);
                                        continue;
                                }
                                //AliInfo("Passed chi2 filter.");
 
-                               if(AliTRDReconstructor::StreamLevel() >= 2){
-                                       Float_t minmax[2] = { -100.0,  100.0 };
-                                       for (Int_t iLayer = 0; iLayer < 4; iLayer++) {
-                                               Float_t max = c[iLayer]->GetZ() + cseed[planes[iLayer]].GetPadLength() * 0.5 + 1.0 - cseed[planes[iLayer]].GetZref(0);
-                                               if (max < minmax[1]) minmax[1] = max;
-                                               Float_t min = c[iLayer]->GetZ()-cseed[planes[iLayer]].GetPadLength() * 0.5 - 1.0 - cseed[planes[iLayer]].GetZref(0);
-                                               if (min > minmax[0]) minmax[0] = min;
-                                       }
-                                       Double_t xpos[4];
-                                       for(Int_t l = 0; l < kNSeedPlanes; l++) xpos[l] = layer[l]->GetX();
-                                       TTreeSRedirector &cstreamer = *fgDebugStreamer;
-                                                       cstreamer << "MakeSeeds1"
-                                               << "isFake=" << isFake
-                                               << "config="   << config
-                                               << "Cl0.="   << c[0]
-                                               << "Cl1.="   << c[1]
-                                               << "Cl2.="   << c[2]
-                                               << "Cl3.="   << c[3]
-                                               << "X0="     << xpos[0] //layer[sLayer]->GetX()
-                                               << "X1="     << xpos[1] //layer[sLayer + 1]->GetX()
-                                               << "X2="     << xpos[2] //layer[sLayer + 2]->GetX()
-                                               << "X3="     << xpos[3] //layer[sLayer + 3]->GetX()
-                                               << "Y2exp="  << cond2[0]
-                                               << "Z2exp="  << cond2[1]
-                                               << "Chi2R="  << chi2[0]
-                                               << "Chi2Z="  << chi2[1]
-                                               << "Seed0.=" << &cseed[planes[0]]
-                                               << "Seed1.=" << &cseed[planes[1]]
-                                               << "Seed2.=" << &cseed[planes[2]]
-                                               << "Seed3.=" << &cseed[planes[3]]
-                                               << "Zmin="   << minmax[0]
-                                               << "Zmax="   << minmax[1]
-                                               << "\n" ;
-                               }               
-                               
                                // try attaching clusters to tracklets
                                Int_t nUsedCl = 0;
                                Int_t nlayers = 0;
@@ -2030,18 +2049,25 @@ Int_t AliTRDtrackerV1::MakeSeeds(AliTRDtrackingChamber **stack, AliTRDseedV1 *ss
                                }
                                if(nlayers < kNSeedPlanes){ 
                                        //AliInfo(Form("Failed updating all seeds %d [%d].", nlayers, kNSeedPlanes));
+                                       AliTRDtrackerDebug::SetCandidateNumber(AliTRDtrackerDebug::GetCandidateNumber() + 1);
                                        continue;
                                }
                                // fit tracklets and cook likelihood
-                               FitRieman(&cseed[0], chi2, &planes[0]);
+                               FitTiltedRieman(&cseed[0], kTRUE);// Update Seeds and calculate Likelihood
+                               chi2[0] = GetChi2Y(&cseed[0]);
+                               chi2[1] = GetChi2Z(&cseed[0]);
+                               //Chi2 definitions in testing stage
+                               //chi2[0] = GetChi2YTest(&cseed[0]);
+                               //chi2[1] = GetChi2ZTest(&cseed[0]);
                                Double_t like = CookLikelihood(&cseed[0], planes, chi2); // to be checked
+
                                if (TMath::Log(1.E-9 + like) < AliTRDReconstructor::RecoParam()->GetTrackLikelihood()){
                                        //AliInfo(Form("Failed likelihood %f[%e].", TMath::Log(1.E-9 + like), like));
+                                       AliTRDtrackerDebug::SetCandidateNumber(AliTRDtrackerDebug::GetCandidateNumber() + 1);
                                        continue;
                                }
                                //AliInfo(Form("Passed likelihood %f[%e].", TMath::Log(1.E-9 + like), like));
 
-
                                // book preliminary results
                                seedQuality[ntracks] = like;
                                fSeedLayer[ntracks]  = config;/*sLayer;*/
@@ -2062,30 +2088,53 @@ Int_t AliTRDtrackerV1::MakeSeeds(AliTRDtrackingChamber **stack, AliTRDseedV1 *ss
                                        cseed[jLayer].SetPadLength(padlength[jLayer]);
 
                                        // fit extrapolated seed
-                                       FitTiltedRieman(cseed, kTRUE);
                                        if ((jLayer == 0) && !(cseed[1].IsOK())) continue;
                                        if ((jLayer == 5) && !(cseed[4].IsOK())) continue;
                                        AliTRDseedV1 tseed = cseed[jLayer];
                                        if(!tseed.AttachClustersIter(chamber, 1000.)) continue;
                                        cseed[jLayer] = tseed;
                                        nusedf += cseed[jLayer].GetNUsed(); // debug value
+                                       FitTiltedRieman(cseed,  kTRUE);
                                }
-                               FitTiltedRieman(cseed, kTRUE);
-                               //AliInfo("Extrapolation done.");
 
-                               if(ImproveSeedQuality(stack, cseed) < 4) continue;
+                               // AliInfo("Extrapolation done.");
+                               // Debug Stream containing all the 6 tracklets
+                               if(AliTRDReconstructor::StreamLevel() >= 2){
+                                       TTreeSRedirector &cstreamer = *fgDebugStreamer;
+                                       TLinearFitter *tiltedRieman = GetTiltedRiemanFitter();
+                                       Int_t eventNumber               = AliTRDtrackerDebug::GetEventNumber();
+                                       Int_t candidateNumber   = AliTRDtrackerDebug::GetCandidateNumber();
+                                       cstreamer << "MakeSeeds1"
+                                               << "EventNumber="               << eventNumber
+                                               << "CandidateNumber="   << candidateNumber
+                                               << "S0.="                               << &cseed[0]
+                                               << "S1.="                               << &cseed[1]
+                                               << "S2.="                               << &cseed[2]
+                                           << "S3.="                           << &cseed[3]
+                                           << "S4.="                           << &cseed[4]
+                                           << "S5.="                           << &cseed[5]
+                                           << "FitterT.="                      << tiltedRieman
+                                           << "\n";
+                               }
+                               
+                               if(ImproveSeedQuality(stack, cseed) < 4){
+                                       AliTRDtrackerDebug::SetCandidateNumber(AliTRDtrackerDebug::GetCandidateNumber() + 1);
+                                       continue;
+                               }
                                //AliInfo("Improve seed quality done.");
 
                                // fit full track and cook likelihoods
-                               Double_t curv = FitRieman(&cseed[0], chi2);
-                               Double_t chi2ZF = chi2[0] / TMath::Max((nlayers - 3.), 1.);
-                               Double_t chi2RF = chi2[1] / TMath::Max((nlayers - 3.), 1.);
+//                             Double_t curv = FitRieman(&cseed[0], chi2);
+//                             Double_t chi2ZF = chi2[0] / TMath::Max((nlayers - 3.), 1.);
+//                             Double_t chi2RF = chi2[1] / TMath::Max((nlayers - 3.), 1.);
 
                                // do the final track fitting (Once with vertex constraint and once without vertex constraint)
                                Double_t chi2Vals[3];
                                chi2Vals[0] = FitTiltedRieman(&cseed[0], kFALSE);
                                chi2Vals[1] = FitTiltedRiemanConstraint(&cseed[0], GetZ());
-                               chi2Vals[2] = chi2ZF;
+                               chi2Vals[2] = GetChi2Z(&cseed[0]) / TMath::Max((nlayers - 3.), 1.);
+                               // Chi2 definitions in testing stage
+                               //chi2Vals[2] = GetChi2ZTest(&cseed[0]);
                                fTrackQuality[ntracks] = CalculateTrackLikelihood(&cseed[0], &chi2Vals[0]);
                                //AliInfo("Hyperplane fit done\n");
 
@@ -2111,33 +2160,39 @@ Int_t AliTRDtrackerV1::MakeSeeds(AliTRDtrackingChamber **stack, AliTRDseedV1 *ss
                                Int_t frequency = outlab[1];
                                for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
                                        cseed[iLayer].SetFreq(frequency);
-                                       cseed[iLayer].SetChi2Z(chi2ZF);
+                                       cseed[iLayer].SetChi2Z(chi2[1]);
                                }
            
                                if(AliTRDReconstructor::StreamLevel() >= 2){
                                        TTreeSRedirector &cstreamer = *fgDebugStreamer;
+                                       Int_t eventNumber = AliTRDtrackerDebug::GetEventNumber();
+                                       Int_t candidateNumber = AliTRDtrackerDebug::GetCandidateNumber();
+                                       TLinearFitter *fitterTC = GetTiltedRiemanFitterConstraint();
+                                       TLinearFitter *fitterT = GetTiltedRiemanFitter();
                                        cstreamer << "MakeSeeds2"
-                                               << "C="       << curv
-                                               << "Chi2TR="  << chi2[0]
-                                               << "Chi2TC="  << chi2[1]
-                                               << "Chi2RF="  << chi2RF
-                                               << "Chi2ZF="  << chi2ZF
-                                               << "Nlayers=" << nlayers
-                                               << "NUsedS="  << nUsedCl
-                                               << "NUsed="   << nusedf
-                                               << "Like="    << like
-                                               << "S0.="     << &cseed[0]
-                                               << "S1.="     << &cseed[1]
-                                               << "S2.="     << &cseed[2]
-                                               << "S3.="     << &cseed[3]
-                                               << "S4.="     << &cseed[4]
-                                               << "S5.="     << &cseed[5]
-                                               << "Label="   << label
-                                               << "Freq="    << frequency
+                                               << "EventNumber="               << eventNumber
+                                               << "CandidateNumber="   << candidateNumber
+                                               << "Chi2TR="                    << chi2Vals[0]
+                                               << "Chi2TC="                    << chi2Vals[1]
+                                               << "Nlayers="                   << nlayers
+                                               << "NUsedS="                    << nUsedCl
+                                               << "NUsed="                             << nusedf
+                                               << "Like="                              << like
+                                               << "S0.="                               << &cseed[0]
+                                               << "S1.="                               << &cseed[1]
+                                               << "S2.="                               << &cseed[2]
+                                               << "S3.="                               << &cseed[3]
+                                               << "S4.="                               << &cseed[4]
+                                               << "S5.="                               << &cseed[5]
+                                               << "Label="                             << label
+                                               << "Freq="                              << frequency
+                                               << "FitterT.="                  << fitterT
+                                               << "FitterTC.="                 << fitterTC
                                                << "\n";
                                }
                                
                                ntracks++;
+                               AliTRDtrackerDebug::SetCandidateNumber(AliTRDtrackerDebug::GetCandidateNumber() + 1);
                                if(ntracks == kMaxTracksStack){
                                        AliWarning(Form("Number of seeds reached maximum allowed (%d) in stack.", kMaxTracksStack));
                                        for(int isl=0; isl<4; isl++) delete layer[isl];
@@ -2189,19 +2244,14 @@ AliTRDtrackV1* AliTRDtrackerV1::MakeTrack(AliTRDseedV1 *seeds, Double_t *params)
     delete track;
     track = 0x0;
   } else {
-//     track->CookdEdx();
-//     track->CookdEdxTimBin(-1);
-//     CookLabel(track, 0.9);
+     track->CookdEdx();
+     track->CookdEdxTimBin(-1);
+     track->CookLabel(.9);
   }
 
   return track;
 }
 
-//____________________________________________________________________
-void AliTRDtrackerV1::CookLabel(AliKalmanTrack */*pt*/, Float_t /*wrong*/) const
-{
-       // to be implemented, preferably at the level of TRD tracklet. !!!!!!!
-}
 
 //____________________________________________________________________
 Int_t AliTRDtrackerV1::ImproveSeedQuality(AliTRDtrackingChamber **stack, AliTRDseedV1 *cseed)
@@ -2222,6 +2272,8 @@ Int_t AliTRDtrackerV1::ImproveSeedQuality(AliTRDtrackingChamber **stack, AliTRDs
   // tracklet seed such that the seed quality (see AliTRDseed::GetQuality())
   // can be maximized. If some optimization is found the old seeds are replaced.
   //
+       // debug level: 7
+       //
        
        // make a local working copy
        AliTRDtrackingChamber *chamber = 0x0;
@@ -2258,6 +2310,24 @@ Int_t AliTRDtrackerV1::ImproveSeedQuality(AliTRDtrackingChamber **stack, AliTRDs
                }
 
                chi2 = FitTiltedRieman(bseed, kTRUE);
+               if(AliTRDReconstructor::StreamLevel() >= 7){
+                       Int_t eventNumber               = AliTRDtrackerDebug::GetEventNumber();
+                       Int_t candidateNumber = AliTRDtrackerDebug::GetCandidateNumber();
+                       TLinearFitter *tiltedRieman = GetTiltedRiemanFitter();
+                       TTreeSRedirector &cstreamer = *fgDebugStreamer;
+                       cstreamer << "ImproveSeedQuality"
+                       << "EventNumber="               << eventNumber
+                       << "CandidateNumber="   << candidateNumber
+                       << "Iteration="                         << iter
+                       << "S0.="                                                       << &bseed[0]
+                       << "S1.="                                                       << &bseed[1]
+                       << "S2.="                                                       << &bseed[2]
+                       << "S3.="                                                       << &bseed[3]
+                       << "S4.="                                                       << &bseed[4]
+                       << "S5.="                                                       << &bseed[5]
+                       << "FitterT.="                          << tiltedRieman
+                       << "\n";
+               }
        } // Loop: iter
        
        // we are sure that at least 2 tracklets are OK !
@@ -2299,12 +2369,16 @@ Double_t AliTRDtrackerV1::CalculateTrackLikelihood(AliTRDseedV1 *tracklets, Doub
        Double_t trackLikelihood     = likeChi2Z * likeChi2TR * likeAF;
 
        if(AliTRDReconstructor::StreamLevel() >= 2){
+               Int_t eventNumber = AliTRDtrackerDebug::GetEventNumber();
+               Int_t candidateNumber = AliTRDtrackerDebug::GetCandidateNumber();
                TTreeSRedirector &cstreamer = *fgDebugStreamer;
                cstreamer << "CalculateTrackLikelihood0"
-                       << "LikeChi2Z="         << likeChi2Z
-                       << "LikeChi2TR="        << likeChi2TR
-                       << "LikeChi2TC="        << likeChi2TC
-                       << "LikeAF="                    << likeAF
+                       << "EventNumber="                       << eventNumber
+                       << "CandidateNumber="   << candidateNumber
+                       << "LikeChi2Z="                         << likeChi2Z
+                       << "LikeChi2TR="                        << likeChi2TR
+                       << "LikeChi2TC="                        << likeChi2TC
+                       << "LikeAF="                                    << likeAF
                        << "TrackLikelihood=" << trackLikelihood
                        << "\n";
        }
@@ -2353,26 +2427,37 @@ Double_t AliTRDtrackerV1::CookLikelihood(AliTRDseedV1 *cseed, Int_t planes[4]
        }
        Double_t likea     = TMath::Exp(-sumda*10.6);
        Double_t likechi2y  = 0.0000000001;
-       if (chi2[1] < 0.5) likechi2y += TMath::Exp(-TMath::Sqrt(chi2[1]) * 7.73);
-       Double_t likechi2z = TMath::Exp(-chi2[0] * 0.088) / TMath::Exp(-chi2[0] * 0.019);
+       if (chi2[0] < 0.5) likechi2y += TMath::Exp(-TMath::Sqrt(chi2[0]) * 7.73);
+       Double_t likechi2z = TMath::Exp(-chi2[1] * 0.088) / TMath::Exp(-chi2[1] * 0.019);
        Int_t enc = Int_t(fgFindable*4.*fgNTimeBins);   // Expected Number Of Clusters, normally 72
        Double_t likeN     = TMath::Exp(-(enc - nclusters) * 0.19);
        
        Double_t like      = likea * likechi2y * likechi2z * likeN;
 
-       //AliInfo(Form("sumda(%f) chi2[0](%f) chi2[1](%f) likea(%f) likechi2y(%f) likechi2z(%f) nclusters(%d) likeN(%f)", sumda, chi2[0], chi2[1], likea, likechi2y, likechi2z, nclusters, likeN));
+//     AliInfo(Form("sumda(%f) chi2[0](%f) chi2[1](%f) likea(%f) likechi2y(%f) likechi2z(%f) nclusters(%d) likeN(%f)", sumda, chi2[0], chi2[1], likea, likechi2y, likechi2z, nclusters, likeN));
        if(AliTRDReconstructor::StreamLevel() >= 2){
+               Int_t eventNumber = AliTRDtrackerDebug::GetEventNumber();
+               Int_t candidateNumber = AliTRDtrackerDebug::GetCandidateNumber();
+               // The Debug Stream contains the seed 
                TTreeSRedirector &cstreamer = *fgDebugStreamer;
                cstreamer << "CookLikelihood"
-                       << "sumda="     << sumda
-                       << "chi0="      << chi2[0]
-                       << "chi1="      << chi2[1]
-                       << "likea="     << likea
-                       << "likechi2y=" << likechi2y
-                       << "likechi2z=" << likechi2z
-                       << "nclusters=" << nclusters
-                       << "likeN="     << likeN
-                       << "like="      << like
+                       << "EventNumber="                       << eventNumber
+                       << "CandidateNumber=" << candidateNumber
+                       << "tracklet0.="                        << &cseed[0]
+                       << "tracklet1.="                        << &cseed[1]
+                       << "tracklet2.="                        << &cseed[2]
+                       << "tracklet3.="                        << &cseed[3]
+                       << "tracklet4.="                        << &cseed[4]
+                       << "tracklet5.="                        << &cseed[5]
+                       << "sumda="                                             << sumda
+                       << "chi0="                                              << chi2[0]
+                       << "chi1="                                              << chi2[1]
+                       << "likea="                                             << likea
+                       << "likechi2y="                         << likechi2y
+                       << "likechi2z="                         << likechi2z
+                       << "nclusters="                         << nclusters
+                       << "likeN="                                             << likeN
+                       << "like="                                              << like
                        << "\n";
        }
 
@@ -2609,6 +2694,58 @@ AliCluster* AliTRDtrackerV1::GetCluster(Int_t idx) const
        return idx >= 0 || idx < ncls ? (AliCluster*)fClusters->UncheckedAt(idx) : 0x0;
 }
 
+//____________________________________________________________________
+Float_t AliTRDtrackerV1::CalculateReferenceX(AliTRDseedV1 *tracklets){
+//
+// Calculates the reference x-position for the tilted Rieman fit defined as middle
+// of the stack (middle between layers 2 and 3). For the calculation all the tracklets
+// are taken into account
+// 
+// Parameters: - Array of tracklets(AliTRDseedV1)
+//
+// Output:             - The reference x-position(Float_t)
+//
+       Int_t nDistances = 0;
+       Float_t meanDistance = 0.;
+       Int_t startIndex = 5;
+       for(Int_t il =5; il > 0; il--){
+               if(tracklets[il].IsOK() && tracklets[il -1].IsOK()){
+                       Float_t xdiff = tracklets[il].GetX0() - tracklets[il -1].GetX0();
+                       meanDistance += xdiff;
+                       nDistances++;
+               }
+               if(tracklets[il].IsOK()) startIndex = il;
+       }
+       if(tracklets[0].IsOK()) startIndex = 0;
+       if(!nDistances){
+               // We should normally never get here
+               Float_t xpos[2]; memset(xpos, 0, sizeof(Float_t) * 2);
+               Int_t iok = 0, idiff = 0;
+               // This attempt is worse and should be avoided:
+               // check for two chambers which are OK and repeat this without taking the mean value
+               // Strategy avoids a division by 0;
+               for(Int_t il = 5; il >= 0; il--){
+                       if(tracklets[il].IsOK()){
+                               xpos[iok] = tracklets[il].GetX0();
+                               iok++;
+                               startIndex = il;
+                       }
+                       if(iok) idiff++;        // to get the right difference;
+                       if(iok > 1) break;
+               }
+               if(iok > 1){
+                       meanDistance = (xpos[0] - xpos[1])/idiff;
+               }
+               else{
+                       // we have do not even have 2 layers which are OK? The we do not need to fit at all
+                       return 331.;
+               }
+       }
+       else{
+               meanDistance /= nDistances;
+       }
+       return tracklets[startIndex].GetX0() + (2.5 - startIndex) * meanDistance - 0.5 * (AliTRDgeometry::AmThick() + AliTRDgeometry::DrThick());
+}
 
 //_____________________________________________________________________________
 Int_t AliTRDtrackerV1::Freq(Int_t n, const Int_t *inlist
@@ -2668,3 +2805,31 @@ Int_t AliTRDtrackerV1::Freq(Int_t n, const Int_t *inlist
   return countPos;
 
 }
+
+//_____________________________________________________________________________
+Float_t AliTRDtrackerV1::GetChi2Y(AliTRDseedV1 *tracklets) const
+{
+//     Chi2 definition on y-direction
+
+       Float_t chi2 = 0;
+       for(Int_t ipl = 0; ipl < kNPlanes; ipl++){
+               if(!tracklets[ipl].IsOK()) continue;
+               Double_t distLayer = tracklets[ipl].GetYfit(0) - tracklets[ipl].GetYref(0); 
+               chi2 += distLayer * distLayer;
+       }
+       return chi2;
+}
+
+//_____________________________________________________________________________
+Float_t AliTRDtrackerV1::GetChi2Z(AliTRDseedV1 *tracklets) const 
+{
+//     Chi2 definition on z-direction
+
+       Float_t chi2 = 0;
+       for(Int_t ipl = 0; ipl < kNPlanes; ipl++){
+               if(!tracklets[ipl].IsOK()) continue;
+               Double_t distLayer = tracklets[ipl].GetMeanz() - tracklets[ipl].GetZref(0); 
+               chi2 += distLayer * distLayer;
+       }
+       return chi2;
+}
index a5afc1e..696cf5d 100644 (file)
@@ -92,10 +92,9 @@ public:
 protected:
   Bool_t         AdjustSector(AliTRDtrackV1 *track); 
        Double_t       BuildSeedingConfigs(AliTRDtrackingChamber **stack, Int_t *configs);
-       static Float_t CalculateChi2Z(AliTRDseedV1 *tracklets, Double_t offset, Double_t slope);
+       static Float_t CalculateChi2Z(AliTRDseedV1 *tracklets, Double_t offset, Double_t slope, Double_t xref);
        Int_t          Clusters2TracksSM(Int_t sector, AliESDEvent *esd);
        Int_t          Clusters2TracksStack(AliTRDtrackingChamber **stack, TClonesArray *esdTrackList);
-       void           CookLabel(AliKalmanTrack *pt, Float_t wrong) const;
        AliTRDseedV1*  GetTracklet(AliTRDtrackV1 *trk, Int_t plane, Int_t &idx);
        Bool_t         GetTrackPoint(Int_t index, AliTrackPoint &p) const;      
        Int_t          MakeSeeds(AliTRDtrackingChamber **stack, AliTRDseedV1 *sseed, Int_t *ipar);
@@ -106,17 +105,22 @@ protected:
 private:
        AliTRDtrackerV1(const AliTRDtrackerV1 &tracker);
        AliTRDtrackerV1 &operator=(const AliTRDtrackerV1 &tracker);
-       Double_t       CookLikelihood(AliTRDseedV1 *cseed, Int_t planes[4], Double_t *chi2);
-       Double_t       CalculateTrackLikelihood(AliTRDseedV1 *tracklets, Double_t *chi2);
-       Int_t          ImproveSeedQuality(AliTRDtrackingChamber **stack, AliTRDseedV1 *tracklet);
+       Double_t        CookLikelihood(AliTRDseedV1 *cseed, Int_t planes[4], Double_t *chi2);
+       Double_t        CalculateTrackLikelihood(AliTRDseedV1 *tracklets, Double_t *chi2);
+       Int_t           ImproveSeedQuality(AliTRDtrackingChamber **stack, AliTRDseedV1 *tracklet);
+       static void     FitLeastSquare(Int_t nPoints, Float_t *x, Float_t *y, Float_t *errors, Float_t *fitparams);
+       static Float_t  CalculateReferenceX(AliTRDseedV1 *tracklets);
+       
+       Float_t     GetChi2Y(AliTRDseedV1 *tracklets) const;
+       Float_t     GetChi2Z(AliTRDseedV1 *tracklets) const;
 
 private:
   AliTRDgeometry          *fGeom;                          // Pointer to TRD geometry
   AliTRDtrackingSector    fTrSec[kTrackingSectors];       // Array of tracking sectors;    
 
-       TClonesArray        *fClusters;                       // List of clusters
-       TClonesArray        *fTracklets;                      // List of tracklets
-       TClonesArray        *fTracks;                         // List of tracks
+  TClonesArray        *fClusters;                       // List of clusters
+  TClonesArray        *fTracklets;                      // List of tracklets
+  TClonesArray        *fTracks;                         // List of tracks
   
   // should go to the recoParam
   static const Double_t    fgkMaxChi2;                     // Max increment in track chi2 
@@ -135,8 +139,9 @@ private:
        static TLinearFitter *fgTiltedRieman;                 //  Fitter for the tilted Rieman fit without vertex constriant
        static TLinearFitter *fgTiltedRiemanConstrained;      //  Fitter for the tilted Rieman fit with vertex constraint       
        static AliRieman     *fgRieman;                       //  Fitter for the untilted Rieman fit
-  static TTreeSRedirector *fgDebugStreamer;             //!Debug streamer
-
+  
+       static TTreeSRedirector *fgDebugStreamer;             //!Debug streamer
+       
        ClassDef(AliTRDtrackerV1, 2)                          //  TRD tracker development class
 
 };