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