]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - HLT/hough/AliL3Hough.cxx
Remove warnings with gcc4.0 (F.Carminati).
[u/mrichter/AliRoot.git] / HLT / hough / AliL3Hough.cxx
index 659c2dedc137ed99a2681c76d808af1d0d6225fe..c1a544f4f31a4a773d9405ac6590c88b86d82f05 100644 (file)
@@ -31,6 +31,7 @@
 #include "AliL3TrackArray.h"
 #include "AliL3HoughTrack.h"
 #include "AliL3DDLDataFileHandler.h"
+#include "AliL3HoughKalmanTrack.h"
 
 #include "TThread.h"
 
@@ -91,6 +92,7 @@ AliL3Hough::AliL3Hough()
   fPeakRatio   = 0.5;
   fInputFile   = 0;
   fInputPtr    = 0;
+  fRawEvent    = 0;
   
   SetTransformerParams();
   SetThreshold();
@@ -199,6 +201,17 @@ void AliL3Hough::CleanUp()
   //cout << "Cleaned class mem " << endl;
 }
 
+void AliL3Hough::Init(Int_t netasegments,Int_t tv,AliRawEvent *rawevent,Float_t zvertex)
+{
+  //Normal constructor
+  fNEtaSegments  = netasegments;
+  fVersion       = tv;
+  fRawEvent      = rawevent;
+  fZVertex       = zvertex;
+
+  Init();
+}
+
 void AliL3Hough::Init(Char_t *path,Bool_t binary,Int_t netasegments,Bool_t bit8,Int_t tv,Char_t *infile,Char_t *ptr,Float_t zvertex)
 {
   //Normal init of the AliL3Hough
@@ -282,35 +295,42 @@ void AliL3Hough::Init(Bool_t doit, Bool_t addhists)
       else
 #ifdef use_aliroot
        {
-         if(!fInputFile) {
-           if(!fInputPtr) {
-             /* In case of reading digits file */
-             fMemHandler[i] = new AliL3FileHandler(kTRUE); //use static index
-             if(!fBinary) {
+         if(!fRawEvent) {
+           if(!fInputFile) {
+             if(!fInputPtr) {
+               /* In case of reading digits file */
+               fMemHandler[i] = new AliL3FileHandler(kTRUE); //use static index
+               if(!fBinary) {
 #if use_newio
-               if(!fRunLoader) {
+                 if(!fRunLoader) {
 #endif
-                 Char_t filename[1024];
-                 sprintf(filename,"%s/digitfile.root",fPath);
-                 fMemHandler[i]->SetAliInput(filename);
+                   Char_t filename[1024];
+                   sprintf(filename,"%s/digitfile.root",fPath);
+                   fMemHandler[i]->SetAliInput(filename);
 #if use_newio
-               }
-               else {
-                 fMemHandler[i]->SetAliInput(fRunLoader);
-               }
+                 }
+                 else {
+                   fMemHandler[i]->SetAliInput(fRunLoader);
+                 }
 #endif
+               }
+             }
+             else {
+               /* In case of reading from DATE */
+               fMemHandler[i] = new AliL3DDLDataFileHandler();
+               fMemHandler[i]->SetReaderInput(fInputPtr,-1);
              }
            }
            else {
-             /* In case of reading from DATE */
+             /* In case of reading rawdata from ROOT file */
              fMemHandler[i] = new AliL3DDLDataFileHandler();
-             fMemHandler[i]->SetReaderInput(fInputPtr,-1);
+             fMemHandler[i]->SetReaderInput(fInputFile);
            }
          }
          else {
-           /* In case of reading rawdata from ROOT file */
+           /* In case of reading rawdata using AliRawEvent */
            fMemHandler[i] = new AliL3DDLDataFileHandler();
-           fMemHandler[i]->SetReaderInput(fInputFile);
+           fMemHandler[i]->SetReaderInput(fRawEvent);
          }
        }
 #else
@@ -319,8 +339,14 @@ void AliL3Hough::Init(Bool_t doit, Bool_t addhists)
     }
 
   fPeakFinder = new AliL3HoughMaxFinder("KappaPhi",50000);
-  fMerger = new AliL3HoughMerger(fNPatches);
-  fInterMerger = new AliL3HoughIntMerger();
+  if(fVersion!=4) {
+    fMerger = new AliL3HoughMerger(fNPatches);
+    fInterMerger = new AliL3HoughIntMerger();
+  }
+  else {
+    fMerger = 0;
+    fInterMerger = 0;
+  }
   fGlobalMerger = 0;
   fBenchmark = new AliL3Benchmark();
 }
@@ -328,6 +354,9 @@ void AliL3Hough::Init(Bool_t doit, Bool_t addhists)
 void AliL3Hough::SetTransformerParams(Float_t ptres,Float_t ptmin,Float_t ptmax,Int_t ny,Int_t patch)
 {
   // Setup the parameters for the Hough Transformer
+  // This includes the bin size and limits for
+  // the parameter space histograms
+
   Int_t mrow;
   Float_t psi=0;
   if(patch==-1)
@@ -390,29 +419,55 @@ void AliL3Hough::SetTransformerParams(Int_t nx,Int_t ny,Float_t ptmin,Int_t /*pa
 {
   // Setup the parameters for the Hough Transformer
 
-
-  Int_t mrow=79;
-  Double_t lineradius = sqrt(pow(AliL3Transform::Row2X(mrow),2) + pow(AliL3Transform::GetMaxY(mrow),2));
-  Double_t alpha1 = AliL3Transform::GetMaxY(mrow)/pow(lineradius,2);
-  Double_t kappa = 1*AliL3Transform::GetBField()*AliL3Transform::GetBFact()/ptmin;
+  Double_t lineradius = 1.0/(AliL3HoughTransformerRow::GetBeta1()*sqrt(1.0+tan(AliL3Transform::Pi()*10/180)*tan(AliL3Transform::Pi()*10/180)));
+  Double_t alpha1 = AliL3HoughTransformerRow::GetBeta1()*tan(AliL3Transform::Pi()*10/180);
+  Double_t kappa = 1*AliL3Transform::GetBField()*AliL3Transform::GetBFact()/(ptmin*0.9);
   Double_t psi = AliL3Transform::Deg2Rad(10) - asin(lineradius*kappa/2);
   //  cout<<"Calculated psi range "<<psi<<" in patch "<<patch<<endl;
-  AliL3HoughTrack track;
-  track.SetTrackParameters(kappa,psi,1);
-  Float_t hit[3];
-  Int_t mrow2 = 158;
-  track.GetCrossingPoint(mrow2,hit);
-  Double_t lineradius2 = sqrt(pow(AliL3Transform::Row2X(mrow2),2) + pow(AliL3Transform::GetMaxY(mrow2),2));
-  Double_t alpha2 = hit[1]/pow(lineradius2,2);
+  Double_t alpha2 = alpha1 - (AliL3HoughTransformerRow::GetBeta1()-AliL3HoughTransformerRow::GetBeta2())*tan(psi);
   //  cout<<"Calculated alphas range "<<alpha1<<" "<<alpha2<<" in patch "<<patch<<endl;
 
   Int_t i=0;
   while(i < 6)
     {
-      fLowPt[i] = 1.15*alpha1;
+      fLowPt[i] = 1.1*alpha1;
       fNBinY[i] = ny;
       fNBinX[i] = nx;
-      fPhi[i] = 1.15*alpha2;
+      fPhi[i] = alpha2;
+      i++;
+    }
+}
+
+void AliL3Hough::CalcTransformerParams(Float_t ptmin)
+{
+  // Setup the parameters for the Row Hough Transformer
+  // Automatically adjusts the number of bins in X and Y in a way
+  // that the size of the hough bin is 2x (in X) and 2.5 (in Y) the
+  // size of the tpc pads
+
+  Double_t lineradius = 1.0/(AliL3HoughTransformerRow::GetBeta1()*sqrt(1.0+tan(AliL3Transform::Pi()*10/180)*tan(AliL3Transform::Pi()*10/180)));
+  Double_t alpha1 = AliL3HoughTransformerRow::GetBeta1()*tan(AliL3Transform::Pi()*10/180);
+  Double_t kappa = 1*AliL3Transform::GetBField()*AliL3Transform::GetBFact()/(ptmin*0.9);
+  Double_t psi = AliL3Transform::Deg2Rad(10) - asin(lineradius*kappa/2);
+  //  cout<<"Calculated psi range "<<psi<<endl;
+  Double_t alpha2 = alpha1 - (AliL3HoughTransformerRow::GetBeta1()-AliL3HoughTransformerRow::GetBeta2())*tan(psi);
+  alpha1 *= 1.1;
+  //  cout<<"Calculated alphas range "<<alpha1<<" "<<alpha2<<endl;
+
+  Double_t sizex = 2.0*AliL3Transform::GetPadPitchWidthLow()*AliL3HoughTransformerRow::GetBeta1()*AliL3HoughTransformerRow::GetBeta1();
+  Double_t sizey = 2.5*AliL3Transform::GetPadPitchWidthUp()*AliL3HoughTransformerRow::GetBeta2()*AliL3HoughTransformerRow::GetBeta2();
+
+  Int_t nx = 2*(Int_t)(alpha1/sizex)+1;
+  Int_t ny = 2*(Int_t)(alpha2/sizey)+1;
+  //  cout<<"Calculated number of bins "<<nx<<" "<<ny<<endl;
+
+  Int_t i=0;
+  while(i < 6)
+    {
+      fLowPt[i] = alpha1;
+      fNBinY[i] = ny;
+      fNBinX[i] = nx;
+      fPhi[i] = alpha2;
       i++;
     }
 }
@@ -768,6 +823,13 @@ void AliL3Hough::AddAllHistogramsRows()
 
 void AliL3Hough::PrepareForNextPatch(Int_t nextpatch)
 {
+  // Prepare the parameter space for the processing of
+  // the next read patch. According to the already
+  // accumulated number of gaps in parameter space
+  // bins, the routine updates the dynamic
+  // pointers used in order to jump rapidly during the
+  // filling of the parameter space.
+
   char buf[256];
   sprintf(buf,"Prepare For Patch %d",nextpatch);
   fBenchmark->Start(buf);
@@ -801,34 +863,36 @@ void AliL3Hough::PrepareForNextPatch(Int_t nextpatch)
       Int_t ymax = hist->GetLastYbin();
       Int_t nxbins = hist->GetNbinsX()+2;
 
-      UChar_t lastyvalue = 0;
-      Int_t endybin = ymin - 1;
-      for(Int_t ybin=ymin; ybin<=ymax; ybin++)
-       {
-         UChar_t lastxvalue = 0;
-         UChar_t maxvalue = 0;
-         Int_t endxbin = xmin - 1;
-         for(Int_t xbin=xmin; xbin<=xmax; xbin++)
-           {
-             Int_t bin = xbin + ybin*nxbins;
-             UChar_t value = 0;
-             if(gapcount[bin] < MAX_N_GAPS) {
-               value = 1;
-               maxvalue = 1;
-               if(tracklastrow[bin] > lastpatchlastrow) {
-                 if(lastpatchlastrow > currentrowcount[bin])
-                   gapcount[bin] += (lastpatchlastrow-currentrowcount[bin]-1);
-               }
-               else {
-                 if(tracklastrow[bin] > currentrowcount[bin])
-                   gapcount[bin] += (tracklastrow[bin]-currentrowcount[bin]-1);
+      if(fLastPatch != -1) {
+       UChar_t lastyvalue = 0;
+       Int_t endybin = ymin - 1;
+       for(Int_t ybin=nextrow[ymin]; ybin<=ymax; ybin = nextrow[++ybin])
+         {
+           UChar_t lastxvalue = 0;
+           UChar_t maxvalue = 0;
+           Int_t endxbin = xmin - 1;
+           for(Int_t xbin=xmin; xbin<=xmax; xbin++)
+             {
+               Int_t bin = xbin + ybin*nxbins;
+               UChar_t value = 0;
+               if(gapcount[bin] < MAX_N_GAPS) {
+                 if(tracklastrow[bin] > lastpatchlastrow) {
+                   if(lastpatchlastrow > currentrowcount[bin])
+                     gapcount[bin] += (lastpatchlastrow-currentrowcount[bin]-1);
+                 }
+                 else {
+                   if(tracklastrow[bin] > currentrowcount[bin])
+                     gapcount[bin] += (tracklastrow[bin]-currentrowcount[bin]-1);
+                 }
+                 if(gapcount[bin] < MAX_N_GAPS) {
+                   value = 1;
+                   maxvalue = 1;
+                   if(trackfirstrow[bin] < nextpatchfirstrow)
+                     currentrowcount[bin] = nextpatchfirstrow;
+                   else
+                     currentrowcount[bin] = trackfirstrow[bin];
+                 }
                }
-               if(trackfirstrow[bin] < nextpatchfirstrow)
-                 currentrowcount[bin] = nextpatchfirstrow;
-               else
-                 currentrowcount[bin] = trackfirstrow[bin];
-             }
-             if(fLastPatch != -1) {
                if(value > 0)
                  {
                    nextbin[xbin + ybin*nxbins] = (UChar_t)xbin;
@@ -846,8 +910,6 @@ void AliL3Hough::PrepareForNextPatch(Int_t nextpatch)
                  }
                lastxvalue = value;
              }
-           }
-         if(fLastPatch != -1) {
            UChar_t *tempnextbin = nextbin + endxbin + 1 + ybin*nxbins;
            memset(tempnextbin,(UChar_t)(xmax+1),xmax-endxbin);
            if(maxvalue > 0)
@@ -862,10 +924,40 @@ void AliL3Hough::PrepareForNextPatch(Int_t nextpatch)
              }
            lastyvalue = maxvalue;
          }
-       }
-      if(fLastPatch != -1) {
        UChar_t *tempnextrow = nextrow + endybin + 1;
-       memset(tempnextrow,(UChar_t)(ymax+1),ymax-endybin);
+       memset(tempnextrow,(UChar_t)(ymax+1),ymax-endybin+1);
+      }
+      else {
+       UChar_t lastyvalue = 0;
+       Int_t endybin = ymin - 1;
+       for(Int_t ybin=ymin; ybin<=ymax; ybin++)
+         {
+           UChar_t maxvalue = 0;
+           for(Int_t xbin=xmin; xbin<=xmax; xbin++)
+             {
+               Int_t bin = xbin + ybin*nxbins;
+               if(gapcount[bin] < MAX_N_GAPS) {
+                 maxvalue = 1;
+                 if(trackfirstrow[bin] < nextpatchfirstrow)
+                   currentrowcount[bin] = nextpatchfirstrow;
+                 else
+                   currentrowcount[bin] = trackfirstrow[bin];
+               }
+             }
+           if(maxvalue > 0)
+             {
+               nextrow[ybin] = (UChar_t)ybin;
+               if(maxvalue > lastyvalue)
+                 {
+                   UChar_t *tempnextrow = nextrow + endybin + 1;
+                   memset(tempnextrow,(UChar_t)ybin,ybin-endybin-1);
+                 }
+               endybin = ybin;
+             }
+           lastyvalue = maxvalue;
+         }
+       UChar_t *tempnextrow = nextrow + endybin + 1;
+       memset(tempnextrow,(UChar_t)(ymax+1),ymax-endybin+1);
       }
     }
 
@@ -915,6 +1007,11 @@ void AliL3Hough::FindTrackCandidatesRow()
   for(Int_t i=0; i<npatches; i++)
     {
       AliL3HoughBaseTransformer *tr = fHoughTransformer[i];
+      AliL3Histogram *h = tr->GetHistogram(0);
+      Float_t deltax = h->GetBinWidthX()*AliL3HoughTransformerRow::GetDAlpha();
+      Float_t deltay = h->GetBinWidthY()*AliL3HoughTransformerRow::GetDAlpha();
+      Float_t deltaeta = (tr->GetEtaMax()-tr->GetEtaMin())/tr->GetNEtaSegments()*AliL3HoughTransformerRow::GetDEta();
+      Float_t zvertex = tr->GetZVertex();
       fTracks[i]->Reset();
       fPeakFinder->Reset();
       
@@ -924,38 +1021,36 @@ void AliL3Hough::FindTrackCandidatesRow()
          if(hist->GetNEntries()==0) continue;
          fPeakFinder->SetHistogram(hist);
          fPeakFinder->SetEtaSlice(j);
-         fPeakFinder->SetTrackLUTs(((AliL3HoughTransformerRow *)tr)->fTrackNRows,((AliL3HoughTransformerRow *)tr)->fTrackFirstRow,((AliL3HoughTransformerRow *)tr)->fTrackLastRow);
+         fPeakFinder->SetTrackLUTs(((AliL3HoughTransformerRow *)tr)->GetTrackNRows(),((AliL3HoughTransformerRow *)tr)->GetTrackFirstRow(),((AliL3HoughTransformerRow *)tr)->GetTrackLastRow(),((AliL3HoughTransformerRow *)tr)->GetNextRow(j));
 #ifdef do_mc
          LOG(AliL3Log::kInformational,"AliL3Hough::FindTrackCandidates()","")
            <<"Starting "<<j<<" etaslice"<<ENDLOG;
 #endif
          fPeakFinder->SetThreshold(fPeakThreshold[i]);
          fPeakFinder->FindAdaptedRowPeaks(1,0,0);//Maxima finder for HoughTransformerRow
-
-         //fPeakFinder->FindMaxima(fPeakThreshold[i]); //Simple maxima finder
        }
   
       for(Int_t k=0; k<fPeakFinder->GetEntries(); k++)
        {
-         if(fPeakFinder->GetWeight(k) < 0) continue;
+         //      if(fPeakFinder->GetWeight(k) < 0) continue;
          AliL3HoughTrack *track = (AliL3HoughTrack*)fTracks[i]->NextTrack();
-         Float_t psi = atan((fPeakFinder->GetXPeak(k)-fPeakFinder->GetYPeak(k))/(AliL3HoughTransformerRow::GetBeta1()-AliL3HoughTransformerRow::GetBeta2()));
-         Float_t kappa = 2.0*(fPeakFinder->GetXPeak(k)*cos(psi)-AliL3HoughTransformerRow::GetBeta1()*sin(psi));
-         //          track->SetTrackParameters(fPeakFinder->GetXPeak(k),fPeakFinder->GetYPeak(k),fPeakFinder->GetWeight(k));
-         track->SetTrackParameters(kappa,psi,fPeakFinder->GetWeight(k));
+         Double_t starteta = tr->GetEta(fPeakFinder->GetStartEta(k),fCurrentSlice);
+         Double_t endeta = tr->GetEta(fPeakFinder->GetEndEta(k),fCurrentSlice);
+         Double_t eta = (starteta+endeta)/2.0;
+         track->SetTrackParametersRow(fPeakFinder->GetXPeak(k),fPeakFinder->GetYPeak(k),eta,fPeakFinder->GetWeight(k));
+         track->SetPterr(deltax); track->SetPsierr(deltay); track->SetTglerr(deltaeta);
          track->SetBinXY(fPeakFinder->GetXPeak(k),fPeakFinder->GetYPeak(k),fPeakFinder->GetXPeakSize(k),fPeakFinder->GetYPeakSize(k));
+         track->SetZ0(zvertex);
          Int_t etaindex = (fPeakFinder->GetStartEta(k)+fPeakFinder->GetEndEta(k))/2;
          track->SetEtaIndex(etaindex);
-         Float_t starteta = tr->GetEta(fPeakFinder->GetStartEta(k),fCurrentSlice);
-         Float_t endeta = tr->GetEta(fPeakFinder->GetEndEta(k),fCurrentSlice);
-         track->SetEta((starteta+endeta)/2.0);
-         track->SetRowRange(AliL3Transform::GetFirstRow(0),AliL3Transform::GetLastRow(5));
+         Int_t rows[2];
+         ((AliL3HoughTransformerRow *)tr)->GetTrackLength(fPeakFinder->GetXPeak(k),fPeakFinder->GetYPeak(k),rows);
+         track->SetRowRange(rows[0],rows[1]);
          track->SetSector(fCurrentSlice);
          track->SetSlice(fCurrentSlice);
 #ifdef do_mc
          Int_t label = tr->GetTrackID(etaindex,fPeakFinder->GetXPeak(k),fPeakFinder->GetYPeak(k));
          track->SetMCid(label);
-         //      cout<<"Track found with label "<<label<<" at "<<fPeakFinder->GetXPeak(k)<<" "<<fPeakFinder->GetYPeak(k)<<" with weight "<<fPeakFinder->GetWeight(k)<<endl; 
 #endif
        }
       LOG(AliL3Log::kInformational,"AliL3Hough::FindTrackCandidates()","")
@@ -1211,6 +1306,42 @@ void AliL3Hough::WriteTracks(Int_t slice,Char_t *path)
     }
 }
 
+#ifdef use_aliroot
+Int_t AliL3Hough::FillESD(AliESD *esd)
+{
+  // Fill the found hough transform tracks
+  // into the ESD. The tracks are stored as
+  // AliESDHLTtrack objects.
+
+  // No TPC PID so far,assuming pions
+  Double_t prob[AliPID::kSPECIES];
+  for(Int_t i=0;i<AliPID::kSPECIES;i++) {
+    if(i==AliPID::kPion) prob[i]=1.0;
+    else prob[i]=0.1;
+  }
+
+  if(!fGlobalTracks) return 0;
+  Int_t nglobaltracks = 0;
+  for(Int_t i=0; i<fGlobalTracks->GetNTracks(); i++)
+    {
+      AliL3HoughTrack *tpt = (AliL3HoughTrack *)fGlobalTracks->GetCheckedTrack(i);
+      if(!tpt) continue; 
+      
+      if(tpt->GetWeight()<0) continue;
+      AliL3HoughKalmanTrack *tpctrack = new AliL3HoughKalmanTrack(*tpt);
+      if(!tpctrack) continue;
+      AliESDtrack *esdtrack2 = new AliESDtrack() ; 
+      esdtrack2->UpdateTrackParams(tpctrack,AliESDtrack::kTPCin);
+      esdtrack2->SetESDpid(prob);
+      esd->AddTrack(esdtrack2);
+      nglobaltracks++;
+      delete esdtrack2;
+      delete tpctrack;
+    }
+  return nglobaltracks;
+}
+#endif
+
 void AliL3Hough::WriteDigits(Char_t *outfile)
 {
   //Write the current data to a new rootfile.
@@ -1237,6 +1368,9 @@ Double_t AliL3Hough::GetCpuTime()
 
 void *AliL3Hough::ProcessInThread(void *args)
 {
+  // Called in case Hough transform tracking
+  // is executed in a thread
+
   AliL3Hough *instance = (AliL3Hough *)args;
   Int_t minslice = instance->GetMinSlice();
   Int_t maxslice = instance->GetMaxSlice();
@@ -1253,6 +1387,10 @@ void *AliL3Hough::ProcessInThread(void *args)
 
 void AliL3Hough::StartProcessInThread(Int_t minslice,Int_t maxslice)
 {
+  // Starts the Hough transform tracking as a
+  // separate thread. Takes as parameters the
+  // range of TPC slices (sectors) to be reconstructed
+
   if(!fThread) {
     char buf[255];
     sprintf(buf,"houghtrans_%d_%d",minslice,maxslice);
@@ -1266,5 +1404,14 @@ void AliL3Hough::StartProcessInThread(Int_t minslice,Int_t maxslice)
 
 Int_t AliL3Hough::WaitForThreadFinish()
 {
+  // Routine is used in case we run the
+  // Hough transform tracking in several
+  // threads and want to sync them before
+  // writing the results to the ESD
+
+#if ROOT_VERSION_CODE < 262403
   return TThread::Join(fThread->GetId());
+#else
+  return fThread->Join(fThread->GetId());
+#endif
 }