]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - TRD/AliTRDseedV1.cxx
- fix cluster ownership problems in time bin container
[u/mrichter/AliRoot.git] / TRD / AliTRDseedV1.cxx
index 7b603e098950addb6e9527a5dc03fe5cab1f0da8..43792bc2278996ce31e336f708ab6c64eb71a164 100644 (file)
 #include "AliTRDtrackerV1.h"
 #include "AliTRDReconstructor.h"
 #include "AliTRDrecoParam.h"
-#include "AliTRDgeometry.h"
 #include "Cal/AliTRDCalPID.h"
 
 ClassImp(AliTRDseedV1)
 
 //____________________________________________________________________
-AliTRDseedV1::AliTRDseedV1(Int_t plane
+AliTRDseedV1::AliTRDseedV1(Int_t det
   :AliTRDseed()
   ,fReconstructor(0x0)
-  ,fPlane(plane)
+  ,fClusterIter(0x0)
+  ,fClusterIdx(0)
+  ,fDet(det)
   ,fMom(0.)
   ,fSnp(0.)
   ,fTgl(0.)
@@ -70,7 +71,9 @@ AliTRDseedV1::AliTRDseedV1(Int_t plane)
 AliTRDseedV1::AliTRDseedV1(const AliTRDseedV1 &ref)
   :AliTRDseed((AliTRDseed&)ref)
   ,fReconstructor(ref.fReconstructor)
-  ,fPlane(ref.fPlane)
+  ,fClusterIter(0x0)
+  ,fClusterIdx(0)
+  ,fDet(ref.fDet)
   ,fMom(ref.fMom)
   ,fSnp(ref.fSnp)
   ,fTgl(ref.fTgl)
@@ -129,7 +132,9 @@ void AliTRDseedV1::Copy(TObject &ref) const
   //AliInfo("");
   AliTRDseedV1 &target = (AliTRDseedV1 &)ref; 
 
-  target.fPlane         = fPlane;
+  target.fClusterIter   = 0x0;
+  target.fClusterIdx    = 0;
+  target.fDet           = fDet;
   target.fMom           = fMom;
   target.fSnp           = fSnp;
   target.fTgl           = fTgl;
@@ -226,7 +231,8 @@ void AliTRDseedV1::CookdEdx(Int_t nslices)
     nclusters[slice]++;
   } // End of loop over clusters
 
-  if(fReconstructor->GetPIDMethod() == AliTRDReconstructor::kLQPID){
+  //if(fReconstructor->GetPIDMethod() == AliTRDReconstructor::kLQPID){
+  if(nslices == AliTRDReconstructor::kLQslices){
   // calculate mean charge per slice (only LQ PID)
     for(int is=0; is<nslices; is++){ 
       if(nclusters[is]) fdEdx[is] /= nclusters[is];
@@ -283,7 +289,7 @@ Double_t* AliTRDseedV1::GetProbability()
   
   // Sets the a priori probabilities
   for(int ispec=0; ispec<AliPID::kSPECIES; ispec++) {
-    fProb[ispec] = pd->GetProbability(ispec, fMom, &fdEdx[0], length, fPlane); 
+    fProb[ispec] = pd->GetProbability(ispec, fMom, &fdEdx[0], length, GetPlane());     
   }
 
   return &fProb[0];
@@ -361,30 +367,11 @@ Bool_t    AliTRDseedV1::AttachClustersIter(AliTRDtrackingChamber *chamber, Float_t
   }
 
   AliTRDchamberTimeBin *layer = 0x0;
-  if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker)>=7 && c){
-    TClonesArray clusters("AliTRDcluster", 24);
-    clusters.SetOwner(kTRUE);
-    AliTRDcluster *cc = 0x0;
-    Int_t det=-1, ncl, ncls = 0;
-    for (Int_t iTime = 0; iTime < AliTRDtrackerV1::GetNTimeBins(); iTime++) {
-      if(!(layer = chamber->GetTB(iTime))) continue;
-      if(!(ncl = Int_t(*layer))) continue;
-      for(int ic=0; ic<ncl; ic++){ 
-        cc = (*layer)[ic];
-        det = cc->GetDetector();
-        new(clusters[ncls++]) AliTRDcluster(*cc);
-      }
-    }
-    AliInfo(Form("N clusters[%d] = %d", fPlane, ncls));
-    
-    Int_t ref = c ? 1 : 0;
-    TTreeSRedirector &cstreamer = *AliTRDtrackerV1::DebugStreamer();
-    cstreamer << "AttachClustersIter"
-      << "det="        << det 
-      << "ref="        << ref 
-      << "clusters.="  << &clusters
+  if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker)>=7){
+    AliTRDtrackingChamber *ch = new AliTRDtrackingChamber(*chamber); 
+    (*AliTRDtrackerV1::DebugStreamer()) << "AttachClustersIter"
+      << "chamber.="   << ch
       << "tracklet.="  << this
-      << "cl.="        << c
       << "\n"; 
   }
 
@@ -429,7 +416,7 @@ Bool_t      AliTRDseedV1::AttachClustersIter(AliTRDtrackingChamber *chamber, Float_t
       fZ[iTime]        = cl->GetZ();
       ncl++;
     }
-    if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker)>=7) AliInfo(Form("iter = %d ncl [%d] = %d", iter, fPlane, ncl));
+    if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker)>=7) AliInfo(Form("iter = %d ncl [%d] = %d", iter, fDet, ncl));
     
     if(ncl>1){ 
       // calculate length of the time bin (calibration aware)
@@ -470,7 +457,7 @@ Bool_t      AliTRDseedV1::AttachClustersIter(AliTRDtrackingChamber *chamber, Float_t
       
       AliTRDseed::Update();
     }
-    if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker)>=7) AliInfo(Form("iter = %d nclFit [%d] = %d", iter, fPlane, fN2));
+    if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker)>=7) AliInfo(Form("iter = %d nclFit [%d] = %d", iter, fDet, fN2));
     
     if(IsOK()){
       tquality = GetQuality(kZcorr);
@@ -481,7 +468,7 @@ Bool_t      AliTRDseedV1::AttachClustersIter(AliTRDtrackingChamber *chamber, Float_t
   } // Loop: iter
   if (!IsOK()) return kFALSE;
 
-  CookLabels();
+  if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker)>=1) CookLabels();
   UpdateUsed();
   return kTRUE;        
 }
@@ -630,7 +617,7 @@ Bool_t      AliTRDseedV1::AttachClusters(AliTRDtrackingChamber *chamber
 }
 
 //____________________________________________________________________
-Bool_t AliTRDseedV1::Fit()
+Bool_t AliTRDseedV1::Fit(Bool_t tilt)
 {
   //
   // Linear fit of the tracklet
@@ -647,60 +634,87 @@ Bool_t AliTRDseedV1::Fit()
   //
 
   const Int_t kClmin = 8;
+  const Float_t q0 = 100.;
+  const Float_t clSigma0 = 2.E-2;    //[cm]
+  const Float_t clSlopeQ = -1.19E-2; //[1/cm]
+
+  // get track direction
+  Double_t y0   = fYref[0];
+  Double_t dydx = fYref[1]; 
+  Double_t z0   = fZref[0];
+  Double_t dzdx = fZref[1];
+  Double_t yt, zt;
+
   const Int_t kNtb = AliTRDtrackerV1::GetNTimeBins();
   AliTRDtrackerV1::AliTRDLeastSquare fitterY, fitterZ;
 
   // convertion factor from square to gauss distribution for sigma
   Double_t convert = 1./TMath::Sqrt(12.);
-
+  
   // book cluster information
-  Double_t xc[knTimebins+1], yc[knTimebins], zc[knTimebins+1], sy[knTimebins], sz[knTimebins+1];
+  Double_t xc[knTimebins], yc[knTimebins], zc[knTimebins], sy[knTimebins], sz[knTimebins];
   Int_t zRow[knTimebins];
-  AliTRDcluster *c = 0x0;
-  Int_t nc = 0;
-  for (Int_t ic=0; ic<kNtb; ic++) {
+
+
+  fN = 0;
+  AliTRDcluster *c=0x0, **jc = &fClusters[0];
+  for (Int_t ic=0; ic<kNtb; ic++, ++jc) {
     zRow[ic] = -1;
     xc[ic]  = -1.;
     yc[ic]  = 999.;
     zc[ic]  = 999.;
     sy[ic]  = 0.;
     sz[ic]  = 0.;
-    if(!(c = fClusters[ic])) continue;
+    if(!(c = (*jc))) continue;
     if(!c->IsInChamber()) continue;
     Float_t w = 1.;
     if(c->GetNPads()>4) w = .5;
     if(c->GetNPads()>5) w = .2;
-    zRow[nc] = c->GetPadRow();
-    xc[nc]   = fX0 - c->GetX();
-    yc[nc]   = c->GetY();
-    zc[nc]   = c->GetZ();
-    sy[nc]   = w; // all clusters have the same sigma
-    sz[nc]   = fPadLength*convert;
-    fitterZ.AddPoint(&xc[nc], zc[nc], sz[nc]);
-    nc++;
+    zRow[fN] = c->GetPadRow();
+    xc[fN]   = fX0 - c->GetX();
+    yc[fN]   = c->GetY();
+    zc[fN]   = c->GetZ();
+
+    // extrapolated y value for the track
+    yt = y0 - xc[fN]*dydx; 
+    // extrapolated z value for the track
+    zt = z0 - xc[fN]*dzdx; 
+    // tilt correction
+    if(tilt) yc[fN] -= fTilt*(zc[fN] - zt); 
+
+    // elaborate cluster error
+    Float_t qr = c->GetQ() - q0;
+    sy[fN]   = qr < 0. ? clSigma0*TMath::Exp(clSlopeQ*qr) : clSigma0;
+
+    fitterY.AddPoint(&xc[fN], yc[fN]-yt, sy[fN]);
+
+    sz[fN]   = fPadLength*convert;
+    fitterZ.AddPoint(&xc[fN], zc[fN], sz[fN]);
+    fN++;
   }
   // to few clusters
-  if (nc < kClmin) return kFALSE; 
-  
+  if (fN < kClmin) return kFALSE; 
+
+  // fit XY plane
+  fitterY.Eval();
+  fYfit[0] = y0+fitterY.GetFunctionParameter(0);
+  fYfit[1] = dydx-fitterY.GetFunctionParameter(1);
 
-  Int_t zN[2*35];
-  Int_t nz = AliTRDtrackerV1::Freq(nc, zRow, zN, kFALSE);
+  // check par row crossing
+  Int_t zN[2*AliTRDseed::knTimebins];
+  Int_t nz = AliTRDtrackerV1::Freq(fN, zRow, zN, kFALSE);
   // more than one pad row crossing
   if(nz>2) return kFALSE; 
-  
-  // estimate reference parameter at average x
-  Double_t y0 = fYref[0];
-  Double_t dydx = fYref[1]; 
-  Double_t dzdx = fZref[1];
-  zc[nc]  = fZref[0];
+
 
   // determine z offset of the fit
+  Float_t zslope = 0.;
   Int_t nchanges = 0, nCross = 0;
   if(nz==2){ // tracklet is crossing pad row
     // Find the break time allowing one chage on pad-rows
     // with maximal number of accepted clusters
     Int_t padRef = zRow[0];
-    for (Int_t ic=1; ic<nc; ic++) {
+    for (Int_t ic=1; ic<fN; ic++) {
       if(zRow[ic] == padRef) continue;
       
       // debug
@@ -710,12 +724,12 @@ Bool_t AliTRDseedV1::Fit()
     
       // evaluate parameters of the crossing point
       Float_t sx = (xc[ic-1] - xc[ic])*convert;
-      xc[nc] = .5 * (xc[ic-1] + xc[ic]);
-      zc[nc] = .5 * (zc[ic-1] + zc[ic]);
-      sz[nc] = TMath::Max(dzdx * sx, .01);
-      dzdx   = zc[ic-1] > zc[ic] ? 1. : -1.;
-      padRef = zRow[ic];
-      nCross = ic;
+      fCross[0] = .5 * (xc[ic-1] + xc[ic]);
+      fCross[2] = .5 * (zc[ic-1] + zc[ic]);
+      fCross[3] = TMath::Max(dzdx * sx, .01);
+      zslope    = zc[ic-1] > zc[ic] ? 1. : -1.;
+      padRef    = zRow[ic];
+      nCross    = ic;
       nchanges++;
     }
   }
@@ -723,44 +737,25 @@ Bool_t AliTRDseedV1::Fit()
   // condition on nCross and reset nchanges TODO
 
   if(nchanges==1){
-    if(dzdx * fZref[1] < 0.){
+    if(dzdx * zslope < 0.){
       AliInfo("tracklet direction does not correspond to the track direction. TODO.");
     }
     SetBit(kRowCross, kTRUE); // mark pad row crossing
-    fCross[0] = xc[nc]; fCross[2] = zc[nc]; fCross[3] = sz[nc]; 
-    fitterZ.AddPoint(&xc[nc], zc[nc], sz[nc]);
+    fitterZ.AddPoint(&fCross[0], fCross[2], fCross[3]);
     fitterZ.Eval();
-    dzdx = fZref[1]; // we don't trust Parameter[1] ??;
-    zc[nc] = fitterZ.GetFunctionParameter(0); 
+    //zc[nc] = fitterZ.GetFunctionParameter(0); 
+    fCross[1] = fYfit[0] - fCross[0] * fYfit[1];
+    fCross[0] = fX0 - fCross[0];
   } else if(nchanges > 1){ // debug
-    AliInfo("ERROR in n changes!!!");
+    AliError("N pad row crossing > 1.");
     return kFALSE;
   }
 
-  
-  // estimate deviation from reference direction
-  dzdx *= fTilt;
-  for (Int_t ic=0; ic<nc; ic++) {
-    yc[ic] -= y0 + xc[ic]*(dydx + dzdx) + fTilt * (zc[ic] - zc[nc]);
-    fitterY.AddPoint(&xc[ic], yc[ic], sy[ic]);
-  }
-  fitterY.Eval();
-  fYfit[0] = y0+fitterY.GetFunctionParameter(0);
-  fYfit[1] = dydx+fitterY.GetFunctionParameter(1);
-  if(nchanges) fCross[1] = fYfit[0] + fCross[0] * fYfit[1];
-
-//     printf("\nnz = %d\n", nz);
-//     for(int ic=0; ic<35; ic++) printf("%d row[%d]\n", ic, zRow[ic]);        
-// 
-//     for(int ic=0; ic<nz; ic++) printf("%d n[%d]\n", ic, zN[ic]);    
+  UpdateUsed();
 
   return kTRUE;
 }
 
-//___________________________________________________________________
-void AliTRDseedV1::Draw(Option_t*)
-{
-}
 
 //___________________________________________________________________
 void AliTRDseedV1::Print(Option_t*) const
@@ -769,35 +764,27 @@ void AliTRDseedV1::Print(Option_t*) const
   // Printing the seedstatus
   //
 
-  printf("Seed status :\n");
-  printf("  fTilt      = %f\n", fTilt);
-  printf("  fPadLength = %f\n", fPadLength);
-  printf("  fX0        = %f\n", fX0);
-  for(int ic=0; ic<AliTRDtrackerV1::GetNTimeBins(); ic++) {
-          const Char_t *isUsable = fUsable[ic]?"Yes":"No";
-    printf("  %d X[%f] Y[%f] Z[%f] Indexes[%d] clusters[%p] usable[%s]\n"
-                , ic
-                , fX[ic]
-                , fY[ic]
-                , fZ[ic]
-                , fIndexes[ic]
-                , ((void*) fClusters[ic])
-                , isUsable);
-        }
+  AliInfo(Form("Tracklet X0[%7.2f] Det[%d]", fX0, fDet));
+  printf("  Tilt[%+6.2f] PadLength[%5.2f]\n", fTilt, fPadLength);
+  AliTRDcluster* const* jc = &fClusters[0];
+  for(int ic=0; ic<AliTRDtrackerV1::GetNTimeBins(); ic++, jc++) {
+    if(!(*jc)) continue;
+    printf("  %2d X[%7.2f] Y[%7.2f] Z[%7.2f] Idx[%d] c[%p] usable[%s]\n", 
+      ic, (*jc)->GetX(), (*jc)->GetY(), (*jc)->GetZ(), 
+      fIndexes[ic], (void*)(*jc), fUsable[ic]?"y":"n");
+  }
 
   printf("  fYref[0] =%f fYref[1] =%f\n", fYref[0], fYref[1]);
   printf("  fZref[0] =%f fZref[1] =%f\n", fZref[0], fZref[1]);
   printf("  fYfit[0] =%f fYfit[1] =%f\n", fYfit[0], fYfit[1]);
-  printf("  fYfitR[0]=%f fYfitR[1]=%f\n", fYfitR[0], fYfitR[1]);
   printf("  fZfit[0] =%f fZfit[1] =%f\n", fZfit[0], fZfit[1]);
-  printf("  fZfitR[0]=%f fZfitR[1]=%f\n", fZfitR[0], fZfitR[1]);
   printf("  fSigmaY =%f\n", fSigmaY);
   printf("  fSigmaY2=%f\n", fSigmaY2);            
   printf("  fMeanz  =%f\n", fMeanz);
   printf("  fZProb  =%f\n", fZProb);
   printf("  fLabels[0]=%d fLabels[1]=%d\n", fLabels[0], fLabels[1]);
   printf("  fN      =%d\n", fN);
-  printf("  fN2     =%d (>8 isOK)\n",fN2);
+  printf("  fN2     =%d (>4 isOK - to be redesigned)\n",fN2);
   printf("  fNUsed  =%d\n", fNUsed);
   printf("  fFreq   =%d\n", fFreq);
   printf("  fNChange=%d\n",  fNChange);