]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - HLT/src/AliL3Evaluate.cxx
Fix for memory leak in fClusters
[u/mrichter/AliRoot.git] / HLT / src / AliL3Evaluate.cxx
index 6da45d54d871ce062237f98eb6c070fd5b9441cf..5ecb2e345d390b8847fc340c800b903182cf92d6 100644 (file)
-//Author:        Anders Strand Vestbo
-//Last Modified: 5.01.2001
+// @(#) $Id$
 
-#include <stdio.h>
+// Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
+//*-- Copyright &copy ALICE HLT Group
+
+#include "AliL3StandardIncludes.h"
 #include <TFile.h>
 #include <TH1.h>
 #include <TParticle.h>
 #include <TTree.h>
 #include <TClonesArray.h>
 
-#include "AliRun.h"
-#include "AliSimDigits.h"
-#include "AliTPC.h"
-#include "AliTPCClustersArray.h"
-#include "AliTPCClustersRow.h"
-#include "AliTPCcluster.h"
-#include "AliTPCParam.h"
+#include <AliRun.h>
+#include <AliSimDigits.h>
+#include <AliTPC.h>
+#include <AliTPCcluster.h>
+#include <AliTPCClustersArray.h>
+#include <AliTPCClustersRow.h>
+#include <AliTPCParam.h>
+#include <AliComplexCluster.h>
+#include <AliStack.h>
+
+#if __GNUC__ == 3
+#include <fstream>
+#include <iosfwd>
+#else
+#include <fstream.h>
+#endif
 
+#include "AliL3Logging.h"
 #include "AliL3Transform.h"
 #include "AliL3SpacePointData.h"
 #include "AliL3Track.h"
 #include "AliL3FileHandler.h"
 #include "AliL3TrackArray.h"
 #include "AliL3Evaluate.h"
-#include "AliL3Logging.h"
 
-//AliL3Evaluate
-//Class for tracking evaluation.
+#if __GNUC__ == 3
+using namespace std;
+#endif
+
+/** \class AliL3Evaluate
+<pre>
+//_____________________________________________________________
+// AliL3Evaluate
+//
+// Evaluation class for tracking. Plots efficiencies etc..
+//
+</pre>
+*/
 
 ClassImp(AliL3Evaluate)
 
 AliL3Evaluate::AliL3Evaluate()
 {
-  fMCFile = NULL;
-  fTracks = NULL;
-  fMCclusterfile = NULL;
+  fTracks = 0;
+  fNFastPoints = 0;
+  fMcIndex = 0;
+  fMcId = 0;
+  fMinSlice=0;
+  fMaxSlice=0;
+  fGoodTracks=0;
+  fNtupleRes=0;
+  fDigitsTree=0;
+  fPtRes=0;
+  fNGoodTracksPt=0;
+  fNFoundTracksPt=0;
+  fNFakeTracksPt=0;
+  fTrackEffPt=0;
+  fFakeTrackEffPt=0;
+  fNGoodTracksEta=0;
+  fNFoundTracksEta=0;
+  fNFakeTracksEta=0;
+  fTrackEffEta=0;
+  fFakeTrackEffEta=0;
+  fMcIndex=0;
+  fMcId=0;
+  fNtuppel=0;
+  fStandardComparison=kTRUE;
 }
 
-AliL3Evaluate::AliL3Evaluate(Char_t *mcfile,Int_t *slice)
+AliL3Evaluate::AliL3Evaluate(Char_t *datapath,Int_t min_clusters,Int_t minhits,Double_t minpt,Double_t maxpt,Int_t *slice)
 {
-  //Normal constructor. Input are the rootfile containing the 
-  //original MC information of the simulated event. 
 
-  fMCFile = new TFile(mcfile,"READ");
-  if(!fMCFile->IsOpen())
+  if(slice)
     {
-      LOG(AliL3Log::kError,"AliL3Evaluation::AliL3Evaluation","File Open")
-       <<"Inputfile "<<mcfile<<" does not exist"<<ENDLOG;
-      return;
+      fMinSlice=slice[0];
+      fMaxSlice=slice[1];
     }
-  
-  fParam = (AliTPCParam*)fMCFile->Get("75x40_100x60");
-  fTransform = new AliL3Transform();
-  
-  fMinSlice = slice[0];
-  fMaxSlice = slice[1];
-
-}
-
-AliL3Evaluate::~AliL3Evaluate()
-{
-  if(fDigitsTree) 
-    fDigitsTree->Delete();
-  
-
-  /*if(fMCFile) 
+  else
     {
-      fMCFile->Close();
-      delete fMCFile;
+      fMinSlice=0;
+      fMaxSlice=35;
     }
-  */
-  if(fTransform)
-    delete fTransform;
-  if(fTracks)
-    delete fTracks;
-  if(fPtRes)
-    delete fPtRes;
-  if(fNGoodTracksPt)
-    delete fNGoodTracksPt;
-  if(fNFoundTracksPt)
-    delete fNFoundTracksPt;
-  if(fNFakeTracksPt)
-    delete fNFakeTracksPt;
-  if(fTrackEffPt)
-    delete fTrackEffPt;
-  if(fFakeTrackEffPt)
-    delete fFakeTrackEffPt;
-  if(fNGoodTracksEta)
-    delete fNGoodTracksEta;
-  if(fNFoundTracksEta)
-    delete fNFoundTracksEta;
-  if(fNFakeTracksEta)
-    delete fNFakeTracksEta;
-  if(fTrackEffEta)
-    delete fTrackEffEta;
-  if(fFakeTrackEffEta)
-    delete fFakeTrackEffEta;
+  sprintf(fPath,"%s",datapath);
+  fMaxFalseClusters = 0.1;
+  fGoodFound = 0;
+  fGoodGen = 0;
+  fMinPointsOnTrack = min_clusters;
+  fMinHitsFromParticle = minhits;
+  fMinGoodPt = minpt;
+  fMaxGoodPt = maxpt;
+  memset(fClusters,0,36*6*sizeof(AliL3SpacePointData*));
+  fTracks=0;
+  fGoodTracks=0;
+  fNtupleRes=0;
+  fDigitsTree=0;
+  fPtRes=0;
+  fNGoodTracksPt=0;
+  fNFoundTracksPt=0;
+  fNFakeTracksPt=0;
+  fTrackEffPt=0;
+  fFakeTrackEffPt=0;
+  fNGoodTracksEta=0;
+  fNFoundTracksEta=0;
+  fNFakeTracksEta=0;
+  fTrackEffEta=0;
+  fFakeTrackEffEta=0;
+  fMcIndex=0;
+  fMcId=0;
+  fNtuppel=0;
+  fStandardComparison=kTRUE;
 }
 
-void AliL3Evaluate::Setup(Char_t *trackfile,Char_t *path)
+void AliL3Evaluate::LoadData(Int_t event,Bool_t sp)
 {
-  //Read in the hit and track information from produced files.
-  
-  Char_t fname[256];
-  AliL3FileHandler *clusterfile[36][5];
+  Char_t fname[1024];
+  AliL3FileHandler *clusterfile[36][6];
   for(Int_t s=fMinSlice; s<=fMaxSlice; s++)
     {
-      for(Int_t p=0; p<5; p++)
+      for(Int_t p=0; p<AliL3Transform::GetNPatches(); p++)
        {
-          fClusters[s][p] = 0;
+         Int_t patch;
+         if(sp==kTRUE)
+           patch=-1;
+         else
+           patch=p;
+         
+         if(fClusters[s][p])
+           delete fClusters[s][p];
+         fClusters[s][p] = 0;
          clusterfile[s][p] = new AliL3FileHandler();
-         sprintf(fname,"%s/points_%d_%d.raw",path,s,p);
+         if(event<0)
+           sprintf(fname,"%s/points_%d_%d.raw",fPath,s,patch);
+         else
+           sprintf(fname,"%s/points_%d_%d_%d.raw",fPath,event,s,patch);
          if(!clusterfile[s][p]->SetBinaryInput(fname))
            {
              LOG(AliL3Log::kError,"AliL3Evaluation::Setup","File Open")
@@ -121,332 +155,467 @@ void AliL3Evaluate::Setup(Char_t *trackfile,Char_t *path)
          fClusters[s][p] = (AliL3SpacePointData*)clusterfile[s][p]->Allocate();
          clusterfile[s][p]->Binary2Memory(fNcl[s][p],fClusters[s][p]);
          clusterfile[s][p]->CloseBinaryInput();
+         if(sp==kTRUE)
+           break;
        }
     }
-  
-  
+   
+  sprintf(fname,"%s/tracks_%d.raw",fPath,event);
   AliL3FileHandler *tfile = new AliL3FileHandler();
-  if(!tfile->SetBinaryInput(trackfile)){
+  if(!tfile->SetBinaryInput(fname)){
     LOG(AliL3Log::kError,"AliL3Evaluation::Setup","File Open")
-    <<"Inputfile "<<trackfile<<" does not exist"<<ENDLOG; 
+      <<"Inputfile "<<fname<<" does not exist"<<ENDLOG; 
     return;
   }
+  if(fTracks)
+    delete fTracks;
   fTracks = new AliL3TrackArray();
   tfile->Binary2TrackArray(fTracks);
   tfile->CloseBinaryInput();
   delete tfile;
+  fTracks->QSort();
 }
 
-void AliL3Evaluate::SetupSlow(Char_t *trackfile,Char_t *path)
-{
-  //Setup for using the slow simulator.
-  
-  fIsSlow = true;
-  Setup(trackfile,path);
-  
-  if(!SetDigitsTree())
-    LOG(AliL3Log::kError,"AliL3Evaluation::SetupSlow","Digits Tree")
-    <<"Error setting up digits tree"<<ENDLOG;
-  if(!SetMCParticleArray())
-    LOG(AliL3Log::kError,"AliL3Evaluation::Setup","Particle array")
-    <<"Error setting up particle array"<<ENDLOG;
-}
 
-void AliL3Evaluate::SetupFast(Char_t *trackfile,Char_t *mcClusterfile,Char_t *path)
+AliL3Evaluate::~AliL3Evaluate()
 {
-  //Setup for using the fast simulator.
-
-  fIsSlow = false;
-  
-  fMCclusterfile = new TFile(mcClusterfile);
-  if(!fMCclusterfile->IsOpen())
-    LOG(AliL3Log::kError,"AliL3Evaluation::SetupFast","File Open")
-      <<"Inputfile "<<mcClusterfile<<" does not exist"<<ENDLOG; 
-
-  Setup(trackfile,path);
-
-  if(!SetMCParticleArray())
-    LOG(AliL3Log::kError,"AliL3Evaluation::SetupFast","Particle array")
-      <<"Error setting up particle array"<<ENDLOG;
+  if(fGoodTracks) delete fGoodTracks;
+  if(fDigitsTree) fDigitsTree->Delete();
+  if(fTracks) delete fTracks;
+  if(fPtRes) delete fPtRes;
+  if(fNGoodTracksPt) delete fNGoodTracksPt;
+  if(fNFoundTracksPt) delete fNFoundTracksPt;
+  if(fNFakeTracksPt) delete fNFakeTracksPt;
+  if(fTrackEffPt) delete fTrackEffPt;
+  if(fFakeTrackEffPt) delete fFakeTrackEffPt;
+  if(fNGoodTracksEta) delete fNGoodTracksEta;
+  if(fNFoundTracksEta) delete fNFoundTracksEta;
+  if(fNFakeTracksEta) delete fNFakeTracksEta;
+  if(fTrackEffEta) delete fTrackEffEta;
+  if(fFakeTrackEffEta) delete fFakeTrackEffEta;
+  if(fMcIndex) delete [] fMcIndex;
+  if(fMcId)    delete [] fMcId;
+  if(fNtuppel) delete fNtuppel;
+  if(fNtupleRes) delete fNtupleRes;
+  for(Int_t s=0; s<=35; s++)
+    for(Int_t p=0; p<6; p++)
+      if(fClusters[s][p]) delete fClusters[s][p];
 }
 
-void AliL3Evaluate::CreateHistos(Int_t nbin,Int_t xlow,Int_t xup)
+void AliL3Evaluate::AssignPIDs()
 {
-  //Create the histograms 
-  
-  fPtRes = new TH1F("fPtRes","Relative Pt resolution",30,-10.,10.); 
-  fNGoodTracksPt = new TH1F("fNGoodTracksPt","Good tracks vs pt",nbin,xlow,xup);    
-  fNFoundTracksPt = new TH1F("fNFoundTracksPt","Found tracks vs pt",nbin,xlow,xup);
-  fNFakeTracksPt = new TH1F("fNFakeTracksPt","Fake tracks vs pt",nbin,xlow,xup);
-  fTrackEffPt = new TH1F("fTrackEffPt","Tracking efficiency vs pt",nbin,xlow,xup);
-  fFakeTrackEffPt = new TH1F("fFakeTrackEffPt","Efficiency for fake tracks vs pt",nbin,xlow,xup);
-  
-  fNGoodTracksEta = new TH1F("fNGoodTracksEta","Good tracks vs eta",20,-50,50);
-  fNFoundTracksEta = new TH1F("fNFoundTracksEta","Found tracks vs eta",20,-50,50);
-  fNFakeTracksEta = new TH1F("fNFakeTracksEta","Fake tracks vs eta",20,-50,50);
-  fTrackEffEta = new TH1F("fTrackEffEta","Tracking efficienct vs eta",20,-50,50);
-  fFakeTrackEffEta = new TH1F("fFakeTrackEffEta","Efficiency for fake tracks vs eta",20,-50,50);
-
+  fTracks->QSort();
+  LOG(AliL3Log::kDebug,"AliL3Evaluate::AssignPIDs","Track Loop")
+    <<"Assigning pid to the found tracks...."<<ENDLOG;
+  for(Int_t i=0; i<fTracks->GetNTracks(); i++)
+    {
+      AliL3Track *track = (AliL3Track*)fTracks->GetCheckedTrack(i);
+      if(!track) continue; 
+      if(track->GetNumberOfPoints() < fMinPointsOnTrack)
+       track->SetPID(0);
+      else {
+       Float_t pid = GetTrackPID(track);
+       track->SetPID(pid);
+      }
+    }
 }
-
-Bool_t AliL3Evaluate::SetDigitsTree()
+  
+void AliL3Evaluate::AssignIDs()
 {
-
-  fMCFile->cd();
-  fDigitsTree = (TTree*)fMCFile->Get("TreeD_75x40_100x60");
-  if(!fDigitsTree) return false;
-  fDigitsTree->GetBranch("Segment")->SetAddress(&fDigits);
-  for(Int_t i=0; i<fDigitsTree->GetEntries(); i++)
+  //Assign MC id to the tracks.
+#ifndef do_mc
+  cerr<<"AliL3Evaluate::AssignIDs() : You need to compile with the do_mc flag!"<<endl;
+  return;
+#endif
+  fGoodFound=0;
+  fTracks->QSort();
+  LOG(AliL3Log::kDebug,"AliL3Evaluate::AssignIDs","Track Loop")
+    <<"Assigning MC id to the found tracks...."<<ENDLOG;
+  for(Int_t i=0; i<fTracks->GetNTracks(); i++)
     {
-      if(!fDigitsTree->GetEvent(i)) continue;
-      Int_t se,ro,slice,slicerow;
-      fParam->AdjustSectorRow(fDigits->GetID(),se,ro);
-      fTransform->Sector2Slice(slice,slicerow,se,ro);
-      fRowid[slice][slicerow] = i;
+      AliL3Track *track = (AliL3Track*)fTracks->GetCheckedTrack(i);
+      if(!track) continue; 
+      if(track->GetNumberOfPoints() < fMinPointsOnTrack) break;
+      
+      fGoodFound++;
+      Int_t tID = GetMCTrackLabel(track);
+      track->SetMCid(tID);
     }
-  
-  return true;
+  //cout<<"Found "<<fGoodFound<<" good tracks "<<endl;
 }
 
-Bool_t AliL3Evaluate::SetMCParticleArray()
+Float_t AliL3Evaluate::GetTrackPID(AliL3Track *track)
 {
-  fMCFile->cd();
-  AliRun *gAlice = (AliRun*)fMCFile->Get("gAlice");
-  if (!gAlice) 
-    {
-      LOG(AliL3Log::kError,"AliL3Evaluate::SetParticleArray","gAlice")
-       <<"AliRun object non existing on file"<<ENDLOG;
-      return false;
+  track->CalculateHelix();
+  // Track dEdx
+  Int_t nc=track->GetNHits();
+  UInt_t *hits = track->GetHitNumbers();
+  Float_t sampleDEdx[159];
+  for (Int_t iHit = 0; iHit < nc; iHit++) {
+    UInt_t hitID = hits[iHit];
+    Int_t iSector = (hitID>>25) & 0x7f;
+    Int_t patch = (hitID>>22) & 0x7;
+    UInt_t position = hitID&0x3fffff;
+    AliL3SpacePointData *points = fClusters[iSector][patch];
+    if(!points) continue; 
+    if(position>=fNcl[iSector][patch]) 
+      {
+       LOG(AliL3Log::kError,"AliL3Evaluate::GetMCTrackLabel","Clusterarray")
+         <<AliL3Log::kDec<<"ERROR"<<ENDLOG;
+       continue;
+      }
+    UChar_t padrow = points[position].fPadRow;
+    Float_t pWidth = AliL3Transform::GetPadPitchWidthLow();
+    if (padrow>63)
+      pWidth = AliL3Transform::GetPadPitchWidthUp(); 
+    Float_t corr=1.; if (padrow>63) corr=0.67;
+    sampleDEdx[iHit] = points[position].fCharge/pWidth*corr;
+    Double_t crossingangle = track->GetCrossingAngle(padrow,iSector);
+    Double_t s = sin(crossingangle);
+    Double_t t = track->GetTgl();
+    sampleDEdx[iHit] *= sqrt((1-s*s)/(1+t*t));
+  }
+
+  /* Cook dEdx */
+  Int_t i;
+  Int_t swap;//stupid sorting
+  do {
+    swap=0;
+    for (i=0; i<nc-1; i++) {
+      if (sampleDEdx[i]<=sampleDEdx[i+1]) continue;
+      Float_t tmp=sampleDEdx[i];
+      sampleDEdx[i]=sampleDEdx[i+1]; sampleDEdx[i+1]=tmp;
+      swap++;
     }
-  
-  gAlice->GetEvent(0);
-  fParticles=gAlice->Particles(); 
-  return true;
-}
+  } while (swap);
 
+  Double_t low=0.05; Double_t up=0.7;
+  Int_t nl=Int_t(low*nc), nu=Int_t(up*nc);
+  Float_t trackDEdx=0;
+  for (i=nl; i<=nu; i++) trackDEdx += sampleDEdx[i];
+  trackDEdx /= (nu-nl+1);
 
-TObjArray *AliL3Evaluate::DefineGoodTracks(Int_t slice,Int_t *padrow,Int_t good_number,Int_t *particle_id)
-{
-  //Loop over MC particles, and mark the good ones
-  //(which the tracker should find...)
+  //  cout<<" PID: "<<nc<<" "<<nl<<" "<<nu<<" "<<trackDEdx<<" "<<track->GetPt()<<endl;
+  return trackDEdx;
+}
 
-  Int_t np=fParticles->GetEntriesFast();
-  AliTPC *TPC = (AliTPC*)gAlice->GetDetector("TPC");
-  TPC->SetParam(fParam);
-  Int_t ver = TPC->IsVersion();
-  LOG(AliL3Log::kInformational,"AliL3Evaluate::DefineGoodTracks","TPC version")
-    <<"TPC version "<<ver<<" found on file"<<ENDLOG;
-  
-  //Int_t nrow_up=TPC->GetParam()->GetNRowUp();
-  //Int_t nrows=TPC->GetParam()->GetNRowLow()+nrow_up;
-  Int_t zero=TPC->GetParam()->GetZeroSup();
+struct S {Int_t lab; Int_t max;};
+Int_t AliL3Evaluate::GetMCTrackLabel(AliL3Track *track){ 
+  //Returns the MCtrackID of the belonging clusters.
+  //If MCLabel < 0, means that track is fake.
+  //Definitions are identical to offline.
+  //Fake track means:
+  // - more than 10 percent of clusters are assigned incorrectly
+  // - more than half of the innermost 10% of clusters were assigned incorrectly.
   
-  //Int_t number_of_rows = padrow[1] - padrow[0] + 1;
   
-  Int_t *good = new Int_t[np];
-  for(Int_t ii=0; ii<np; ii++)
-    good[ii] = 0;
-
-  if(ver==1)
+#ifdef do_mc
+  Int_t num_of_clusters = track->GetNumberOfPoints();
+  S *s=new S[num_of_clusters];
+  Int_t i;
+  for (i=0; i<num_of_clusters; i++) s[i].lab=s[i].max=0;
+  UInt_t *hitnum = track->GetHitNumbers();  
+  UInt_t id;
+    
+  Int_t lab=123456789;
+  for (i=0; i<num_of_clusters; i++) 
     {
-      if(fIsSlow)
-       LOG(AliL3Log::kError,"AliL3Evaluate::DefineGoodTracks","TPC version")
-         <<"TPC version "<<ver<<" does not match."<<ENDLOG;
-      fMCclusterfile->cd();
-      AliTPCClustersArray carray;
-      carray.Setup(fParam);
-      carray.SetClusterType("AliTPCcluster");
-      Bool_t clusterok = carray.ConnectTree("Segment Tree");
-      if(!clusterok) 
-       LOG(AliL3Log::kError,"AliL3Evaluate::DefineGoodTracks","Cluster Array")
-         <<"Error loading clusters from rootfile"<<ENDLOG;
-            
-      for(Int_t i=0; i<carray.GetTree()->GetEntries(); i++)  
+      //Tricks to get the clusters belonging to this track:
+      id = hitnum[i];
+      Int_t slice = (id>>25) & 0x7f;
+      Int_t patch = (id>>22) & 0x7;
+      UInt_t pos = id&0x3fffff;              
+      
+      AliL3SpacePointData *points = fClusters[slice][patch];
+      if(!points) continue; 
+      if(pos>=fNcl[slice][patch]) 
        {
-         Int_t sec,row,sl,lr;
-         AliSegmentID *s = carray.LoadEntry(i);
-         fParam->AdjustSectorRow(s->GetID(),sec,row);
-         fTransform->Sector2Slice(sl,lr,sec,row);
-         
-         if(sl != slice) {carray.ClearRow(sec,row); continue;}
-         if(lr < padrow[0]) {carray.ClearRow(sec,row); continue;}
-         if(lr > padrow[1]) {carray.ClearRow(sec,row); continue;}
-         AliTPCClustersRow *cRow = carray.GetRow(sec,row);
-         for(Int_t j=0; j<cRow->GetArray()->GetEntriesFast(); j++)
-           {
-             AliTPCcluster *cluster=(AliTPCcluster*)(*cRow)[j];
-             Int_t lab=cluster->GetLabel(0);
-             if(lab<0) continue;
-             lab=TMath::Abs(lab);
-             good[lab]++;
-           }
-         if(carray.GetRow(sec,row)) 
-           carray.ClearRow(sec,row);
+         LOG(AliL3Log::kError,"AliL3Evaluate::GetMCTrackLabel","Clusterarray")
+           <<AliL3Log::kDec<<"ERROR"<<ENDLOG;
+         continue;
        }
+      
+      //Get the label of the cluster:
+      lab=points[pos].fTrackID[0];
+
+      Int_t j;
+      for (j=0; j<num_of_clusters; j++)
+        if (s[j].lab==lab || s[j].max==0) break;
+      s[j].lab=lab;
+      s[j].max++;
     }
-  else if(ver==2)
+  
+  Int_t max=0;
+  for (i=0; i<num_of_clusters; i++) 
+    if (s[i].max>max) {max=s[i].max; lab=s[i].lab;}
+  
+  if(lab == -1)
+    return -1; //If most clusters is -1, this is a noise track.
+  if(lab < 0)
+    cerr<<"AliL3Evaluate::GetMCTrackLabel : Track label negative :"<<lab<<endl;
+  
+  delete[] s;
+  
+  for (i=0; i<num_of_clusters; i++) 
     {
-      if(!fIsSlow)
-       LOG(AliL3Log::kError,"AliL3Evaluate::DefineGoodTracks","TPC version")
-         <<"TPC version "<<ver<<" does not match."<<ENDLOG;
-      Int_t *count = new Int_t[np]; //np number of particles.
-      Int_t i;
-      for (i=0; i<np; i++) count[i]=0;
-      for (i=padrow[0]; i<=padrow[1]; i++) {
-       Int_t index = fRowid[slice][i];
-       if (!fDigitsTree->GetEvent(index)) continue;
-       Int_t sec,row;
-       fParam->AdjustSectorRow(fDigits->GetID(),sec,row);
-       fDigits->First();
-       while (fDigits->Next()) {
-         Int_t it=fDigits->CurrentRow(), ip=fDigits->CurrentColumn();
-         Short_t dig = fDigits->GetDigit(it,ip);
-         Int_t idx0=fDigits->GetTrackID(it,ip,0); 
-         Int_t idx1=fDigits->GetTrackID(it,ip,1);
-         Int_t idx2=fDigits->GetTrackID(it,ip,2);
-         if (idx0>=0 && dig>=zero) count[idx0]+=1;
-         if (idx1>=0 && dig>=zero) count[idx1]+=1;
-         if (idx2>=0 && dig>=zero) count[idx2]+=1;
-       }
-       for (Int_t j=0; j<np; j++) {
-         if (count[j]>1) {//at least two digits at this padrow 
-           good[j]++;
-         }
-         count[j]=0;
+      id = hitnum[i];
+      Int_t slice = (id>>25) & 0x7f;
+      Int_t patch = (id>>22) & 0x7;
+      UInt_t pos = id&0x3fffff;              
+      
+      AliL3SpacePointData *points = fClusters[slice][patch];
+      if(!points) continue; 
+      if(pos>=fNcl[slice][patch]) 
+       {
+         LOG(AliL3Log::kError,"AliL3Evaluate::GetMCTrackLabel","Clusterarray")
+           <<AliL3Log::kDec<<"ERROR"<<ENDLOG;
+         continue;
        }
-      }
-      delete[] count;
+      
+      if (abs(points[pos].fTrackID[1]) == lab || 
+         abs(points[pos].fTrackID[2]) == lab ) max++;
     }
-  else 
+  
+  
+  //Check if more than 10% of the clusters were assigned incorrectly:
+  if (1.-Float_t(max)/num_of_clusters > fMaxFalseClusters) 
     {
-      LOG(AliL3Log::kError,"AliL3Evaluation::FillEffHistos","TPC version")
-       <<"No valid TPC version found"<<ENDLOG;
-      return 0;
+      return -lab;
     }
-  
-  Float_t torad=TMath::Pi()/180;
-  
-  Float_t phi_min = slice*20 - 10;
-  Float_t phi_max = slice*20 + 10;
-  TObjArray *good_part = new TObjArray();
-  
-  for(Int_t i=0; i<fParticles->GetEntriesFast(); i++)
-   {
-     TParticle *p = (TParticle*)fParticles->UncheckedAt(i);
-     if(p->GetFirstMother()>0) continue; //secondary particle
-     if(good[i] < good_number) {continue;}
-     
-     Double_t ptg=p->Pt(),pxg=p->Px(),pyg=p->Py(),pzg=p->Pz();
-     Double_t phi_part = TMath::ATan2(pyg,pxg);
-     if (phi_part < 0) phi_part += 2*TMath::Pi();
-     
-
-     if(phi_part < phi_min*torad || phi_part > phi_max*torad) {continue;}
-     if(ptg<0.100) continue;
-     if(fabs(pzg/ptg)>0.999) {continue;}
-     Int_t entries = good_part->GetEntriesFast();
-     good_part->AddLast(p);
-     particle_id[entries] = i;
-   }
-  delete [] good;
-  LOG(AliL3Log::kInformational,"AliL3Evaluate::DefineGoodTracks","NPart")
-  <<AliL3Log::kDec<<"Found "<<good_part->GetEntriesFast()
-  <<" good Particles (Tracks) out of "<<fParticles->GetEntriesFast()<<ENDLOG;
-
-  return good_part;
+  else //Check if at least half of the 10% innermost clusters are assigned correctly.
+    {
+      Int_t tail=Int_t(0.10*num_of_clusters);
+      max=0;
+      for (i=1; i<=tail; i++) 
+       {
+         id = hitnum[num_of_clusters - i];
+         Int_t slice = (id>>25) & 0x7f;
+         Int_t patch = (id>>22) & 0x7;
+         UInt_t pos = id&0x3fffff;           
+         
+         AliL3SpacePointData *points = fClusters[slice][patch];
+         if(lab == abs(points[pos].fTrackID[0]) ||
+            lab == abs(points[pos].fTrackID[1]) ||
+            lab == abs(points[pos].fTrackID[2])) max++;
+       }
+      if (max < Int_t(0.5*tail)) return -lab;
+    }
+
+  return lab;
+#else //If we are running with mc_ids or not
+  return 0;
+#endif
+
 }
 
-void AliL3Evaluate::EvaluatePatch(Int_t slice,Int_t patch,Int_t min_points,Int_t good_number)
+void AliL3Evaluate::GetFastClusterIDs(Char_t *path)
 {
-  //Make efficiency plots for tracking on patch level (before any merging).
-  
-  Int_t row[5][2] = {{ 0, 45},{46,77},{78,109},{110,141},{142,173}};
-  Int_t *particle_id = new Int_t[fParticles->GetEntriesFast()];
-  TObjArray *good_particles = DefineGoodTracks(slice,row[patch],good_number,particle_id);
-  SetMinPoints(min_points);
-  AssignIDs();
-  CreateHistos();
-  FillEffHistos(good_particles,particle_id);
-  CalcEffHistos();
-  delete good_particles;
-  delete [] particle_id;
+  //Get the MC id of space points in case of using the fast simulator. 
+  char fname[256];
+  sprintf(fname,"%s/point_mc.dat",path);
+  FILE *infile = fopen(fname,"r");
+  if(!infile) return;
+  Int_t hitid,hitmc,i;
+  
+  for(i=0; ; i++)
+    if(fscanf(infile,"%d %d",&hitid,&hitmc)==EOF) break;
+  rewind(infile);
+  fNFastPoints = i;
+  fMcId = new Int_t[fNFastPoints];
+  fMcIndex = new UInt_t[fNFastPoints];
+  
+  for(i=0; i<fNFastPoints; i++)
+    {
+      if(fscanf(infile,"%d %d",&hitid,&hitmc)==EOF) break;
+      fMcId[i] = hitmc;
+      fMcIndex[i] = hitid;
+    }
+  fclose(infile);
 }
 
-void AliL3Evaluate::EvaluateSlice(Int_t slice,Int_t min_points,Int_t good_number)
+void AliL3Evaluate::CreateHistos(Int_t nbin,Float_t xlow,Float_t xup)
 {
-  //Make efficiency plots for tracking on a slice (after merging).
-  //min_points = minimum points on track to be considered for evaluation
-  //good_number = minimum hits (padrows) produced by simulated track for consideration.
-
-  Int_t row[2] = {0,173};
-  Int_t *particle_id = new Int_t[fParticles->GetEntriesFast()];
-  TObjArray *good_particles = DefineGoodTracks(slice,row,good_number,particle_id);
-  SetMinPoints(min_points);
-  AssignIDs();
-  CreateHistos();
-  FillEffHistos(good_particles,particle_id);
-  CalcEffHistos();
-  delete good_particles;
-  delete [] particle_id;
+  //Create the histograms 
+  
+  LOG(AliL3Log::kInformational,"AliL3Evaluate::CreateHistos","Allocating")
+    <<"Creating histograms..."<<ENDLOG;
+  
+  fNtuppel = new TNtuple("fNtuppel","Pt resolution","pt_gen:pt_found:nHits");
+  fNtuppel->SetDirectory(0);
+  fPtRes = new TH1F("fPtRes","Relative Pt resolution",30,-10.,10.); 
+  fNGoodTracksPt = new TH1F("fNGoodTracksPt","Good tracks vs pt",nbin,xlow,xup);    
+  fNFoundTracksPt = new TH1F("fNFoundTracksPt","Found tracks vs pt",nbin,xlow,xup);
+  fNFakeTracksPt = new TH1F("fNFakeTracksPt","Fake tracks vs pt",nbin,xlow,xup);
+  fTrackEffPt = new TH1F("fTrackEffPt","Tracking efficiency vs pt",nbin,xlow,xup);
+  fFakeTrackEffPt = new TH1F("fFakeTrackEffPt","Efficiency for fake tracks vs pt",nbin,xlow,xup);
+  
+  fNGoodTracksEta = new TH1F("fNGoodTracksEta","Good tracks vs eta",20,-50,50);
+  fNFoundTracksEta = new TH1F("fNFoundTracksEta","Found tracks vs eta",20,-50,50);
+  fNFakeTracksEta = new TH1F("fNFakeTracksEta","Fake tracks vs eta",20,-50,50);
+  fTrackEffEta = new TH1F("fTrackEffEta","Tracking efficienct vs eta",20,-50,50);
+  fFakeTrackEffEta = new TH1F("fFakeTrackEffEta","Efficiency for fake tracks vs eta",20,-50,50);
+
 }
 
-void AliL3Evaluate::EvaluateGlobal(Int_t min_points,Int_t good_number)
+void AliL3Evaluate::GetGoodParticles(Char_t *path,Int_t event,Int_t *padrowrange)
 {
-  //Make efficiency plots for tracking on several slices.
-  
-  Int_t row[2] = {0,173};
-  Int_t *particle_id = new Int_t[fParticles->GetEntriesFast()];
-  SetMinPoints(min_points);
-  AssignIDs();
-  CreateHistos();
-  for(Int_t slice=fMinSlice;slice<=fMaxSlice;slice++){
-    TObjArray *good_particles = DefineGoodTracks(slice,row,good_number,particle_id);
-    FillEffHistos(good_particles,particle_id);
-    delete good_particles;
+  //Read the good particles from file. This file should already have been
+  //generated by macro AliTPCComparison.C.
+  
+  Char_t filename[1024];
+  if(event<0 && !padrowrange)
+    sprintf(filename,"%s/good_tracks_tpc",path);
+  else if(event>=0 && !padrowrange)
+    sprintf(filename,"%s/good_tracks_tpc_%d",path,event);
+  else
+    sprintf(filename,"%s/good_tracks_tpc_%d_%d_%d",path,event,padrowrange[0],padrowrange[1]);
+  ifstream in(filename);
+  if(!in)
+    {
+      cerr<<"AliL3Evaluate::GetGoodParticles : Problems opening file :"<<filename<<endl;
+      return;
+    }
+  Int_t MaxTracks=20000;
+  if(fGoodTracks)
+    delete [] fGoodTracks;
+  fGoodGen=0;
+  fGoodTracks = new GoodTrack[MaxTracks];
+  
+  if(fStandardComparison){
+    while (in>>fGoodTracks[fGoodGen].label>>fGoodTracks[fGoodGen].code>>
+          fGoodTracks[fGoodGen].px>>fGoodTracks[fGoodGen].py>>fGoodTracks[fGoodGen].pz>>
+          fGoodTracks[fGoodGen].x>>fGoodTracks[fGoodGen].y>>fGoodTracks[fGoodGen].z)
+      {
+        fGoodTracks[fGoodGen].nhits=-1;
+        fGoodTracks[fGoodGen].sector=-1; 
+       fGoodGen++;
+       if (fGoodGen==MaxTracks) 
+         {
+           cerr<<"AliL3Evaluate::GetGoodParticles : Too many good tracks !\n";
+           break;
+         }
+      }
+  } else {
+    while (in>>fGoodTracks[fGoodGen].label>>fGoodTracks[fGoodGen].code>>
+          fGoodTracks[fGoodGen].px>>fGoodTracks[fGoodGen].py>>fGoodTracks[fGoodGen].pz>>
+          fGoodTracks[fGoodGen].x>>fGoodTracks[fGoodGen].y >>fGoodTracks[fGoodGen].z>>fGoodTracks[fGoodGen].nhits>>fGoodTracks[fGoodGen].sector) 
+      {
+       fGoodGen++;
+       if (fGoodGen==MaxTracks) 
+         {
+           cerr<<"AliL3Evaluate::GetGoodParticles : Too many good tracks !\n";
+           break;
+         }
+      }
   }
-  CalcEffHistos();
 }
 
-void AliL3Evaluate::FillEffHistos(TObjArray *good_particles,Int_t *particle_id)
-{  
-  //Fill the efficiency histograms.
+//has to be modified for fakes.
 
-  for(Int_t i=0; i<good_particles->GetEntriesFast(); i++)
+void AliL3Evaluate::FillEffHistos()
+{  
+  if(!fGoodTracks)
     {
-      TParticle *p = (TParticle*)good_particles->UncheckedAt(i);
-      Double_t ptg=p->Pt(),pzg=p->Pz();
+      cerr<<"AliL3Evaluate::FillEffHistos : No good tracks"<<endl;
+      return;
+    }
+  //cout<<"Comparing "<<fGoodGen<<" good tracks ..."<<endl;
+  for(Int_t i=0; i<fGoodGen; i++)
+    {
+      //cout<<"Checking particle "<<i<<endl;
+      if(!fStandardComparison) 
+       if(fGoodTracks[i].nhits < fMinHitsFromParticle) continue;
+      Float_t ptg = TMath::Sqrt(fGoodTracks[i].px*fGoodTracks[i].px + fGoodTracks[i].py*fGoodTracks[i].py);
+      if(ptg < fMinGoodPt || ptg > fMaxGoodPt) continue;
+      Float_t pzg=fGoodTracks[i].pz;
       Float_t dipangle=TMath::ATan2(pzg,ptg)*180./TMath::Pi();
+      
+      //If we are only considering tracks on one side of the TPC:
+      if(fMaxSlice <= 17)
+       if(dipangle < 0)
+         continue;
+
       fNGoodTracksPt->Fill(ptg);
       fNGoodTracksEta->Fill(dipangle);
       Int_t found = 0;
+      
       for(Int_t k=0; k<fTracks->GetNTracks(); k++)
        {
          AliL3Track *track = fTracks->GetCheckedTrack(k);
          if(!track) continue;
          Int_t nHits = track->GetNumberOfPoints();
          if(nHits < fMinPointsOnTrack) break;
-         
          Int_t tracklabel;
          tracklabel = track->GetMCid();
          
-         if(TMath::Abs(tracklabel) != particle_id[i]) continue;
+         if(TMath::Abs(tracklabel) != fGoodTracks[i].label) continue;
          found=1;
-         if(tracklabel == particle_id[i]) {fNFoundTracksPt->Fill(ptg); fNFoundTracksEta->Fill(dipangle);}
-         else {fNFakeTracksPt->Fill(ptg); fNFakeTracksEta->Fill(dipangle);}
          Float_t pt=track->GetPt();
-         fPtRes->Fill((pt-ptg)/ptg*100.);
+         if(tracklabel == fGoodTracks[i].label) 
+           {
+             fNFoundTracksPt->Fill(ptg); 
+             fNFoundTracksEta->Fill(dipangle);
+             fNtuppel->Fill(ptg,pt,nHits);
+             fPtRes->Fill((pt-ptg)/ptg*100.);
+           }
+         else 
+           {
+             fNFakeTracksPt->Fill(ptg); 
+             fNFakeTracksEta->Fill(dipangle);
+           }
+         //fPtRes->Fill((pt-ptg)/ptg*100.);
+         //fNtuppel->Fill(ptg,pt,nHits);
          break;
          
        }
+      //if(!found)
+      //cout<<"Track "<<fGoodTracks[i].label<<" was not found"<<endl;
     }
 }
 
-void AliL3Evaluate::CalcEffHistos(){  
+void AliL3Evaluate::FillEffHistosNAIVE()
+{  
+  //Fill the efficiency histograms.
+  
+  cout<<endl<<"Note: Doing NAIVE evaluation "<<endl;
+  for(Int_t i=0; i<fGoodGen; i++)
+    {
+      if(!fStandardComparison) 
+       if(fGoodTracks[i].nhits < fMinHitsFromParticle) continue;
+      Double_t ptg=TMath::Sqrt(fGoodTracks[i].px*fGoodTracks[i].px + fGoodTracks[i].py*fGoodTracks[i].py);
+      if(ptg < fMinGoodPt || ptg > fMaxGoodPt) continue;
+      Double_t pzg=fGoodTracks[i].pz;
+      Float_t dipangle=TMath::ATan2(pzg,ptg)*180./TMath::Pi();
+      //printf("filling particle with pt %f and dipangle %f\n",ptg,dipangle);
+      fNGoodTracksPt->Fill(ptg);
+      fNGoodTracksEta->Fill(dipangle);
+      
+    }
   
+  for(Int_t k=0; k<fTracks->GetNTracks(); k++)
+    {
+      AliL3Track *track = fTracks->GetCheckedTrack(k);
+      if(!track) continue;
+      Int_t nHits = track->GetNumberOfPoints();
+      if(nHits < fMinPointsOnTrack) break;
+      if(track->GetPt()<fMinGoodPt || track->GetPt() > fMaxGoodPt) continue;
+      if(fabs(track->GetPseudoRapidity())>0.9) continue;
+
+      fNFoundTracksPt->Fill(track->GetPt()); fNFoundTracksEta->Fill(track->GetPseudoRapidity());
+      //Float_t pt=track->GetPt();
+      //fPtRes->Fill((pt-ptg)/ptg*100.);
+      //fNtuppel->Fill(ptg,pt,nHits);
+            
+    }
+}
+
+void AliL3Evaluate::CalcEffHistos()
+{  
+
   Stat_t ngood=fNGoodTracksPt->GetEntries();
   Stat_t nfound=fNFoundTracksPt->GetEntries();
   Stat_t nfake=fNFakeTracksPt->GetEntries();
-  
+
   LOG(AliL3Log::kInformational,"AliL3Evaluate::FillEffHistos","Efficiency")
     <<AliL3Log::kDec<<"There was "<<ngood<<" generated good tracks"<<ENDLOG;
   LOG(AliL3Log::kInformational,"AliL3Evaluate::FillEffHistos","Efficiency")
@@ -455,7 +624,9 @@ void AliL3Evaluate::CalcEffHistos(){
     <<AliL3Log::kDec<<"Integral efficiency is about "<<nfound/ngood*100<<ENDLOG;
   LOG(AliL3Log::kInformational,"AliL3Evaluate::FillEffHistos","Efficiency")
     <<AliL3Log::kDec<<"Fake tracks relative is about "<<nfake/ngood*100<<ENDLOG;
-  
+  //LOG(AliL3Log::kInformational,"AliL3Evaluate::FillEffHistos","Efficiency")
+  //<<AliL3Log::kDec<<"Naive efficiency "<<(Double_t)fGoodFound/(Double_t)fGoodGen<<ENDLOG;
+
   fNFoundTracksPt->Sumw2(); fNGoodTracksPt->Sumw2();
   fTrackEffPt->Divide(fNFoundTracksPt,fNGoodTracksPt,1,1.,"b");
   fFakeTrackEffPt->Divide(fNFakeTracksPt,fNGoodTracksPt,1,1.,"b");
@@ -465,7 +636,7 @@ void AliL3Evaluate::CalcEffHistos(){
   fFakeTrackEffPt->SetFillStyle(3013);
   fTrackEffPt->SetLineColor(4);
   fFakeTrackEffPt->SetFillColor(2);
-  
+
   fNFoundTracksEta->Sumw2(); fNGoodTracksEta->Sumw2();
   fTrackEffEta->Divide(fNFoundTracksEta,fNGoodTracksEta,1,1.,"b");
   fFakeTrackEffEta->Divide(fNFakeTracksEta,fNGoodTracksEta,1,1.,"b");
@@ -478,268 +649,498 @@ void AliL3Evaluate::CalcEffHistos(){
        
 }
 
-void AliL3Evaluate::AssignIDs()
+void AliL3Evaluate::Write2File(Char_t *outputfile)
 {
-  //Assign MC id to the tracks.
-  
-  UInt_t *index=0,tmp_ind=0;
-  Int_t *pID=0,npoints=0;
+  //Write histograms to file:
   
-  if(!fIsSlow)
+  TFile *of = TFile::Open(outputfile,"RECREATE");
+  if(!of->IsOpen())
     {
-      pID = GetFastIDs(tmp_ind,npoints);
-      index = (UInt_t*)tmp_ind;
+      LOG(AliL3Log::kError,"AliL3Evaluate::Write2File","File Open")
+       <<"Problems opening rootfile"<<ENDLOG;
+      return;
     }
   
-  fTracks->QSort();
-  LOG(AliL3Log::kDebug,"AliL3Evaluate::AssignIDs","Track Loop")
-    <<"Assigning MC id to the found tracks...."<<ENDLOG;
-  for(Int_t i=0; i<fTracks->GetNTracks(); i++)
-    {
-      AliL3Track *track = (AliL3Track*)fTracks->GetCheckedTrack(i);
-      if(!track) 
-       {
-         LOG(AliL3Log::kWarning,"AliL3Evaluate::AssignIDs","Track Loop")
-           <<AliL3Log::kDec<<"No track in track array, index "<<i<<ENDLOG;
-         continue;
-       }
-      if(track->GetNumberOfPoints() < fMinPointsOnTrack) break;
-      Int_t tID;
-      if(!fIsSlow)
-       tID = GetMCTrackLabel(track,index,pID,npoints);
-      else
-       tID = GetMCTrackLabel(track);
-      track->SetMCid(tID);
-      printf("track %i id %d\n",i,tID);
-    }
+  of->cd();
+  fNtuppel->Write();
+  fPtRes->Write();
+  fNGoodTracksPt->Write();
+  fNFoundTracksPt->Write();
+  fNFakeTracksPt->Write();
+  fTrackEffPt->Write();
+  fFakeTrackEffPt->Write();
+  fNGoodTracksEta->Write();
+  fNFoundTracksEta->Write();
+  fNFakeTracksEta->Write();
+  fTrackEffEta->Write();
+  fFakeTrackEffEta->Write();
   
-  if(pID)   delete [] pID;
-  if(index) delete [] index;
+  of->Close();
+
 }
 
+TNtuple *AliL3Evaluate::GetNtuple()
+{
+  if(!fNtupleRes)
+    {
+      fNtupleRes = new TNtuple("ntuppel","Residuals","residual_trans:residual_long:zHit:pt:dipangle:beta:padrow:nHits");
+      fNtupleRes->SetDirectory(0);//Bug in older version of root.
+    }
+  return fNtupleRes;
+}
 
-struct S {Int_t lab; Int_t max;};
-Int_t AliL3Evaluate::GetMCTrackLabel(AliL3Track *track,UInt_t *index,Int_t *pID,Int_t npoints) 
+void AliL3Evaluate::CalculateResiduals()
 {
-  //Returns the MCtrackID of the belonging clusters.
-  //If MCLabel < 0, means that track is fake.
-  //Fake track means that more than 10 percent of clusters are assigned incorrectly.
-  
-  Int_t num_of_clusters = track->GetNumberOfPoints();
-  S *s=new S[num_of_clusters];
-  Int_t i;
-  for (i=0; i<num_of_clusters; i++) s[i].lab=s[i].max=0;
-  UInt_t *hitnum = track->GetHitNumbers();  
-  UInt_t id;
 
-  Int_t **trackID;
-    
-  if(fIsSlow)
-    trackID = GetClusterIDs(track);
-  else
-    trackID = GetClusterIDs(track,index,pID,npoints);
+  TNtuple *ntuppel = GetNtuple();
   
-
-  Int_t lab=123456789;
-  for (i=0; i<num_of_clusters; i++) 
+  for(Int_t i=0; i<fTracks->GetNTracks(); i++)
     {
-      //Tricks to get the clusters belonging to this track:
-      id = hitnum[i];
-      Int_t slice = (id>>25) & 0x7f;
-      Int_t patch = (id>>22) & 0x7;
-      UInt_t pos = id&0x3fffff;              
       
-      AliL3SpacePointData *points = fClusters[slice][patch];
+      AliL3Track *track = (AliL3Track*)fTracks->GetCheckedTrack(i);
+      if(!track) continue;
+      if(track->GetNHits() < fMinPointsOnTrack) break;
       
-      if(!points) 
-       continue;
+      track->CalculateHelix();
+      UInt_t *hitnum = track->GetHitNumbers();
+      UInt_t id;
       
-      if(pos>=fNcl[slice][patch]) 
+      Float_t xyz[3];
+      Int_t padrow;
+      for(Int_t j=0; j<track->GetNumberOfPoints()-1; j++)
        {
-         LOG(AliL3Log::kError,"AliL3Evaluate::GetMCTrackLabel","Clusterarray")
-           <<AliL3Log::kDec<<"ERROR"<<ENDLOG;
-         continue;
+         id = hitnum[j];
+         Int_t slice = (id>>25) & 0x7f;
+         Int_t patch = (id>>22) & 0x7;
+         UInt_t pos = id&0x3fffff;           
+
+         //if(slice<18) continue;
+         
+         AliL3SpacePointData *points = fClusters[slice][patch];
+         
+         if(!points) 
+           {
+             LOG(AliL3Log::kError,"AliL3Evaluate::CalculateResiduals","Clusterarray")
+               <<"No points at slice "<<slice<<" patch "<<patch<<" pos "<<pos<<ENDLOG;
+             continue;
+           }
+         if(pos>=fNcl[slice][patch]) 
+           {
+             LOG(AliL3Log::kError,"AliL3Evaluate::CalculateResiduals","Clusterarray")
+               <<AliL3Log::kDec<<"ERROR"<<ENDLOG;
+             continue;
+           }
+         
+         xyz[0] = points[pos].fX;
+         xyz[1] = points[pos].fY;
+         xyz[2] = points[pos].fZ;
+         padrow = points[pos].fPadRow;
+         //AliL3Transform::Global2Local(xyz,slice,kTRUE);
+         AliL3Transform::Global2LocHLT(xyz,slice);
+         
+         Float_t angle = 0;
+         AliL3Transform::Local2GlobalAngle(&angle,slice);
+         if(!track->CalculateReferencePoint(angle,AliL3Transform::Row2X(padrow)))
+           {
+             LOG(AliL3Log::kError,"AliL3Evaluate::CalculateResiduals","Crossing point")
+               <<"Track does not crossing padrow "<<padrow<<" in slice "<<slice<<ENDLOG;
+             continue;
+           }
+         
+         Float_t xyz_cross[3] = {track->GetPointX(),track->GetPointY(),track->GetPointZ()};
+         //AliL3Transform::Global2Local(xyz_cross,slice,kTRUE);          
+         AliL3Transform::Global2LocHLT(xyz_cross,slice);
+         
+         Double_t beta = track->GetCrossingAngle(padrow,slice);
+         
+         Double_t yres = xyz_cross[1] - xyz[1];
+         Double_t zres = xyz_cross[2] - xyz[2];
+         Double_t dipangle = atan(track->GetTgl());
+         ntuppel->Fill(yres,zres,xyz_cross[2],track->GetPt(),dipangle,beta,padrow,track->GetNumberOfPoints());
+         
        }
-      
-      //Get the label of the cluster:
-      //printf("label %d %d %d\n",trackID[i][0],trackID[i][1],trackID[i][2]);
-      lab=abs(trackID[i][0]);
-      Int_t j;
-      for (j=0; j<num_of_clusters; j++)
-        if (s[j].lab==lab || s[j].max==0) break;
-      s[j].lab=lab;
-      s[j].max++;
     }
+}
+
+enum tagprimary {kPrimaryCharged = 0x4000};
+void AliL3Evaluate::EvaluatePoints(Char_t *rootfile,Char_t *exactfile,Char_t *tofile,Int_t nevent,Bool_t offline,Bool_t sp)
+{
+  //Compare points to the exact crossing points of track and padrows.
+  //The input file to this function, contains the exact clusters calculated
+  //in AliTPC::Hits2ExactClusters.
   
-  Int_t max=0;
-  for (i=0; i<num_of_clusters; i++) 
-    if (s[i].max>max) {max=s[i].max; lab=s[i].lab;}
+#ifndef do_mc
+  cerr<<"AliL3Evaluate::EvaluatePoints : Compile with do_mc flag!"<<endl;
+  return;
+#else
+  cout<<"Evaluating points"<<endl;
+  TNtuple *ntuppel = new TNtuple("ntuppel_res","Cluster properties",
+                                "slice:padrow:charge:resy:resz:zHit:pt:beta:sigmaY2:sigmaZ2:psigmaY2:psigmaZ2");
+  ntuppel->SetDirectory(0);
   
-  delete[] s;
+  TNtuple *ntuppel2 = new TNtuple("ntuppel_eff","Efficiency","slice:padrow:nfound:ngen");
+  ntuppel2->SetDirectory(0);
+
+  TFile *exfile = TFile::Open(rootfile);
+  if(!exfile)
+    {
+      cerr<<"Error opening rootfile "<<rootfile<<endl;
+      return;
+    }
+  gAlice = (AliRun*)exfile->Get("gAlice");
+  if (!gAlice) 
+    {
+      LOG(AliL3Log::kError,"AliL3Evaluate::InitMC","gAlice")
+       <<"AliRun object non existing on file"<<ENDLOG;
+      return;
+    }
+
+  AliTPCParam *param = (AliTPCParam*)exfile->Get(AliL3Transform::GetParamName());
   
-  for (i=0; i<num_of_clusters; i++) 
+  TFile *exact = TFile::Open(exactfile);
+  if(!exact)
     {
-      id = hitnum[i];
-      Int_t slice = (id>>25) & 0x7f;
-      Int_t patch = (id>>22) & 0x7;
-      UInt_t pos = id&0x3fffff;              
-      
-      AliL3SpacePointData *points = fClusters[slice][patch];
+      cerr<<"AliL3Evaluate::EvaluatePoints : Problems opening file :"<<exactfile<<endl;
+      return;
+    }
+  
+  AliStack *astack=gAlice->Stack();
+
+  AliTPCClustersArray *arr=0;
+  for(Int_t event=0; event<nevent; event++)
+    {
+      LoadData(event,sp);   
+      exfile->cd();
+      if(arr)
+       delete arr;
+      Int_t nparticles = gAlice->GetEvent(event);
+      Int_t nprimaries = 0;//FindPrimaries(nparticles);
+      cout<<"Event "<<event<<" had "<<nparticles<<" particles and "<<nprimaries<<" primaries"<<endl;
+      exact->cd();
       
-      if(!points) 
-       continue;
+      //Get the exact clusters from file:
+      AliTPCClustersArray *arr = new AliTPCClustersArray;
+      arr->Setup(param);
+      arr->SetClusterType("AliComplexCluster");
+      char treeName[500];
+      sprintf(treeName,"TreeCExact_%s_%d",param->GetTitle(),event);
+      Bool_t clusterok = arr->ConnectTree(treeName);//Segment Tree (for offline clusters)
+      if(!clusterok) {printf("AliL3Evaluate::EvaluatePoints : Error in clusterloading\n"); return;}
       
-      if(pos>=fNcl[slice][patch]) 
+      //cout<<"Entering loop with "<<(Int_t)arr->GetTree()->GetEntries()<<endl;
+      for(Int_t i=0; i<arr->GetTree()->GetEntries(); i++)
        {
-         LOG(AliL3Log::kError,"AliL3Evaluate::GetMCTrackLabel","Clusterarray")
-           <<AliL3Log::kDec<<"ERROR"<<ENDLOG;
-         continue;
+         //Get the exact clusters for this row:
+         Int_t cursec,currow;
+         AliSegmentID *s = arr->LoadEntry(i);
+         param->AdjustSectorRow(s->GetID(),cursec,currow);
+         
+         AliTPCClustersRow *ro = (AliTPCClustersRow *)arr->GetRow(cursec,currow);
+         TClonesArray *clusters = ro->GetArray();
+         int num_of_offline=clusters->GetEntriesFast();
+         
+         //Get the found clusters:
+         Int_t slice,padrow;
+         AliL3Transform::Sector2Slice(slice,padrow,cursec,currow);
+         if(slice < fMinSlice) continue;
+         if(slice > fMaxSlice) break;
+         
+         Int_t patch = AliL3Transform::GetPatch(padrow);
+         if(sp)
+           patch=0;
+         AliL3SpacePointData *points = fClusters[slice][patch];
+         if(!points)
+           continue;
+         
+         //cout<<"Slice "<<slice<<" padrow "<<padrow<<" has "<<num_of_offline<<" clusters "<<endl;
+         Int_t clustercount=0;
+         Int_t crosscount=0;
+         for(Int_t m=0; m<num_of_offline; m++)
+           {
+             AliComplexCluster *cluster = (AliComplexCluster *)clusters->UncheckedAt(m);
+#ifdef use_newio
+             Int_t mcId = cluster->GetTrack(0);
+#else
+             Int_t mcId = cluster->fTracks[0];
+#endif       
+             if(mcId <0) continue;
+       
+#ifdef use_newio      
+             if(cluster->GetY() < 1 || cluster->GetY() > AliL3Transform::GetNPads(padrow) - 2 ||
+                cluster->GetX() < 1 || cluster->GetX() > AliL3Transform::GetNTimeBins() - 2)
+               continue;
+#else
+             if(cluster->fY < 1 || cluster->fY > AliL3Transform::GetNPads(padrow) - 2 ||
+                cluster->fX < 1 || cluster->fX > AliL3Transform::GetNTimeBins() - 2)
+               continue;
+#endif       
+             Float_t xyz_ex[3];
+             
+#ifdef use_newio
+             AliL3Transform::Raw2Local(xyz_ex,cursec,currow,cluster->GetY(),cluster->GetX());
+#else        
+             AliL3Transform::Raw2Local(xyz_ex,cursec,currow,cluster->fY,cluster->fX);
+#endif       
+             //In function AliTPC::Hits2ExactClusters the time offset is not included,
+             //so we have to substract it again here.
+             if(slice<18)
+               xyz_ex[2]-=AliL3Transform::GetZOffset();
+             else
+               xyz_ex[2]+=AliL3Transform::GetZOffset();
+             
+             //Outside our cone:
+             if(param->GetPadRowRadii(cursec,currow)<230./250.*fabs(xyz_ex[2]))
+               continue;
+             
+             TParticle *part = astack->Particle(mcId);
+             crosscount++;
+             
+             if(part->Pt() < fMinGoodPt) continue;
+             
+             //Dont take secondaries, because in width calculation we assume primaries:
+             //if(!(part->TestBit(kPrimaryCharged))) continue;
+             if(part->GetFirstMother()>=0) continue;
+             
+             Int_t tempcount=0;
+             for(UInt_t c=0; c<fNcl[slice][patch]; c++)
+               {
+                 if((Int_t)points[c].fPadRow!=padrow) continue;
+                 Float_t xyz_cl[3] = {points[c].fX,points[c].fY,points[c].fZ};
+                 
+                 if(!offline)
+                   AliL3Transform::Global2Local(xyz_cl,cursec);
+                 tempcount++;
+                 
+                 if(points[c].fTrackID[0] != mcId &&
+                    points[c].fTrackID[1] != mcId &&
+                    points[c].fTrackID[2] != mcId)
+                   continue;
+                 
+                 //Residuals:
+                 Float_t resy = xyz_cl[1] - xyz_ex[1];
+                 Float_t resz = xyz_cl[2] - xyz_ex[2];
+                 
+                 //Cluster shape
+                 Int_t charge = (Int_t)points[c].fCharge;
+                 Float_t beta = GetCrossingAngle(part,slice,padrow,xyz_ex);
+                 Double_t tanl = xyz_ex[2]/sqrt(xyz_ex[0]*xyz_ex[0]+xyz_ex[1]*xyz_ex[1]);
+                 Float_t psigmaY2 = AliL3Transform::GetParSigmaY2(padrow,xyz_ex[2],beta);
+                 Float_t psigmaZ2 = AliL3Transform::GetParSigmaZ2(padrow,xyz_ex[2],tanl);
+                 Float_t sigmaY2 = points[c].fSigmaY2;
+                 Float_t sigmaZ2 = points[c].fSigmaZ2;
+                 ntuppel->Fill(slice,padrow,charge,resy,resz,xyz_ex[2],part->Pt(),beta,sigmaY2,sigmaZ2,psigmaY2,psigmaZ2);
+               }
+             clustercount=tempcount;
+           }
+         ntuppel2->Fill(slice,padrow,clustercount,crosscount);
+         arr->ClearRow(cursec,currow);
        }
-      
-      if (abs(trackID[i][1]) == lab || 
-         abs(trackID[i][2]) == lab ) max++;
-    }
-  
-  //check if more than 10% of the clusters are incorrectly assigned (fake track):
-  if (1.-Float_t(max)/num_of_clusters > 0.10) 
-    {
-      return -lab;
     }
+  exfile->Close();
+  exact->Close();
+
+  TFile *ofile = TFile::Open(tofile,"RECREATE");
+  ntuppel->Write();
+  ntuppel2->Write();
+  ofile->Close();
   
-  delete [] trackID;
-  return lab;
+#endif
 }
 
-
-Int_t **AliL3Evaluate::GetClusterIDs(AliL3Track *track,UInt_t *index,Int_t *pID,Int_t npoints)
+void AliL3Evaluate::GetCFeff(Char_t *path,Char_t *outfile,Int_t nevent,Bool_t sp)
 {
-  //Return the MC information of all clusters belonging to track.
+  //Evaluate the cluster finder efficiency.
   
-  Int_t num_of_clusters = track->GetNumberOfPoints();
-  Int_t **trackID = new Int_t*[num_of_clusters];
+#ifndef do_mc
+  cerr<<"AliL3Evaluate::GetCFeff : Compile with do_mc flag"<<endl;
+  return;
+#else
+  TNtuple *ntuppel = new TNtuple("ntuppel","Cluster finder efficiency","slice:row:ncrossings:nclusters");
+  ntuppel->SetDirectory(0);
   
-  UInt_t *hitnum = track->GetHitNumbers();  
-  UInt_t id;
+  Char_t filename[1024];
+  sprintf(filename,"%s/alirunfile.root",path);
+  TFile *rfile = TFile::Open(filename);
+  gAlice = (AliRun*)rfile->Get("gAlice");
 
+  AliStack *astack=gAlice->Stack();
   
-  Float_t xyz[3];
-  Int_t sector,padrow;
-  for(Int_t i=0; i<num_of_clusters; i++)
-    {
-      id = hitnum[i];
-      Int_t slice = (id>>25) & 0x7f;
-      Int_t patch = (id>>22) & 0x7;
-      UInt_t pos = id&0x3fffff;              
+  AliTPCParam *param = (AliTPCParam*)rfile->Get(AliL3Transform::GetParamName());
       
-      AliL3SpacePointData *points = fClusters[slice][patch];
+  Int_t zero=param->GetZeroSup();
+
+  sprintf(filename,"%s/digitfile.root",path);
+  TFile *dfile = TFile::Open(filename);
+  
+  for(Int_t event=0; event<nevent; event++)
+    {
+      LoadData(event,sp);
+      rfile->cd();
+      gAlice->GetEvent(event);
+      Int_t np = astack->GetNtrack();
+      cout<<"Processing event "<<event<<" with "<<np<<" particles "<<endl;
+      dfile->cd();
+      sprintf(filename,"TreeD_75x40_100x60_150x60_%d",event);
+      TTree *TD=(TTree*)gDirectory->Get(filename);
+      AliSimDigits da, *digits=&da;
+      TD->GetBranch("Segment")->SetAddress(&digits);
       
-      if(!points) 
-       {
-         LOG(AliL3Log::kError,"AliL3Evaluate::GetClusterIDs","Clusterarray")
-           <<"No points at slice "<<slice<<" patch "<<patch<<" pos "<<pos<<ENDLOG;
-         continue;
-       }
-      if(pos>=fNcl[slice][patch]) 
-       {
-         LOG(AliL3Log::kError,"AliL3Evaluate::GetClusterIDs","Clusterarray")
-           <<AliL3Log::kDec<<"ERROR"<<ENDLOG;
-         continue;
-       }
+      Int_t crossed=0,recs=0;
+      Int_t *count = new Int_t[np]; //np number of particles.
+      Int_t i;
+      Float_t xyz[3];
+      for (i=0; i<np; i++) count[i]=0;
       
-      xyz[0] = points[pos].fX;
-      xyz[1] = points[pos].fY;
-      xyz[2] = points[pos].fZ;
-      //sector = points[pos].fSector;
-      padrow = points[pos].fPadRow;
-      Int_t se,ro;
-      fTransform->Slice2Sector(slice,padrow,se,ro);
-      fTransform->Global2Raw(xyz,se,ro);
       
-      if(fIsSlow)
-       {
-         Int_t p = fRowid[slice][padrow];
-         
-         if(!fDigitsTree->GetEvent(p)) 
-           LOG(AliL3Log::kError,"AliL3Evaluate::GetClusterIDs","Digits Tree")
-             <<"Error reading digits tree"<<ENDLOG;
-         
-         trackID[i] = new Int_t[3];
-         trackID[i][0] = fDigits->GetTrackID((Int_t)xyz[2],(Int_t)xyz[1],0);
-         trackID[i][1] = fDigits->GetTrackID((Int_t)xyz[2],(Int_t)xyz[1],1);
-         trackID[i][2] = fDigits->GetTrackID((Int_t)xyz[2],(Int_t)xyz[1],2);
-       }
-      else
+      Int_t sec,row,sl,sr;
+      for(Int_t i=0; i<(Int_t)TD->GetEntries(); i++)
        {
-         Int_t tmp_pid=0;
-         for(Int_t ii=0; ii<npoints; ii++)
+         crossed=recs=0;
+         if (!TD->GetEvent(i)) continue;
+         param->AdjustSectorRow(digits->GetID(),sec,row);
+         AliL3Transform::Sector2Slice(sl,sr,sec,row);
+         if(sl < fMinSlice) continue;
+         if(sl > fMaxSlice) break;
+         cout<<"Processing slice "<<sl<<" row "<<sr<<endl;
+         digits->First();
+         do {
+           Int_t it=digits->CurrentRow(), ip=digits->CurrentColumn();
+           Short_t dig = digits->GetDigit(it,ip);
+           
+           if(dig<=param->GetZeroSup()) continue;
+           AliL3Transform::Raw2Local(xyz,sec,row,ip,it);
+           if(param->GetPadRowRadii(sec,row)<230./250.*fabs(xyz[2]))
+             continue;
+           
+           Int_t idx0=digits->GetTrackID(it,ip,0); 
+           Int_t idx1=digits->GetTrackID(it,ip,1);
+           Int_t idx2=digits->GetTrackID(it,ip,2);
+           
+           if (idx0>=0 && dig>=zero) count[idx0]+=1;
+           if (idx1>=0 && dig>=zero) count[idx1]+=1;
+           if (idx2>=0 && dig>=zero) count[idx2]+=1;
+         } while (digits->Next());
+         for (Int_t j=0; j<np; j++) 
+           {
+             TParticle *part = astack->Particle(j);
+             if(part->Pt() < fMinGoodPt) continue;
+             if(part->GetFirstMother() >= 0) continue;
+             if (count[j]>1) //at least two digits at this padrow 
+               {
+                 crossed++;
+                 count[j]=0;
+               }
+           }
+
+         Int_t patch = AliL3Transform::GetPatch(sr);
+         if(sp==kTRUE)
+           patch=0;
+         AliL3SpacePointData *points = fClusters[sl][patch];
+         if(!points)
+           continue;
+         for(UInt_t k=0; k<fNcl[sl][patch]; k++)
            {
-             tmp_pid = pID[ii];
-             if(index[ii] == id) break;
+             if(points[k].fPadRow!=sr) continue;
+             recs++;
            }
-         trackID[i] = new Int_t[3];
-         trackID[i][0] = tmp_pid;
-         trackID[i][1] = -1;
-         trackID[i][2] = -1;
+         ntuppel->Fill(sl,sr,crossed,recs);
        }
+      
+      TD->Delete();
+      delete[] count;
     }
-  return trackID;
+  TFile *file = TFile::Open(outfile,"RECREATE");
+  ntuppel->Write();
+  file->Close();
+  
+  rfile->Close();
+  dfile->Close();
+#endif
 }
 
-Int_t *AliL3Evaluate::GetFastIDs(UInt_t &tmp_ind,Int_t &npoints)
+Float_t AliL3Evaluate::GetCrossingAngle(TParticle *part,Int_t slice,Int_t /*padrow*/,Float_t *xyz)
 {
-  //Get the MC id of space points in case of using the fast simulator. 
+  //Calculate the padrow crossing angle of the particle
+  
+  Double_t kappa = AliL3Transform::GetBField()*AliL3Transform::GetBFact()/part->Pt();
+  
+  Double_t radius = 1/fabs(kappa);
+  if(part->GetPdgCode() > 0) kappa = -kappa;
 
-  FILE *infile = fopen("point_mc.dat","r");
-  if(!infile) return 0;
-  Int_t hitid,hitmc,i;
+  Float_t angl[1] = {part->Phi()};
   
-  for(i=0; ; i++)
-    if(fscanf(infile,"%d %d",&hitid,&hitmc)==EOF) break;
-  npoints = i;
-  rewind(infile);
-  Int_t *pID = new Int_t[npoints];
-  UInt_t *ind = new UInt_t[npoints];
-  tmp_ind = (UInt_t)ind;
+  AliL3Transform::Global2LocalAngle(angl,slice);
   
-  for(i=0; i<npoints; i++)
-    {
-      if(fscanf(infile,"%d %d",&hitid,&hitmc)==EOF) break;
-      pID[i] = hitmc;
-      ind[i] = hitid;
-    }
-  fclose(infile);
+  Double_t charge = -1.*kappa;
 
-  return pID;
+  Double_t trackPhi0 = angl[0] + charge*0.5*AliL3Transform::Pi()/fabs(charge);
+  
+  Double_t x0=0;
+  Double_t y0=0;
+  Double_t xc = x0 - radius * cos(trackPhi0);
+  Double_t yc = y0 - radius * sin(trackPhi0);
 
+  Double_t tangent[2];
+  tangent[0] = -1.*(xyz[1] - yc)/radius;
+  tangent[1] = (xyz[0] - xc)/radius;
+  
+  Double_t perp_padrow[2] = {1,0}; //locally in slice
+  
+  Double_t cos_beta = fabs(tangent[0]*perp_padrow[0] + tangent[1]*perp_padrow[1]);
+  if(cos_beta > 1) cos_beta=1;
+  return acos(cos_beta);
 }
 
-void AliL3Evaluate::Write2File(Char_t *outputfile)
+Int_t AliL3Evaluate::FindPrimaries(Int_t nparticles)
 {
-  //Write histograms to file:
-  
-  TFile *of = new TFile(outputfile,"RECREATE");
-  if(!of->IsOpen())
-    {
-      LOG(AliL3Log::kError,"AliL3Evaluate::Write2File","File Open")
-       <<"Problems opening rootfile"<<ENDLOG;
-      return;
+  // cuts:
+  Double_t vertcut = 0.001;
+  Double_t decacut = 3.;
+  Double_t timecut = 0.;
+  Int_t nprch1=0;
+  AliStack *astack=gAlice->Stack();
+  TParticle * part = astack->Particle(0);
+  Double_t xori = part->Vx();
+  Double_t yori = part->Vy();
+  Double_t zori = part->Vz();
+  for(Int_t iprim = 0; iprim<nparticles; iprim++){   //loop on  tracks
+    
+    part = astack->Particle(iprim);
+    const char * xxx=strstr(part->GetName(),"XXX");
+    if(xxx)continue;
+    
+    TParticlePDG *ppdg = part->GetPDG();
+    if(TMath::Abs(ppdg->Charge())!=3)continue;  // only charged (no quarks)
+    
+    Double_t dist=TMath::Sqrt((part->Vx()-xori)*(part->Vx()-xori)+(part->Vy()-yori)*(part->Vy()-yori)+(part->Vz()-zori)*(part->Vz()-zori));
+    if(dist>vertcut)continue;  // cut on the vertex
+    
+    if(part->T()>timecut)continue;
+    
+    Double_t ptot=TMath::Sqrt(part->Px()*part->Px()+part->Py()*part->Py()+part->Pz()*part->Pz());
+    if(ptot==(TMath::Abs(part->Pz())))continue; // no beam particles
+    
+    Bool_t prmch = kTRUE;   // candidate primary track
+    Int_t fidau=part->GetFirstDaughter();  // cut on daughters
+    Int_t lasdau=0;
+    Int_t ndau=0;
+    if(fidau>=0){
+      lasdau=part->GetLastDaughter();
+      ndau=lasdau-fidau+1;
     }
-  
-  of->cd();
-  fPtRes->Write();
-  fNGoodTracksPt->Write();
-  fNFoundTracksPt->Write();
-  fNFakeTracksPt->Write();
-  fTrackEffPt->Write();
-  fFakeTrackEffPt->Write();
-  fNGoodTracksEta->Write();
-  fNFoundTracksEta->Write();
-  fNFakeTracksEta->Write();
-  fTrackEffEta->Write();
-  fFakeTrackEffEta->Write();
-  
-  of->Close();
-  delete of;
+    if(ndau>0){
+      for(Int_t j=fidau;j<=lasdau;j++){
+        TParticle *dau=astack->Particle(j);
+        Double_t distd=TMath::Sqrt((dau->Vx()-xori)*(dau->Vx()-xori)+(dau->Vy()-yori)*(dau->Vy()-yori)+(dau->Vz()-zori)*(dau->Vz()-zori));
+        if(distd<decacut)prmch=kFALSE;  // eliminate if the decay is near the vertex
+      }
+    }
+    
+    if(prmch){
+      nprch1++;
+      part->SetBit(kPrimaryCharged);
+    }
+  }
 
+  return nprch1;
 }
-