HLT/src/AliL3Evaluate.cxx

   1 //Author:        Anders Strand Vestbo
   2 //Last Modified: 5.01.2001
   3
   4 #include <stdio.h>
   5 #include <TFile.h>
   6 #include <TH1.h>
   7 #include <TParticle.h>
   8 #include <TTree.h>
   9 #include <TClonesArray.h>
  10
  11 #include "AliRun.h"
  12 #include "AliSimDigits.h"
  13 #include "AliTPC.h"
  14 #include "AliTPCClustersArray.h"
  15 #include "AliTPCClustersRow.h"
  16 #include "AliTPCcluster.h"
  17 #include "AliTPCParam.h"
  18
  19 #include "AliL3Defs.h"
  20 #include "AliL3Transform.h"
  21 #include "AliL3SpacePointData.h"
  22 #include "AliL3Track.h"
  23 #include "AliL3FileHandler.h"
  24 #include "AliL3TrackArray.h"
  25 #include "AliL3Evaluate.h"
  26 #include "AliL3Logging.h"
  27
  28 //AliL3Evaluate
  29 //Class for tracking evaluation.
  30
  31 ClassImp(AliL3Evaluate)
  32
  33 AliL3Evaluate::AliL3Evaluate()
  34 {
  35   fMCFile = NULL;
  36   fTracks = NULL;
  37   fMCclusterfile = NULL;
  38   fNFastPoints = 0;
  39   fMcIndex = 0;
  40   fMcId = 0;
  41   fMinSlice=0;
  42   fMaxSlice=0;
  43 }
  44
  45 AliL3Evaluate::AliL3Evaluate(Char_t *mcfile,Int_t *slice)
  46 {
  47   //Normal constructor. Input are the rootfile containing the
  48   //original MC information of the simulated event.
  49
  50   fMCFile = new TFile(mcfile,"READ");
  51   if(!fMCFile->IsOpen())
  52     {
  53       LOG(AliL3Log::kError,"AliL3Evaluation::AliL3Evaluation","File Open")
  54         <<"Inputfile "<<mcfile<<" does not exist"<<ENDLOG;
  55       return;
  56     }
  57
  58   fParam = (AliTPCParam*)fMCFile->Get("75x40_100x60");
  59   fTransform = new AliL3Transform();
  60
  61   fMinSlice = slice[0];
  62   fMaxSlice = slice[1];
  63
  64 }
  65
  66 AliL3Evaluate::AliL3Evaluate(Int_t *slice)
  67 {
  68
  69   fMinSlice = slice[0];
  70   fMaxSlice = slice[1];
  71   fTransform = new AliL3Transform();
  72
  73 }
  74
  75 AliL3Evaluate::~AliL3Evaluate()
  76 {
  77   if(fDigitsTree) fDigitsTree->Delete();
  78   if(fMCFile) {
  79     fMCFile->Close();
  80     delete fMCFile;
  81   }
  82   if(fTransform) delete fTransform;
  83   if(fTracks) delete fTracks;
  84   if(fPtRes) delete fPtRes;
  85   if(fNGoodTracksPt) delete fNGoodTracksPt;
  86   if(fNFoundTracksPt) delete fNFoundTracksPt;
  87   if(fNFakeTracksPt) delete fNFakeTracksPt;
  88   if(fTrackEffPt) delete fTrackEffPt;
  89   if(fFakeTrackEffPt) delete fFakeTrackEffPt;
  90   if(fNGoodTracksEta) delete fNGoodTracksEta;
  91   if(fNFoundTracksEta) delete fNFoundTracksEta;
  92   if(fNFakeTracksEta) delete fNFakeTracksEta;
  93   if(fTrackEffEta) delete fTrackEffEta;
  94   if(fFakeTrackEffEta) delete fFakeTrackEffEta;
  95   if(fMcIndex) delete [] fMcIndex;
  96   if(fMcId)    delete [] fMcId;
  97   if(fNtuppel) delete fNtuppel;
  98 }
  99
 100 void AliL3Evaluate::Setup(Char_t *trackfile,Char_t *path)
 101 {
 102   //Read in the hit and track information from produced files.
 103
 104   Char_t fname[256];
 105   AliL3FileHandler *clusterfile[36][5];
 106   for(Int_t s=fMinSlice; s<=fMaxSlice; s++)
 107     {
 108       for(Int_t p=0; p<5; p++)
 109         {
 110           fClusters[s][p] = 0;
 111           clusterfile[s][p] = new AliL3FileHandler();
 112           sprintf(fname,"%s/points_%d_%d.raw",path,s,p);
 113           if(!clusterfile[s][p]->SetBinaryInput(fname))
 114             {
 115               LOG(AliL3Log::kError,"AliL3Evaluation::Setup","File Open")
 116                 <<"Inputfile "<<fname<<" does not exist"<<ENDLOG;
 117               delete clusterfile[s][p];
 118               clusterfile[s][p] = 0;
 119               continue;
 120             }
 121           fClusters[s][p] = (AliL3SpacePointData*)clusterfile[s][p]->Allocate();
 122           clusterfile[s][p]->Binary2Memory(fNcl[s][p],fClusters[s][p]);
 123           clusterfile[s][p]->CloseBinaryInput();
 124         }
 125     }
 126
 127   AliL3FileHandler *tfile = new AliL3FileHandler();
 128   if(!tfile->SetBinaryInput(trackfile)){
 129     LOG(AliL3Log::kError,"AliL3Evaluation::Setup","File Open")
 130     <<"Inputfile "<<trackfile<<" does not exist"<<ENDLOG;
 131     return;
 132   }
 133   fTracks = new AliL3TrackArray();
 134   tfile->Binary2TrackArray(fTracks);
 135   tfile->CloseBinaryInput();
 136   delete tfile;
 137 }
 138
 139 void AliL3Evaluate::SetupSlow(Char_t *trackfile,Char_t *path)
 140 {
 141   //Setup for using the slow simulator.
 142
 143   fIsSlow = true;
 144   Setup(trackfile,path);
 145
 146   if(!SetDigitsTree())
 147     LOG(AliL3Log::kError,"AliL3Evaluation::SetupSlow","Digits Tree")
 148     <<"Error setting up digits tree"<<ENDLOG;
 149   if(!SetMCParticleArray())
 150     LOG(AliL3Log::kError,"AliL3Evaluation::Setup","Particle array")
 151     <<"Error setting up particle array"<<ENDLOG;
 152 }
 153
 154 void AliL3Evaluate::SetupFast(Char_t *trackfile,Char_t *mcClusterfile,Char_t *path)
 155 {
 156   //Setup for using the fast simulator.
 157
 158   fIsSlow = false;
 159   GetFastClusterIDs(path);
 160
 161   fMCclusterfile = new TFile(mcClusterfile);
 162   if(!fMCclusterfile->IsOpen())
 163     LOG(AliL3Log::kError,"AliL3Evaluation::SetupFast","File Open")
 164       <<"Inputfile "<<mcClusterfile<<" does not exist"<<ENDLOG;
 165
 166   Setup(trackfile,path);
 167
 168   if(!SetMCParticleArray())
 169     LOG(AliL3Log::kError,"AliL3Evaluation::SetupFast","Particle array")
 170       <<"Error setting up particle array"<<ENDLOG;
 171 }
 172
 173 Bool_t AliL3Evaluate::SetDigitsTree()
 174 {
 175
 176   fMCFile->cd();
 177   fDigitsTree = (TTree*)fMCFile->Get("TreeD_75x40_100x60");
 178   if(!fDigitsTree) return false;
 179   fDigitsTree->GetBranch("Segment")->SetAddress(&fDigits);
 180   for(Int_t i=0; i<fDigitsTree->GetEntries(); i++)
 181     {
 182       if(!fDigitsTree->GetEvent(i)) continue;
 183       Int_t se,ro,slice,slicerow;
 184       fParam->AdjustSectorRow(fDigits->GetID(),se,ro);
 185       fTransform->Sector2Slice(slice,slicerow,se,ro);
 186       fRowid[slice][slicerow] = i;
 187     }
 188
 189   return true;
 190 }
 191
 192 Bool_t AliL3Evaluate::SetMCParticleArray()
 193 {
 194   fMCFile->cd();
 195   AliRun *gAlice = (AliRun*)fMCFile->Get("gAlice");
 196   if (!gAlice)
 197     {
 198       LOG(AliL3Log::kError,"AliL3Evaluate::SetParticleArray","gAlice")
 199         <<"AliRun object non existing on file"<<ENDLOG;
 200       return false;
 201     }
 202
 203   gAlice->GetEvent(0);
 204   fParticles=gAlice->Particles();
 205   return true;
 206 }
 207
 208
 209 TObjArray *AliL3Evaluate::DefineGoodTracks(Int_t slice,Int_t *padrow,Int_t good_number,Int_t *particle_id)
 210 {
 211   //Loop over MC particles, and mark the good ones
 212   //(which the tracker should find...)
 213
 214   Int_t np=fParticles->GetEntriesFast();
 215   AliTPC *TPC = (AliTPC*)gAlice->GetDetector("TPC");
 216   TPC->SetParam(fParam);
 217   Int_t ver = TPC->IsVersion();
 218   LOG(AliL3Log::kInformational,"AliL3Evaluate::DefineGoodTracks","TPC version")
 219     <<"TPC version "<<ver<<" found on file"<<ENDLOG;
 220
 221   //Int_t nrow_up=TPC->GetParam()->GetNRowUp();
 222   //Int_t nrows=TPC->GetParam()->GetNRowLow()+nrow_up;
 223   Int_t zero=TPC->GetParam()->GetZeroSup();
 224
 225   //Int_t number_of_rows = padrow[1] - padrow[0] + 1;
 226
 227   Int_t *good = new Int_t[np];
 228   for(Int_t ii=0; ii<np; ii++)
 229     good[ii] = 0;
 230
 231   if(ver==1)
 232     {
 233       if(fIsSlow)
 234         LOG(AliL3Log::kError,"AliL3Evaluate::DefineGoodTracks","TPC version")
 235           <<"TPC version "<<ver<<" does not match."<<ENDLOG;
 236       fMCclusterfile->cd();
 237       AliTPCClustersArray carray;
 238       carray.Setup(fParam);
 239       carray.SetClusterType("AliTPCcluster");
 240       Bool_t clusterok = carray.ConnectTree("Segment Tree");
 241       if(!clusterok)
 242         LOG(AliL3Log::kError,"AliL3Evaluate::DefineGoodTracks","Cluster Array")
 243           <<"Error loading clusters from rootfile"<<ENDLOG;
 244
 245       for(Int_t i=0; i<carray.GetTree()->GetEntries(); i++)
 246         {
 247           Int_t sec,row,sl,lr;
 248           AliSegmentID *s = carray.LoadEntry(i);
 249           fParam->AdjustSectorRow(s->GetID(),sec,row);
 250           fTransform->Sector2Slice(sl,lr,sec,row);
 251
 252           if(sl != slice) {carray.ClearRow(sec,row); continue;}
 253           if(lr < padrow[0]) {carray.ClearRow(sec,row); continue;}
 254           if(lr > padrow[1]) {carray.ClearRow(sec,row); continue;}
 255           AliTPCClustersRow *cRow = carray.GetRow(sec,row);
 256           for(Int_t j=0; j<cRow->GetArray()->GetEntriesFast(); j++)
 257             {
 258               AliTPCcluster *cluster=(AliTPCcluster*)(*cRow)[j];
 259               Int_t lab=cluster->GetLabel(0);
 260               if(lab<0) continue;
 261               lab=TMath::Abs(lab);
 262               good[lab]++;
 263             }
 264           if(carray.GetRow(sec,row))
 265             carray.ClearRow(sec,row);
 266         }
 267     }
 268   else if(ver==2)
 269     {
 270       if(!fIsSlow)
 271         LOG(AliL3Log::kError,"AliL3Evaluate::DefineGoodTracks","TPC version")
 272           <<"TPC version "<<ver<<" does not match."<<ENDLOG;
 273       Int_t *count = new Int_t[np]; //np number of particles.
 274       Int_t i;
 275       for (i=0; i<np; i++) count[i]=0;
 276       for (i=padrow[0]; i<=padrow[1]; i++) {
 277         Int_t index = fRowid[slice][i];
 278         if (!fDigitsTree->GetEvent(index)) continue;
 279         Int_t sec,row;
 280         fParam->AdjustSectorRow(fDigits->GetID(),sec,row);
 281         fDigits->First();
 282         while (fDigits->Next()) {
 283           Int_t it=fDigits->CurrentRow(), ip=fDigits->CurrentColumn();
 284           Short_t dig = fDigits->GetDigit(it,ip);
 285           Int_t idx0=fDigits->GetTrackID(it,ip,0);
 286           Int_t idx1=fDigits->GetTrackID(it,ip,1);
 287           Int_t idx2=fDigits->GetTrackID(it,ip,2);
 288           if (idx0>=0 && dig>=zero) count[idx0]+=1;
 289           if (idx1>=0 && dig>=zero) count[idx1]+=1;
 290           if (idx2>=0 && dig>=zero) count[idx2]+=1;
 291         }
 292         for (Int_t j=0; j<np; j++) {
 293           if (count[j]>1) {//at least two digits at this padrow
 294             good[j]++;
 295           }
 296           count[j]=0;
 297         }
 298       }
 299       delete[] count;
 300     }
 301   else
 302     {
 303       LOG(AliL3Log::kError,"AliL3Evaluation::FillEffHistos","TPC version")
 304         <<"No valid TPC version found"<<ENDLOG;
 305       return 0;
 306     }
 307
 308   Float_t torad=TMath::Pi()/180;
 309
 310   Float_t phi_min = slice*20 - 10;
 311   Float_t phi_max = slice*20 + 10;
 312   TObjArray *good_part = new TObjArray();
 313
 314   for(Int_t i=0; i<fParticles->GetEntriesFast(); i++)
 315    {
 316      TParticle *p = (TParticle*)fParticles->UncheckedAt(i);
 317      if(p->GetFirstMother()>0) continue; //secondary particle
 318      if(good[i] < good_number) {continue;}
 319
 320      Double_t ptg=p->Pt(),pxg=p->Px(),pyg=p->Py(),pzg=p->Pz();
 321      Double_t phi_part = TMath::ATan2(pyg,pxg);
 322      if (phi_part < 0) phi_part += 2*TMath::Pi();
 323
 324
 325      if(phi_part < phi_min*torad || phi_part > phi_max*torad) {continue;}
 326      if(ptg<0.100) continue;
 327      if(fabs(pzg/ptg)>0.999) {continue;}
 328      Int_t entries = good_part->GetEntriesFast();
 329      good_part->AddLast(p);
 330      particle_id[entries] = i;
 331    }
 332   delete [] good;
 333   LOG(AliL3Log::kInformational,"AliL3Evaluate::DefineGoodTracks","NPart")
 334   <<AliL3Log::kDec<<"Found "<<good_part->GetEntriesFast()
 335   <<" good Particles (Tracks) out of "<<fParticles->GetEntriesFast()<<ENDLOG;
 336
 337   return good_part;
 338 }
 339
 340 void AliL3Evaluate::EvaluatePatch(Int_t slice,Int_t patch,Int_t min_points,Int_t good_number)
 341 {
 342   //Make efficiency plots for tracking on patch level (before any merging).
 343
 344   Int_t row[5][2] = {{ 0, 45},{46,77},{78,109},{110,141},{142,173}};
 345   Int_t *particle_id = new Int_t[fParticles->GetEntriesFast()];
 346   TObjArray *good_particles = DefineGoodTracks(slice,row[patch],good_number,particle_id);
 347   SetMinPoints(min_points);
 348   AssignIDs();
 349   CreateHistos();
 350   FillEffHistos(good_particles,particle_id);
 351   CalcEffHistos();
 352   delete good_particles;
 353   delete [] particle_id;
 354 }
 355
 356 void AliL3Evaluate::EvaluateSlice(Int_t slice,Int_t min_points,Int_t good_number)
 357 {
 358   //Make efficiency plots for tracking on a slice (after merging).
 359   //min_points = minimum points on track to be considered for evaluation
 360   //good_number = minimum hits (padrows) produced by simulated track for consideration.
 361
 362   Int_t row[2] = {0,173};
 363   Int_t *particle_id = new Int_t[fParticles->GetEntriesFast()];
 364   TObjArray *good_particles = DefineGoodTracks(slice,row,good_number,particle_id);
 365   SetMinPoints(min_points);
 366   AssignIDs();
 367   CreateHistos();
 368   FillEffHistos(good_particles,particle_id);
 369   CalcEffHistos();
 370   delete good_particles;
 371   delete [] particle_id;
 372 }
 373
 374 void AliL3Evaluate::EvaluateGlobal(Int_t min_points,Int_t good_number)
 375 {
 376   //Make efficiency plots for tracking on several slices.
 377
 378   Int_t row[2] = {0,173};
 379   Int_t *particle_id = new Int_t[fParticles->GetEntriesFast()];
 380   SetMinPoints(min_points);
 381   AssignIDs();
 382   CreateHistos();
 383   for(Int_t slice=fMinSlice;slice<=fMaxSlice;slice++){
 384     TObjArray *good_particles = DefineGoodTracks(slice,row,good_number,particle_id);
 385     FillEffHistos(good_particles,particle_id);
 386     delete good_particles;
 387   }
 388   CalcEffHistos();
 389 }
 390
 391 void AliL3Evaluate::AssignIDs()
 392 {
 393   //Assign MC id to the tracks.
 394
 395   fTracks->QSort();
 396   LOG(AliL3Log::kDebug,"AliL3Evaluate::AssignIDs","Track Loop")
 397     <<"Assigning MC id to the found tracks...."<<ENDLOG;
 398   for(Int_t i=0; i<fTracks->GetNTracks(); i++)
 399     {
 400       AliL3Track *track = (AliL3Track*)fTracks->GetCheckedTrack(i);
 401       if(!track) continue;
 402       if(track->GetNumberOfPoints() < fMinPointsOnTrack) break;
 403       Int_t tID = GetMCTrackLabel(track);
 404       track->SetMCid(tID);
 405       printf("track %i id %d\n",i,tID);
 406     }
 407 }
 408
 409
 410 struct S {Int_t lab; Int_t max;};
 411 Int_t AliL3Evaluate::GetMCTrackLabel(AliL3Track *track){
 412   //Returns the MCtrackID of the belonging clusters.
 413   //If MCLabel < 0, means that track is fake.
 414   //Fake track means that more than 10 percent of clusters are assigned incorrectly.
 415
 416   Int_t num_of_clusters = track->GetNumberOfPoints();
 417   S *s=new S[num_of_clusters];
 418   Int_t i;
 419   for (i=0; i<num_of_clusters; i++) s[i].lab=s[i].max=0;
 420   UInt_t *hitnum = track->GetHitNumbers();
 421   UInt_t id;
 422
 423   Int_t **  trackID = GetClusterIDs(track);
 424
 425   Int_t lab=123456789;
 426   for (i=0; i<num_of_clusters; i++)
 427     {
 428       //Tricks to get the clusters belonging to this track:
 429       id = hitnum[i];
 430       Int_t slice = (id>>25) & 0x7f;
 431       Int_t patch = (id>>22) & 0x7;
 432       UInt_t pos = id&0x3fffff;
 433
 434       AliL3SpacePointData *points = fClusters[slice][patch];
 435       if(!points) continue;
 436       if(pos>=fNcl[slice][patch])
 437         {
 438           LOG(AliL3Log::kError,"AliL3Evaluate::GetMCTrackLabel","Clusterarray")
 439             <<AliL3Log::kDec<<"ERROR"<<ENDLOG;
 440           continue;
 441         }
 442
 443       //Get the label of the cluster:
 444       //printf("label %d %d %d\n",trackID[i][0],trackID[i][1],trackID[i][2]);
 445       lab=abs(trackID[i][0]);
 446       Int_t j;
 447       for (j=0; j<num_of_clusters; j++)
 448         if (s[j].lab==lab || s[j].max==0) break;
 449       s[j].lab=lab;
 450       s[j].max++;
 451     }
 452
 453   Int_t max=0;
 454   for (i=0; i<num_of_clusters; i++)
 455     if (s[i].max>max) {max=s[i].max; lab=s[i].lab;}
 456
 457   delete[] s;
 458
 459   for (i=0; i<num_of_clusters; i++)
 460     {
 461       id = hitnum[i];
 462       Int_t slice = (id>>25) & 0x7f;
 463       Int_t patch = (id>>22) & 0x7;
 464       UInt_t pos = id&0x3fffff;
 465
 466       AliL3SpacePointData *points = fClusters[slice][patch];
 467       if(!points) continue;
 468       if(pos>=fNcl[slice][patch])
 469         {
 470           LOG(AliL3Log::kError,"AliL3Evaluate::GetMCTrackLabel","Clusterarray")
 471             <<AliL3Log::kDec<<"ERROR"<<ENDLOG;
 472           continue;
 473         }
 474
 475       if (abs(trackID[i][1]) == lab ||
 476           abs(trackID[i][2]) == lab ) max++;
 477     }
 478
 479   //check if more than 10% of the clusters are incorrectly assigned (fake track):
 480   if (1.-Float_t(max)/num_of_clusters > 0.10)
 481     {
 482       return -lab;
 483     }
 484
 485   delete [] trackID;
 486   return lab;
 487 }
 488
 489
 490 Int_t **AliL3Evaluate::GetClusterIDs(AliL3Track *track)
 491 {
 492   //Return the MC information of all clusters belonging to track.
 493
 494   Int_t num_of_clusters = track->GetNumberOfPoints();
 495   Int_t **trackID = new Int_t*[num_of_clusters];
 496
 497   UInt_t *hitnum = track->GetHitNumbers();
 498   UInt_t id;
 499
 500   Float_t xyz[3];
 501   Int_t padrow;
 502   for(Int_t i=0; i<num_of_clusters; i++)
 503     {
 504       id = hitnum[i];
 505       Int_t slice = (id>>25) & 0x7f;
 506       Int_t patch = (id>>22) & 0x7;
 507       UInt_t pos = id&0x3fffff;
 508
 509       AliL3SpacePointData *points = fClusters[slice][patch];
 510
 511       if(!points)
 512         {
 513           LOG(AliL3Log::kError,"AliL3Evaluate::GetClusterIDs","Clusterarray")
 514             <<"No points at slice "<<slice<<" patch "<<patch<<" pos "<<pos<<ENDLOG;
 515           continue;
 516         }
 517       if(pos>=fNcl[slice][patch])
 518         {
 519           LOG(AliL3Log::kError,"AliL3Evaluate::GetClusterIDs","Clusterarray")
 520             <<AliL3Log::kDec<<"ERROR"<<ENDLOG;
 521           continue;
 522         }
 523
 524       xyz[0] = points[pos].fX;
 525       xyz[1] = points[pos].fY;
 526       xyz[2] = points[pos].fZ;
 527       //sector = points[pos].fSector;
 528       padrow = points[pos].fPadRow;
 529       Int_t se,ro;
 530       fTransform->Slice2Sector(slice,padrow,se,ro);
 531       fTransform->Global2Raw(xyz,se,ro);
 532
 533       if(fIsSlow)
 534         {
 535           Int_t p = fRowid[slice][padrow];
 536
 537           if(!fDigitsTree->GetEvent(p))
 538             LOG(AliL3Log::kError,"AliL3Evaluate::GetClusterIDs","Digits Tree")
 539               <<"Error reading digits tree"<<ENDLOG;
 540
 541           trackID[i] = new Int_t[3];
 542           trackID[i][0] = fDigits->GetTrackID((Int_t)xyz[2],(Int_t)xyz[1],0);
 543           trackID[i][1] = fDigits->GetTrackID((Int_t)xyz[2],(Int_t)xyz[1],1);
 544           trackID[i][2] = fDigits->GetTrackID((Int_t)xyz[2],(Int_t)xyz[1],2);
 545         }
 546       else
 547         {
 548           Int_t tmp_pid=0;
 549           for(Int_t ii=0; ii<fNFastPoints; ii++)
 550             {
 551               tmp_pid = fMcId[ii];
 552               if(fMcIndex[ii] == id) break;
 553             }
 554           trackID[i] = new Int_t[3];
 555           trackID[i][0] = tmp_pid;
 556           trackID[i][1] = -1;
 557           trackID[i][2] = -1;
 558         }
 559     }
 560   return trackID;
 561 }
 562
 563 void AliL3Evaluate::GetFastClusterIDs(Char_t *path)
 564 {
 565   //Get the MC id of space points in case of using the fast simulator.
 566   char fname[256];
 567   sprintf(fname,"%s/point_mc.dat",path);
 568   FILE *infile = fopen(fname,"r");
 569   if(!infile) return;
 570   Int_t hitid,hitmc,i;
 571
 572   for(i=0; ; i++)
 573     if(fscanf(infile,"%d %d",&hitid,&hitmc)==EOF) break;
 574   rewind(infile);
 575   fNFastPoints = i;
 576   fMcId = new Int_t[fNFastPoints];
 577   fMcIndex = new UInt_t[fNFastPoints];
 578
 579   for(i=0; i<fNFastPoints; i++)
 580     {
 581       if(fscanf(infile,"%d %d",&hitid,&hitmc)==EOF) break;
 582       fMcId[i] = hitmc;
 583       fMcIndex[i] = hitid;
 584     }
 585   fclose(infile);
 586 }
 587
 588 void AliL3Evaluate::CreateHistos(Int_t nbin,Int_t xlow,Int_t xup)
 589 {
 590   //Create the histograms
 591
 592   fNtuppel = new TNtuple("fNtuppel","Pt resolution","pt_gen:pt_found:nHits");
 593
 594   fPtRes = new TH1F("fPtRes","Relative Pt resolution",30,-10.,10.);
 595   fNGoodTracksPt = new TH1F("fNGoodTracksPt","Good tracks vs pt",nbin,xlow,xup);
 596   fNFoundTracksPt = new TH1F("fNFoundTracksPt","Found tracks vs pt",nbin,xlow,xup);
 597   fNFakeTracksPt = new TH1F("fNFakeTracksPt","Fake tracks vs pt",nbin,xlow,xup);
 598   fTrackEffPt = new TH1F("fTrackEffPt","Tracking efficiency vs pt",nbin,xlow,xup);
 599   fFakeTrackEffPt = new TH1F("fFakeTrackEffPt","Efficiency for fake tracks vs pt",nbin,xlow,xup);
 600
 601   fNGoodTracksEta = new TH1F("fNGoodTracksEta","Good tracks vs eta",20,-50,50);
 602   fNFoundTracksEta = new TH1F("fNFoundTracksEta","Found tracks vs eta",20,-50,50);
 603   fNFakeTracksEta = new TH1F("fNFakeTracksEta","Fake tracks vs eta",20,-50,50);
 604   fTrackEffEta = new TH1F("fTrackEffEta","Tracking efficienct vs eta",20,-50,50);
 605   fFakeTrackEffEta = new TH1F("fFakeTrackEffEta","Efficiency for fake tracks vs eta",20,-50,50);
 606 }
 607
 608 void AliL3Evaluate::FillEffHistos(TObjArray *good_particles,Int_t *particle_id)
 609 {
 610   //Fill the efficiency histograms.
 611
 612   for(Int_t i=0; i<good_particles->GetEntriesFast(); i++)
 613     {
 614       TParticle *p = (TParticle*)good_particles->UncheckedAt(i);
 615       Double_t ptg=p->Pt(),pzg=p->Pz();
 616       Float_t dipangle=TMath::ATan2(pzg,ptg)*180./TMath::Pi();
 617       fNGoodTracksPt->Fill(ptg);
 618       fNGoodTracksEta->Fill(dipangle);
 619       Int_t found = 0;
 620       for(Int_t k=0; k<fTracks->GetNTracks(); k++)
 621         {
 622           AliL3Track *track = fTracks->GetCheckedTrack(k);
 623           if(!track) continue;
 624           Int_t nHits = track->GetNumberOfPoints();
 625           if(nHits < fMinPointsOnTrack) break;
 626
 627           Int_t tracklabel;
 628           tracklabel = track->GetMCid();
 629
 630           if(TMath::Abs(tracklabel) != particle_id[i]) continue;
 631           found=1;
 632           if(tracklabel == particle_id[i]) {fNFoundTracksPt->Fill(ptg); fNFoundTracksEta->Fill(dipangle);}
 633           else {fNFakeTracksPt->Fill(ptg); fNFakeTracksEta->Fill(dipangle);}
 634           Float_t pt=track->GetPt();
 635           fPtRes->Fill((pt-ptg)/ptg*100.);
 636           fNtuppel->Fill(ptg,pt,nHits);
 637           break;
 638
 639         }
 640     }
 641 }
 642
 643 void AliL3Evaluate::CalcEffHistos(){
 644
 645   Stat_t ngood=fNGoodTracksPt->GetEntries();
 646   Stat_t nfound=fNFoundTracksPt->GetEntries();
 647   Stat_t nfake=fNFakeTracksPt->GetEntries();
 648
 649   LOG(AliL3Log::kInformational,"AliL3Evaluate::FillEffHistos","Efficiency")
 650     <<AliL3Log::kDec<<"There was "<<ngood<<" generated good tracks"<<ENDLOG;
 651   LOG(AliL3Log::kInformational,"AliL3Evaluate::FillEffHistos","Efficiency")
 652     <<AliL3Log::kDec<<"Found "<<nfound<<" tracks"<<ENDLOG;
 653   LOG(AliL3Log::kInformational,"AliL3Evaluate::FillEffHistos","Efficiency")
 654     <<AliL3Log::kDec<<"Integral efficiency is about "<<nfound/ngood*100<<ENDLOG;
 655   LOG(AliL3Log::kInformational,"AliL3Evaluate::FillEffHistos","Efficiency")
 656     <<AliL3Log::kDec<<"Fake tracks relative is about "<<nfake/ngood*100<<ENDLOG;
 657
 658   fNFoundTracksPt->Sumw2(); fNGoodTracksPt->Sumw2();
 659   fTrackEffPt->Divide(fNFoundTracksPt,fNGoodTracksPt,1,1.,"b");
 660   fFakeTrackEffPt->Divide(fNFakeTracksPt,fNGoodTracksPt,1,1.,"b");
 661   fTrackEffPt->SetMaximum(1.4);
 662   fTrackEffPt->SetXTitle("P_{T} [GeV]");
 663   fTrackEffPt->SetLineWidth(2);
 664   fFakeTrackEffPt->SetFillStyle(3013);
 665   fTrackEffPt->SetLineColor(4);
 666   fFakeTrackEffPt->SetFillColor(2);
 667
 668   fNFoundTracksEta->Sumw2(); fNGoodTracksEta->Sumw2();
 669   fTrackEffEta->Divide(fNFoundTracksEta,fNGoodTracksEta,1,1.,"b");
 670   fFakeTrackEffEta->Divide(fNFakeTracksEta,fNGoodTracksEta,1,1.,"b");
 671   fTrackEffEta->SetMaximum(1.4);
 672   fTrackEffEta->SetXTitle("#lambda [degrees]");
 673   fTrackEffEta->SetLineWidth(2);
 674   fFakeTrackEffEta->SetFillStyle(3013);
 675   fTrackEffEta->SetLineColor(4);
 676   fFakeTrackEffEta->SetFillColor(2);
 677
 678 }
 679
 680 void AliL3Evaluate::Write2File(Char_t *outputfile)
 681 {
 682   //Write histograms to file:
 683
 684   TFile *of = new TFile(outputfile,"RECREATE");
 685   if(!of->IsOpen())
 686     {
 687       LOG(AliL3Log::kError,"AliL3Evaluate::Write2File","File Open")
 688         <<"Problems opening rootfile"<<ENDLOG;
 689       return;
 690     }
 691
 692   of->cd();
 693   fNtuppel->Write();
 694   fPtRes->Write();
 695   fNGoodTracksPt->Write();
 696   fNFoundTracksPt->Write();
 697   fNFakeTracksPt->Write();
 698   fTrackEffPt->Write();
 699   fFakeTrackEffPt->Write();
 700   fNGoodTracksEta->Write();
 701   fNFoundTracksEta->Write();
 702   fNFakeTracksEta->Write();
 703   fTrackEffEta->Write();
 704   fFakeTrackEffEta->Write();
 705
 706   of->Close();
 707   delete of;
 708
 709 }
 710
 711 TNtuple *AliL3Evaluate::CalculateResiduals()
 712 {
 713
 714   TNtuple *ntuppel=new TNtuple("ntuppel","Residuals","residual_trans:residual_long:zHit:pt:dipangle:beta:padrow:nHits");
 715
 716   for(int f=fMinSlice; f<=fMaxSlice; f++)
 717     {
 718       AliL3FileHandler *tfile = new AliL3FileHandler();
 719       char fname[256];
 720       sprintf(fname,"tracks_tr_%d_0.raw",f);
 721       if(!tfile->SetBinaryInput(fname)){
 722         LOG(AliL3Log::kError,"AliL3Evaluation::Setup","File Open")
 723           <<"Inputfile "<<fname<<" does not exist"<<ENDLOG;
 724         return 0;
 725       }
 726       fTracks = new AliL3TrackArray();
 727       tfile->Binary2TrackArray(fTracks);
 728       tfile->CloseBinaryInput();
 729       delete tfile;
 730
 731
 732       for(Int_t i=0; i<fTracks->GetNTracks(); i++)
 733         {
 734
 735           AliL3Track *track = (AliL3Track*)fTracks->GetCheckedTrack(i);
 736           if(!track) continue;
 737
 738           track->CalculateHelix();
 739           UInt_t *hitnum = track->GetHitNumbers();
 740           UInt_t id;
 741
 742           Float_t xyz[3];
 743           Int_t padrow;
 744           for(Int_t j=0; j<track->GetNumberOfPoints(); j++)
 745             {
 746               id = hitnum[j];
 747               Int_t slice = (id>>25) & 0x7f;
 748               Int_t patch = (id>>22) & 0x7;
 749               UInt_t pos = id&0x3fffff;
 750
 751               AliL3SpacePointData *points = fClusters[slice][patch];
 752
 753               if(!points)
 754                 {
 755                   LOG(AliL3Log::kError,"AliL3Evaluate::CalculateResiduals","Clusterarray")
 756                     <<"No points at slice "<<slice<<" patch "<<patch<<" pos "<<pos<<ENDLOG;
 757                   continue;
 758                 }
 759               if(pos>=fNcl[slice][patch])
 760                 {
 761                   LOG(AliL3Log::kError,"AliL3Evaluate::CalculateResiduals","Clusterarray")
 762                     <<AliL3Log::kDec<<"ERROR"<<ENDLOG;
 763                   continue;
 764                 }
 765
 766               xyz[0] = points[pos].fX;
 767               xyz[1] = points[pos].fY;
 768               xyz[2] = points[pos].fZ;
 769               padrow = points[pos].fPadRow;
 770               fTransform->Global2Local(xyz,slice);
 771
 772               Float_t xyz_cross[3];
 773               track->GetCrossingPoint(padrow,xyz_cross);
 774               Double_t beta = track->GetCrossingAngle(padrow);
 775
 776               Double_t yres = xyz_cross[1] - xyz[1];
 777               Double_t zres = xyz_cross[2] - xyz[2];
 778
 779               Double_t dipangle = atan(track->GetTgl());
 780               ntuppel->Fill(yres,zres,xyz_cross[2],track->GetPt(),dipangle,beta,padrow,track->GetNumberOfPoints());
 781
 782             }
 783         }
 784       if(fTracks)
 785         delete fTracks;
 786     }
 787   return ntuppel;
 788 }
 789