HLT/src/AliL3Evaluate.cxx

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