Checking in for the weekend. Compressing/uncompressing works. Restoring data - buildi...
authorvestbo <vestbo@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 15 Mar 2002 15:50:01 +0000 (15:50 +0000)
committervestbo <vestbo@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 15 Mar 2002 15:50:01 +0000 (15:50 +0000)
HLT/comp/AliL3Compress.cxx
HLT/comp/AliL3Compress.h
HLT/comp/AliL3ModelTrack.cxx
HLT/comp/AliL3ModelTrack.h
HLT/comp/AliL3Modeller.cxx
HLT/comp/AliL3Modeller.h

index 2b285f8..1e6fa8e 100644 (file)
@@ -6,23 +6,61 @@
 #include <stdio.h>
 #include <stream.h>
 #include <stdlib.h>
+#include <TH1.h>
+#include <TH2.h>
+#include <TRandom.h>
 
 #include "AliL3Compress.h"
 #include "AliL3TrackArray.h"
 #include "AliL3ModelTrack.h"
+#include "AliL3Transform.h"
+#include "AliL3MemHandler.h"
 #include "bitio.h"
 
+//_____________________________________________________________
+//
+//  AliL3Compress
+//
+// Class for compressing and uncompressing data.
+
 ClassImp(AliL3Compress)
 
 AliL3Compress::AliL3Compress()
 {
   fTracks=0;
+  SetBitNumbers(7,7,10,4);
+  fSlice =0;
+  fPatch=0;
+  fDigits=0;
+  fDPt=0;
+}
+
+AliL3Compress::AliL3Compress(Int_t slice,Int_t patch,Int_t pad,Int_t time,Int_t charge,Int_t shape)
+{
+  fSlice=slice;
+  fPatch=patch;
+  SetBitNumbers(pad,time,charge,shape);
+  fTracks=0;
+  fDigits=0;
+  fDPt=0;
 }
 
 AliL3Compress::~AliL3Compress()
 {
   if(fTracks)
     delete fTracks;
+  if(fDigits)
+    delete [] fDigits;
+  if(fDPt)
+    delete [] fDPt;
+}
+
+void AliL3Compress::SetBitNumbers(Int_t pad,Int_t time,Int_t charge,Int_t shape)
+{
+  fNumPadBits=pad;
+  fNumTimeBits=time;
+  fNumChargeBits=charge;
+  fNumShapeBits=shape;
 }
 
 void AliL3Compress::WriteFile(AliL3TrackArray *tracks,Char_t *filename)
@@ -38,6 +76,8 @@ void AliL3Compress::WriteFile(AliL3TrackArray *tracks,Char_t *filename)
     {
       AliL3ModelTrack *track = (AliL3ModelTrack*)tracks->GetCheckedTrack(i);
       if(!track) continue;
+            
+      track->FillModel();
       model = track->GetModel();
       if(model->fNClusters==0) continue;
       clusters = track->GetClusters();
@@ -45,7 +85,9 @@ void AliL3Compress::WriteFile(AliL3TrackArray *tracks,Char_t *filename)
       if(fwrite(model,sizeof(AliL3TrackModel),1,file)!=1) break;
       //cout<<"Writing "<<(int)model->fNClusters<<" clusters to file"<<endl;
       if(fwrite(clusters,model->fNClusters*sizeof(AliL3ClusterModel),1,file)!=1) break;
+      //track->Print();
       count++;
+      
     }
   cout<<"Wrote "<<count<<" tracks "<<endl;
   fclose(file);
@@ -54,21 +96,33 @@ void AliL3Compress::WriteFile(AliL3TrackArray *tracks,Char_t *filename)
 void AliL3Compress::ReadFile(Char_t *filename)
 {
   FILE *file = fopen(filename,"r");
-  
+  if(!file)
+    {
+      cerr<<"Cannot open file "<<filename<<endl;
+      return;
+    }
+
   if(fTracks)
     delete fTracks;
   fTracks = new AliL3TrackArray("AliL3ModelTrack");
   
+  cout<<"Reading file "<<filename<<endl;
   while(!feof(file))
     {
       AliL3ModelTrack *track = (AliL3ModelTrack*)fTracks->NextTrack();
-      track->Init(0,0);
+      track->Init(fSlice,fPatch);
       AliL3TrackModel *model = track->GetModel();
       AliL3ClusterModel *clusters = track->GetClusters();
+      //cout<<"Reading model "<<(int)model<<endl;
       if(fread(model,sizeof(AliL3TrackModel),1,file)!=1) break;
+      //cout<<"Reading clusters "<<(int)clusters<<endl;
       if(fread(clusters,(model->fNClusters)*sizeof(AliL3ClusterModel),1,file)!=1) break;
+      //cout<<"Filling track"<<endl;
+      track->FillTrack();
+      //track->Print();
     }
-  
+
+  fTracks->RemoveLast();
   cout<<"Read "<<fTracks->GetNTracks()<<" tracks from file"<<endl;
   fclose(file);
 }
@@ -82,7 +136,10 @@ void AliL3Compress::CompressFile(Char_t *infile,Char_t *outfile)
   AliL3TrackModel track;
   AliL3ClusterModel cluster;
   Int_t temp;
-
+  Short_t power;
+  
+  Int_t timeo,pado,chargeo,shapeo;
+  timeo=pado=chargeo=shapeo=0;
   while(!feof(input))
     {
       if(fread(&track,sizeof(AliL3TrackModel),1,input)!=1) break;
@@ -92,37 +149,90 @@ void AliL3Compress::CompressFile(Char_t *infile,Char_t *outfile)
          if(putc(output->rack,output->file )!=output->rack)
            cerr<<"AliL3Compress::ComressFile : Error writing to bitfile"<<endl;
          output->mask=0x80;
+         output->rack=0;
        }
+      
+      //Write track parameters:
       fwrite(&track,sizeof(AliL3TrackModel),1,output->file);
       for(Int_t i=0; i<track.fNClusters; i++)
        {
          if(fread(&cluster,sizeof(AliL3ClusterModel),1,input)!=1) break;
-         Int_t flag = (Int_t)cluster.fPresent;
-         OutputBit(output,flag);
-         if(!flag) continue;
+         
+         //Write empty flag:
+         temp = (Int_t)cluster.fPresent;
+         OutputBit(output,temp);
+         if(!temp) continue;
+         
+         //Write time information:
          temp = (Int_t)cluster.fDTime;
          if(temp<0)
            OutputBit(output,0);
          else
            OutputBit(output,1);
-         OutputBits(output,abs(temp),8);
+         power = 1<<fNumTimeBits;
+         if(abs(temp)>=power)
+           {
+             timeo++;
+             temp=power - 1;
+           }
+         OutputBits(output,abs(temp),fNumTimeBits);
+         
+         //Write pad information:
          temp = (Int_t)cluster.fDPad;
          if(temp<0)
            OutputBit(output,0);
          else
            OutputBit(output,1);
-         OutputBits(output,abs(temp),8);
+         power = 1<<fNumPadBits;
+         if(abs(temp)>=power)
+           {
+             pado++;
+             temp=power - 1;
+           }
+         OutputBits(output,abs(temp),fNumPadBits);
+         
+         //Write charge information:
          temp = (Int_t)cluster.fDCharge;
-         OutputBits(output,temp,10);
+         if(temp<0)
+           OutputBit(output,0);
+         else
+           OutputBit(output,1);
+         power = 1<<fNumChargeBits;
+         if(abs(temp)>=power)
+           {
+             chargeo++;
+             temp=power - 1;
+           }
+         OutputBits(output,abs(temp),fNumChargeBits);
          
-         //Short_t temp=(Short_t)cluster.fDTime;
-         // cout<<"flag "<<(int)flag<<" dtime "<<(int)cluster.fDTime<<" dpad "<<(int)cluster.fDPad<<" charge "<<cluster.fDCharge<<endl;
+         //Write shape information:
+         temp = (Int_t)cluster.fDSigmaY2;
+         power = 1<<fNumShapeBits;
+         if(abs(temp) >= power)
+           {
+             shapeo++;
+             temp = power - 1;
+           }
+         OutputBits(output,abs(temp),fNumShapeBits);
+         
+         temp = (Int_t)cluster.fDSigmaZ2;
+         if(abs(temp) >= power)
+           {
+             shapeo++;
+             temp=power - 1;
+           }
+         OutputBits(output,abs(temp),fNumShapeBits);
        }
-      
     }
   
   fclose(input);
   CloseOutputBitFile(output);
+
+  cout<<endl<<"There was following number of overflows: "<<endl
+      <<"Pad "<<pado<<endl
+      <<"Time "<<timeo<<endl
+      <<"Charge "<<chargeo<<endl
+      <<"Shape "<<shapeo<<endl;
 }
 
 void AliL3Compress::ExpandFile(Char_t *infile,Char_t *outfile)
@@ -132,17 +242,22 @@ void AliL3Compress::ExpandFile(Char_t *infile,Char_t *outfile)
 
   AliL3TrackModel trackmodel;
   AliL3ClusterModel *clusters=0;
+  Int_t count=0;
   
+  clusters = new AliL3ClusterModel[(NumRows[fPatch])];
   while(!feof(input->file))
     {
+      input->mask=0x80;//make sure we read a new byte from file.
       
+      //Read and write track:
       if(fread(&trackmodel,sizeof(AliL3TrackModel),1,input->file)!=1) break;
       fwrite(&trackmodel,sizeof(AliL3TrackModel),1,output);
-      input->mask=0x80;//make sure we read a new byte from file.
-      clusters = new AliL3ClusterModel[(trackmodel.fNClusters)];
-      for(Int_t i=0; i<trackmodel.fNClusters; i++)
+      
+      for(Int_t i=0; i<NumRows[fPatch]; i++)
        {
          Int_t temp,sign;
+         
+         //Read empty flag:
          temp = InputBit(input);
          if(!temp) 
            {
@@ -150,22 +265,312 @@ void AliL3Compress::ExpandFile(Char_t *infile,Char_t *outfile)
              continue;
            }
          clusters[i].fPresent=kTRUE;
+         
+         //Read time information:
          sign=InputBit(input);
-         temp = InputBits(input,8);
+         temp = InputBits(input,fNumTimeBits);
          if(!sign)
            temp*=-1;
          clusters[i].fDTime = temp;
+         
+         //Read pad information:
          sign=InputBit(input);
-         temp = InputBits(input,8);
+         temp = InputBits(input,fNumPadBits);
          if(!sign)
            temp*=-1;
          clusters[i].fDPad = temp;
-         temp=InputBits(input,10);
+         
+         //Read charge information:
+         sign = InputBit(input);
+         temp=InputBits(input,fNumChargeBits);
+         if(!sign)
+           temp*=-1;
          clusters[i].fDCharge = temp;
+         
+         //Read shape information:
+         temp = InputBits(input,fNumShapeBits);
+         clusters[i].fDSigmaY2 = temp;
+         
+         temp = InputBits(input,fNumShapeBits);
+         clusters[i].fDSigmaZ2 = temp;
        }
+      
+
+      count++;
       fwrite(clusters,(trackmodel.fNClusters)*sizeof(AliL3ClusterModel),1,output);
-      delete [] clusters;
+      
     }
+  
+  delete [] clusters;
   fclose(output);
   CloseInputBitFile(input);
 }
+
+void AliL3Compress::CreateDigitArray(Int_t maxnumber)
+{
+  fNUsed=0;
+  fNDigits = 0;
+  fMaxDigits=maxnumber;
+  if(fDigits) delete [] fDigits;
+  fDigits = new AliL3RandomDigitData[maxnumber];
+  if(fDPt) delete [] fDPt;
+  fDPt = new AliL3RandomDigitData*[maxnumber];
+}
+
+void AliL3Compress::RestoreData(Char_t *uncompfile)
+{
+  
+  //Read the uncompressed file:
+  ReadFile(uncompfile);
+  
+  CreateDigitArray(100000);
+  
+  Float_t pad,time,sigmaY2,sigmaZ2;
+  Int_t charge;
+  for(Int_t j=NRows[fPatch][0]; j<=NRows[fPatch][1]; j++)
+    {
+      for(Int_t i=0; i<fTracks->GetNTracks(); i++)
+       {
+         AliL3ModelTrack *track = (AliL3ModelTrack*)fTracks->GetCheckedTrack(i);
+         if(!track) continue;
+         if(!track->GetPad(j,pad) || 
+            !track->GetTime(j,time) || 
+            !track->GetClusterCharge(j,charge) ||
+            !track->GetXYWidth(j,sigmaY2) || 
+            !track->GetZWidth(j,sigmaZ2))
+           continue;
+         
+         CreateDigits(j,pad,time,charge,sigmaY2,sigmaZ2);
+       }
+    }
+  
+  QSort(fDPt,0,fNDigits);
+}
+
+void AliL3Compress::PrintDigits()
+{
+  Int_t pad,time,charge,row;
+  for(Int_t i=0; i<fNDigits; i++)
+    {
+      row = fDPt[i]->fRow;
+      pad = fDPt[i]->fPad;
+      time = fDPt[i]->fTime;
+      charge = fDPt[i]->fCharge;
+      if(i>0 && row != fDPt[i-1]->fRow)
+       cout<<"---Padrow "<<row<<"---"<<endl;
+      cout<<"Pad "<<pad<<" time "<<time<<" charge "<<charge<<endl;
+    }
+}
+
+void AliL3Compress::WriteRestoredData(Char_t *remainfile,Char_t *restoredfile)
+{
+  //Get the remaining raw data array:
+  AliL3MemHandler *mem = new AliL3MemHandler();
+  mem->SetBinaryInput(remainfile);
+  UInt_t numdigits;
+  AliL3DigitRowData *origRow = mem->CompBinary2Memory(numdigits);
+  mem->CloseBinaryInput();
+  
+  //Allocate memory for the merged data:
+  UInt_t size = mem->GetAllocatedSize() + fNDigits*sizeof(AliL3DigitData);
+  Byte_t *data = new Byte_t[size];
+  memset(data,0,size);
+  AliL3DigitRowData *tempRow = (AliL3DigitRowData*)data;
+
+  Int_t ndigits,action,charge;
+  UShort_t pad,time;
+  for(Int_t i=NRows[fPatch][0]; i<=NRows[fPatch][1]; i++)
+    {
+      tempRow->fRow = i;
+      ndigits=0;
+      AliL3DigitData *origDig = origRow->fDigitData;
+      AliL3DigitData *tempDig = tempRow->fDigitData;
+      
+      while(1)
+       {
+         for(UInt_t j=0; j<origRow->fNDigit; j++)
+           {
+             pad = origDig[j].fPad;
+             time = origDig[j].fTime;
+             charge = origDig[j].fCharge;
+             while((action=ComparePoints(i,pad,time)) == 1)
+               {
+                 tempDig[ndigits].fPad = fDPt[fNUsed]->fPad;
+                 tempDig[ndigits].fTime = fDPt[fNUsed]->fTime;
+                 tempDig[ndigits].fCharge = fDPt[fNUsed]->fCharge;
+                 ndigits++;
+                 fNUsed++;
+               }
+             if(action == 0)
+               {
+                 tempDig[ndigits].fPad = pad;
+                 tempDig[ndigits].fTime = time;
+                 tempDig[ndigits].fCharge = charge;
+                 ndigits++;
+               }
+           }
+         if(fNUsed >= fNDigits) break;
+         if(fDPt[fNUsed]->fRow != i) //we are on a new row
+           break;
+         tempDig[ndigits].fPad = fDPt[fNUsed]->fPad;
+         tempDig[ndigits].fTime = fDPt[fNUsed]->fTime;
+         tempDig[ndigits].fCharge = fDPt[fNUsed]->fCharge;
+         ndigits++;
+         fNUsed++;
+       }
+      tempRow->fNDigit = ndigits;
+      Int_t size = sizeof(AliL3DigitData)*tempRow->fNDigit + sizeof(AliL3DigitRowData);
+      Byte_t *byte_pt = (Byte_t*)tempRow;
+      byte_pt += size;
+      tempRow = (AliL3DigitRowData*)byte_pt;
+      mem->UpdateRowPointer(origRow);
+    }
+  
+  mem->Free();
+  mem->SetBinaryOutput(restoredfile);
+  mem->Memory2CompBinary((UInt_t)NumRows[fPatch],(AliL3DigitRowData*)data);
+  mem->CloseBinaryOutput();
+  
+  delete [] data;
+  delete mem;
+}
+
+void AliL3Compress::CreateDigits(Int_t row,Float_t pad,Float_t time,Int_t charge,Float_t sigmaY2,Float_t sigmaZ2)
+{
+  //Create raw data out of the cluster.
+  
+  AliL3Transform *tr = new AliL3Transform();
+  TRandom *random = new TRandom();
+  
+  Int_t entries=10000;
+  TH1F *hist1 = new TH1F("hist1","",tr->GetNPads(row),0,tr->GetNPads(row)-1);
+  TH1F *hist2 = new TH1F("hist2","",tr->GetNTimeBins(),0,tr->GetNTimeBins()-1);
+  TH2F *hist3 = new TH2F("hist3","",tr->GetNPads(row),0,tr->GetNPads(row)-1,tr->GetNTimeBins(),0,tr->GetNTimeBins()-1);
+  
+  //Convert back the sigmas:
+  Float_t padw,timew;
+  if(fPatch < 3)
+    padw = tr->GetPadPitchWidthLow();
+  else
+    padw = tr->GetPadPitchWidthUp();
+  timew = tr->GetZWidth();
+
+  if(fPatch < 3)
+    sigmaY2 = sigmaY2/2.07;
+  sigmaY2 = sigmaY2/0.108;
+  sigmaY2 = sigmaY2/(padw*padw);
+  sigmaY2 = sigmaY2 - 1./12;
+  
+  if(fPatch < 3)
+    sigmaZ2 = sigmaZ2/1.77;
+  sigmaZ2 = sigmaZ2/0.169;
+  sigmaZ2 = sigmaZ2/(timew*timew);
+  sigmaZ2 = sigmaZ2 - 1./12;
+  
+  if(sigmaY2 <= 0 || sigmaZ2 <= 0)
+    {
+      cerr<<"AliL3Compress::CreateDigits() : Wrong sigmas : "<<sigmaY2<<" "<<sigmaZ2<<endl;
+      return;
+    }
+  
+  //Create the distributions in pad and time:
+  for(Int_t i=0; i<entries; i++)
+    {
+      hist1->Fill(random->Gaus(pad,sqrt(sigmaY2)));
+      hist2->Fill(random->Gaus(time,sqrt(sigmaZ2)));
+    }
+    
+  //Create the cluster:
+  Int_t bin1,bin2;
+  Double_t content1,content2,dpad,dtime;
+  for(Int_t i=0; i<hist1->GetEntries(); i++)
+    {
+      bin1 = hist1->GetBin(i);
+      content1 = hist1->GetBinContent(bin1);
+      if((Int_t)content1==0) continue;
+      content1 = charge*content1/entries;
+      dpad = hist1->GetBinCenter(bin1);
+      for(Int_t j=0; j<hist2->GetEntries(); j++)
+       {
+         bin2 = hist2->GetBin(j);
+         content2 = hist2->GetBinContent(bin2);
+         if((Int_t)content2==0) continue;
+         content2 = content1*content2/entries;
+         dtime = hist2->GetBinCenter(bin2);
+         hist3->Fill(dpad,dtime,content2);
+       }
+    }
+  
+  //Fill it into the digit array:
+  for(Int_t i=0; i<hist3->GetNbinsX(); i++)
+    {
+      for(Int_t j=0; j<hist3->GetNbinsY(); j++)
+       {
+         bin1 = hist3->GetBin(i,j);
+         content1 = hist3->GetBinContent(bin1);
+         if((Int_t)content1 < 3) continue;
+         if(content1 >= 1024)
+           content1 = 1023;
+         if(fNDigits >= fMaxDigits)
+           {
+             cerr<<"AliL3Compress::CreateDigits() : Array index out of range : "<<fNDigits<<endl;
+             return;
+           }
+         fDigits[fNDigits].fCharge=(Int_t)content1;
+         fDigits[fNDigits].fRow = row;
+         fDigits[fNDigits].fPad = (Int_t)hist3->GetXaxis()->GetBinCenter(i);
+         fDigits[fNDigits].fTime = (Int_t)hist3->GetYaxis()->GetBinCenter(j);
+         fDPt[fNDigits] = &fDigits[fNDigits];
+         fNDigits++;
+       }
+    }
+  
+  delete random;
+  delete hist1;
+  delete hist2;
+  delete hist3;
+  delete tr;
+}
+
+void AliL3Compress::QSort(AliL3RandomDigitData **a, Int_t first, Int_t last)
+{
+  
+  // Sort array of AliL3RandomDigitData pointers using a quicksort algorithm.
+  // Uses CompareDigits() to compare objects.
+  // Thanks to Root!
+  
+  static AliL3RandomDigitData *tmp;
+  static int i;           // "static" to save stack space
+  int j;
+  
+  while (last - first > 1) {
+    i = first;
+    j = last;
+    for (;;) {
+      while (++i < last && CompareDigits(a[i], a[first]) < 0)
+       ;
+      while (--j > first && CompareDigits(a[j], a[first]) > 0)
+       ;
+      if (i >= j)
+       break;
+      
+      tmp  = a[i];
+      a[i] = a[j];
+      a[j] = tmp;
+    }
+    if (j == first) {
+      ++first;
+      continue;
+    }
+    tmp = a[first];
+    a[first] = a[j];
+    a[j] = tmp;
+    if (j - first < last - (j + 1)) {
+      QSort(a, first, j);
+      first = j + 1;   // QSort(j + 1, last);
+    } else {
+      QSort(a, j + 1, last);
+      last = j;        // QSort(first, j);
+    }
+  }
+}
index ddca173..e5a6e1d 100644 (file)
@@ -3,6 +3,7 @@
 
 #include "AliL3RootTypes.h"
 #include "AliL3Models.h"
+#include "AliL3DigitData.h"
 
 class AliL3TrackArray;
 
@@ -10,21 +11,67 @@ class AliL3Compress {
   
  private:
   AliL3TrackArray *fTracks; //!
+  AliL3RandomDigitData *fDigits; //!
+  AliL3RandomDigitData **fDPt; //!
+  Int_t fNDigits;
+  Int_t fNUsed;
+  Int_t fMaxDigits;
   
+  Int_t fNumPadBits;
+  Int_t fNumTimeBits;
+  Int_t fNumChargeBits;
+  Int_t fNumShapeBits;
+  Int_t fSlice;
+  Int_t fPatch;
   
+  void CreateDigitArray(Int_t maxnumber);
+  void CreateDigits(Int_t row,Float_t pad,Float_t time,Int_t charge,Float_t ywidth,Float_t zwidth);
+  void QSort(AliL3RandomDigitData **a, Int_t first, Int_t last);
+  Int_t ComparePoints(Int_t row,UShort_t pad,UShort_t time);
+  Int_t CompareDigits(AliL3RandomDigitData *a,AliL3RandomDigitData *b);
+
  public:
   AliL3Compress();
+  AliL3Compress(Int_t slice,Int_t patch,Int_t pad=7,Int_t time=7,Int_t charge=10,Int_t shape=4);
   virtual ~AliL3Compress();
   
+  void SetBitNumbers(Int_t pad,Int_t time,Int_t charge,Int_t shape);
   void WriteFile(AliL3TrackArray *tracks,Char_t *filename);
   void ReadFile(Char_t *filename);
   void CompressFile(Char_t *infile,Char_t *outfile);
   void ExpandFile(Char_t *infile,Char_t *outfile);
-  
+  void RestoreData(Char_t *uncompfile);
+  void WriteRestoredData(Char_t *remainfile,Char_t *restoredfile);
+  void PrintDigits();
+
   AliL3TrackArray *GetTracks() {return fTracks;}
   
   ClassDef(AliL3Compress,1) 
 
 };
 
+inline Int_t  AliL3Compress::ComparePoints(Int_t row,UShort_t pad,UShort_t time)
+{
+  
+  if(fDPt[fNUsed]->fRow != row) return 0;
+  
+  if(fDPt[fNUsed]->fPad < pad) return 1;
+  if(fDPt[fNUsed]->fPad == pad && fDPt[fNUsed]->fTime < time) return 1;
+  
+  return 0;
+
+}
+
+inline Int_t AliL3Compress::CompareDigits(AliL3RandomDigitData *a,AliL3RandomDigitData *b)
+{
+  if(a->fRow < b->fRow) return -1;
+    
+  if(a->fPad==b->fPad && a->fTime == b->fTime && a->fRow == b->fRow) return 0;
+  
+  if(a->fPad<b->fPad && a->fRow == b->fRow) return -1;
+  if(a->fPad==b->fPad && a->fTime<b->fTime && a->fRow == b->fRow) return -1;
+  
+  return 1;
+}
+
 #endif
index 6a47461..19a4ff1 100644 (file)
@@ -8,7 +8,6 @@
 #include <math.h>
 
 #include "AliL3ModelTrack.h"
-#include "AliL3Defs.h"
 #include "AliL3Transform.h"
 
 ClassImp(AliL3ModelTrack)
@@ -17,7 +16,7 @@ AliL3ModelTrack::AliL3ModelTrack()
 {
   fNClusters = 0;
   fClusters = 0;
-  fOverlap = -1;
+  fOverlap = 0;
   fPad=0;
   fTime=0;
   fClusterCharge=0;
@@ -33,6 +32,8 @@ AliL3ModelTrack::~AliL3ModelTrack()
     delete [] fPad;
   if(fTime)
     delete [] fTime;
+  if(fOverlap)
+    delete [] fOverlap;
   if(fTrackModel)
     delete fTrackModel;
 }
@@ -40,75 +41,194 @@ AliL3ModelTrack::~AliL3ModelTrack()
 void AliL3ModelTrack::Init(Int_t slice,Int_t patch)
 {
   fNClusters = 0;
+  fSlice=slice;
+  fPatch=patch;
   Int_t nrows = NumRows[patch];
   fClusters = new AliL3ClusterModel[nrows];
-  memset((void*)fClusters,0,nrows*sizeof(AliL3ClusterModel));
-  
-  fPad = new Float_t[NRowsSlice];
-  fTime = new Float_t[NRowsSlice];
+  fPad = new Float_t[nrows];
+  fTime = new Float_t[nrows];
   fTrackModel = new AliL3TrackModel;
-  memset(fTrackModel,0,sizeof(AliL3TrackModel));
+  fOverlap = new Int_t[nrows];
   
-  fClusterCharge = 100;
+  memset(fClusters,0,nrows*sizeof(AliL3ClusterModel));
+  memset(fPad,0,nrows*sizeof(Float_t));
+  memset(fTime,0,nrows*sizeof(Float_t));
+  memset(fTrackModel,0,sizeof(AliL3TrackModel));
+  for(Int_t i=0; i<nrows; i++)
+    fOverlap[i]=-1;
+
+  fClusterCharge = 260;
   AliL3Transform transform;
   
-  fXYResidualQ = 0.01/transform.GetPadPitchWidth(patch);
-  fZResidualQ = 0.01/transform.GetPadPitchWidth(patch);
+  fXYResidualQ = 0.1/transform.GetPadPitchWidth(patch);
+  fZResidualQ = 0.1/transform.GetPadPitchWidth(patch);
+  
+
+  fXYWidthQ = 0.01;
+  fZWidthQ = 0.01;
 }
 
 
-void AliL3ModelTrack::SetCluster(Float_t fpad,Float_t ftime,Float_t charge,Float_t sigmaY2,Float_t sigmaZ2)
+void AliL3ModelTrack::SetCluster(Int_t row,Float_t fpad,Float_t ftime,Float_t charge,Float_t sigmaY2,Float_t sigmaZ2)
 {
-
-  AliL3ClusterModel *cl = &fClusters[fNClusters];
+  Int_t index = row - NRows[fPatch][0];
+  if(index != fNClusters)
+    cout<<"AliL3ModelTrack::SetCluster() : Mismatch ; index: "<<index<<" nclusters "<<fNClusters<<endl;
+  
+  if(index < 0 || index > NumRows[fPatch])
+    {
+      cerr<<"AliL3ModelTrack::SetCluster() : Wrong index: "<<index<<" row "<<row<<endl;
+      return;
+    }
+  AliL3ClusterModel *cl = GetClusterModel(row);
   if(!charge)
     cl->fPresent = kFALSE;
   else
     {
       cl->fPresent = kTRUE;
-      cl->fDTime = (ftime - GetTimeHit(fNClusters))/fXYResidualQ;
-      cl->fDPad = (fpad - GetPadHit(fNClusters))/fZResidualQ;
-      cl->fDCharge = charge;
-      cl->fDSigmaY2 = sigmaY2;
-      cl->fDSigmaZ2 = sigmaZ2;
+      cl->fDTime = (ftime - GetTimeHit(row))/fXYResidualQ;
+      cl->fDPad = (fpad - GetPadHit(row))/fZResidualQ;
+      cl->fDCharge = charge - fClusterCharge;
+      cl->fDSigmaY2 = sigmaY2/fXYWidthQ;
+      cl->fDSigmaZ2 = sigmaZ2/fZWidthQ;
     }
-  //cout<<"Pad "<<fpad<<" dtime "<<ftime<<" charge "<<charge<<" sigmaY2 "<<sigmaY2<<" sigmaZ2 "<<sigmaZ2<<endl;
+  
   fNClusters++;
 }
 
 
-
 void AliL3ModelTrack::FillModel()
 {
-  //fTrackModel = new AliL3TrackModel;
+  //Fill the track structure
+  
+  if(!fTrackModel)
+    {
+      cerr<<"AliL3ModelTrack::FillModel() : No trackmodel "<<endl;
+      return;
+    }
   fTrackModel->fKappa = GetKappa();
   fTrackModel->fFirstPointX = GetFirstPointX();
   fTrackModel->fFirstPointY = GetFirstPointY();
   fTrackModel->fFirstPointZ = GetFirstPointZ();
   fTrackModel->fTgl = GetTgl();
   fTrackModel->fPsi = GetPsi();
-  fTrackModel->fLength = GetLength();
+  fTrackModel->fLength = (Short_t)GetLength();
   fTrackModel->fClusterCharge = fClusterCharge;
   fTrackModel->fNClusters = fNClusters;
 
 }
 
-void AliL3ModelTrack::Print()
+void AliL3ModelTrack::FillTrack()
 {
-  //Print info
+  //Fill the track parameters from the structure.
+  
+  if(!fTrackModel)
+    {
+      cerr<<"AliL3ModelTrack::FillTrack() : No data!!"<<endl;
+      return;
+    }
+  SetKappa(fTrackModel->fKappa);
+  SetCharge((-1*(Int_t)copysign(1.,GetKappa())));
+  SetFirstPoint(fTrackModel->fFirstPointX,fTrackModel->fFirstPointY,fTrackModel->fFirstPointZ);
+  SetTgl(fTrackModel->fTgl);
+  SetPsi(fTrackModel->fPsi);
+  SetLength(fTrackModel->fLength);
+  fClusterCharge=fTrackModel->fClusterCharge;
+  fNClusters = fTrackModel->fNClusters;
+  SetPt((BFACT*BField)/fabs(GetKappa()));
+  
+  CalculateHelix();
+  
+  AliL3Transform transform;
+  Float_t hit[3];
+  Int_t sector,row;
+  for(Int_t i=NRows[fPatch][0]; i<=NRows[fPatch][1]; i++)
+    {
+      AliL3ClusterModel *cl = GetClusterModel(i);
+      if(!cl) continue;
+      GetCrossingPoint(i,hit);
+      transform.Slice2Sector(fSlice,i,sector,row);
+      transform.Local2Raw(hit,sector,row);
+      SetPadHit(i,hit[1]);
+      SetTimeHit(i,hit[2]);
+    }
+}
 
-  cout<<"---------------------"<<endl;
-  for(Int_t i=0; i<fNClusters; i++)
+void AliL3ModelTrack::SetPadHit(Int_t row,Float_t pad)
+{
+  Int_t index = row-NRows[fPatch][0];
+  if(index < 0 || index > NumRows[fPatch])
     {
-      AliL3ClusterModel *cl = &fClusters[i];
-      if(!cl->fPresent)
-       cout<<i<<" Empty"<<endl;
-      else
-       {
-         cout<<i<<" Dpad "<<cl->fDPad<<" Dtime "<<cl->fDTime<<" Dcharge "<<cl->fDCharge;
-         cout<<" Padcrossing "<<GetPadHit(i)<<" Timecrossing "<<GetTimeHit(i)<<endl;
-       }
+      cerr<<"AliL3ModelTrack::SetPadHit() : Wrong index: "<<index<<endl;
+      return;
     }
+  fPad[index]=pad;
+  
+}
+
+void AliL3ModelTrack::SetTimeHit(Int_t row,Float_t time)
+{
+  Int_t index = row-NRows[fPatch][0];
+  if(index < 0 || index > NumRows[fPatch])
+    {
+      cerr<<"AliL3ModelTrack::SetTimeHit() : Wrong index: "<<index<<endl;
+      return;
+    }
+  fTime[index]=time;
+}
+
+void AliL3ModelTrack::SetOverlap(Int_t row,Int_t id)
+{
+  Int_t index = row-NRows[fPatch][0];
+  if(index < 0 || index > NumRows[fPatch])
+    {
+      cerr<<"AliL3ModelTrack::SetOverlap() : Wrong index: "<<index<<endl;
+      return;
+    }
+  fOverlap[index]=id;
+}
+
+Bool_t AliL3ModelTrack::GetPad(Int_t row,Float_t &pad)
+{
+  //(ftime - GetTimeHit(fNClusters))/fXYResidualQ;
+  //(fpad - GetPadHit(fNClusters))/fZResidualQ;
+
+  AliL3ClusterModel *cl = GetClusterModel(row);
+  pad = cl->fDPad*fXYResidualQ + GetPadHit(row);
+
+  return (Bool_t)cl->fPresent;
+}
+
+Bool_t AliL3ModelTrack::GetTime(Int_t row,Float_t &time)
+{
+  AliL3ClusterModel *cl = GetClusterModel(row);
+  time = cl->fDTime*fZResidualQ + GetTimeHit(row);
+
+  return (Bool_t)cl->fPresent;
+}
+
+Bool_t AliL3ModelTrack::GetClusterCharge(Int_t row,Int_t &charge)
+{
+  AliL3ClusterModel *cl = GetClusterModel(row);
+  charge = (Int_t)cl->fDCharge + fClusterCharge;
+  
+  return (Bool_t)cl->fPresent;
+}
+
+Bool_t AliL3ModelTrack::GetXYWidth(Int_t row,Float_t &width)
+{
+  AliL3ClusterModel *cl = GetClusterModel(row);
+  width = cl->fDSigmaY2*fXYWidthQ;
+  
+  return (Bool_t)cl->fPresent;
+}
+
+Bool_t AliL3ModelTrack::GetZWidth(Int_t row,Float_t &width)
+{
+  AliL3ClusterModel *cl = GetClusterModel(row);
+  width = cl->fDSigmaZ2*fZWidthQ;
+  
+  return (Bool_t)cl->fPresent;
 }
 
 Bool_t AliL3ModelTrack::GetPadResidual(Int_t row,Float_t &res)
@@ -123,6 +243,78 @@ Bool_t AliL3ModelTrack::GetTimeResidual(Int_t row,Float_t &res)
   return fClusters[row].fPresent;
 }
 
+Float_t AliL3ModelTrack::GetPadHit(Int_t row)
+{
+  Int_t index = row-NRows[fPatch][0];
+  if(index < 0 || index > NumRows[fPatch])
+    {
+      cerr<<"AliL3ModelTrack::GetPadHit() : Wrong index: "<<index<<" row "<<row<<endl;
+      return 0;
+    }
+  return fPad[index];
+}
+
+Float_t AliL3ModelTrack::GetTimeHit(Int_t row)
+{
+  Int_t index = row-NRows[fPatch][0];
+  if(index < 0 || index > NumRows[fPatch])
+    {
+      cerr<<"AliL3ModelTrack::GetTimeHit() : Wrong index: "<<index<<" row "<<row<<endl;
+      return 0;
+    }
+  return fTime[index];
+}
+
+Int_t AliL3ModelTrack::GetOverlap(Int_t row)
+{
+  Int_t index = row-NRows[fPatch][0];
+  if(index < 0 || index > NumRows[fPatch])
+    {
+      cerr<<"AliL3ModelTrack::GetOverlap() : Wrong index: "<<index<<endl;
+      return 0;
+    }
+  return fOverlap[index];
+}
+
+AliL3ClusterModel *AliL3ModelTrack::GetClusterModel(Int_t row)
+{
+  if(!fClusters) return 0; 
+  Int_t index = row-NRows[fPatch][0];
+  if(index < 0 || index > NumRows[fPatch])
+    {
+      cerr<<"AliL3ModelTrack::GetClusterModel() : Wrong index: "<<index<<endl;
+      return 0;
+    }
+  return &fClusters[index];
+}
+
+void AliL3ModelTrack::Print()
+{
+  //Print info
+
+  cout<<"---------------------"<<endl;
+  cout<<"First point "<<GetFirstPointX()<<" "<<GetFirstPointY()<<" "<<GetFirstPointZ()<<endl;
+  cout<<"Last point "<<GetLastPointX()<<" "<<GetLastPointY()<<" "<<GetLastPointZ()<<endl;
+  cout<<"Pt "<<GetPt()<<" kappa "<<GetKappa()<<" tgl "<<GetTgl()<<" psi "<<GetPsi()<<" charge "<<GetCharge()<<endl;
+  cout<<"Center "<<GetCenterX()<<" "<<GetCenterY()<<endl<<endl;
+  cout<<"NHits "<<GetNClusters()<<endl;
+  cout<<"Clusters:"<<endl;
+
+  for(Int_t i=NRows[fPatch][0]; i<=NRows[fPatch][1]; i++)
+    {
+      AliL3ClusterModel *cl = GetClusterModel(i);
+      
+      if(!cl->fPresent)
+       cout<<i<<" Empty"<<" Padcrossing "<<GetPadHit(i)<<" Timecrossing "<<GetTimeHit(i)<<" ";
+      else
+       {
+         cout<<i<<" Dpad "<<cl->fDPad<<" Dtime "<<cl->fDTime<<" Dcharge "<<cl->fDCharge;
+         cout<<" Padcrossing "<<GetPadHit(i)<<" Timecrossing "<<GetTimeHit(i)<<" ";
+       }
+      cout<<"Overlapping index "<<GetOverlap(i)<<endl;
+    }
+}
+
 //----------Code below taken from AliTPCTracker.cxx-----------------------
 //Functions that give the expected cluster errors based on track parameters.
 Double_t AliL3ModelTrack::GetParSigmaY2(Double_t r)//, Double_t tgl, Double_t pt)
index e6baf13..4ee05d8 100644 (file)
@@ -3,6 +3,7 @@
 
 #include "AliL3Track.h"
 #include "AliL3Models.h"
+#include "AliL3Defs.h"
 
 class AliL3ModelTrack : public AliL3Track {
 
@@ -12,9 +13,15 @@ class AliL3ModelTrack : public AliL3Track {
   AliL3ClusterModel *fClusters; //!
   AliL3TrackModel *fTrackModel; //!
   Short_t fNClusters;
-  Int_t fOverlap;
+  Int_t *fOverlap; //!
   Float_t fXYResidualQ; //Quantization steps.
   Float_t fZResidualQ;
+  Float_t fXYResolution;
+  Float_t fZResolution;
+  Float_t fXYWidthQ;
+  Float_t fZWidthQ;
+  Int_t fSlice;
+  Int_t fPatch;
   
   //Crossing points with padrows
   Float_t *fPad; //!
@@ -25,20 +32,28 @@ class AliL3ModelTrack : public AliL3Track {
   virtual ~AliL3ModelTrack();
   
   void Init(Int_t slice,Int_t patch);
-  void SetCluster(Float_t dpad,Float_t dtime,Float_t charge,Float_t sigmaY2,Float_t sigmaZ2);
+  void SetCluster(Int_t row,Float_t dpad,Float_t dtime,Float_t charge,Float_t sigmaY2,Float_t sigmaZ2);
   void FillModel();
+  void FillTrack();
   void Print();
   
-  void SetPadHit(Int_t row,Float_t f) {fPad[row]=f;}
-  void SetTimeHit(Int_t row,Float_t f) {fTime[row]=f;}
-  void SetOverlap(Int_t i) {fOverlap=i;}
+  void SetPadHit(Int_t row,Float_t f);
+  void SetTimeHit(Int_t row,Float_t f);
+  void SetOverlap(Int_t row,Int_t id);
+  void SetXYResolution(Float_t f) {fXYResolution=f;}
+  void SetZResolution(Float_t f) {fZResolution=f;}
   
   AliL3ClusterModel *GetClusters() {return fClusters;}
-  AliL3ClusterModel *GetClusterModel(Int_t i) {if(!fClusters) return 0; return &fClusters[i];}
   AliL3TrackModel *GetModel() {return fTrackModel;}
-  Int_t GetOverlap() {return fOverlap;}
-  Float_t GetPadHit(Int_t row) {return fPad[row];}
-  Float_t GetTimeHit(Int_t row) {return fTime[row];}
+  AliL3ClusterModel *GetClusterModel(Int_t row);
+  Int_t GetOverlap(Int_t row);
+  Float_t GetPadHit(Int_t row);
+  Float_t GetTimeHit(Int_t row);
+  Bool_t GetPad(Int_t row,Float_t &pad);
+  Bool_t GetTime(Int_t row,Float_t &time);
+  Bool_t GetClusterCharge(Int_t row,Int_t &charge);
+  Bool_t GetXYWidth(Int_t row,Float_t &width);
+  Bool_t GetZWidth(Int_t row,Float_t &width);
   Bool_t GetPadResidual(Int_t row,Float_t &res);
   Bool_t GetTimeResidual(Int_t row,Float_t &res);
   Int_t GetNClusters() {return fNClusters;}
index d13a740..a5d2f16 100644 (file)
@@ -1,6 +1,6 @@
 //$Id$
 
-// Author: Anders Vestbo <mailto:vestbo$fi.uib.no>
+// Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
 //*-- Copyright &copy ASV
 
 #include <stream.h>
@@ -43,6 +43,8 @@ void AliL3Modeller::Init(Int_t slice,Int_t patch,Char_t *path)
   fPatch = patch;
   fPadOverlap=4;
   fTimeOverlap=4;
+  fTrackThreshold=10;
+  
   fTransform = new AliL3Transform();
   fTracks = new AliL3TrackArray("AliL3ModelTrack");
 
@@ -83,7 +85,7 @@ void AliL3Modeller::Init(Int_t slice,Int_t patch,Char_t *path)
   SetInputData(digits);
 }
 
-void AliL3Modeller::Process()
+void AliL3Modeller::FindClusters()
 {
   
   if(!fTracks)
@@ -102,10 +104,11 @@ void AliL3Modeller::Process()
 
   Int_t ntimes = fTransform->GetNTimeBins()+1;
   Int_t npads = fTransform->GetNPads(NRows[fPatch][1])+1;//Max num of pads.
-  Digit *row = new Digit[(ntimes)*(npads)];
+  Int_t bounds = ntimes*npads;
+  Digit *row = new Digit[bounds];
   
   Int_t seq_charge;
-  Int_t pad,time;
+  Int_t pad,time,index;
   Short_t charge;
   Cluster cluster;
 
@@ -120,22 +123,24 @@ void AliL3Modeller::Process()
          charge = digPt[j].fCharge;
          row[ntimes*pad+time].fCharge = charge;
          row[ntimes*pad+time].fUsed = kFALSE;
+         //cout<<"Row "<<i<<" pad "<<pad<<" time "<<time<<" charge "<<charge<<endl;
        }
       
       for(Int_t k=0; k<fTracks->GetNTracks(); k++)
        {
          AliL3ModelTrack *track = (AliL3ModelTrack*)fTracks->GetCheckedTrack(k);
          if(!track) continue;
-         if(track->GetOverlap()>=0) continue;//Track is overlapping
-         if(track->GetPadHit(i)<0 || track->GetTimeHit(i)<0)
+         //if(track->GetOverlap(i)>=0) continue;//Track is overlapping
+         
+         if(track->GetPadHit(i)<0 || track->GetTimeHit(i)<0 || track->GetOverlap(i)>=0)
            {
-             track->SetCluster(0,0,0,0,0); //The track has left the patch.
+             track->SetCluster(i,0,0,0,0,0); //The track has left the patch.
              continue;
            }
          
          Int_t hitpad = (Int_t)rint(track->GetPadHit(i));
          Int_t hittime = (Int_t)rint(track->GetTimeHit(i));
-         //cout<<"Checking track with pad "<<hitpad<<" time "<<hittime<<endl;
+         //cout<<"Checking track on row "<<i<<" with pad "<<hitpad<<" time "<<hittime<<endl;
          pad = hitpad;
          time = hittime;
          Int_t padsign=-1;
@@ -145,16 +150,51 @@ void AliL3Modeller::Process()
          
          while(1)//Process this padrow
            {
+             if(pad < 0 || pad >= fTransform->GetNPads(i)) 
+               {
+                 //cout<<"Pad<0 on row "<<i<<endl;
+                 FillCluster(track,&cluster,i);
+                 break;
+               }
              seq_charge=0;
              timesign=-1;
              time = hittime;
              while(1) //Process sequence on this pad:
                {
-                 charge = row[ntimes*pad+time].fCharge;
+                 if(time < 0) break;
+                 index = ntimes*pad + time;
+                 if(index < 0 || index >= bounds)
+                   {
+                     cerr<<"AliL3Modeller::FindClusters : Index out of range : "<<index
+                       <<" on row "<<i<<" pad "<<pad<<" time "<<time<<endl;
+                     break;
+                   }
+                 charge = row[index].fCharge;
                  if(charge==0 && timesign==-1)
-                   {time=hittime+1; timesign=1; continue;}
+                   {
+                     if(seq_charge==0 && abs(time-hittime) <= fTimeOverlap/2)
+                       {
+                         time--;
+                         continue;
+                       }
+                     else
+                       {
+                         time = hittime+1;
+                         timesign=1;
+                         continue;
+                       }
+                   }
                  else if(charge==0 && timesign==1)
-                   break;
+                   {
+                     if(seq_charge==0 && abs(time-hittime) <= fTimeOverlap/2)
+                       {
+                         time++;
+                         continue;
+                       }
+                     else
+                       break;
+                   }
+                 
                  //cout<<"Doing pad "<<pad<<" time "<<time<<" charge "<<charge<<endl;
                  
                  seq_charge += charge;
@@ -175,7 +215,7 @@ void AliL3Modeller::Process()
                {
                  if(padsign==-1) 
                    {
-                     if(cluster.fCharge==0 && abs(pad-hitpad) < fPadOverlap/2)
+                     if(cluster.fCharge==0 && abs(pad-hitpad) <= fPadOverlap/2)
                        {
                          pad--; //In this case, we haven't found anything yet, 
                        }        //so we will try to expand our search within the natural boundaries.
@@ -186,30 +226,63 @@ void AliL3Modeller::Process()
                        }
                      continue;
                    }
-                 else if(padsign==1 && cluster.fCharge==0 && abs(pad-hitpad) < fPadOverlap/2)
+                 
+                 else if(padsign==1)
                    {
-                     pad++;
-                     continue;
+                     if(cluster.fCharge==0 && abs(pad-hitpad) <= fPadOverlap/2)
+                       {
+                         pad++;     //In this case, we haven't found anything yet, 
+                         continue;  //so we will try to expand our search within the natural boundaries.
+                       }
+                     else //We are out of range, or cluster if finished.
+                       {
+                         //cout<<"Outof range; charge "<<cluster.fCharge<<" paddiff "<<abs(pad-hitpad)<<endl;
+                         FillCluster(track,&cluster,i);
+                         break;
+                       }
                    }
                  else //Nothing more in this cluster
                    {
-                     Float_t fcharge = (Float_t)cluster.fCharge;
-                     Float_t fpad = ((Float_t)cluster.fPad/fcharge);
-                     Float_t ftime = ((Float_t)cluster.fTime/fcharge);
-                     Float_t sigmaY2,sigmaZ2;
-                     CalcClusterWidth(&cluster,sigmaY2,sigmaZ2);
-                     //cout<<"row "<<i<<" charge "<<fcharge<<endl;
-                     track->SetCluster(fpad,ftime,fcharge,sigmaY2,sigmaZ2);
+                     //cout<<"Filling final cluster"<<endl;
+                     FillCluster(track,&cluster,i);
                      break;
                    } 
                }
-             // pad += padsign;
            }
+         //cout<<"done"<<endl;
        }
+      
       FillZeros(rowPt,row);
       fMemHandler->UpdateRowPointer(rowPt);
     }
   delete [] row;
+  cout<<"done processing"<<endl;
+  
+  
+  //Debug:
+  for(Int_t i=0; i<fTracks->GetNTracks(); i++)
+    {
+      AliL3ModelTrack *track = (AliL3ModelTrack*)fTracks->GetCheckedTrack(i);
+      if(!track) continue;
+      if(track->GetNClusters() != NumRows[fPatch])
+       cerr<<endl<<"Mismatching hitcounts; nclusters: "<<track->GetNClusters()<<" nrows "<<NumRows[fPatch]<<endl<<endl;
+    }
+  
+}
+
+void AliL3Modeller::FillCluster(AliL3ModelTrack *track,Cluster *cluster,Int_t row)
+{
+  if(cluster->fCharge==0)
+    {
+      track->SetCluster(row,0,0,0,0,0);
+      return;
+    }
+  Float_t fcharge = (Float_t)cluster->fCharge;
+  Float_t fpad = ((Float_t)cluster->fPad/fcharge);
+  Float_t ftime = ((Float_t)cluster->fTime/fcharge);
+  Float_t sigmaY2,sigmaZ2;
+  CalcClusterWidth(cluster,sigmaY2,sigmaZ2);
+  track->SetCluster(row,fpad,ftime,fcharge,sigmaY2,sigmaZ2);
   
 }
 
@@ -237,7 +310,7 @@ void AliL3Modeller::WriteRemaining(Char_t *output)
   rowPt = (AliL3DigitRowData*)fRowData;
   Int_t digitcount=0;
   Int_t ndigits[(NumRows[fPatch])];
-  for(Int_t i=NRows[fPatch][0]; i<NRows[fPatch][1]; i++)
+  for(Int_t i=NRows[fPatch][0]; i<=NRows[fPatch][1]; i++)
     {
       AliL3DigitData *digPt = (AliL3DigitData*)rowPt->fDigitData;
       ndigits[(i-NRows[fPatch][0])]=0;
@@ -257,7 +330,7 @@ void AliL3Modeller::WriteRemaining(Char_t *output)
   AliL3DigitRowData *tempPt = (AliL3DigitRowData*)data;
   rowPt = (AliL3DigitRowData*)fRowData;
   
-  for(Int_t i=NRows[fPatch][0]; i<NRows[fPatch][1]; i++)
+  for(Int_t i=NRows[fPatch][0]; i<=NRows[fPatch][1]; i++)
     {
       Int_t localcount=0;
       tempPt->fRow = i;
@@ -293,6 +366,7 @@ void AliL3Modeller::WriteRemaining(Char_t *output)
   mem->Memory2CompBinary((UInt_t)NumRows[fPatch],(AliL3DigitRowData*)data);
   mem->CloseBinaryOutput();
   delete mem;
+  delete [] data;
 }
 
 
@@ -305,22 +379,45 @@ void AliL3Modeller::CalculateCrossingPoints()
       return;
     }
   Float_t hit[3];
+  
+  
+  //Remove tracks which are no good:
+  for(Int_t i=0; i<fTracks->GetNTracks(); i++)
+    {
+      AliL3ModelTrack *track = (AliL3ModelTrack*)fTracks->GetCheckedTrack(i);
+      if(!track) continue;
+      if(track->GetFirstPointX() > fTransform->Row2X(NRows[fPatch][1]) || track->GetPt()<0.1)
+       fTracks->Remove(i);
+    }
+  
+  Int_t sector,row;
   for(Int_t i=NRows[fPatch][1]; i>=NRows[fPatch][0]; i--)
     {
       for(Int_t j=0; j<fTracks->GetNTracks(); j++)
        {
          AliL3ModelTrack *track = (AliL3ModelTrack*)fTracks->GetCheckedTrack(j);
          if(!track) continue;
-         if(track->GetNHits() < 100)
-           fTracks->Remove(j);
+                 
+         if(track->GetNHits() < fTrackThreshold)
+           {
+             fTracks->Remove(j);
+             continue;
+           }
+         
          if(!track->GetCrossingPoint(i,hit)) 
            {
-             cerr<<"AliL3Modeller::CalculateCrossingPoints : Track does not intersect line "<<endl;
+             cerr<<"AliL3Modeller::CalculateCrossingPoints : Track "<<j<<" does not intersect row "<<i<<" :"<<endl<<
+               "First point "<<track->GetFirstPointX()<<
+               " nhits "<<track->GetNHits()<<endl;//" tgl "<<track->GetTgl()<<" psi "<<track->GetPsi()<<" charge "<<track->GetCharge()<<endl;
+               //"Center "<<track->GetCenterX()<<" "<<track->GetCenterY()<<endl<<endl<<
+               //"--------"<<endl;
+             fTracks->Remove(j);
              continue;
            }
          //cout<<" x "<<track->GetPointX()<<" y "<<track->GetPointY()<<" z "<<track->GetPointZ()<<endl;
          
-         fTransform->Local2Raw(hit,fSlice,i);
+         fTransform->Slice2Sector(fSlice,i,sector,row);
+         fTransform->Local2Raw(hit,sector,row);
          if(hit[1]<0 || hit[1]>fTransform->GetNPads(i) ||
             hit[2]<0 || hit[2]>fTransform->GetNTimeBins())
            {//Track is leaving the patch, so flag the track hits (<0)
@@ -339,14 +436,14 @@ void AliL3Modeller::CalculateCrossingPoints()
        }
     }
   fTracks->Compress();
-  //cout<<"And there are "<<fTracks->GetNTracks()<<" tracks remaining"<<endl;
+  cout<<"And there are "<<fTracks->GetNTracks()<<" tracks remaining"<<endl;
 }
 
 void AliL3Modeller::CheckForOverlaps()
 {
   //Flag the tracks that overlap
   
-  cout<<"Checking for overlaps"<<endl;
+  cout<<"Checking for overlaps...";
   
   for(Int_t i=0; i<fTracks->GetNTracks(); i++)
     {
@@ -356,21 +453,21 @@ void AliL3Modeller::CheckForOverlaps()
        {
          AliL3ModelTrack *track2 = (AliL3ModelTrack*)fTracks->GetCheckedTrack(j);
          if(!track2) continue;
-         for(Int_t k=NRows[fPatch][0]; k<NRows[fPatch][1]; k++)
+         for(Int_t k=NRows[fPatch][0]; k<=NRows[fPatch][1]; k++)
            {
              if(track1->GetPadHit(k)<0 || track1->GetTimeHit(k)<0 ||
                 track2->GetPadHit(k)<0 || track2->GetTimeHit(k)<0)
                continue;
-             if(fabs(track1->GetPadHit(k)-track2->GetPadHit(k)) < fPadOverlap &&
-                fabs(track1->GetTimeHit(k)-track2->GetTimeHit(k)) < fTimeOverlap)
+             if(fabs(track1->GetPadHit(k)-track2->GetPadHit(k)) <= fPadOverlap &&
+                fabs(track1->GetTimeHit(k)-track2->GetTimeHit(k)) <= fTimeOverlap)
                {
-                 track1->SetOverlap(j);
-                 track2->SetOverlap(i);
+                 track2->SetOverlap(k,i);
+                 track1->SetOverlap(k,j);
                }
            }
        }
     }
-  
+  cout<<"done"<<endl;
 }
 
 
index 963f80c..9c16a37 100644 (file)
@@ -8,6 +8,7 @@ class AliL3TrackArray;
 class AliL3MemHandler;
 class AliL3DigitRowData;
 class AliL3Transform;
+class AliL3ModelTrack;
 
 struct Cluster {
   UInt_t fCharge;
@@ -29,17 +30,19 @@ class AliL3Modeller {
   AliL3TrackArray *fTracks; //!
   AliL3MemHandler *fMemHandler; //!
   AliL3DigitRowData *fRowData;//!
-
+  
   AliL3Transform *fTransform; //!
   Int_t fNClusters;
   Int_t fMaxClusters;
   
   Float_t fPadOverlap;
   Float_t fTimeOverlap;
+  Int_t fTrackThreshold; //minimum weigth track need in order to be included.(=Nhits/weight)
   
   Int_t fSlice;
   Int_t fPatch;
   
+  void FillCluster(AliL3ModelTrack *track,Cluster *cluster,Int_t row);
   void CalcClusterWidth(Cluster *cl,Float_t &sigmaY2,Float_t &sigmaZ2);
   void FillZeros(AliL3DigitRowData *digPt,Digit *row);
 
@@ -49,7 +52,7 @@ class AliL3Modeller {
   virtual ~AliL3Modeller();
   
   void Init(Int_t slice,Int_t patch,Char_t *path);
-  void Process();
+  void FindClusters();
   void CheckForOverlaps();
   void CalculateCrossingPoints();
   void WriteRemaining(Char_t *output);