HLT/src/AliL3Evaluate.cxx

   1 //$Id$
   2
   3 // Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
   4 //*-- Copyright &copy ASV
   5
   6 #include "AliL3StandardIncludes.h"
   7 #include <TFile.h>
   8 #include <TH1.h>
   9 #include <TParticle.h>
  10 #include <TTree.h>
  11 #include <TClonesArray.h>
  12
  13 #include <AliRun.h>
  14 #include <AliSimDigits.h>
  15 #include <AliTPC.h>
  16 #include <AliTPCcluster.h>
  17 #include <AliTPCClustersArray.h>
  18 #include <AliTPCClustersRow.h>
  19 #include <AliTPCParam.h>
  20 #include <AliComplexCluster.h>
  21
  22 #if GCCVERSION == 3
  23 #include <fstream>
  24 #include <iosfwd>
  25 #else
  26 #include <fstream.h>
  27 #endif
  28
  29 #include "AliL3Logging.h"
  30 #include "AliL3Transform.h"
  31 #include "AliL3SpacePointData.h"
  32 #include "AliL3Track.h"
  33 #include "AliL3FileHandler.h"
  34 #include "AliL3TrackArray.h"
  35 #include "AliL3Evaluate.h"
  36
  37 #if GCCVERSION == 3
  38 using namespace std;
  39 #endif
  40
  41 //_____________________________________________________________
  42 // AliL3Evaluate
  43 //
  44 // Evaluation class for tracking. Plots efficiencies etc..
  45 //
  46
  47
  48 ClassImp(AliL3Evaluate)
  49
  50 AliL3Evaluate::AliL3Evaluate()
  51 {
  52   fTracks = 0;
  53   fNFastPoints = 0;
  54   fMcIndex = 0;
  55   fMcId = 0;
  56   fMinSlice=0;
  57   fMaxSlice=0;
  58   fGoodTracks=0;
  59 }
  60
  61 AliL3Evaluate::AliL3Evaluate(Char_t *datapath,Int_t min_clusters,Int_t minhits,Double_t minpt,Int_t *slice)
  62 {
  63
  64   if(slice)
  65     {
  66       fMinSlice=slice[0];
  67       fMaxSlice=slice[1];
  68     }
  69   else
  70     {
  71       fMinSlice=0;
  72       fMaxSlice=35;
  73     }
  74   fMaxFalseClusters = 0.1;
  75   fGoodFound = 0;
  76   fGoodGen = 0;
  77   fMinPointsOnTrack = min_clusters;
  78   fMinHitsFromParticle = minhits;
  79   fMinGoodPt = minpt;
  80   Char_t fname[1024];
  81   AliL3FileHandler *clusterfile[36][6];
  82   for(Int_t s=fMinSlice; s<=fMaxSlice; s++)
  83     {
  84       for(Int_t p=0; p<AliL3Transform::GetNPatches(); p++)
  85         {
  86           fClusters[s][p] = 0;
  87           clusterfile[s][p] = new AliL3FileHandler();
  88           sprintf(fname,"%s/points_%d_%d.raw",datapath,s,p);
  89           if(!clusterfile[s][p]->SetBinaryInput(fname))
  90             {
  91               LOG(AliL3Log::kError,"AliL3Evaluation::Setup","File Open")
  92                 <<"Inputfile "<<fname<<" does not exist"<<ENDLOG;
  93               delete clusterfile[s][p];
  94               clusterfile[s][p] = 0;
  95               continue;
  96             }
  97           fClusters[s][p] = (AliL3SpacePointData*)clusterfile[s][p]->Allocate();
  98           clusterfile[s][p]->Binary2Memory(fNcl[s][p],fClusters[s][p]);
  99           clusterfile[s][p]->CloseBinaryInput();
 100         }
 101     }
 102
 103   sprintf(fname,"%s/tracks.raw",datapath);
 104   AliL3FileHandler *tfile = new AliL3FileHandler();
 105   if(!tfile->SetBinaryInput(fname)){
 106     LOG(AliL3Log::kError,"AliL3Evaluation::Setup","File Open")
 107       <<"Inputfile "<<fname<<" does not exist"<<ENDLOG;
 108     return;
 109   }
 110   fTracks = new AliL3TrackArray();
 111   tfile->Binary2TrackArray(fTracks);
 112   tfile->CloseBinaryInput();
 113   delete tfile;
 114 }
 115
 116
 117 AliL3Evaluate::~AliL3Evaluate()
 118 {
 119   if(fGoodTracks) delete fGoodTracks;
 120   if(fDigitsTree) fDigitsTree->Delete();
 121   if(fTracks) delete fTracks;
 122   if(fPtRes) delete fPtRes;
 123   if(fNGoodTracksPt) delete fNGoodTracksPt;
 124   if(fNFoundTracksPt) delete fNFoundTracksPt;
 125   if(fNFakeTracksPt) delete fNFakeTracksPt;
 126   if(fTrackEffPt) delete fTrackEffPt;
 127   if(fFakeTrackEffPt) delete fFakeTrackEffPt;
 128   if(fNGoodTracksEta) delete fNGoodTracksEta;
 129   if(fNFoundTracksEta) delete fNFoundTracksEta;
 130   if(fNFakeTracksEta) delete fNFakeTracksEta;
 131   if(fTrackEffEta) delete fTrackEffEta;
 132   if(fFakeTrackEffEta) delete fFakeTrackEffEta;
 133   if(fMcIndex) delete [] fMcIndex;
 134   if(fMcId)    delete [] fMcId;
 135   if(fNtuppel) delete fNtuppel;
 136 }
 137
 138
 139 void AliL3Evaluate::AssignIDs()
 140 {
 141   //Assign MC id to the tracks.
 142 #ifndef do_mc
 143   cerr<<"AliL3Evaluate::AssignIDs() : You need to compile with the do_mc flag!"<<endl;
 144   return;
 145 #endif
 146
 147   fTracks->QSort();
 148   LOG(AliL3Log::kDebug,"AliL3Evaluate::AssignIDs","Track Loop")
 149     <<"Assigning MC id to the found tracks...."<<ENDLOG;
 150   for(Int_t i=0; i<fTracks->GetNTracks(); i++)
 151     {
 152       AliL3Track *track = (AliL3Track*)fTracks->GetCheckedTrack(i);
 153       if(!track) continue;
 154       if(track->GetNumberOfPoints() < fMinPointsOnTrack) break;
 155
 156       fGoodFound++;
 157       Int_t tID = GetMCTrackLabel(track);
 158       track->SetMCid(tID);
 159     }
 160   cout<<"Found "<<fGoodFound<<" good tracks "<<endl;
 161 }
 162
 163
 164 struct S {Int_t lab; Int_t max;};
 165 Int_t AliL3Evaluate::GetMCTrackLabel(AliL3Track *track){
 166   //Returns the MCtrackID of the belonging clusters.
 167   //If MCLabel < 0, means that track is fake.
 168   //Fake track means that more than 10 percent of clusters are assigned incorrectly.
 169
 170 #ifdef do_mc
 171   Int_t num_of_clusters = track->GetNumberOfPoints();
 172   S *s=new S[num_of_clusters];
 173   Int_t i;
 174   for (i=0; i<num_of_clusters; i++) s[i].lab=s[i].max=0;
 175   UInt_t *hitnum = track->GetHitNumbers();
 176   UInt_t id;
 177
 178   Int_t lab=123456789;
 179   for (i=0; i<num_of_clusters; i++)
 180     {
 181       //Tricks to get the clusters belonging to this track:
 182       id = hitnum[i];
 183       Int_t slice = (id>>25) & 0x7f;
 184       Int_t patch = (id>>22) & 0x7;
 185       UInt_t pos = id&0x3fffff;
 186
 187       AliL3SpacePointData *points = fClusters[slice][patch];
 188       if(!points) continue;
 189       if(pos>=fNcl[slice][patch])
 190         {
 191           LOG(AliL3Log::kError,"AliL3Evaluate::GetMCTrackLabel","Clusterarray")
 192             <<AliL3Log::kDec<<"ERROR"<<ENDLOG;
 193           continue;
 194         }
 195
 196       //Get the label of the cluster:
 197
 198       lab=abs(points[pos].fTrackID[0]);
 199       if(lab < 0)
 200         cout<<"Track had negative id : "<<lab<<" padrow "<<(Int_t)points[pos].fPadRow<<" nhits "<<num_of_clusters<<" pt "<<track->GetPt()<<endl;
 201
 202       Int_t j;
 203       for (j=0; j<num_of_clusters; j++)
 204         if (s[j].lab==lab || s[j].max==0) break;
 205       s[j].lab=lab;
 206       s[j].max++;
 207     }
 208
 209   Int_t max=0;
 210   for (i=0; i<num_of_clusters; i++)
 211     if (s[i].max>max) {max=s[i].max; lab=s[i].lab;}
 212
 213   delete[] s;
 214
 215   for (i=0; i<num_of_clusters; i++)
 216     {
 217       id = hitnum[i];
 218       Int_t slice = (id>>25) & 0x7f;
 219       Int_t patch = (id>>22) & 0x7;
 220       UInt_t pos = id&0x3fffff;
 221
 222       AliL3SpacePointData *points = fClusters[slice][patch];
 223       if(!points) continue;
 224       if(pos>=fNcl[slice][patch])
 225         {
 226           LOG(AliL3Log::kError,"AliL3Evaluate::GetMCTrackLabel","Clusterarray")
 227             <<AliL3Log::kDec<<"ERROR"<<ENDLOG;
 228           continue;
 229         }
 230
 231       if (abs(points[pos].fTrackID[1]) == lab ||
 232           abs(points[pos].fTrackID[2]) == lab ) max++;
 233     }
 234
 235   if (1.-Float_t(max)/num_of_clusters > fMaxFalseClusters)
 236     {
 237       return -lab;
 238     }
 239
 240   return lab;
 241 #else //If we are running with mc_ids or not
 242   return 0;
 243 #endif
 244
 245 }
 246
 247 void AliL3Evaluate::GetFastClusterIDs(Char_t *path)
 248 {
 249   //Get the MC id of space points in case of using the fast simulator.
 250   char fname[256];
 251   sprintf(fname,"%s/point_mc.dat",path);
 252   FILE *infile = fopen(fname,"r");
 253   if(!infile) return;
 254   Int_t hitid,hitmc,i;
 255
 256   for(i=0; ; i++)
 257     if(fscanf(infile,"%d %d",&hitid,&hitmc)==EOF) break;
 258   rewind(infile);
 259   fNFastPoints = i;
 260   fMcId = new Int_t[fNFastPoints];
 261   fMcIndex = new UInt_t[fNFastPoints];
 262
 263   for(i=0; i<fNFastPoints; i++)
 264     {
 265       if(fscanf(infile,"%d %d",&hitid,&hitmc)==EOF) break;
 266       fMcId[i] = hitmc;
 267       fMcIndex[i] = hitid;
 268     }
 269   fclose(infile);
 270 }
 271
 272 void AliL3Evaluate::CreateHistos(Int_t nbin,Float_t xlow,Float_t xup)
 273 {
 274   //Create the histograms
 275
 276   LOG(AliL3Log::kInformational,"AliL3Evaluate::CreateHistos","Allocating")
 277     <<"Creating histograms..."<<ENDLOG;
 278
 279   fNtuppel = new TNtuple("fNtuppel","Pt resolution","pt_gen:pt_found:nHits");
 280
 281   fPtRes = new TH1F("fPtRes","Relative Pt resolution",30,-10.,10.);
 282   fNGoodTracksPt = new TH1F("fNGoodTracksPt","Good tracks vs pt",nbin,xlow,xup);
 283   fNFoundTracksPt = new TH1F("fNFoundTracksPt","Found tracks vs pt",nbin,xlow,xup);
 284   fNFakeTracksPt = new TH1F("fNFakeTracksPt","Fake tracks vs pt",nbin,xlow,xup);
 285   fTrackEffPt = new TH1F("fTrackEffPt","Tracking efficiency vs pt",nbin,xlow,xup);
 286   fFakeTrackEffPt = new TH1F("fFakeTrackEffPt","Efficiency for fake tracks vs pt",nbin,xlow,xup);
 287
 288   fNGoodTracksEta = new TH1F("fNGoodTracksEta","Good tracks vs eta",20,-50,50);
 289   fNFoundTracksEta = new TH1F("fNFoundTracksEta","Found tracks vs eta",20,-50,50);
 290   fNFakeTracksEta = new TH1F("fNFakeTracksEta","Fake tracks vs eta",20,-50,50);
 291   fTrackEffEta = new TH1F("fTrackEffEta","Tracking efficienct vs eta",20,-50,50);
 292   fFakeTrackEffEta = new TH1F("fFakeTrackEffEta","Efficiency for fake tracks vs eta",20,-50,50);
 293
 294 }
 295
 296 void AliL3Evaluate::GetGoodParticles(Char_t *path)
 297 {
 298   //Read the good particles from file. This file should already have been
 299   //generated by macro AliTPCComparison.C.
 300
 301   Char_t filename[1024];
 302   sprintf(filename,"%s/good_tracks_tpc",path);
 303   ifstream in(filename);
 304   if(!in)
 305     {
 306       cerr<<"AliL3Evaluate::GetGoodParticles : Problems opening file :"<<filename<<endl;
 307       return;
 308     }
 309   Int_t MaxTracks=20000;
 310   fGoodTracks = new GoodTrack[MaxTracks];
 311
 312   while (in>>fGoodTracks[fGoodGen].label>>fGoodTracks[fGoodGen].code>>
 313          fGoodTracks[fGoodGen].px>>fGoodTracks[fGoodGen].py>>fGoodTracks[fGoodGen].pz>>
 314          fGoodTracks[fGoodGen].x>>fGoodTracks[fGoodGen].y >>fGoodTracks[fGoodGen].z>>fGoodTracks[fGoodGen].nhits)
 315
 316     {
 317       fGoodGen++;
 318       if (fGoodGen==MaxTracks)
 319         {
 320           cerr<<"AliL3Evaluate::GetGoodParticles : Too many good tracks !\n";
 321           break;
 322         }
 323     }
 324
 325 }
 326
 327 void AliL3Evaluate::FillEffHistos()
 328 {
 329   if(!fGoodTracks)
 330     {
 331       cerr<<"AliL3Evaluate::FillEffHistos : No good tracks"<<endl;
 332       return;
 333     }
 334   cout<<"Comparing..."<<endl;
 335   for(Int_t i=0; i<fGoodGen; i++)
 336     {
 337       //cout<<"Checking particle "<<i<<endl;
 338       if(fGoodTracks[i].nhits < fMinHitsFromParticle) continue;
 339       Float_t ptg = TMath::Sqrt(fGoodTracks[i].px*fGoodTracks[i].px + fGoodTracks[i].py*fGoodTracks[i].py);
 340       Float_t pzg=fGoodTracks[i].pz;
 341       Float_t dipangle=TMath::ATan2(pzg,ptg)*180./TMath::Pi();
 342
 343       fNGoodTracksPt->Fill(ptg);
 344       fNGoodTracksEta->Fill(dipangle);
 345       Int_t found = 0;
 346
 347       for(Int_t k=0; k<fTracks->GetNTracks(); k++)
 348         {
 349           AliL3Track *track = fTracks->GetCheckedTrack(k);
 350           if(!track) continue;
 351           Int_t nHits = track->GetNumberOfPoints();
 352           if(nHits < fMinPointsOnTrack) break;
 353           Int_t tracklabel;
 354           tracklabel = track->GetMCid();
 355
 356           if(TMath::Abs(tracklabel) != fGoodTracks[i].label) continue;
 357           found=1;
 358           if(tracklabel == fGoodTracks[i].label) {fNFoundTracksPt->Fill(ptg); fNFoundTracksEta->Fill(dipangle);}
 359           else {fNFakeTracksPt->Fill(ptg); fNFakeTracksEta->Fill(dipangle);}
 360           Float_t pt=track->GetPt();
 361           fPtRes->Fill((pt-ptg)/ptg*100.);
 362           //fNtuppel->Fill(ptg,pt,nHits);
 363           break;
 364
 365         }
 366     }
 367 }
 368
 369 void AliL3Evaluate::FillEffHistosNAIVE()
 370 {
 371   //Fill the efficiency histograms.
 372
 373   cout<<endl<<"Note: Doing NAIVE evaluation "<<endl;
 374   for(Int_t i=0; i<fGoodGen; i++)
 375     {
 376       Double_t ptg=TMath::Sqrt(fGoodTracks[i].px*fGoodTracks[i].px + fGoodTracks[i].py*fGoodTracks[i].py);
 377       Double_t pzg=fGoodTracks[i].pz;
 378       Float_t dipangle=TMath::ATan2(pzg,ptg)*180./TMath::Pi();
 379       //printf("filling particle with pt %f and dipangle %f\n",ptg,dipangle);
 380       fNGoodTracksPt->Fill(ptg);
 381       fNGoodTracksEta->Fill(dipangle);
 382
 383     }
 384
 385   for(Int_t k=0; k<fTracks->GetNTracks(); k++)
 386     {
 387       AliL3Track *track = fTracks->GetCheckedTrack(k);
 388       if(!track) continue;
 389       Int_t nHits = track->GetNumberOfPoints();
 390       if(nHits < fMinPointsOnTrack) break;
 391       if(track->GetPt()<fMinGoodPt) continue;
 392       if(fabs(track->GetPseudoRapidity())>0.9) continue;
 393
 394       fNFoundTracksPt->Fill(track->GetPt()); fNFoundTracksEta->Fill(track->GetPseudoRapidity());
 395       //Float_t pt=track->GetPt();
 396       //fPtRes->Fill((pt-ptg)/ptg*100.);
 397       //fNtuppel->Fill(ptg,pt,nHits);
 398
 399     }
 400 }
 401
 402 void AliL3Evaluate::CalcEffHistos(){
 403
 404   Stat_t ngood=fNGoodTracksPt->GetEntries();
 405   Stat_t nfound=fNFoundTracksPt->GetEntries();
 406   Stat_t nfake=fNFakeTracksPt->GetEntries();
 407
 408   LOG(AliL3Log::kInformational,"AliL3Evaluate::FillEffHistos","Efficiency")
 409     <<AliL3Log::kDec<<"There was "<<ngood<<" generated good tracks"<<ENDLOG;
 410   LOG(AliL3Log::kInformational,"AliL3Evaluate::FillEffHistos","Efficiency")
 411     <<AliL3Log::kDec<<"Found "<<nfound<<" tracks"<<ENDLOG;
 412   LOG(AliL3Log::kInformational,"AliL3Evaluate::FillEffHistos","Efficiency")
 413     <<AliL3Log::kDec<<"Integral efficiency is about "<<nfound/ngood*100<<ENDLOG;
 414   LOG(AliL3Log::kInformational,"AliL3Evaluate::FillEffHistos","Efficiency")
 415     <<AliL3Log::kDec<<"Fake tracks relative is about "<<nfake/ngood*100<<ENDLOG;
 416
 417   LOG(AliL3Log::kInformational,"AliL3Evaluate::FillEffHistos","Efficiency")
 418     <<AliL3Log::kDec<<"Naive efficiency "<<(Double_t)fGoodFound/(Double_t)fGoodGen<<ENDLOG;
 419
 420
 421   fNFoundTracksPt->Sumw2(); fNGoodTracksPt->Sumw2();
 422   fTrackEffPt->Divide(fNFoundTracksPt,fNGoodTracksPt,1,1.,"b");
 423   fFakeTrackEffPt->Divide(fNFakeTracksPt,fNGoodTracksPt,1,1.,"b");
 424   fTrackEffPt->SetMaximum(1.4);
 425   fTrackEffPt->SetXTitle("P_{T} [GeV]");
 426   fTrackEffPt->SetLineWidth(2);
 427   fFakeTrackEffPt->SetFillStyle(3013);
 428   fTrackEffPt->SetLineColor(4);
 429   fFakeTrackEffPt->SetFillColor(2);
 430
 431   fNFoundTracksEta->Sumw2(); fNGoodTracksEta->Sumw2();
 432   fTrackEffEta->Divide(fNFoundTracksEta,fNGoodTracksEta,1,1.,"b");
 433   fFakeTrackEffEta->Divide(fNFakeTracksEta,fNGoodTracksEta,1,1.,"b");
 434   fTrackEffEta->SetMaximum(1.4);
 435   fTrackEffEta->SetXTitle("#lambda [degrees]");
 436   fTrackEffEta->SetLineWidth(2);
 437   fFakeTrackEffEta->SetFillStyle(3013);
 438   fTrackEffEta->SetLineColor(4);
 439   fFakeTrackEffEta->SetFillColor(2);
 440
 441 }
 442
 443 void AliL3Evaluate::Write2File(Char_t *outputfile)
 444 {
 445   //Write histograms to file:
 446
 447   TFile *of = TFile::Open(outputfile,"RECREATE");
 448   if(!of->IsOpen())
 449     {
 450       LOG(AliL3Log::kError,"AliL3Evaluate::Write2File","File Open")
 451         <<"Problems opening rootfile"<<ENDLOG;
 452       return;
 453     }
 454
 455   of->cd();
 456   fNtuppel->Write();
 457   fPtRes->Write();
 458   fNGoodTracksPt->Write();
 459   fNFoundTracksPt->Write();
 460   fNFakeTracksPt->Write();
 461   fTrackEffPt->Write();
 462   fFakeTrackEffPt->Write();
 463   fNGoodTracksEta->Write();
 464   fNFoundTracksEta->Write();
 465   fNFakeTracksEta->Write();
 466   fTrackEffEta->Write();
 467   fFakeTrackEffEta->Write();
 468
 469   of->Close();
 470
 471 }
 472
 473 TNtuple *AliL3Evaluate::CalculateResiduals()
 474 {
 475
 476   TNtuple *ntuppel=new TNtuple("ntuppel","Residuals","residual_trans:residual_long:zHit:pt:dipangle:beta:padrow:nHits");
 477
 478   for(int f=fMinSlice; f<=fMaxSlice; f++)
 479     {
 480       AliL3FileHandler *tfile = new AliL3FileHandler();
 481       char fname[256];
 482       sprintf(fname,"tracks_tr_%d_0.raw",f);
 483       if(!tfile->SetBinaryInput(fname)){
 484         LOG(AliL3Log::kError,"AliL3Evaluation::Setup","File Open")
 485           <<"Inputfile "<<fname<<" does not exist"<<ENDLOG;
 486         return 0;
 487       }
 488       fTracks = new AliL3TrackArray();
 489       tfile->Binary2TrackArray(fTracks);
 490       tfile->CloseBinaryInput();
 491       delete tfile;
 492       printf("Looking in slice %d\n",f);
 493       for(Int_t i=0; i<fTracks->GetNTracks(); i++)
 494         {
 495
 496           AliL3Track *track = (AliL3Track*)fTracks->GetCheckedTrack(i);
 497           if(!track) continue;
 498
 499           track->CalculateHelix();
 500           UInt_t *hitnum = track->GetHitNumbers();
 501           UInt_t id;
 502
 503           Float_t xyz[3];
 504           Int_t padrow;
 505           for(Int_t j=0; j<track->GetNumberOfPoints()-1; j++)
 506             {
 507               id = hitnum[j];
 508               Int_t slice = (id>>25) & 0x7f;
 509               Int_t patch = (id>>22) & 0x7;
 510               UInt_t pos = id&0x3fffff;
 511
 512               //if(slice!=1) continue;
 513
 514               AliL3SpacePointData *points = fClusters[slice][patch];
 515
 516               if(!points)
 517                 {
 518                   LOG(AliL3Log::kError,"AliL3Evaluate::CalculateResiduals","Clusterarray")
 519                     <<"No points at slice "<<slice<<" patch "<<patch<<" pos "<<pos<<ENDLOG;
 520                   continue;
 521                 }
 522               if(pos>=fNcl[slice][patch])
 523                 {
 524                   LOG(AliL3Log::kError,"AliL3Evaluate::CalculateResiduals","Clusterarray")
 525                     <<AliL3Log::kDec<<"ERROR"<<ENDLOG;
 526                   continue;
 527                 }
 528
 529               xyz[0] = points[pos].fX;
 530               xyz[1] = points[pos].fY;
 531               xyz[2] = points[pos].fZ;
 532               padrow = points[pos].fPadRow;
 533               //AliL3Transform::Global2Local(xyz,slice);
 534
 535               Float_t xyz_cross[3];
 536               track->GetCrossingPoint(padrow,xyz_cross);
 537               Double_t beta = track->GetCrossingAngle(padrow);
 538
 539               Double_t yres = xyz_cross[1] - xyz[1];
 540               Double_t zres = xyz_cross[2] - xyz[2];
 541
 542               Double_t dipangle = atan(track->GetTgl());
 543               ntuppel->Fill(yres,zres,xyz_cross[2],track->GetPt(),dipangle,beta,padrow,track->GetNumberOfPoints());
 544
 545             }
 546         }
 547       if(fTracks)
 548         delete fTracks;
 549     }
 550   return ntuppel;
 551 }
 552
 553 TNtuple *AliL3Evaluate::EvaluatePoints(Char_t *rootfile)
 554 {
 555   //Compare points to the exact crossing points of track and padrows.
 556   //The input file to this function, contains the exact clusters calculated
 557   //in AliTPC::Hits2ExactClusters.
 558
 559   cout<<"Evaluating points"<<endl;
 560   TNtuple *ntuppel = new TNtuple("ntuppel","residuals","slice:padrow:resy:resz:zHit:pt");
 561
 562   TFile *exfile = TFile::Open(rootfile);
 563   if(!exfile)
 564     {
 565       cerr<<"Error opening rootfile "<<rootfile<<endl;
 566       return 0;
 567     }
 568
 569   AliTPCParam *param = (AliTPCParam*)exfile->Get(AliL3Transform::GetParamName());
 570
 571   //Get the exact clusters from file:
 572   AliTPCClustersArray *arr = new AliTPCClustersArray;
 573   arr->Setup(param);
 574   arr->SetClusterType("AliComplexCluster");
 575   char treeName[500];
 576   sprintf(treeName,"TreeCExact_%s",param->GetTitle());
 577   Bool_t clusterok = arr->ConnectTree(treeName);//Segment Tree (for offline clusters)
 578   if(!clusterok) {printf("AliL3Evaluate::EvaluatePoints : Error in clusterloading\n"); return 0;}
 579
 580   cout<<"Entering loop with "<<(Int_t)arr->GetTree()->GetEntries()<<endl;
 581   for(Int_t i=0; i<arr->GetTree()->GetEntries(); i++)
 582     {
 583       //Get the exact clusters for this row:
 584       Int_t cursec,currow;
 585       AliSegmentID *s = arr->LoadEntry(i);
 586       param->AdjustSectorRow(s->GetID(),cursec,currow);
 587
 588       AliTPCClustersRow *ro = (AliTPCClustersRow *)arr->GetRow(cursec,currow);
 589       TClonesArray *clusters = ro->GetArray();
 590       int num_of_offline=clusters->GetEntriesFast();
 591
 592       //Get the found clusters:
 593       Int_t slice,padrow;
 594       AliL3Transform::Sector2Slice(slice,padrow,cursec,currow);
 595       if(slice<fMinSlice || slice>fMaxSlice) continue;
 596       AliL3SpacePointData *points = fClusters[slice][0];
 597
 598       Int_t index = fRowid[slice][padrow];
 599       if(!fDigitsTree->GetEvent(index))
 600         printf("AliL3Evaluate::EvaluatePoints : ERROR IN DIGITSTREE\n");
 601       printf("Checking slice %d padrow %d with %d clusters\n",slice,padrow,num_of_offline);
 602
 603       for(UInt_t c=0; c<fNcl[slice][0]; c++)
 604         {
 605           if(points[c].fPadRow!=padrow) continue;
 606           for(Int_t m=0; m<num_of_offline; m++)
 607             {
 608               AliComplexCluster *cluster = (AliComplexCluster *)clusters->UncheckedAt(m);
 609               Int_t mcId = cluster->fTracks[0];
 610               //Int_t mcId = cluster->GetLabel(0);
 611               if(mcId <0) continue;
 612               TParticle *part = gAlice->Particle(mcId);
 613
 614               Float_t xyz_cl[3] = {points[c].fX,points[c].fY,points[c].fZ};
 615               Float_t xyz_ex[3];
 616               AliL3Transform::Global2Raw(xyz_cl,cursec,currow);
 617               if(fDigits->GetTrackID((Int_t)rint(xyz_cl[2]),(Int_t)rint(xyz_cl[1]),0)!=mcId &&
 618                  fDigits->GetTrackID((Int_t)rint(xyz_cl[2]),(Int_t)rint(xyz_cl[1]),1)!=mcId &&
 619                  fDigits->GetTrackID((Int_t)rint(xyz_cl[2]),(Int_t)rint(xyz_cl[1]),2)!=mcId)
 620                 continue;
 621               AliL3Transform::Raw2Local(xyz_ex,cursec,currow,cluster->fY,cluster->fX);
 622               AliL3Transform::Raw2Local(xyz_cl,cursec,currow,xyz_cl[1],xyz_cl[2]);
 623               Float_t resy = xyz_cl[1] - xyz_ex[1];//cluster->GetY()
 624               Float_t resz = xyz_cl[2] - xyz_ex[2];//cluster->GetZ()
 625
 626               ntuppel->Fill(slice,padrow,resy,resz,xyz_ex[2],part->Pt());
 627             }
 628         }
 629       arr->ClearRow(cursec,currow);
 630     }
 631
 632   return ntuppel;
 633 }
 634
 635 void AliL3Evaluate::GetCFeff(Char_t *outfile)
 636 {
 637
 638   TNtuple *ntuppel = new TNtuple("ntuppel","Cluster finder efficiency","row:ncrossings:nclusters");
 639
 640   AliTPC *TPC = (AliTPC*)gAlice->GetDetector("TPC");
 641
 642   TPC->SetParam(fParam);
 643
 644   Int_t ver = TPC->IsVersion();
 645   LOG(AliL3Log::kInformational,"AliL3Evaluate::GetCFeff","TPC version")
 646     <<"TPC version "<<ver<<" found on file"<<ENDLOG;
 647
 648   Int_t zero=TPC->GetParam()->GetZeroSup();
 649
 650   Int_t np = gAlice->GetNtrack();
 651
 652   Int_t crossed,recs;
 653   Int_t *count = new Int_t[np]; //np number of particles.
 654   Int_t i;
 655   Float_t xyz[3];
 656   for (i=0; i<np; i++) count[i]=0;
 657   for(Int_t sl=fMinSlice; sl<=fMaxSlice; sl++)
 658     {
 659       for (i=0; i<=175; i++)
 660         {
 661           crossed=0;
 662           recs=0;
 663           Int_t index = fRowid[sl][i];
 664           if (!fDigitsTree->GetEvent(index)) continue;
 665           Int_t sec,row;
 666           fParam->AdjustSectorRow(fDigits->GetID(),sec,row);
 667           fDigits->First();
 668           do {
 669             Int_t it=fDigits->CurrentRow(), ip=fDigits->CurrentColumn();
 670             Short_t dig = fDigits->GetDigit(it,ip);
 671
 672             if(dig<=fParam->GetZeroSup()) continue;
 673             if(it < fParam->GetMaxTBin()-1 && it > 0)
 674               if(fDigits->GetDigit(it+1,ip) <= fParam->GetZeroSup()
 675                  && fDigits->GetDigit(it-1,ip) <= fParam->GetZeroSup())
 676                 continue;
 677
 678             AliL3Transform::Raw2Local(xyz,sec,row,ip,it);
 679             if(fParam->GetPadRowRadii(sec,row)<230./250.*fabs(xyz[2]))
 680               continue;
 681
 682
 683             Int_t idx0=fDigits->GetTrackID(it,ip,0);
 684             Int_t idx1=fDigits->GetTrackID(it,ip,1);
 685             Int_t idx2=fDigits->GetTrackID(it,ip,2);
 686
 687             if (idx0>=0 && dig>=zero) count[idx0]+=1;
 688             if (idx1>=0 && dig>=zero) count[idx1]+=1;
 689             if (idx2>=0 && dig>=zero) count[idx2]+=1;
 690           } while (fDigits->Next());
 691
 692           for (Int_t j=0; j<np; j++)
 693             {
 694               if (count[j]>1) //at least two digits at this padrow
 695                 {
 696                   crossed++;
 697                   count[j]=0;
 698                 }
 699             }
 700           AliL3SpacePointData *points = fClusters[sl][0];
 701           for(UInt_t k=0; k<fNcl[sl][0]; k++)
 702             {
 703               if(points[k].fPadRow!=i) continue;
 704               recs++;
 705             }
 706           ntuppel->Fill(i,crossed,recs);
 707         }
 708
 709     }
 710   delete[] count;
 711
 712   TFile *file = TFile::Open(outfile,"RECREATE");
 713   ntuppel->Write();
 714   file->Close();
 715
 716 }
 717