]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - EVE/EveDet/AliEveTRDData.cxx
AliTPCkalmanFit.h AliTPCkalmanFit.cxx -
[u/mrichter/AliRoot.git] / EVE / EveDet / AliEveTRDData.cxx
index 60d50d618b7581906fb6c508975be1c6ca118ba5..9f58f749494d4b32d1c8cc6b828948be5b1bd986 100644 (file)
@@ -9,23 +9,32 @@
 
 #include "TVector.h"
 #include "TLinearFitter.h"
+
 #include "TEveTrans.h"
+#include "TEveManager.h"
+
+#include "EveBase/AliEveEventManager.h"
 
 #include "AliEveTRDData.h"
 #include "AliEveTRDModuleImp.h"
+#include "AliEveTRDLoader.h"
+#include "AliEveTRDLoaderImp.h"
 
 #include "AliLog.h"
 #include "AliPID.h"
 #include "AliTrackPointArray.h"
+#include "AliRieman.h"
 
 #include "AliTRDhit.h"
 #include "AliTRDcluster.h"
 #include "AliTRDseedV1.h"
+#include "AliTRDtrackletMCM.h"
+#include "AliTRDmcmSim.h"
 #include "AliTRDtrackV1.h"
 #include "AliTRDtrackerV1.h"
 #include "AliTRDpadPlane.h"
 #include "AliTRDdigitsManager.h"
-#include "AliTRDdataArrayDigits.h"
+#include "AliTRDarrayADC.h"
 #include "AliTRDSignalIndex.h"
 #include "AliTRDgeometry.h"
 #include "AliTRDtransform.h"
@@ -89,13 +98,13 @@ void AliEveTRDDigits::ComputeRepresentation()
     for (Int_t ic = 0; ic < ncols; ic++) {
       dy = pp->GetColSize(ic);
       for (Int_t it = 0; it < ntbs; it++) {
-        q = fData.GetDataUnchecked(ir, ic, it);
+        q = fData.GetData(ir, ic, it);
         if (q < threshold) continue;
 
         Double_t x[6] = {0., 0., Double_t(q), 0., 0., 0.}; 
-        Int_t  roc[3] = {ir, ic, 0}; 
-        Bool_t    out = kTRUE;
-        transform.Transform(&x[0], &roc[0], UInt_t(it), out, 0);
+        //Int_t  roc[3] = {ir, ic, 0}; 
+        //Bool_t    out = kTRUE;
+        //transform.Transform(&x[0], &roc[0], UInt_t(it), out, 0);
 
         scale = q < 512 ? q/512. : 1.;
         color  = 50+int(scale*50.);
@@ -107,7 +116,7 @@ void AliEveTRDDigits::ComputeRepresentation()
       }  // end time loop
     }  // end col loop
   }  // end row loop
-  fData.Compress(1);
+  fData.Compress();
   
   // rotate to global coordinates
   //RefitPlex();
@@ -121,8 +130,8 @@ void AliEveTRDDigits::SetData(AliTRDdigitsManager *digits)
   // Set data source.
 
   Int_t det = fParent->GetID();
-  AliTRDdataArrayDigits *data = digits->GetDigits(det);
-  if(!data->HasData()) return;
+  AliTRDarrayADC *data = digits->GetDigits(det);
+  if(!data->GetDim()) return;
   data->Expand();
 
   AliTRDSignalIndex *indexes = digits->GetIndexes(det);
@@ -142,12 +151,12 @@ void AliEveTRDDigits::SetData(AliTRDdigitsManager *digits)
       }
       adc = data->GetData(row, col, time);
       if(adc <= 1) continue;
-      fData.SetDataUnchecked(row, col, time, adc);
+      fData.SetData(row, col, time, adc);
       //fIndex->AddIndexTBin(row,col,time);
       //printf("\tr[%d] c[%d] t[%d] ADC[%d]\n", row, col, time, adc);
     } 
   }
-  fData.Compress(1);
+  fData.Compress();
 }
 
 
@@ -179,7 +188,7 @@ AliEveTRDHits::AliEveTRDHits() : TEvePointSet("hits", 20)
 {
   // Constructor.
   SetMarkerSize(.1);
-  SetMarkerColor(2);
+  SetMarkerColor(kGreen);
   SetOwnIds(kTRUE);
 }
 
@@ -194,13 +203,13 @@ void AliEveTRDHits::PointSelected(Int_t n)
 {
   // Handle an individual point selection from GL.
 
-  AliTRDhit *h = dynamic_cast<AliTRDhit*>(GetPointId(n));
-  printf("\nDetector             : %d\n", h->GetDetector());
-  printf("Region of production : %c\n", h->FromAmplification() ? 'A' : 'D');
-  printf("TR photon            : %s\n", h->FromTRphoton() ? "Yes" : "No");
-  printf("Charge               : %d\n", h->GetCharge());
-  printf("MC track label       : %d\n", h->GetTrack());
-  printf("Time from collision  : %f\n", h->GetTime());
+  AliTRDhit *h = 0x0;
+  if(!(h = dynamic_cast<AliTRDhit*>(GetPointId(n)))) return;
+  printf("Id[%3d] Det[%3d] Reg[%c] TR[%c] Q[%3d] MC[%d] t[%f]\n", 
+    n, h->GetDetector(), 
+    h->FromAmplification() ? 'A' : 'D', 
+    h->FromTRphoton() ? 'y' : 'n', 
+    h->GetCharge(), h->GetTrack(), h->GetTime());
 }
 
 
@@ -254,6 +263,76 @@ void AliEveTRDClusters::Print(Option_t *o) const
   }
 }
 
+//______________________________________________________________________________
+void AliEveTRDClusters::Load(Char_t *w, Bool_t stk) const
+{
+  Int_t typ = -1;
+  if(strcmp(w, "hit")==0) typ = 0;
+  else if(strcmp(w, "dig")==0) typ = 1;
+  else if(strcmp(w, "cls")==0) typ = 2;
+  else if(strcmp(w, "all")==0) typ = 3;
+  else{
+    AliInfo("The following arguments are accepted:");
+    AliInfo("   \"hit\" : loading of MC hits");
+    AliInfo("   \"dig\" : loading of digits");
+    AliInfo("   \"cls\" : loading of reconstructed clusters");
+    AliInfo("   \"all\" : loading of MC hits+digits+clusters");
+    return;
+  }
+
+  AliTRDcluster *c = 0x0;
+  Int_t n = 0;
+  while((n = GetN() && !(c = dynamic_cast<AliTRDcluster*>(GetPointId(n))))) n++;
+  if(!c) return;
+
+  Int_t det = c->GetDetector();
+  AliEveTRDLoader *loader = 0x0;
+  switch(typ){
+  case 0:  
+    loader = new AliEveTRDLoader("Hits");
+    if(!loader->Open("TRD.Hits.root")){ 
+      delete loader;
+      return;
+    }
+    loader->SetDataType(AliEveTRDLoader::kTRDHits);
+    break;
+  case 1:
+    loader = new AliEveTRDLoader("Digits");
+    if(!loader->Open("TRD.Digits.root")){ 
+      delete loader;
+      return;
+    }
+    loader->SetDataType(AliEveTRDLoader::kTRDDigits);
+    break;
+  case 2:
+    loader = new AliEveTRDLoader("Clusters");
+    if(!loader->Open("TRD.RecPoints.root")){ 
+      delete loader;
+      return;
+    }
+    loader->SetDataType(AliEveTRDLoader::kTRDClusters);
+    break;
+  case 3:
+    loader = new AliEveTRDLoaderSim("MC");
+    if(!loader->Open("galice.root")){ 
+      delete loader;
+      return;
+    }
+    loader->SetDataType(AliEveTRDLoader::kTRDHits | AliEveTRDLoader::kTRDDigits | AliEveTRDLoader::kTRDClusters);
+    break;
+  default: return;
+  }
+
+  loader->AddChambers(AliTRDgeometry::GetSector(det),AliTRDgeometry::GetStack(det), stk ? -1 : AliTRDgeometry::GetLayer(det));
+  // load first event
+  loader->GoToEvent(AliEveEventManager::GetCurrent()->GetEventId());
+  
+  // register loader with alieve
+  gEve->AddElement(loader);
+  //loader->SpawnEditor();
+  gEve->Redraw3D();
+}
+
 ///////////////////////////////////////////////////////////
 /////////////   AliEveTRDTracklet         /////////////////////
 ///////////////////////////////////////////////////////////
@@ -266,46 +345,111 @@ AliEveTRDTracklet::AliEveTRDTracklet(AliTRDseedV1 *trklt):TEveLine()
   SetName("tracklet");
   
   SetUserData(trklt);
-  Int_t det = -1, sec;
+  Float_t dx;
+  Float_t x0   = trklt->GetX0();
+  Float_t y0   = trklt->GetYref(0);
+  Float_t z0   = trklt->GetZref(0);
+  Float_t dydx = trklt->GetYref(1);
+  Float_t dzdx = trklt->GetZref(1);
+  Float_t tilt = trklt->GetTilt();
   Float_t g[3];
   AliTRDcluster *c = 0x0;
-  for(Int_t ic=0; ic<AliTRDseed::knTimebins; ic++){
+  for(Int_t ic=0; ic<AliTRDseedV1::kNclusters; ic++){
     if(!(c = trklt->GetClusters(ic))) continue;
     if(!fClusters) AddElement(fClusters = new AliEveTRDClusters());
-    det = c->GetDetector();
+    dx = x0 - c->GetX();
+    //Float_t yt = y0 - dx*dydx;
+    Float_t zt = z0 - dx*dzdx;
+    // backup yc - for testing purposes
+    Float_t yc = c->GetY(); 
+    c->SetY(yc-tilt*(c->GetZ()-zt));
     c->GetGlobalXYZ(g); 
     Int_t id = fClusters->SetNextPoint(g[0], g[1], g[2]);    
     //Int_t id = fClusters->SetNextPoint(c->GetX(), c->GetY(), c->GetZ());    
+    c->SetY(yc);
     fClusters->SetPointId(id, new AliTRDcluster(*c));
   } 
-  if(fClusters) fClusters->SetTitle(Form("N[%d]", trklt->GetN2()));
-
+  if(fClusters){
+    fClusters->SetTitle(Form("N[%d]", trklt->GetN2()));
+    fClusters->SetMarkerColor(kMagenta);
+  }
 
-  SetTitle(Form("Det[%d] Plane[%d] P[%7.3f]", det, trklt->GetPlane(), trklt->GetMomentum()));
+  SetTitle(Form("Det[%d] Plane[%d] P[%7.3f]", trklt->GetDetector(), trklt->GetPlane(), trklt->GetMomentum()));
   SetLineColor(kRed);
   //SetOwnIds(kTRUE);
   
-  sec = det/30;
+  // init tracklet line
+  Int_t sec = AliTRDgeometry::GetSector(trklt->GetDetector());
   Double_t alpha = AliTRDgeometry::GetAlpha() * (sec<9 ? sec + .5 : sec - 17.5); 
-  Double_t x0 = trklt->GetX0(), 
-    y0f = trklt->GetYfit(0), 
-    ysf = trklt->GetYfit(1),
-    z0r = trklt->GetZref(0), 
-    zsr = trklt->GetZref(1);
-  Double_t xg =  x0 * TMath::Cos(alpha) - y0f * TMath::Sin(alpha); 
-  Double_t yg = x0 * TMath::Sin(alpha) + y0f * TMath::Cos(alpha);
-  SetPoint(0, xg, yg, z0r);
-  //SetPoint(0, x0, y0f, z0r);
-
-
-  //SetPointId(0, new AliTRDseedV1(*trackletObj));
-  Double_t x1 = x0-3.5, 
-    y1f = y0f - ysf*3.5,
-    z1r = z0r - zsr*3.5; 
-  xg =  x1 * TMath::Cos(alpha) - y1f * TMath::Sin(alpha); 
-  yg = x1 * TMath::Sin(alpha) + y1f * TMath::Cos(alpha);
-  SetPoint(1, xg, yg, z1r);
-  //SetPoint(1, x1, y1f, z1r);
+
+  //trklt->Fit(kTRUE);
+  y0   = trklt->GetYfit(0);
+  dydx = trklt->GetYfit(1);
+  Double_t xg =  x0 * TMath::Cos(alpha) - y0 * TMath::Sin(alpha); 
+  Double_t yg = x0 * TMath::Sin(alpha) + y0 * TMath::Cos(alpha);
+  SetPoint(0, xg, yg, z0); 
+  //SetPoint(0, x0, y0, z0);
+
+
+  dx = .5*AliTRDgeometry::CamHght()+AliTRDgeometry::CdrHght();
+  x0 -= dx; 
+  y0 -= dydx*dx,
+  z0 -= dzdx*dx; 
+  xg = x0 * TMath::Cos(alpha) - y0 * TMath::Sin(alpha); 
+  yg = x0 * TMath::Sin(alpha) + y0 * TMath::Cos(alpha);
+  SetPoint(1, xg, yg, z0);
+  //SetPoint(1, x0, y0, z0);
+}
+
+AliEveTRDTracklet::AliEveTRDTracklet(AliTRDtrackletMCM *tracklet) : 
+  TEveLine(),
+  fClusters(0x0)
+{
+  SetName("TRD tracklet");
+  AliTRDtrackletMCM *trkl = new AliTRDtrackletMCM(*tracklet);
+  SetUserData(trkl);
+
+  SetTitle(Form("Det: %i, ROB: %i, MCM: %i, Label: %i\n0x%08x", 
+                trkl->GetDetector(), trkl->GetROB(), trkl->GetMCM(), trkl->GetLabel(),
+                trkl->GetTrackletWord()));
+  SetLineColor(kRed);
+
+  AliTRDgeometry *geo = new AliTRDgeometry();
+
+  Float_t length = 10;
+  Double_t x[3];
+  Double_t p[3];
+  x[0] = trkl->GetX();
+  x[1] = trkl->GetY();
+  x[2] = trkl->GetZ();
+  
+  geo->RotateBack(trkl->GetDetector(), x, p);
+  SetPoint(0, p[0], p[1], p[2]);
+
+  x[0] -= length;
+  x[1] += length * trkl->GetdYdX();
+  x[2] *= x[0] / (x[0] + length);
+  geo->RotateBack(trkl->GetDetector(), x, p);
+  SetPoint(1, p[0], p[1], p[2]);
+}
+
+AliEveTRDTracklet::~AliEveTRDTracklet() 
+{
+  AliTRDtrackletMCM *trkl = dynamic_cast<AliTRDtrackletMCM*> ((AliTRDtrackletBase*) GetUserData());
+  delete trkl;
+}
+
+void AliEveTRDTracklet::ShowMCM(Option_t *opt) const
+{
+  AliTRDtrackletMCM *trkl = dynamic_cast<AliTRDtrackletMCM*> ((AliTRDtrackletBase*) GetUserData());
+  if (!trkl)
+    return;
+  printf("Det: %3i, ROB: %i, MCM: %2i\n", trkl->GetDetector(), trkl->GetROB(), trkl->GetMCM());
+  AliTRDmcmSim *mcm = new AliTRDmcmSim();
+  AliRunLoader *rl = AliEveEventManager::AssertRunLoader();
+  mcm->LoadMCM(rl, trkl->GetDetector(), trkl->GetROB(), trkl->GetMCM());
+  mcm->Tracklet();
+  mcm->Draw(opt);
 }
 
 //______________________________________________________________________________
@@ -327,6 +471,7 @@ AliEveTRDTrack::AliEveTRDTrack(AliTRDtrackV1 *trk)
   ,fESDStatus(0)
   ,fAlpha(0.)
   ,fPoints(0x0)
+  ,fRim(0x0)
 {
   // Constructor.
   SetUserData(trk);
@@ -334,13 +479,27 @@ AliEveTRDTrack::AliEveTRDTrack(AliTRDtrackV1 *trk)
 
   AliTRDtrackerV1::SetNTimeBins(24);
 
+  fRim = new AliRieman(trk->GetNumberOfClusters());
   AliTRDseedV1 *tracklet = 0x0;
   for(Int_t il=0; il<AliTRDgeometry::kNlayer; il++){
     if(!(tracklet = trk->GetTracklet(il))) continue;
     if(!tracklet->IsOK()) continue;
     AddElement(new AliEveTRDTracklet(tracklet));
-  }
 
+    AliTRDcluster *c = 0x0;
+//     tracklet->ResetClusterIter(kFALSE);
+//     while((c = tracklet->PrevCluster())){
+    for(Int_t ic=AliTRDseedV1::kNtb; ic--;){
+      if(!(c=tracklet->GetClusters(ic))) continue;
+      Float_t xc = c->GetX();
+      Float_t yc = c->GetY();
+      Float_t zc = c->GetZ();
+      Float_t zt = tracklet->GetZref(0) - (tracklet->GetX0()-xc)*tracklet->GetZref(1); 
+      yc -= tracklet->GetTilt()*(zc-zt);
+      fRim->AddPoint(xc, yc, zc, .05, 2.3);
+    }
+  }
+  if(trk->GetNumberOfTracklets()>1) fRim->Update();
   SetStatus(fTrackState);
 }
 
@@ -351,6 +510,13 @@ AliEveTRDTrack::~AliEveTRDTrack()
   //delete dynamic_cast<AliTRDtrackV1*>(GetUserData());
 }
 
+//______________________________________________________________________________
+void AliEveTRDTrack::Print(Option_t *o) const
+{
+  AliTRDtrackV1 *track = (AliTRDtrackV1*)GetUserData();
+  if(!track) return;
+  track->Print(o);
+}
 
 //______________________________________________________________________________
 void AliEveTRDTrack::SetStatus(UChar_t s)
@@ -365,15 +531,23 @@ void AliEveTRDTrack::SetStatus(UChar_t s)
   if(!fPoints){ 
     fPoints = new AliTrackPoint[nc];
 
-    AliTRDcluster *c = trk->GetCluster(0);
-    Double_t x = c->GetX();
-    Int_t sec = c->GetDetector()/30;
-    fAlpha = AliTRDgeometry::GetAlpha() * (sec<9 ? sec + .5 : sec - 17.5); 
+    // define the radial span of the track in the TRD
+    Double_t xmin = -1., xmax = -1.;
+    Int_t det = 0;
+    AliTRDseedV1 *trklt = 0x0;
+    for(Int_t ily=AliTRDgeometry::kNlayer; ily--;){
+      if(!(trklt = trk->GetTracklet(ily))) continue;
+      if(xmin<0.) xmin = trklt->GetX0() - AliTRDgeometry::CamHght() - AliTRDgeometry::CdrHght();
+      if(xmax<0.) xmax = trklt->GetX0();
+      det = trklt->GetDetector();
+    }
+    Int_t sec = det/AliTRDgeometry::kNdets;
+    fAlpha = AliTRDgeometry::GetAlpha() * (sec<9 ? sec + .5 : sec - 17.5); //trk->GetAlpha()
 
-    Double_t dx = (trk->GetCluster(trk->GetNumberOfClusters()-1)->GetX()-x)/nc;
+    Double_t dx =(xmax - xmin)/nc;
     for(Int_t ip=0; ip<nc; ip++){
-      fPoints[ip].SetXYZ(x, 0., 0.);
-      x+=dx;
+      fPoints[ip].SetXYZ(xmin, 0., 0.);
+      xmin+=dx;
     }
     BUILD = kTRUE;
   }
@@ -389,12 +563,18 @@ void AliEveTRDTrack::SetStatus(UChar_t s)
         if(trk->GetNumberOfTracklets() >=4) AliTRDtrackerV1::FitKalman(trk, 0x0, kFALSE, nc, fPoints);
       } else { 
         //printf("Rieman track\n");
-        if(trk->GetNumberOfTracklets() >=4) AliTRDtrackerV1::FitRiemanTilt(trk, 0x0, kTRUE, nc, fPoints);
+        // if(trk->GetNumberOfTracklets() >=4) AliTRDtrackerV1::FitRiemanTilt(trk, 0x0, kTRUE, nc, fPoints);
+
+/*        Float_t x = 0.;
+        for(Int_t ip = nc; ip--;){
+          x = fPoints[ip].GetX();
+          fPoints[ip].SetXYZ(x, fRim->GetYat(x), fRim->GetZat(x));
+        }*/
       }
     }
   
     Float_t global[3];
-    for(Int_t ip=0; ip<nc; ip++){
+    for(Int_t ip=0; ip<0/*nc*/; ip++){
       fPoints[ip].Rotate(-fAlpha).GetXYZ(global);
       SetPoint(ip, global[0], global[1], global[2]);
     }
@@ -413,7 +593,7 @@ void AliEveTRDTrack::SetStatus(UChar_t s)
         SetLineColor(kMagenta);
       }
     } else {
-      if(TESTBIT(s, kPID) == AliTRDReconstructor::kLQPID){
+      if(TESTBIT(s, kPID) == AliTRDpidUtil::kLQ){
         //printf("PID color kLQPID\n");
         //trk->GetReconstructor()->SetOption("!nn");
       } else {
@@ -442,7 +622,7 @@ void AliEveTRDTrack::SetStatus(UChar_t s)
     SetLineWidth(2);
   }
   
-  Char_t *model = "line";
+  const Char_t *model = "line";
   if(!TESTBIT(s, kTrackCosmics)){
     if(TESTBIT(s, kTrackModel)) model = "kalman";
     else model = "rieman";