HLT/src/AliL3Evaluate.cxx

   1 //Author:        Anders Strand Vestbo
   2 //Last Modified: 5.01.2001
   3
   4 #include <stdio.h>
   5 #include <math.h>
   6 #include <fstream.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 "AliTPCClustersArray.h"
  17 #include "AliTPCClustersRow.h"
  18 #include "AliTPCcluster.h"
  19 #include "AliTPCParam.h"
  20 #include "AliComplexCluster.h"
  21
  22 #include "AliL3Defs.h"
  23 #include "AliL3Transform.h"
  24 #include "AliL3SpacePointData.h"
  25 #include "AliL3Track.h"
  26 #include "AliL3FileHandler.h"
  27 #include "AliL3TrackArray.h"
  28 #include "AliL3Evaluate.h"
  29 #include "AliL3Logging.h"
  30
  31 //AliL3Evaluate
  32 //Class for tracking evaluation.
  33
  34 ClassImp(AliL3Evaluate)
  35
  36 AliL3Evaluate::AliL3Evaluate()
  37 {
  38   fDigitsFile = NULL;
  39   fTracks = NULL;
  40   fMCclusterfile = NULL;
  41   fNFastPoints = 0;
  42   fMcIndex = 0;
  43   fMcId = 0;
  44   fMinSlice=0;
  45   fMaxSlice=0;
  46   fTransform = new AliL3Transform();
  47 }
  48
  49 AliL3Evaluate::AliL3Evaluate(Char_t *mcfile,Int_t *slice)
  50 {
  51   //Normal constructor. Input are the rootfile containing the
  52   //original MC information of the simulated event.
  53
  54   fEventFile = new TFile(mcfile,"READ");
  55
  56   fParam = (AliTPCParam*)fEventFile->Get("75x40_100x60");
  57   fTransform = new AliL3Transform();
  58
  59   fMinSlice = slice[0];
  60   fMaxSlice = slice[1];
  61   fMinGoodPt = 0.1;
  62   fNoOverlap = kFALSE;
  63 }
  64
  65 AliL3Evaluate::AliL3Evaluate(Int_t *slice)
  66 {
  67   //ctor to use if you do not need any rootfile.
  68
  69
  70   fMinSlice = slice[0];
  71   fMaxSlice = slice[1];
  72   fTransform = new AliL3Transform();
  73
  74 }
  75
  76 AliL3Evaluate::~AliL3Evaluate()
  77 {
  78   if(fDigitsTree) fDigitsTree->Delete();
  79   if(fDigitsFile) {
  80     fDigitsFile->Close();
  81     delete fDigitsFile;
  82   }
  83   if(fTransform) delete fTransform;
  84   if(fTracks) delete fTracks;
  85   if(fPtRes) delete fPtRes;
  86   if(fNGoodTracksPt) delete fNGoodTracksPt;
  87   if(fNFoundTracksPt) delete fNFoundTracksPt;
  88   if(fNFakeTracksPt) delete fNFakeTracksPt;
  89   if(fTrackEffPt) delete fTrackEffPt;
  90   if(fFakeTrackEffPt) delete fFakeTrackEffPt;
  91   if(fNGoodTracksEta) delete fNGoodTracksEta;
  92   if(fNFoundTracksEta) delete fNFoundTracksEta;
  93   if(fNFakeTracksEta) delete fNFakeTracksEta;
  94   if(fTrackEffEta) delete fTrackEffEta;
  95   if(fFakeTrackEffEta) delete fFakeTrackEffEta;
  96   if(fMcIndex) delete [] fMcIndex;
  97   if(fMcId)    delete [] fMcId;
  98   if(fNtuppel) delete fNtuppel;
  99 }
 100
 101 void AliL3Evaluate::Setup(Char_t *trackfile,Char_t *path)
 102 {
 103   //Read in the hit and track information from produced files.
 104
 105   fGoodFound = 0;
 106   fGoodGen = 0;
 107   Char_t fname[256];
 108   AliL3FileHandler *clusterfile[36][6];
 109   for(Int_t s=fMinSlice; s<=fMaxSlice; s++)
 110     {
 111       for(Int_t p=0; p<6; p++)
 112         {
 113           fClusters[s][p] = 0;
 114           clusterfile[s][p] = new AliL3FileHandler();
 115           sprintf(fname,"%s/points_%d_%d.raw",path,s,p);
 116           if(!clusterfile[s][p]->SetBinaryInput(fname))
 117             {
 118               LOG(AliL3Log::kError,"AliL3Evaluation::Setup","File Open")
 119                 <<"Inputfile "<<fname<<" does not exist"<<ENDLOG;
 120               delete clusterfile[s][p];
 121               clusterfile[s][p] = 0;
 122               continue;
 123             }
 124           fClusters[s][p] = (AliL3SpacePointData*)clusterfile[s][p]->Allocate();
 125           clusterfile[s][p]->Binary2Memory(fNcl[s][p],fClusters[s][p]);
 126           clusterfile[s][p]->CloseBinaryInput();
 127         }
 128     }
 129
 130
 131   AliL3FileHandler *tfile = new AliL3FileHandler();
 132   if(!tfile->SetBinaryInput(trackfile)){
 133     LOG(AliL3Log::kError,"AliL3Evaluation::Setup","File Open")
 134     <<"Inputfile "<<trackfile<<" does not exist"<<ENDLOG;
 135     return;
 136   }
 137   fTracks = new AliL3TrackArray();
 138   tfile->Binary2TrackArray(fTracks);
 139   tfile->CloseBinaryInput();
 140   delete tfile;
 141 }
 142
 143 void AliL3Evaluate::SetupSlow(Char_t *trackfile,Char_t *digitsfile,Char_t *path)
 144 {
 145   //Setup for using the slow simulator.
 146
 147   fDigitsFile = new TFile(digitsfile,"READ");
 148   if(!fDigitsFile->IsOpen())
 149     {
 150       LOG(AliL3Log::kError,"AliL3Evaluate::SetupSlow","File Open")
 151         <<"Inputfile "<<digitsfile<<" does not exist"<<ENDLOG;
 152       return;
 153     }
 154   fIsSlow = true;
 155   Setup(trackfile,path);
 156
 157   if(!InitMC())
 158     LOG(AliL3Log::kError,"AliL3Evaluation::SetupSlow","Digits Tree")
 159     <<"Error setting up digits tree"<<ENDLOG;
 160
 161 }
 162
 163 void AliL3Evaluate::SetupFast(Char_t *trackfile,Char_t *mcClusterfile,Char_t *path)
 164 {
 165   //Setup for using the fast simulator.
 166
 167   fIsSlow = false;
 168   GetFastClusterIDs(path);
 169
 170   fMCclusterfile = new TFile(mcClusterfile);
 171   if(!fMCclusterfile->IsOpen())
 172     LOG(AliL3Log::kError,"AliL3Evaluation::SetupFast","File Open")
 173       <<"Inputfile "<<mcClusterfile<<" does not exist"<<ENDLOG;
 174
 175   Setup(trackfile,path);
 176   InitMC();
 177 }
 178
 179 Bool_t AliL3Evaluate::InitMC()
 180 {
 181   if(fIsSlow)
 182     {
 183       fDigitsFile->cd();
 184       fDigitsTree = (TTree*)fDigitsFile->Get("TreeD_75x40_100x60_0");
 185       if(!fDigitsTree)
 186         {
 187           LOG(AliL3Log::kError,"AliL3Evaluate::InitMC","Digits Tree")
 188             <<AliL3Log::kHex<<"Error getting digitstree "<<(Int_t)fDigitsTree<<ENDLOG;
 189           return false;
 190         }
 191       fDigitsTree->GetBranch("Segment")->SetAddress(&fDigits);
 192       for(Int_t i=0; i<fDigitsTree->GetEntries(); i++)
 193         {
 194           if(!fDigitsTree->GetEvent(i)) continue;
 195           Int_t se,ro,slice,slicerow;
 196           fParam->AdjustSectorRow(fDigits->GetID(),se,ro);
 197           fTransform->Sector2Slice(slice,slicerow,se,ro);
 198           fRowid[slice][slicerow] = i;
 199         }
 200     }
 201
 202   fEventFile->cd();
 203   gAlice = (AliRun*)fEventFile->Get("gAlice");
 204   if (!gAlice)
 205     {
 206       LOG(AliL3Log::kError,"AliL3Evaluate::InitMC","gAlice")
 207         <<"AliRun object non existing on file"<<ENDLOG;
 208       return false;
 209     }
 210
 211   gAlice->GetEvent(0);
 212
 213   return true;
 214
 215 }
 216
 217
 218
 219 void AliL3Evaluate::DefineGoodTracks(Int_t *slice,Int_t *padrow,Int_t good_number,Char_t *fname)
 220 {
 221   //Loop over MC particles, and mark the good ones
 222   //(which the tracker should find...)
 223
 224
 225   ifstream in(fname);
 226   if(in)
 227     {
 228       LOG(AliL3Log::kInformational,"AliL3Evaluate::DefineGoodTracks","infile")
 229         <<"Reading good particles from file "<<fname<<ENDLOG;
 230       while (in>>fGoodTracks[fGoodGen].label>>fGoodTracks[fGoodGen].code>>
 231              fGoodTracks[fGoodGen].px>>fGoodTracks[fGoodGen].py>>fGoodTracks[fGoodGen].pz>>
 232              fGoodTracks[fGoodGen].pt)
 233         {
 234           fGoodGen++;
 235           if (fGoodGen==15000)
 236             {
 237               LOG(AliL3Log::kWarning,"AliL3Evaluate::DefineGoodTracks","fGoodGen")
 238                 <<"Too many good tracks"<<ENDLOG;
 239               break;
 240             }
 241         }
 242       if (!in.eof()) cerr<<"Read error (good_tracks_tpc) !\n";
 243       return;
 244     }
 245
 246
 247   AliTPC *TPC = (AliTPC*)gAlice->GetDetector("TPC");
 248
 249   TPC->SetParam(fParam);
 250
 251   Int_t ver = TPC->IsVersion();
 252   LOG(AliL3Log::kInformational,"AliL3Evaluate::DefineGoodTracks","TPC version")
 253     <<"TPC version "<<ver<<" found on file"<<ENDLOG;
 254
 255   Int_t zero=TPC->GetParam()->GetZeroSup();
 256
 257   Int_t np = gAlice->GetNtrack();
 258   Int_t *good = new Int_t[np];
 259   for(Int_t ii=0; ii<np; ii++)
 260     good[ii] = 0;
 261
 262   if(ver==1)
 263     {
 264       if(fIsSlow)
 265         LOG(AliL3Log::kError,"AliL3Evaluate::DefineGoodTracks","TPC version")
 266           <<"TPC version "<<ver<<" does not match."<<ENDLOG;
 267       fMCclusterfile->cd();
 268       AliTPCClustersArray carray;
 269       carray.Setup(fParam);
 270       carray.SetClusterType("AliTPCcluster");
 271       Char_t cname[100];
 272       Int_t eventn = 0;
 273       sprintf(cname,"TreeC_TPC_%d",eventn);
 274       Bool_t clusterok = carray.ConnectTree(cname);
 275       if(!clusterok)
 276         LOG(AliL3Log::kError,"AliL3Evaluate::DefineGoodTracks","Cluster Array")
 277           <<"Error loading clusters from rootfile"<<ENDLOG;
 278
 279       for(Int_t i=0; i<carray.GetTree()->GetEntries(); i++)
 280         {
 281           Int_t sec,row,sl,lr;
 282           AliSegmentID *s = carray.LoadEntry(i);
 283           fParam->AdjustSectorRow(s->GetID(),sec,row);
 284           fTransform->Sector2Slice(sl,lr,sec,row);
 285
 286           if(sl != slice[0]) {carray.ClearRow(sec,row); continue;}
 287           if(lr < padrow[0]) {carray.ClearRow(sec,row); continue;}
 288           if(lr > padrow[1]) {carray.ClearRow(sec,row); continue;}
 289           AliTPCClustersRow *cRow = carray.GetRow(sec,row);
 290           for(Int_t j=0; j<cRow->GetArray()->GetEntriesFast(); j++)
 291             {
 292               AliTPCcluster *cluster=(AliTPCcluster*)(*cRow)[j];
 293               Int_t lab=cluster->GetLabel(0);
 294               if(lab<0) continue;
 295               lab=TMath::Abs(lab);
 296               good[lab]++;
 297             }
 298           if(carray.GetRow(sec,row))
 299             carray.ClearRow(sec,row);
 300         }
 301     }
 302   else if(ver==2)
 303     {
 304       if(!fIsSlow)
 305         LOG(AliL3Log::kError,"AliL3Evaluate::DefineGoodTracks","TPC version")
 306           <<"TPC version "<<ver<<" does not match."<<ENDLOG;
 307       Int_t *count = new Int_t[np]; //np number of particles.
 308       Int_t i;
 309       for (i=0; i<np; i++) count[i]=0;
 310       for(Int_t sl=slice[0]; sl<=slice[1]; sl++)
 311         {
 312           for (i=padrow[0]; i<=padrow[1]; i++) {
 313             Int_t index = fRowid[sl][i];
 314             if (!fDigitsTree->GetEvent(index)) continue;
 315             Int_t sec,row;
 316             fParam->AdjustSectorRow(fDigits->GetID(),sec,row);
 317             fDigits->First();
 318             do {
 319               Int_t it=fDigits->CurrentRow(), ip=fDigits->CurrentColumn();
 320               Short_t dig = fDigits->GetDigit(it,ip);
 321               Int_t idx0=fDigits->GetTrackID(it,ip,0);
 322               Int_t idx1=fDigits->GetTrackID(it,ip,1);
 323               Int_t idx2=fDigits->GetTrackID(it,ip,2);
 324
 325               if (idx0>=0 && dig>=zero) count[idx0]+=1;
 326               if (idx1>=0 && dig>=zero) count[idx1]+=1;
 327               if (idx2>=0 && dig>=zero) count[idx2]+=1;
 328             } while (fDigits->Next());
 329
 330             for (Int_t j=0; j<np; j++)
 331               {
 332                 if (count[j]>1) //at least two digits at this padrow
 333                   good[j]++;
 334
 335                 count[j]=0;
 336               }
 337           }
 338         }
 339       delete[] count;
 340     }
 341
 342   else
 343     {
 344       LOG(AliL3Log::kError,"AliL3Evaluation::FillEffHistos","TPC version")
 345         <<"No valid TPC version found"<<ENDLOG;
 346       return;
 347     }
 348
 349
 350   //The following code has been taken from offlinemacro->AliTPCComparison.C
 351
 352   TTree *TH=gAlice->TreeH();
 353   Int_t npart=(Int_t)TH->GetEntries();
 354   Int_t max = 15000;
 355   while (npart--) {
 356     AliTPChit *hit0=0;
 357
 358     TPC->ResetHits();
 359     TH->GetEvent(npart);
 360     AliTPChit * hit = (AliTPChit*) TPC->FirstHit(-1);
 361     while (hit){
 362       if (hit->fQ==0.) break;
 363       hit =  (AliTPChit*) TPC->NextHit();
 364     }
 365     if (hit) {
 366       hit0 = new AliTPChit(*hit); //Make copy of hit
 367       hit = hit0;
 368     }
 369     else continue;
 370     AliTPChit *hit1=(AliTPChit*)TPC->NextHit();
 371     if (hit1==0) continue;
 372     if (hit1->fQ != 0.) continue;
 373     Int_t i=hit->Track();
 374     TParticle *p = (TParticle*)gAlice->Particle(i);
 375
 376     //printf("Checking particle %d with code %d\n",i,p->GetPdgCode());
 377     if (p->GetFirstMother()>=0) continue;  //secondary particle
 378     if (good[i] < good_number) continue;
 379     if (p->Pt()<fMinGoodPt) continue;
 380     if (TMath::Abs(p->Pz()/p->Pt())>0.999) continue;
 381     printf("Found good particle %d, nHits %d\n",i,good[i]);
 382
 383     fGoodTracks[fGoodGen].label=i;
 384     fGoodTracks[fGoodGen].code=p->GetPdgCode();
 385     //**** px py pz - in global coordinate system, x y z - in local !
 386     fGoodTracks[fGoodGen].px=hit->X(); fGoodTracks[fGoodGen].py=hit->Y(); fGoodTracks[fGoodGen].pz=hit->Z();
 387     fGoodTracks[fGoodGen].pt = p->Pt();
 388     fGoodGen++;
 389     //nt++;
 390     if (hit0) delete hit0;
 391     if (fGoodGen==max) {cerr<<"Too many good tracks !n"; break;}
 392   }
 393
 394   delete [] good;
 395
 396   LOG(AliL3Log::kInformational,"AliL3Evaluate::DefineGoodParticles","Eval")
 397     <<AliL3Log::kDec<<"Found "<<fGoodGen<<"good tracks"<<ENDLOG;
 398
 399   ofstream out(fname);
 400   if(out)
 401     {
 402       for (Int_t ngd=0; ngd<fGoodGen; ngd++)
 403         out<<fGoodTracks[ngd].label<<' '<<fGoodTracks[ngd].code<<' '<<
 404           fGoodTracks[ngd].px<<' '<<fGoodTracks[ngd].py<<' '<<fGoodTracks[ngd].pz<<' '<<
 405           fGoodTracks[ngd].pt<<endl;
 406
 407     }
 408   else
 409     LOG(AliL3Log::kError,"AliL3Evaluate::DefineGoodParticles","file")
 410       <<"Can not open file containing good tracks"<<ENDLOG;
 411
 412   out.close();
 413
 414 }
 415
 416 void AliL3Evaluate::EvaluatePatch(Int_t slice,Int_t patch,Int_t min_points,Int_t good_number)
 417 {
 418   //Make efficiency plots for tracking on patch level (before any merging).
 419
 420   Int_t row[6][2] = {{ 0, 45},{46,77},{78,109},{110,141},{142,175}};
 421   Int_t sl[2] ={slice,slice};
 422   DefineGoodTracks(sl,row[patch],good_number);
 423   SetMinPoints(min_points);
 424   AssignIDs();
 425   CreateHistos();
 426   FillEffHistos();
 427   CalcEffHistos();
 428 }
 429
 430 void AliL3Evaluate::EvaluateSlice(Int_t slice,Int_t min_points,Int_t good_number)
 431 {
 432   //Make efficiency plots for tracking on a slice (after merging).
 433   //min_points = minimum points on track to be considered for evaluation
 434   //good_number = minimum hits (padrows) produced by simulated track for consideration.
 435
 436   Int_t row[2] = {0,175};
 437   Int_t sl[2] ={slice,slice};
 438   DefineGoodTracks(sl,row,good_number);
 439
 440   SetMinPoints(min_points);
 441
 442   AssignIDs();
 443   CreateHistos();
 444   FillEffHistos();
 445   CalcEffHistos();
 446 }
 447
 448 void AliL3Evaluate::EvaluateGlobal(Int_t min_points,Int_t good_number,Char_t *fname)
 449 {
 450   //Make efficiency plots for tracking on several slices.
 451
 452   Int_t row[2] = {0,175};
 453   SetMinPoints(min_points);
 454   Int_t slice[2] = {fMinSlice,fMaxSlice};
 455   DefineGoodTracks(slice,row,good_number,fname);
 456   AssignIDs();
 457   CreateHistos(20,0,5);
 458   printf("filling histos\n");
 459   FillEffHistos();
 460   printf("done fillling\n");
 461   CalcEffHistos();
 462 }
 463
 464 void AliL3Evaluate::AssignIDs()
 465 {
 466   //Assign MC id to the tracks.
 467
 468   fTracks->QSort();
 469   LOG(AliL3Log::kDebug,"AliL3Evaluate::AssignIDs","Track Loop")
 470     <<"Assigning MC id to the found tracks...."<<ENDLOG;
 471   for(Int_t i=0; i<fTracks->GetNTracks(); i++)
 472     {
 473       AliL3Track *track = (AliL3Track*)fTracks->GetCheckedTrack(i);
 474       if(!track) continue;
 475       if(track->GetNumberOfPoints() < fMinPointsOnTrack) break;
 476
 477       fGoodFound++;
 478       Int_t tID = GetMCTrackLabel(track);
 479       track->SetMCid(tID);
 480       printf("track %i id %d nHits %d\n",i,tID,track->GetNumberOfPoints());
 481     }
 482 }
 483
 484
 485 struct S {Int_t lab; Int_t max;};
 486 Int_t AliL3Evaluate::GetMCTrackLabel(AliL3Track *track){
 487   //Returns the MCtrackID of the belonging clusters.
 488   //If MCLabel < 0, means that track is fake.
 489   //Fake track means that more than 10 percent of clusters are assigned incorrectly.
 490
 491   Int_t num_of_clusters = track->GetNumberOfPoints();
 492   S *s=new S[num_of_clusters];
 493   Int_t i;
 494   for (i=0; i<num_of_clusters; i++) s[i].lab=s[i].max=0;
 495   UInt_t *hitnum = track->GetHitNumbers();
 496   UInt_t id;
 497
 498   Int_t **trackID = GetClusterIDs(track);
 499
 500   Int_t lab=123456789;
 501   for (i=0; i<num_of_clusters; i++)
 502     {
 503       //Tricks to get the clusters belonging to this track:
 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       if(!points) continue;
 511       if(pos>=fNcl[slice][patch])
 512         {
 513           LOG(AliL3Log::kError,"AliL3Evaluate::GetMCTrackLabel","Clusterarray")
 514             <<AliL3Log::kDec<<"ERROR"<<ENDLOG;
 515           continue;
 516         }
 517
 518       //Get the label of the cluster:
 519       //printf("label %d %d %d\n",trackID[i][0],trackID[i][1],trackID[i][2]);
 520       lab=abs(trackID[i][0]);
 521       Int_t j;
 522       for (j=0; j<num_of_clusters; j++)
 523         if (s[j].lab==lab || s[j].max==0) break;
 524       s[j].lab=lab;
 525       s[j].max++;
 526     }
 527
 528   Int_t max=0;
 529   for (i=0; i<num_of_clusters; i++)
 530     if (s[i].max>max) {max=s[i].max; lab=s[i].lab;}
 531
 532   delete[] s;
 533
 534   for (i=0; i<num_of_clusters; i++)
 535     {
 536       id = hitnum[i];
 537       Int_t slice = (id>>25) & 0x7f;
 538       Int_t patch = (id>>22) & 0x7;
 539       UInt_t pos = id&0x3fffff;
 540
 541       AliL3SpacePointData *points = fClusters[slice][patch];
 542       if(!points) continue;
 543       if(pos>=fNcl[slice][patch])
 544         {
 545           LOG(AliL3Log::kError,"AliL3Evaluate::GetMCTrackLabel","Clusterarray")
 546             <<AliL3Log::kDec<<"ERROR"<<ENDLOG;
 547           continue;
 548         }
 549
 550       if (abs(trackID[i][1]) == lab ||
 551           abs(trackID[i][2]) == lab ) max++;
 552     }
 553
 554   //check if more than 10% of the clusters are incorrectly assigned (fake track):
 555   if (1.-Float_t(max)/num_of_clusters > 0.10)
 556     {
 557       return -lab;
 558     }
 559
 560   for(Int_t j=0; j<track->GetNumberOfPoints(); j++)
 561     delete [] trackID[j];
 562   delete [] trackID;
 563   return lab;
 564 }
 565
 566
 567 Int_t **AliL3Evaluate::GetClusterIDs(AliL3Track *track)
 568 {
 569   //Return the MC information of all clusters belonging to track.
 570
 571   Int_t num_of_clusters = track->GetNumberOfPoints();
 572   Int_t **trackID = new Int_t*[num_of_clusters];
 573
 574   UInt_t *hitnum = track->GetHitNumbers();
 575   UInt_t id;
 576
 577   Float_t xyz[3];
 578   Int_t padrow;
 579   for(Int_t i=0; i<num_of_clusters; i++)
 580     {
 581       id = hitnum[i];
 582       Int_t slice = (id>>25) & 0x7f;
 583       Int_t patch = (id>>22) & 0x7;
 584       UInt_t pos = id&0x3fffff;
 585
 586       AliL3SpacePointData *points = fClusters[slice][patch];
 587
 588       if(!points)
 589         {
 590           LOG(AliL3Log::kError,"AliL3Evaluate::GetClusterIDs","Clusterarray")
 591             <<"No points at slice "<<slice<<" patch "<<patch<<" pos "<<pos<<ENDLOG;
 592           continue;
 593         }
 594       if(pos>=fNcl[slice][patch])
 595         {
 596           LOG(AliL3Log::kError,"AliL3Evaluate::GetClusterIDs","Clusterarray")
 597             <<AliL3Log::kDec<<"ERROR"<<ENDLOG;
 598           continue;
 599         }
 600
 601       xyz[0] = points[pos].fX;
 602       xyz[1] = points[pos].fY;
 603       xyz[2] = points[pos].fZ;
 604       //sector = points[pos].fSector;
 605       padrow = points[pos].fPadRow;
 606       Int_t se,ro;
 607       fTransform->Slice2Sector(slice,padrow,se,ro);
 608       fTransform->Global2Raw(xyz,se,ro);
 609
 610       if(fIsSlow)
 611         {
 612           Int_t p = fRowid[slice][padrow];
 613           if(!fDigitsTree->GetEvent(p))
 614             LOG(AliL3Log::kError,"AliL3Evaluate::GetClusterIDs","Digits Tree")
 615               <<"Error reading digits tree"<<ENDLOG;
 616
 617           trackID[i] = new Int_t[3];
 618           trackID[i][0] = fDigits->GetTrackID((Int_t)rint(xyz[2]),(Int_t)rint(xyz[1]),0);
 619           trackID[i][1] = fDigits->GetTrackID((Int_t)rint(xyz[2]),(Int_t)rint(xyz[1]),1);
 620           trackID[i][2] = fDigits->GetTrackID((Int_t)rint(xyz[2]),(Int_t)rint(xyz[1]),2);
 621           //if(trackID[i][0]==6 || trackID[i][0]==32)
 622           //printf("trackID %d, padrow %d pad %d time %d\n",fDigits->GetTrackID((Int_t)rint(xyz[2]),(Int_t)rint(xyz[1]),0),padrow,(int)rint(xyz[1]),(int)rint(xyz[2]));
 623           /*if(trackID[i][0]<0)
 624             {
 625
 626             printf("trackID %d, padrow %d pad %d time %d\n",trackID[i][0],padrow,(int)rint(xyz[1]),(int)rint(xyz[2]));
 627             printf("on the side %d %d %d %d\n",fDigits->GetTrackID(((int)rint(xyz[2])-1),((int)rint(xyz[1])),0),fDigits->GetTrackID(((int)rint(xyz[2])+1),((int)rint(xyz[1])),0),fDigits->GetTrackID(((int)rint(xyz[2])),((int)rint(xyz[1])-1),0),fDigits->GetTrackID(((int)rint(xyz[2])),((int)rint(xyz[1])+1),0));
 628             }*/
 629         }
 630       else
 631         {
 632           Int_t tmp_pid=0;
 633           for(Int_t ii=0; ii<fNFastPoints; ii++)
 634             {
 635               tmp_pid = fMcId[ii];
 636               if(fMcIndex[ii] == id) break;
 637             }
 638           trackID[i] = new Int_t[3];
 639           trackID[i][0] = tmp_pid;
 640           trackID[i][1] = -1;
 641           trackID[i][2] = -1;
 642         }
 643     }
 644
 645   return trackID;
 646 }
 647
 648 void AliL3Evaluate::GetFastClusterIDs(Char_t *path)
 649 {
 650   //Get the MC id of space points in case of using the fast simulator.
 651   char fname[256];
 652   sprintf(fname,"%s/point_mc.dat",path);
 653   FILE *infile = fopen(fname,"r");
 654   if(!infile) return;
 655   Int_t hitid,hitmc,i;
 656
 657   for(i=0; ; i++)
 658     if(fscanf(infile,"%d %d",&hitid,&hitmc)==EOF) break;
 659   rewind(infile);
 660   fNFastPoints = i;
 661   fMcId = new Int_t[fNFastPoints];
 662   fMcIndex = new UInt_t[fNFastPoints];
 663
 664   for(i=0; i<fNFastPoints; i++)
 665     {
 666       if(fscanf(infile,"%d %d",&hitid,&hitmc)==EOF) break;
 667       fMcId[i] = hitmc;
 668       fMcIndex[i] = hitid;
 669     }
 670   fclose(infile);
 671 }
 672
 673 void AliL3Evaluate::CreateHistos(Int_t nbin,Int_t xlow,Int_t xup)
 674 {
 675   //Create the histograms
 676
 677   LOG(AliL3Log::kInformational,"AliL3Evaluate::CreateHistos","Allocating")
 678     <<"Creating histograms..."<<ENDLOG;
 679
 680   fNtuppel = new TNtuple("fNtuppel","Pt resolution","pt_gen:pt_found:nHits");
 681
 682   fPtRes = new TH1F("fPtRes","Relative Pt resolution",30,-10.,10.);
 683   fNGoodTracksPt = new TH1F("fNGoodTracksPt","Good tracks vs pt",nbin,xlow,xup);
 684   fNFoundTracksPt = new TH1F("fNFoundTracksPt","Found tracks vs pt",nbin,xlow,xup);
 685   fNFakeTracksPt = new TH1F("fNFakeTracksPt","Fake tracks vs pt",nbin,xlow,xup);
 686   fTrackEffPt = new TH1F("fTrackEffPt","Tracking efficiency vs pt",nbin,xlow,xup);
 687   fFakeTrackEffPt = new TH1F("fFakeTrackEffPt","Efficiency for fake tracks vs pt",nbin,xlow,xup);
 688
 689   fNGoodTracksEta = new TH1F("fNGoodTracksEta","Good tracks vs eta",20,-50,50);
 690   fNFoundTracksEta = new TH1F("fNFoundTracksEta","Found tracks vs eta",20,-50,50);
 691   fNFakeTracksEta = new TH1F("fNFakeTracksEta","Fake tracks vs eta",20,-50,50);
 692   fTrackEffEta = new TH1F("fTrackEffEta","Tracking efficienct vs eta",20,-50,50);
 693   fFakeTrackEffEta = new TH1F("fFakeTrackEffEta","Efficiency for fake tracks vs eta",20,-50,50);
 694   printf("finished creating histos\n");
 695 }
 696
 697 void AliL3Evaluate::FillEffHistos()
 698 {
 699   //Fill the efficiency histograms.
 700
 701   for(Int_t i=0; i<fGoodGen; i++)
 702     {
 703       Double_t ptg=fGoodTracks[i].pt,pzg=fGoodTracks[i].pz;
 704       Float_t dipangle=TMath::ATan2(pzg,ptg)*180./TMath::Pi();
 705       printf("filling particle with pt %f and dipangle %f\n",ptg,dipangle);
 706       fNGoodTracksPt->Fill(ptg);
 707       fNGoodTracksEta->Fill(dipangle);
 708       Int_t found = 0;
 709       for(Int_t k=0; k<fTracks->GetNTracks(); k++)
 710         {
 711           AliL3Track *track = fTracks->GetCheckedTrack(k);
 712           if(!track) continue;
 713           Int_t nHits = track->GetNumberOfPoints();
 714           if(nHits < fMinPointsOnTrack) break;
 715
 716           Int_t tracklabel;
 717           tracklabel = track->GetMCid();
 718           printf("evaluating track id %d\n",tracklabel);
 719           if(TMath::Abs(tracklabel) != fGoodTracks[i].label) continue;
 720           found=1;
 721           if(tracklabel == fGoodTracks[i].label) {fNFoundTracksPt->Fill(ptg); fNFoundTracksEta->Fill(dipangle);}
 722           else {fNFakeTracksPt->Fill(ptg); fNFakeTracksEta->Fill(dipangle);}
 723           Float_t pt=track->GetPt();
 724           fPtRes->Fill((pt-ptg)/ptg*100.);
 725           fNtuppel->Fill(ptg,pt,nHits);
 726           break;
 727
 728         }
 729     }
 730 }
 731
 732 void AliL3Evaluate::CalcEffHistos(){
 733
 734   Stat_t ngood=fNGoodTracksPt->GetEntries();
 735   Stat_t nfound=fNFoundTracksPt->GetEntries();
 736   Stat_t nfake=fNFakeTracksPt->GetEntries();
 737
 738   LOG(AliL3Log::kInformational,"AliL3Evaluate::FillEffHistos","Efficiency")
 739     <<AliL3Log::kDec<<"There was "<<ngood<<" generated good tracks"<<ENDLOG;
 740   LOG(AliL3Log::kInformational,"AliL3Evaluate::FillEffHistos","Efficiency")
 741     <<AliL3Log::kDec<<"Found "<<nfound<<" tracks"<<ENDLOG;
 742   LOG(AliL3Log::kInformational,"AliL3Evaluate::FillEffHistos","Efficiency")
 743     <<AliL3Log::kDec<<"Integral efficiency is about "<<nfound/ngood*100<<ENDLOG;
 744   LOG(AliL3Log::kInformational,"AliL3Evaluate::FillEffHistos","Efficiency")
 745     <<AliL3Log::kDec<<"Fake tracks relative is about "<<nfake/ngood*100<<ENDLOG;
 746
 747   fNFoundTracksPt->Sumw2(); fNGoodTracksPt->Sumw2();
 748   fTrackEffPt->Divide(fNFoundTracksPt,fNGoodTracksPt,1,1.,"b");
 749   fFakeTrackEffPt->Divide(fNFakeTracksPt,fNGoodTracksPt,1,1.,"b");
 750   fTrackEffPt->SetMaximum(1.4);
 751   fTrackEffPt->SetXTitle("P_{T} [GeV]");
 752   fTrackEffPt->SetLineWidth(2);
 753   fFakeTrackEffPt->SetFillStyle(3013);
 754   fTrackEffPt->SetLineColor(4);
 755   fFakeTrackEffPt->SetFillColor(2);
 756
 757   fNFoundTracksEta->Sumw2(); fNGoodTracksEta->Sumw2();
 758   fTrackEffEta->Divide(fNFoundTracksEta,fNGoodTracksEta,1,1.,"b");
 759   fFakeTrackEffEta->Divide(fNFakeTracksEta,fNGoodTracksEta,1,1.,"b");
 760   fTrackEffEta->SetMaximum(1.4);
 761   fTrackEffEta->SetXTitle("#lambda [degrees]");
 762   fTrackEffEta->SetLineWidth(2);
 763   fFakeTrackEffEta->SetFillStyle(3013);
 764   fTrackEffEta->SetLineColor(4);
 765   fFakeTrackEffEta->SetFillColor(2);
 766
 767 }
 768
 769 void AliL3Evaluate::Write2File(Char_t *outputfile)
 770 {
 771   //Write histograms to file:
 772
 773   TFile *of = new TFile(outputfile,"RECREATE");
 774   if(!of->IsOpen())
 775     {
 776       LOG(AliL3Log::kError,"AliL3Evaluate::Write2File","File Open")
 777         <<"Problems opening rootfile"<<ENDLOG;
 778       return;
 779     }
 780
 781   of->cd();
 782   fNtuppel->Write();
 783   fPtRes->Write();
 784   fNGoodTracksPt->Write();
 785   fNFoundTracksPt->Write();
 786   fNFakeTracksPt->Write();
 787   fTrackEffPt->Write();
 788   fFakeTrackEffPt->Write();
 789   fNGoodTracksEta->Write();
 790   fNFoundTracksEta->Write();
 791   fNFakeTracksEta->Write();
 792   fTrackEffEta->Write();
 793   fFakeTrackEffEta->Write();
 794
 795   of->Close();
 796   delete of;
 797
 798 }
 799
 800 TNtuple *AliL3Evaluate::CalculateResiduals()
 801 {
 802
 803   TNtuple *ntuppel=new TNtuple("ntuppel","Residuals","residual_trans:residual_long:zHit:pt:dipangle:beta:padrow:nHits");
 804
 805   for(int f=fMinSlice; f<=fMaxSlice; f++)
 806     {
 807       AliL3FileHandler *tfile = new AliL3FileHandler();
 808       char fname[256];
 809       sprintf(fname,"tracks_tr_%d_0.raw",f);
 810       if(!tfile->SetBinaryInput(fname)){
 811         LOG(AliL3Log::kError,"AliL3Evaluation::Setup","File Open")
 812           <<"Inputfile "<<fname<<" does not exist"<<ENDLOG;
 813         return 0;
 814       }
 815       fTracks = new AliL3TrackArray();
 816       tfile->Binary2TrackArray(fTracks);
 817       tfile->CloseBinaryInput();
 818       delete tfile;
 819       printf("Looking in slice %d\n",f);
 820       for(Int_t i=0; i<fTracks->GetNTracks(); i++)
 821         {
 822
 823           AliL3Track *track = (AliL3Track*)fTracks->GetCheckedTrack(i);
 824           if(!track) continue;
 825
 826           track->CalculateHelix();
 827           UInt_t *hitnum = track->GetHitNumbers();
 828           UInt_t id;
 829
 830           Float_t xyz[3];
 831           Int_t padrow;
 832           for(Int_t j=0; j<track->GetNumberOfPoints()-1; j++)
 833             {
 834               id = hitnum[j];
 835               Int_t slice = (id>>25) & 0x7f;
 836               Int_t patch = (id>>22) & 0x7;
 837               UInt_t pos = id&0x3fffff;
 838
 839               //if(slice!=1) continue;
 840
 841               AliL3SpacePointData *points = fClusters[slice][patch];
 842
 843               if(!points)
 844                 {
 845                   LOG(AliL3Log::kError,"AliL3Evaluate::CalculateResiduals","Clusterarray")
 846                     <<"No points at slice "<<slice<<" patch "<<patch<<" pos "<<pos<<ENDLOG;
 847                   continue;
 848                 }
 849               if(pos>=fNcl[slice][patch])
 850                 {
 851                   LOG(AliL3Log::kError,"AliL3Evaluate::CalculateResiduals","Clusterarray")
 852                     <<AliL3Log::kDec<<"ERROR"<<ENDLOG;
 853                   continue;
 854                 }
 855
 856               xyz[0] = points[pos].fX;
 857               xyz[1] = points[pos].fY;
 858               xyz[2] = points[pos].fZ;
 859               padrow = points[pos].fPadRow;
 860               //fTransform->Global2Local(xyz,slice);
 861
 862               Float_t xyz_cross[3];
 863               track->GetCrossingPoint(padrow,xyz_cross);
 864               Double_t beta = track->GetCrossingAngle(padrow);
 865
 866               Double_t yres = xyz_cross[1] - xyz[1];
 867               Double_t zres = xyz_cross[2] - xyz[2];
 868
 869               Double_t dipangle = atan(track->GetTgl());
 870               ntuppel->Fill(yres,zres,xyz_cross[2],track->GetPt(),dipangle,beta,padrow,track->GetNumberOfPoints());
 871
 872             }
 873         }
 874       if(fTracks)
 875         delete fTracks;
 876     }
 877   return ntuppel;
 878 }
 879
 880 TNtuple *AliL3Evaluate::EvaluatePoints(Char_t *rootfile)
 881 {
 882   //Compare points to the exact crossing points of track and padrows.
 883   //The input file to this function, contains the exact clusters calculated
 884   //in AliTPC::Hits2ExactClusters.
 885
 886
 887   TNtuple *ntuppel = new TNtuple("ntuppel","residuals","slice:padrow:resy:resz:zHit:pt");
 888
 889   TFile *exfile = new TFile(rootfile);
 890
 891   AliTPCParam *param = (AliTPCParam*)exfile->Get("75x40_100x60");
 892
 893   //Get the exact clusters from file:
 894   AliTPCClustersArray *arr = new AliTPCClustersArray;
 895   arr->Setup(param);
 896   arr->SetClusterType("AliComplexCluster");
 897   char treeName[500];
 898   sprintf(treeName,"TreeCExact_%s",param->GetTitle());
 899   Bool_t clusterok = arr->ConnectTree(treeName);//Segment Tree (for offline clusters)
 900   if(!clusterok) {printf("AliL3Evaluate::EvaluatePoints : Error in clusterloading\n"); return 0;}
 901
 902   for(Int_t i=0; i<arr->GetTree()->GetEntries(); i++)
 903     {
 904       //Get the exact clusters for this row:
 905       Int_t cursec,currow;
 906       AliSegmentID *s = arr->LoadEntry(i);
 907       param->AdjustSectorRow(s->GetID(),cursec,currow);
 908
 909       AliTPCClustersRow *ro = (AliTPCClustersRow *)arr->GetRow(cursec,currow);
 910       TClonesArray *clusters = ro->GetArray();
 911       int num_of_offline=clusters->GetEntriesFast();
 912
 913       //Get the found clusters:
 914       Int_t slice,padrow;
 915       fTransform->Sector2Slice(slice,padrow,cursec,currow);
 916       if(slice<fMinSlice || slice>fMaxSlice) continue;
 917       AliL3SpacePointData *points = fClusters[slice][0];
 918
 919       Int_t index = fRowid[slice][padrow];
 920       if(!fDigitsTree->GetEvent(index))
 921         printf("AliL3Evaluate::EvaluatePoints : ERROR IN DIGITSTREE\n");
 922       printf("Checking slice %d padrow %d with %d clusters\n",slice,padrow,num_of_offline);
 923
 924       for(UInt_t c=0; c<fNcl[slice][0]; c++)
 925         {
 926           if(points[c].fPadRow!=padrow) continue;
 927           for(Int_t m=0; m<num_of_offline; m++)
 928             {
 929               AliComplexCluster *cluster = (AliComplexCluster *)clusters->UncheckedAt(m);
 930               Int_t mcId = cluster->fTracks[0];
 931               //Int_t mcId = cluster->GetLabel(0);
 932               if(mcId <0) continue;
 933               TParticle *part = gAlice->Particle(mcId);
 934
 935               Float_t xyz_cl[3] = {points[c].fX,points[c].fY,points[c].fZ};
 936               Float_t xyz_ex[3];
 937               fTransform->Global2Raw(xyz_cl,cursec,currow);
 938               if(fDigits->GetTrackID((Int_t)rint(xyz_cl[2]),(Int_t)rint(xyz_cl[1]),0)!=mcId &&
 939                  fDigits->GetTrackID((Int_t)rint(xyz_cl[2]),(Int_t)rint(xyz_cl[1]),1)!=mcId &&
 940                  fDigits->GetTrackID((Int_t)rint(xyz_cl[2]),(Int_t)rint(xyz_cl[1]),2)!=mcId)
 941                 continue;
 942               fTransform->Raw2Local(xyz_ex,cursec,currow,cluster->fY,cluster->fX);
 943               fTransform->Raw2Local(xyz_cl,cursec,currow,xyz_cl[1],xyz_cl[2]);
 944               Float_t resy = xyz_cl[1] - xyz_ex[1];//cluster->GetY()
 945               Float_t resz = xyz_cl[2] - xyz_ex[2];//cluster->GetZ()
 946
 947               ntuppel->Fill(slice,padrow,resy,resz,xyz_ex[2],part->Pt());
 948             }
 949         }
 950       arr->ClearRow(cursec,currow);
 951     }
 952
 953   return ntuppel;
 954 }
 955 /*
 956 TNtuple *AliL3Evaluate::EvaluateGEANT()
 957 {
 958   TNtuple *ntuppel = new TNtuple("ntuppel","residuals","resy:ptgen:padrow:zHit:slice");
 959
 960   Int_t good_number=175;
 961   Int_t row[2] = {0,175};
 962
 963   Float_t xyz_cross[3];
 964   Float_t xyz_cl[3];
 965
 966   for(Int_t slice=fMinSlice; slice<=fMaxSlice; slice++)
 967     {
 968       Int_t *particle_id = new Int_t[fParticles->GetEntriesFast()];
 969       TObjArray *good_particles = DefineGoodTracks(slice,row,good_number,particle_id);
 970
 971       if(fMCclusterfile)
 972         {
 973           fMCclusterfile->Close();
 974           delete fMCclusterfile;
 975           fMCclusterfile = new TFile("/nfs/david/subatom/alice/data/V3.04/fast/clusters/hg_8k_v0_s1-3_e0_cl.root");
 976           fMCclusterfile->cd();
 977         }
 978       AliTPCClustersArray arr;
 979       arr.Setup(fParam);
 980       arr.SetClusterType("AliTPCcluster");
 981       Bool_t clusterok = arr.ConnectTree("Segment Tree");
 982       if(!clusterok)
 983         LOG(AliL3Log::kError,"AliL3Evaluate::DefineGoodTracks","Cluster Array")
 984           <<"Error loading clusters from rootfile"<<ENDLOG;
 985
 986       for(Int_t ci=0; ci<arr.GetTree()->GetEntries(); ci++)
 987         {
 988           Int_t sec,row,sl,lr;
 989           AliSegmentID *s = arr.LoadEntry(ci);
 990           fParam->AdjustSectorRow(s->GetID(),sec,row);
 991           fTransform->Sector2Slice(sl,lr,sec,row);
 992           if(sl != slice) {arr.ClearRow(sec,row); continue;}
 993           //if(lr!=pr) {carray.ClearRow(sec,row); continue;}
 994           AliTPCClustersRow *cRow = arr.GetRow(sec,row);
 995           printf("evaluating slice %d row %d\n",slice,lr);
 996           for(Int_t j=0; j<cRow->GetArray()->GetEntriesFast(); j++)
 997             {
 998
 999               AliTPCcluster *cluster=(AliTPCcluster*)(*cRow)[j];
1000               xyz_cl[1] = cluster->GetY();
1001               xyz_cl[2] = cluster->GetZ();
1002               Int_t lab=cluster->GetLabel(0);
1003               if(lab<0) continue;
1004
1005               for(Int_t pa=0; pa<good_particles->GetEntriesFast(); pa++)
1006                 {
1007                   if(particle_id[pa]!=lab) continue;
1008                   TParticle *part = (TParticle*)good_particles->UncheckedAt(pa);
1009                   GetParticleCrossingPoint(part,slice,lr,xyz_cross);
1010                   Double_t yres = xyz_cl[1] - xyz_cross[1];
1011                   ntuppel->Fill(yres,part->Pt(),lr,xyz_cl[2],slice);
1012                 }
1013             }
1014           if(arr.GetRow(sec,row))
1015             arr.ClearRow(sec,row);
1016         }
1017       delete [] particle_id;
1018       delete good_particles;
1019     }
1020   return ntuppel;
1021 }
1022 */
1023 Bool_t AliL3Evaluate::GetParticleCrossingPoint(TParticle *part,Int_t slice,Int_t padrow,Float_t *xyz)
1024 {
1025   //Calcluate the crossing point between a generated particle and given padrow.
1026
1027
1028   Double_t kappa = 0.2*0.0029980/part->Pt();
1029
1030   Double_t radius = 1/fabs(kappa);
1031   if(part->GetPdgCode() > 0) kappa = -kappa;
1032
1033   Int_t charg;
1034   if(kappa>0)
1035     charg=-1;
1036   else
1037     charg=1;
1038
1039   Float_t angl[1] = {part->Phi()};
1040
1041   fTransform->Global2LocalAngle(angl,slice);
1042
1043   Double_t charge = -1.*kappa;
1044   Double_t trackPhi0 = angl[0] + charge*0.5*Pi/fabs(charge);
1045
1046   Double_t x0=0;
1047   Double_t y0=0;
1048   Double_t xc = x0 - radius * cos(trackPhi0);
1049   Double_t yc = y0 - radius * sin(trackPhi0);
1050
1051   //printf("radius %f xc %f yc %f\n",radius,xc,yc);
1052
1053   Double_t xHit = fTransform->Row2X(padrow);
1054   xyz[0] = xHit;
1055   Double_t aa = (xHit - xc)*(xHit - xc);
1056   Double_t r2 = radius*radius;
1057   if(aa > r2)
1058     return false;
1059
1060   Double_t aa2 = sqrt(r2 - aa);
1061   Double_t y1 = yc + aa2;
1062   Double_t y2 = yc - aa2;
1063   xyz[1] = y1;
1064   if(fabs(y2) < fabs(y1)) xyz[1] = y2;
1065
1066   Double_t tgl = part->Pz()/part->Pt();
1067
1068   Double_t yHit = xyz[1];
1069   Double_t angle1 = atan2((yHit - yc),(xHit - xc));
1070   if(angle1 < 0) angle1 += 2.*Pi;
1071   Double_t angle2 = atan2((0 - yc),(0 - xc));
1072   if(angle2 < 0) angle2 += 2.*Pi;
1073   Double_t diff_angle = angle1 - angle2;
1074   diff_angle = fmod(diff_angle,2*Pi);
1075   if((charg*diff_angle) > 0) diff_angle = diff_angle - charg*2.*Pi;
1076   Double_t s_tot = fabs(diff_angle)*radius;
1077   Double_t zHit = 0 + s_tot*tgl;
1078   xyz[2] = zHit;
1079
1080   return true;
1081 }
1082