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