HLT/src/AliL3Evaluate.cxx

   1 // @(#) $Id$
   2
   3 // Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
   4 //*-- Copyright &copy ALICE HLT Group
   5
   6 #include "AliL3StandardIncludes.h"
   7 #include <TFile.h>
   8 #include <TH1.h>
   9 #include <TParticle.h>
  10 #include <TTree.h>
  11 #include <TClonesArray.h>
  12 #include <TNtupleD.h>
  13
  14 #include <AliRun.h>
  15 #include <AliSimDigits.h>
  16 #include <AliTPC.h>
  17 #include <AliTPCcluster.h>
  18 #include <AliTPCClustersArray.h>
  19 #include <AliTPCClustersRow.h>
  20 #include <AliTPCParam.h>
  21 #include <AliComplexCluster.h>
  22
  23 #if GCCVERSION == 3
  24 #include <fstream>
  25 #include <iosfwd>
  26 #else
  27 #include <fstream.h>
  28 #endif
  29
  30 #include "AliL3Logging.h"
  31 #include "AliL3Transform.h"
  32 #include "AliL3SpacePointData.h"
  33 #include "AliL3Track.h"
  34 #include "AliL3FileHandler.h"
  35 #include "AliL3TrackArray.h"
  36 #include "AliL3Evaluate.h"
  37
  38 #if GCCVERSION == 3
  39 using namespace std;
  40 #endif
  41
  42 /** \class AliL3Evaluate
  43 <pre>
  44 //_____________________________________________________________
  45 // AliL3Evaluate
  46 //
  47 // Evaluation class for tracking. Plots efficiencies etc..
  48 //
  49 </pre>
  50 */
  51
  52 ClassImp(AliL3Evaluate)
  53
  54 AliL3Evaluate::AliL3Evaluate()
  55 {
  56   fTracks = 0;
  57   fNFastPoints = 0;
  58   fMcIndex = 0;
  59   fMcId = 0;
  60   fMinSlice=0;
  61   fMaxSlice=0;
  62   fGoodTracks=0;
  63 }
  64
  65 AliL3Evaluate::AliL3Evaluate(Char_t *datapath,Int_t min_clusters,Int_t minhits,Double_t minpt,Int_t *slice)
  66 {
  67
  68   if(slice)
  69     {
  70       fMinSlice=slice[0];
  71       fMaxSlice=slice[1];
  72     }
  73   else
  74     {
  75       fMinSlice=0;
  76       fMaxSlice=35;
  77     }
  78   sprintf(fPath,"%s",datapath);
  79   fMaxFalseClusters = 0.1;
  80   fGoodFound = 0;
  81   fGoodGen = 0;
  82   fMinPointsOnTrack = min_clusters;
  83   fMinHitsFromParticle = minhits;
  84   fMinGoodPt = minpt;
  85   memset(fClusters,0,36*6*sizeof(AliL3SpacePointData*));
  86   fTracks=0;
  87   fGoodTracks=0;
  88 }
  89
  90 void AliL3Evaluate::LoadData(Int_t event,Bool_t sp)
  91 {
  92   Char_t fname[1024];
  93   AliL3FileHandler *clusterfile[36][6];
  94   for(Int_t s=fMinSlice; s<=fMaxSlice; s++)
  95     {
  96       for(Int_t p=0; p<AliL3Transform::GetNPatches(); p++)
  97         {
  98           Int_t patch;
  99           if(sp==kTRUE)
 100             patch=-1;
 101           else
 102             patch=p;
 103
 104           if(fClusters[s][p])
 105             delete fClusters[s][p];
 106           fClusters[s][p] = 0;
 107           clusterfile[s][p] = new AliL3FileHandler();
 108           if(event<0)
 109             sprintf(fname,"%s/points_%d_%d.raw",fPath,s,patch);
 110           else
 111             sprintf(fname,"%s/points_%d_%d_%d.raw",fPath,event,s,patch);
 112           if(!clusterfile[s][p]->SetBinaryInput(fname))
 113             {
 114               LOG(AliL3Log::kError,"AliL3Evaluation::Setup","File Open")
 115                 <<"Inputfile "<<fname<<" does not exist"<<ENDLOG;
 116               delete clusterfile[s][p];
 117               clusterfile[s][p] = 0;
 118               continue;
 119             }
 120           fClusters[s][p] = (AliL3SpacePointData*)clusterfile[s][p]->Allocate();
 121           clusterfile[s][p]->Binary2Memory(fNcl[s][p],fClusters[s][p]);
 122           clusterfile[s][p]->CloseBinaryInput();
 123           if(sp==kTRUE)
 124             break;
 125         }
 126     }
 127
 128   sprintf(fname,"%s/tracks_%d.raw",fPath,event);
 129   AliL3FileHandler *tfile = new AliL3FileHandler();
 130   if(!tfile->SetBinaryInput(fname)){
 131     LOG(AliL3Log::kError,"AliL3Evaluation::Setup","File Open")
 132       <<"Inputfile "<<fname<<" does not exist"<<ENDLOG;
 133     return;
 134   }
 135   if(fTracks)
 136     delete fTracks;
 137   fTracks = new AliL3TrackArray();
 138   tfile->Binary2TrackArray(fTracks);
 139   tfile->CloseBinaryInput();
 140   delete tfile;
 141 }
 142
 143
 144 AliL3Evaluate::~AliL3Evaluate()
 145 {
 146   if(fGoodTracks) delete fGoodTracks;
 147   if(fDigitsTree) fDigitsTree->Delete();
 148   if(fTracks) delete fTracks;
 149   if(fPtRes) delete fPtRes;
 150   if(fNGoodTracksPt) delete fNGoodTracksPt;
 151   if(fNFoundTracksPt) delete fNFoundTracksPt;
 152   if(fNFakeTracksPt) delete fNFakeTracksPt;
 153   if(fTrackEffPt) delete fTrackEffPt;
 154   if(fFakeTrackEffPt) delete fFakeTrackEffPt;
 155   if(fNGoodTracksEta) delete fNGoodTracksEta;
 156   if(fNFoundTracksEta) delete fNFoundTracksEta;
 157   if(fNFakeTracksEta) delete fNFakeTracksEta;
 158   if(fTrackEffEta) delete fTrackEffEta;
 159   if(fFakeTrackEffEta) delete fFakeTrackEffEta;
 160   if(fMcIndex) delete [] fMcIndex;
 161   if(fMcId)    delete [] fMcId;
 162   if(fNtuppel) delete fNtuppel;
 163 }
 164
 165
 166 void AliL3Evaluate::AssignIDs()
 167 {
 168   //Assign MC id to the tracks.
 169 #ifndef do_mc
 170   cerr<<"AliL3Evaluate::AssignIDs() : You need to compile with the do_mc flag!"<<endl;
 171   return;
 172 #endif
 173   fGoodFound=0;
 174   fTracks->QSort();
 175   LOG(AliL3Log::kDebug,"AliL3Evaluate::AssignIDs","Track Loop")
 176     <<"Assigning MC id to the found tracks...."<<ENDLOG;
 177   for(Int_t i=0; i<fTracks->GetNTracks(); i++)
 178     {
 179       AliL3Track *track = (AliL3Track*)fTracks->GetCheckedTrack(i);
 180       if(!track) continue;
 181       if(track->GetNumberOfPoints() < fMinPointsOnTrack) break;
 182
 183       fGoodFound++;
 184       Int_t tID = GetMCTrackLabel(track);
 185       track->SetMCid(tID);
 186     }
 187   //cout<<"Found "<<fGoodFound<<" good tracks "<<endl;
 188 }
 189
 190
 191 struct S {Int_t lab; Int_t max;};
 192 Int_t AliL3Evaluate::GetMCTrackLabel(AliL3Track *track){
 193   //Returns the MCtrackID of the belonging clusters.
 194   //If MCLabel < 0, means that track is fake.
 195   //Definitions are identical to offline.
 196   //Fake track means:
 197   // - more than 10 percent of clusters are assigned incorrectly
 198   // - more than half of the innermost 10% of clusters were assigned incorrectly.
 199
 200
 201 #ifdef do_mc
 202   Int_t num_of_clusters = track->GetNumberOfPoints();
 203   S *s=new S[num_of_clusters];
 204   Int_t i;
 205   for (i=0; i<num_of_clusters; i++) s[i].lab=s[i].max=0;
 206   UInt_t *hitnum = track->GetHitNumbers();
 207   UInt_t id;
 208
 209   Int_t lab=123456789;
 210   for (i=0; i<num_of_clusters; i++)
 211     {
 212       //Tricks to get the clusters belonging to this track:
 213       id = hitnum[i];
 214       Int_t slice = (id>>25) & 0x7f;
 215       Int_t patch = (id>>22) & 0x7;
 216       UInt_t pos = id&0x3fffff;
 217
 218       AliL3SpacePointData *points = fClusters[slice][patch];
 219       if(!points) continue;
 220       if(pos>=fNcl[slice][patch])
 221         {
 222           LOG(AliL3Log::kError,"AliL3Evaluate::GetMCTrackLabel","Clusterarray")
 223             <<AliL3Log::kDec<<"ERROR"<<ENDLOG;
 224           continue;
 225         }
 226
 227       //Get the label of the cluster:
 228       lab=points[pos].fTrackID[0];
 229
 230       Int_t j;
 231       for (j=0; j<num_of_clusters; j++)
 232         if (s[j].lab==lab || s[j].max==0) break;
 233       s[j].lab=lab;
 234       s[j].max++;
 235     }
 236
 237   Int_t max=0;
 238   for (i=0; i<num_of_clusters; i++)
 239     if (s[i].max>max) {max=s[i].max; lab=s[i].lab;}
 240
 241   if(lab == -1)
 242     return -1; //If most clusters is -1, this is a noise track.
 243   if(lab < 0)
 244     cerr<<"AliL3Evaluate::GetMCTrackLabel : Track label negative :"<<lab<<endl;
 245
 246   delete[] s;
 247
 248   for (i=0; i<num_of_clusters; i++)
 249     {
 250       id = hitnum[i];
 251       Int_t slice = (id>>25) & 0x7f;
 252       Int_t patch = (id>>22) & 0x7;
 253       UInt_t pos = id&0x3fffff;
 254
 255       AliL3SpacePointData *points = fClusters[slice][patch];
 256       if(!points) continue;
 257       if(pos>=fNcl[slice][patch])
 258         {
 259           LOG(AliL3Log::kError,"AliL3Evaluate::GetMCTrackLabel","Clusterarray")
 260             <<AliL3Log::kDec<<"ERROR"<<ENDLOG;
 261           continue;
 262         }
 263
 264       if (abs(points[pos].fTrackID[1]) == lab ||
 265           abs(points[pos].fTrackID[2]) == lab ) max++;
 266     }
 267
 268
 269   //Check if more than 10% of the clusters were assigned incorrectly:
 270   if (1.-Float_t(max)/num_of_clusters > fMaxFalseClusters)
 271     {
 272       return -lab;
 273     }
 274   else //Check if at least half of the 10% innermost clusters are assigned correctly.
 275     {
 276       Int_t tail=Int_t(0.10*num_of_clusters);
 277       max=0;
 278       for (i=1; i<=tail; i++)
 279         {
 280           id = hitnum[num_of_clusters - i];
 281           Int_t slice = (id>>25) & 0x7f;
 282           Int_t patch = (id>>22) & 0x7;
 283           UInt_t pos = id&0x3fffff;
 284
 285           AliL3SpacePointData *points = fClusters[slice][patch];
 286           if(lab == abs(points[pos].fTrackID[0]) ||
 287              lab == abs(points[pos].fTrackID[1]) ||
 288              lab == abs(points[pos].fTrackID[2])) max++;
 289         }
 290       if (max < Int_t(0.5*tail)) return -lab;
 291     }
 292
 293   return lab;
 294 #else //If we are running with mc_ids or not
 295   return 0;
 296 #endif
 297
 298 }
 299
 300 void AliL3Evaluate::GetFastClusterIDs(Char_t *path)
 301 {
 302   //Get the MC id of space points in case of using the fast simulator.
 303   char fname[256];
 304   sprintf(fname,"%s/point_mc.dat",path);
 305   FILE *infile = fopen(fname,"r");
 306   if(!infile) return;
 307   Int_t hitid,hitmc,i;
 308
 309   for(i=0; ; i++)
 310     if(fscanf(infile,"%d %d",&hitid,&hitmc)==EOF) break;
 311   rewind(infile);
 312   fNFastPoints = i;
 313   fMcId = new Int_t[fNFastPoints];
 314   fMcIndex = new UInt_t[fNFastPoints];
 315
 316   for(i=0; i<fNFastPoints; i++)
 317     {
 318       if(fscanf(infile,"%d %d",&hitid,&hitmc)==EOF) break;
 319       fMcId[i] = hitmc;
 320       fMcIndex[i] = hitid;
 321     }
 322   fclose(infile);
 323 }
 324
 325 void AliL3Evaluate::CreateHistos(Int_t nbin,Float_t xlow,Float_t xup)
 326 {
 327   //Create the histograms
 328
 329   LOG(AliL3Log::kInformational,"AliL3Evaluate::CreateHistos","Allocating")
 330     <<"Creating histograms..."<<ENDLOG;
 331
 332   fNtuppel = new TNtuple("fNtuppel","Pt resolution","pt_gen:pt_found:nHits");
 333   fNtuppel->SetDirectory(0);
 334   fPtRes = new TH1F("fPtRes","Relative Pt resolution",30,-10.,10.);
 335   fNGoodTracksPt = new TH1F("fNGoodTracksPt","Good tracks vs pt",nbin,xlow,xup);
 336   fNFoundTracksPt = new TH1F("fNFoundTracksPt","Found tracks vs pt",nbin,xlow,xup);
 337   fNFakeTracksPt = new TH1F("fNFakeTracksPt","Fake tracks vs pt",nbin,xlow,xup);
 338   fTrackEffPt = new TH1F("fTrackEffPt","Tracking efficiency vs pt",nbin,xlow,xup);
 339   fFakeTrackEffPt = new TH1F("fFakeTrackEffPt","Efficiency for fake tracks vs pt",nbin,xlow,xup);
 340
 341   fNGoodTracksEta = new TH1F("fNGoodTracksEta","Good tracks vs eta",20,-50,50);
 342   fNFoundTracksEta = new TH1F("fNFoundTracksEta","Found tracks vs eta",20,-50,50);
 343   fNFakeTracksEta = new TH1F("fNFakeTracksEta","Fake tracks vs eta",20,-50,50);
 344   fTrackEffEta = new TH1F("fTrackEffEta","Tracking efficienct vs eta",20,-50,50);
 345   fFakeTrackEffEta = new TH1F("fFakeTrackEffEta","Efficiency for fake tracks vs eta",20,-50,50);
 346
 347 }
 348
 349 void AliL3Evaluate::GetGoodParticles(Char_t *path,Int_t event,Int_t *padrowrange)
 350 {
 351   //Read the good particles from file. This file should already have been
 352   //generated by macro AliTPCComparison.C.
 353
 354   Char_t filename[1024];
 355   if(event<0 && !padrowrange)
 356     sprintf(filename,"%s/good_tracks_tpc",path);
 357   else if(event>=0 && !padrowrange)
 358     sprintf(filename,"%s/good_tracks_tpc_%d",path,event);
 359   else
 360     sprintf(filename,"%s/good_tracks_tpc_%d_%d_%d",path,event,padrowrange[0],padrowrange[1]);
 361   ifstream in(filename);
 362   if(!in)
 363     {
 364       cerr<<"AliL3Evaluate::GetGoodParticles : Problems opening file :"<<filename<<endl;
 365       return;
 366     }
 367   Int_t MaxTracks=20000;
 368   if(fGoodTracks)
 369     delete [] fGoodTracks;
 370   fGoodGen=0;
 371   fGoodTracks = new GoodTrack[MaxTracks];
 372
 373   while (in>>fGoodTracks[fGoodGen].label>>fGoodTracks[fGoodGen].code>>
 374          fGoodTracks[fGoodGen].px>>fGoodTracks[fGoodGen].py>>fGoodTracks[fGoodGen].pz>>
 375          fGoodTracks[fGoodGen].x>>fGoodTracks[fGoodGen].y >>fGoodTracks[fGoodGen].z>>fGoodTracks[fGoodGen].nhits>>fGoodTracks[fGoodGen].sector)
 376
 377     {
 378       fGoodGen++;
 379       if (fGoodGen==MaxTracks)
 380         {
 381           cerr<<"AliL3Evaluate::GetGoodParticles : Too many good tracks !\n";
 382           break;
 383         }
 384     }
 385 }
 386
 387
 388 void AliL3Evaluate::FillEffHistos()
 389 {
 390   if(!fGoodTracks)
 391     {
 392       cerr<<"AliL3Evaluate::FillEffHistos : No good tracks"<<endl;
 393       return;
 394     }
 395   //cout<<"Comparing "<<fGoodGen<<" good tracks ..."<<endl;
 396   for(Int_t i=0; i<fGoodGen; i++)
 397     {
 398       //cout<<"Checking particle "<<i<<endl;
 399       if(fGoodTracks[i].nhits < fMinHitsFromParticle) continue;
 400       Float_t ptg = TMath::Sqrt(fGoodTracks[i].px*fGoodTracks[i].px + fGoodTracks[i].py*fGoodTracks[i].py);
 401       if(ptg < fMinGoodPt) continue;
 402       Float_t pzg=fGoodTracks[i].pz;
 403       Float_t dipangle=TMath::ATan2(pzg,ptg)*180./TMath::Pi();
 404
 405       fNGoodTracksPt->Fill(ptg);
 406       fNGoodTracksEta->Fill(dipangle);
 407       Int_t found = 0;
 408
 409       for(Int_t k=0; k<fTracks->GetNTracks(); k++)
 410         {
 411           AliL3Track *track = fTracks->GetCheckedTrack(k);
 412           if(!track) continue;
 413           Int_t nHits = track->GetNumberOfPoints();
 414           if(nHits < fMinPointsOnTrack) break;
 415           Int_t tracklabel;
 416           tracklabel = track->GetMCid();
 417
 418           if(TMath::Abs(tracklabel) != fGoodTracks[i].label) continue;
 419           found=1;
 420           Float_t pt=track->GetPt();
 421           if(tracklabel == fGoodTracks[i].label)
 422             {
 423               fNFoundTracksPt->Fill(ptg);
 424               fNFoundTracksEta->Fill(dipangle);
 425               fNtuppel->Fill(ptg,pt,nHits);
 426               fPtRes->Fill((pt-ptg)/ptg*100.);
 427             }
 428           else
 429             {
 430               fNFakeTracksPt->Fill(ptg);
 431               fNFakeTracksEta->Fill(dipangle);
 432             }
 433           //fPtRes->Fill((pt-ptg)/ptg*100.);
 434           //fNtuppel->Fill(ptg,pt,nHits);
 435           break;
 436
 437         }
 438       //if(!found)
 439       //        cout<<"Track "<<fGoodTracks[i].label<<" was not found"<<endl;
 440     }
 441 }
 442
 443 void AliL3Evaluate::FillEffHistosNAIVE()
 444 {
 445   //Fill the efficiency histograms.
 446
 447   cout<<endl<<"Note: Doing NAIVE evaluation "<<endl;
 448   for(Int_t i=0; i<fGoodGen; i++)
 449     {
 450       Double_t ptg=TMath::Sqrt(fGoodTracks[i].px*fGoodTracks[i].px + fGoodTracks[i].py*fGoodTracks[i].py);
 451       Double_t pzg=fGoodTracks[i].pz;
 452       Float_t dipangle=TMath::ATan2(pzg,ptg)*180./TMath::Pi();
 453       //printf("filling particle with pt %f and dipangle %f\n",ptg,dipangle);
 454       fNGoodTracksPt->Fill(ptg);
 455       fNGoodTracksEta->Fill(dipangle);
 456
 457     }
 458
 459   for(Int_t k=0; k<fTracks->GetNTracks(); k++)
 460     {
 461       AliL3Track *track = fTracks->GetCheckedTrack(k);
 462       if(!track) continue;
 463       Int_t nHits = track->GetNumberOfPoints();
 464       if(nHits < fMinPointsOnTrack) break;
 465       if(track->GetPt()<fMinGoodPt) continue;
 466       if(fabs(track->GetPseudoRapidity())>0.9) continue;
 467
 468       fNFoundTracksPt->Fill(track->GetPt()); fNFoundTracksEta->Fill(track->GetPseudoRapidity());
 469       //Float_t pt=track->GetPt();
 470       //fPtRes->Fill((pt-ptg)/ptg*100.);
 471       //fNtuppel->Fill(ptg,pt,nHits);
 472
 473     }
 474 }
 475
 476 void AliL3Evaluate::CalcEffHistos()
 477 {
 478
 479   Stat_t ngood=fNGoodTracksPt->GetEntries();
 480   Stat_t nfound=fNFoundTracksPt->GetEntries();
 481   Stat_t nfake=fNFakeTracksPt->GetEntries();
 482
 483   LOG(AliL3Log::kInformational,"AliL3Evaluate::FillEffHistos","Efficiency")
 484     <<AliL3Log::kDec<<"There was "<<ngood<<" generated good tracks"<<ENDLOG;
 485   LOG(AliL3Log::kInformational,"AliL3Evaluate::FillEffHistos","Efficiency")
 486     <<AliL3Log::kDec<<"Found "<<nfound<<" tracks"<<ENDLOG;
 487   LOG(AliL3Log::kInformational,"AliL3Evaluate::FillEffHistos","Efficiency")
 488     <<AliL3Log::kDec<<"Integral efficiency is about "<<nfound/ngood*100<<ENDLOG;
 489   LOG(AliL3Log::kInformational,"AliL3Evaluate::FillEffHistos","Efficiency")
 490     <<AliL3Log::kDec<<"Fake tracks relative is about "<<nfake/ngood*100<<ENDLOG;
 491   //LOG(AliL3Log::kInformational,"AliL3Evaluate::FillEffHistos","Efficiency")
 492   //<<AliL3Log::kDec<<"Naive efficiency "<<(Double_t)fGoodFound/(Double_t)fGoodGen<<ENDLOG;
 493
 494   fNFoundTracksPt->Sumw2(); fNGoodTracksPt->Sumw2();
 495   fTrackEffPt->Divide(fNFoundTracksPt,fNGoodTracksPt,1,1.,"b");
 496   fFakeTrackEffPt->Divide(fNFakeTracksPt,fNGoodTracksPt,1,1.,"b");
 497   fTrackEffPt->SetMaximum(1.4);
 498   fTrackEffPt->SetXTitle("P_{T} [GeV]");
 499   fTrackEffPt->SetLineWidth(2);
 500   fFakeTrackEffPt->SetFillStyle(3013);
 501   fTrackEffPt->SetLineColor(4);
 502   fFakeTrackEffPt->SetFillColor(2);
 503
 504   fNFoundTracksEta->Sumw2(); fNGoodTracksEta->Sumw2();
 505   fTrackEffEta->Divide(fNFoundTracksEta,fNGoodTracksEta,1,1.,"b");
 506   fFakeTrackEffEta->Divide(fNFakeTracksEta,fNGoodTracksEta,1,1.,"b");
 507   fTrackEffEta->SetMaximum(1.4);
 508   fTrackEffEta->SetXTitle("#lambda [degrees]");
 509   fTrackEffEta->SetLineWidth(2);
 510   fFakeTrackEffEta->SetFillStyle(3013);
 511   fTrackEffEta->SetLineColor(4);
 512   fFakeTrackEffEta->SetFillColor(2);
 513
 514 }
 515
 516 void AliL3Evaluate::Write2File(Char_t *outputfile)
 517 {
 518   //Write histograms to file:
 519
 520   TFile *of = TFile::Open(outputfile,"RECREATE");
 521   if(!of->IsOpen())
 522     {
 523       LOG(AliL3Log::kError,"AliL3Evaluate::Write2File","File Open")
 524         <<"Problems opening rootfile"<<ENDLOG;
 525       return;
 526     }
 527
 528   of->cd();
 529   fNtuppel->Write();
 530   fPtRes->Write();
 531   fNGoodTracksPt->Write();
 532   fNFoundTracksPt->Write();
 533   fNFakeTracksPt->Write();
 534   fTrackEffPt->Write();
 535   fFakeTrackEffPt->Write();
 536   fNGoodTracksEta->Write();
 537   fNFoundTracksEta->Write();
 538   fNFakeTracksEta->Write();
 539   fTrackEffEta->Write();
 540   fFakeTrackEffEta->Write();
 541
 542   of->Close();
 543
 544 }
 545
 546 TNtupleD *AliL3Evaluate::CalculateResiduals(Char_t *datapath)
 547 {
 548
 549   TNtupleD *ntuppel=new TNtupleD("ntuppel","Residuals","residual_trans:residual_long:zHit:pt:dipangle:beta:padrow:nHits");
 550   ntuppel->SetDirectory(0);
 551   if(fTracks)
 552     delete fTracks;
 553
 554   AliL3FileHandler *tfile = new AliL3FileHandler();
 555   char fname[256];
 556   sprintf(fname,"%s/tracks_0.raw",datapath);
 557   if(!tfile->SetBinaryInput(fname)){
 558     LOG(AliL3Log::kError,"AliL3Evaluation::Setup","File Open")
 559       <<"Inputfile "<<fname<<" does not exist"<<ENDLOG;
 560     return 0;
 561   }
 562   fTracks = new AliL3TrackArray();
 563   tfile->Binary2TrackArray(fTracks);
 564   tfile->CloseBinaryInput();
 565   delete tfile;
 566
 567   for(Int_t i=0; i<fTracks->GetNTracks(); i++)
 568     {
 569
 570       AliL3Track *track = (AliL3Track*)fTracks->GetCheckedTrack(i);
 571       if(!track) continue;
 572       if(track->GetNHits() < fMinPointsOnTrack) continue;
 573
 574       track->CalculateHelix();
 575       UInt_t *hitnum = track->GetHitNumbers();
 576       UInt_t id;
 577
 578       Float_t xyz[3];
 579       Int_t padrow;
 580       for(Int_t j=0; j<track->GetNumberOfPoints()-1; j++)
 581         {
 582           id = hitnum[j];
 583           Int_t slice = (id>>25) & 0x7f;
 584           Int_t patch = (id>>22) & 0x7;
 585           UInt_t pos = id&0x3fffff;
 586
 587           //if(slice!=1) continue;
 588
 589           AliL3SpacePointData *points = fClusters[slice][patch];
 590
 591           if(!points)
 592             {
 593               LOG(AliL3Log::kError,"AliL3Evaluate::CalculateResiduals","Clusterarray")
 594                 <<"No points at slice "<<slice<<" patch "<<patch<<" pos "<<pos<<ENDLOG;
 595               continue;
 596             }
 597           if(pos>=fNcl[slice][patch])
 598             {
 599               LOG(AliL3Log::kError,"AliL3Evaluate::CalculateResiduals","Clusterarray")
 600                 <<AliL3Log::kDec<<"ERROR"<<ENDLOG;
 601               continue;
 602             }
 603
 604           xyz[0] = points[pos].fX;
 605           xyz[1] = points[pos].fY;
 606           xyz[2] = points[pos].fZ;
 607           padrow = points[pos].fPadRow;
 608           AliL3Transform::Global2Local(xyz,slice,kTRUE);
 609
 610           Float_t angle = 0;
 611           AliL3Transform::Local2GlobalAngle(&angle,slice);
 612           track->CalculateReferencePoint(angle,AliL3Transform::Row2X(padrow));
 613           Float_t xyz_cross[3] = {track->GetPointX(),track->GetPointY(),track->GetPointZ()};
 614           AliL3Transform::Global2Local(xyz_cross,slice,kTRUE);
 615
 616           Double_t beta = track->GetCrossingAngle(padrow,slice);
 617
 618           Double_t yres = xyz_cross[1] - xyz[1];
 619           Double_t zres = xyz_cross[2] - xyz[2];
 620
 621           Double_t dipangle = atan(track->GetTgl());
 622           ntuppel->Fill(yres,zres,xyz_cross[2],track->GetPt(),dipangle,beta,padrow,track->GetNumberOfPoints());
 623
 624         }
 625     }
 626   if(fTracks)
 627     delete fTracks;
 628
 629   return ntuppel;
 630 }
 631
 632 enum tagprimary {kPrimaryCharged = 0x4000};
 633 void AliL3Evaluate::EvaluatePoints(Char_t *rootfile,Char_t *exactfile,Char_t *tofile,Int_t nevent,Bool_t offline,Bool_t sp)
 634 {
 635   //Compare points to the exact crossing points of track and padrows.
 636   //The input file to this function, contains the exact clusters calculated
 637   //in AliTPC::Hits2ExactClusters.
 638
 639 #ifndef do_mc
 640   cerr<<"AliL3Evaluate::EvaluatePoints : Compile with do_mc flag!"<<endl;
 641   return;
 642 #else
 643   cout<<"Evaluating points"<<endl;
 644   TNtuple *ntuppel = new TNtuple("ntuppel_res","Cluster properties",
 645                                  "slice:padrow:charge:resy:resz:zHit:pt:beta:sigmaY2:sigmaZ2:psigmaY2:psigmaZ2");
 646   ntuppel->SetDirectory(0);
 647
 648   TNtuple *ntuppel2 = new TNtuple("ntuppel_eff","Efficiency","slice:padrow:nfound:ngen");
 649   ntuppel2->SetDirectory(0);
 650
 651   TFile *exfile = TFile::Open(rootfile);
 652   if(!exfile)
 653     {
 654       cerr<<"Error opening rootfile "<<rootfile<<endl;
 655       return;
 656     }
 657   gAlice = (AliRun*)exfile->Get("gAlice");
 658   if (!gAlice)
 659     {
 660       LOG(AliL3Log::kError,"AliL3Evaluate::InitMC","gAlice")
 661         <<"AliRun object non existing on file"<<ENDLOG;
 662       return;
 663     }
 664
 665   AliTPCParam *param = (AliTPCParam*)exfile->Get(AliL3Transform::GetParamName());
 666
 667   TFile *exact = TFile::Open(exactfile);
 668   if(!exact)
 669     {
 670       cerr<<"AliL3Evaluate::EvaluatePoints : Problems opening file :"<<exactfile<<endl;
 671       return;
 672     }
 673
 674   AliTPCClustersArray *arr=0;
 675   for(Int_t event=0; event<nevent; event++)
 676     {
 677       LoadData(event,sp);
 678       exfile->cd();
 679       if(arr)
 680         delete arr;
 681       Int_t nparticles = gAlice->GetEvent(event);
 682       Int_t nprimaries = 0;//FindPrimaries(nparticles);
 683       cout<<"Event "<<event<<" had "<<nparticles<<" particles and "<<nprimaries<<" primaries"<<endl;
 684       exact->cd();
 685
 686       //Get the exact clusters from file:
 687       AliTPCClustersArray *arr = new AliTPCClustersArray;
 688       arr->Setup(param);
 689       arr->SetClusterType("AliComplexCluster");
 690       char treeName[500];
 691       sprintf(treeName,"TreeCExact_%s_%d",param->GetTitle(),event);
 692       Bool_t clusterok = arr->ConnectTree(treeName);//Segment Tree (for offline clusters)
 693       if(!clusterok) {printf("AliL3Evaluate::EvaluatePoints : Error in clusterloading\n"); return;}
 694
 695       //cout<<"Entering loop with "<<(Int_t)arr->GetTree()->GetEntries()<<endl;
 696       for(Int_t i=0; i<arr->GetTree()->GetEntries(); i++)
 697         {
 698           //Get the exact clusters for this row:
 699           Int_t cursec,currow;
 700           AliSegmentID *s = arr->LoadEntry(i);
 701           param->AdjustSectorRow(s->GetID(),cursec,currow);
 702
 703           AliTPCClustersRow *ro = (AliTPCClustersRow *)arr->GetRow(cursec,currow);
 704           TClonesArray *clusters = ro->GetArray();
 705           int num_of_offline=clusters->GetEntriesFast();
 706
 707           //Get the found clusters:
 708           Int_t slice,padrow;
 709           AliL3Transform::Sector2Slice(slice,padrow,cursec,currow);
 710           if(slice < fMinSlice) continue;
 711           if(slice > fMaxSlice) break;
 712
 713           Int_t patch = AliL3Transform::GetPatch(padrow);
 714           if(sp)
 715             patch=0;
 716           AliL3SpacePointData *points = fClusters[slice][patch];
 717           if(!points)
 718             continue;
 719
 720           //cout<<"Slice "<<slice<<" padrow "<<padrow<<" has "<<num_of_offline<<" clusters "<<endl;
 721           Int_t clustercount=0;
 722           Int_t crosscount=0;
 723           for(Int_t m=0; m<num_of_offline; m++)
 724             {
 725               AliComplexCluster *cluster = (AliComplexCluster *)clusters->UncheckedAt(m);
 726               Int_t mcId = cluster->fTracks[0];
 727
 728               if(mcId <0) continue;
 729
 730               if(cluster->fY < 1 || cluster->fY > AliL3Transform::GetNPads(padrow) - 2 ||
 731                  cluster->fX < 1 || cluster->fX > AliL3Transform::GetNTimeBins() - 2)
 732                 continue;
 733
 734               Float_t xyz_ex[3];
 735
 736               AliL3Transform::Raw2Local(xyz_ex,cursec,currow,cluster->fY,cluster->fX);
 737
 738               //In function AliTPC::Hits2ExactClusters the time offset is not included,
 739               //so we have to substract it again here.
 740               if(slice<18)
 741                 xyz_ex[2]-=AliL3Transform::GetZOffset();
 742               else
 743                 xyz_ex[2]+=AliL3Transform::GetZOffset();
 744
 745               //Outside our cone:
 746               if(param->GetPadRowRadii(cursec,currow)<230./250.*fabs(xyz_ex[2]))
 747                 continue;
 748
 749               TParticle *part = gAlice->Particle(mcId);
 750               crosscount++;
 751
 752               if(part->Pt() < fMinGoodPt) continue;
 753
 754               //Dont take secondaries, because in width calculation we assume primaries:
 755               //if(!(part->TestBit(kPrimaryCharged))) continue;
 756               if(part->GetFirstMother()>=0) continue;
 757
 758               Int_t tempcount=0;
 759               for(UInt_t c=0; c<fNcl[slice][patch]; c++)
 760                 {
 761                   if((Int_t)points[c].fPadRow!=padrow) continue;
 762                   Float_t xyz_cl[3] = {points[c].fX,points[c].fY,points[c].fZ};
 763
 764                   if(!offline)
 765                     AliL3Transform::Global2Local(xyz_cl,cursec);
 766                   tempcount++;
 767
 768                   if(points[c].fTrackID[0] != mcId &&
 769                      points[c].fTrackID[1] != mcId &&
 770                      points[c].fTrackID[2] != mcId)
 771                     continue;
 772
 773                   //Residuals:
 774                   Float_t resy = xyz_cl[1] - xyz_ex[1];
 775                   Float_t resz = xyz_cl[2] - xyz_ex[2];
 776
 777                   //Cluster shape
 778                   Int_t charge = (Int_t)points[c].fCharge;
 779                   Float_t beta = GetCrossingAngle(part,slice,padrow,xyz_ex);
 780                   Double_t tanl = xyz_ex[2]/sqrt(xyz_ex[0]*xyz_ex[0]+xyz_ex[1]*xyz_ex[1]);
 781                   Float_t psigmaY2 = AliL3Transform::GetParSigmaY2(padrow,xyz_ex[2],beta);
 782                   Float_t psigmaZ2 = AliL3Transform::GetParSigmaZ2(padrow,xyz_ex[2],tanl);
 783                   Float_t sigmaY2 = points[c].fSigmaY2;
 784                   Float_t sigmaZ2 = points[c].fSigmaZ2;
 785                   ntuppel->Fill(slice,padrow,charge,resy,resz,xyz_ex[2],part->Pt(),beta,sigmaY2,sigmaZ2,psigmaY2,psigmaZ2);
 786                 }
 787               clustercount=tempcount;
 788             }
 789           ntuppel2->Fill(slice,padrow,clustercount,crosscount);
 790           arr->ClearRow(cursec,currow);
 791         }
 792     }
 793   exfile->Close();
 794   exact->Close();
 795
 796   TFile *ofile = TFile::Open(tofile,"RECREATE");
 797   ntuppel->Write();
 798   ntuppel2->Write();
 799   ofile->Close();
 800
 801 #endif
 802 }
 803
 804 void AliL3Evaluate::GetCFeff(Char_t *path,Char_t *outfile,Int_t nevent,Bool_t sp)
 805 {
 806   //Evaluate the cluster finder efficiency.
 807
 808 #ifndef do_mc
 809   cerr<<"AliL3Evaluate::GetCFeff : Compile with do_mc flag"<<endl;
 810   return;
 811 #else
 812   TNtuple *ntuppel = new TNtuple("ntuppel","Cluster finder efficiency","slice:row:ncrossings:nclusters");
 813   ntuppel->SetDirectory(0);
 814
 815   Char_t filename[1024];
 816   sprintf(filename,"%s/alirunfile.root",path);
 817   TFile *rfile = TFile::Open(filename);
 818   gAlice = (AliRun*)rfile->Get("gAlice");
 819
 820   AliTPCParam *param = (AliTPCParam*)rfile->Get(AliL3Transform::GetParamName());
 821
 822   Int_t zero=param->GetZeroSup();
 823
 824   sprintf(filename,"%s/digitfile.root",path);
 825   TFile *dfile = TFile::Open(filename);
 826
 827   for(Int_t event=0; event<nevent; event++)
 828     {
 829       LoadData(event,sp);
 830       rfile->cd();
 831       gAlice->GetEvent(event);
 832       Int_t np = gAlice->GetNtrack();
 833       cout<<"Processing event "<<event<<" with "<<np<<" particles "<<endl;
 834       dfile->cd();
 835       sprintf(filename,"TreeD_75x40_100x60_150x60_%d",event);
 836       TTree *TD=(TTree*)gDirectory->Get(filename);
 837       AliSimDigits da, *digits=&da;
 838       TD->GetBranch("Segment")->SetAddress(&digits);
 839
 840       Int_t crossed=0,recs=0;
 841       Int_t *count = new Int_t[np]; //np number of particles.
 842       Int_t i;
 843       Float_t xyz[3];
 844       for (i=0; i<np; i++) count[i]=0;
 845
 846
 847       Int_t sec,row,sl,sr;
 848       for(Int_t i=0; i<(Int_t)TD->GetEntries(); i++)
 849         {
 850           crossed=recs=0;
 851           if (!TD->GetEvent(i)) continue;
 852           param->AdjustSectorRow(digits->GetID(),sec,row);
 853           AliL3Transform::Sector2Slice(sl,sr,sec,row);
 854           if(sl < fMinSlice) continue;
 855           if(sl > fMaxSlice) break;
 856           cout<<"Processing slice "<<sl<<" row "<<sr<<endl;
 857           digits->First();
 858           do {
 859             Int_t it=digits->CurrentRow(), ip=digits->CurrentColumn();
 860             Short_t dig = digits->GetDigit(it,ip);
 861
 862             if(dig<=param->GetZeroSup()) continue;
 863             AliL3Transform::Raw2Local(xyz,sec,row,ip,it);
 864             if(param->GetPadRowRadii(sec,row)<230./250.*fabs(xyz[2]))
 865               continue;
 866
 867             Int_t idx0=digits->GetTrackID(it,ip,0);
 868             Int_t idx1=digits->GetTrackID(it,ip,1);
 869             Int_t idx2=digits->GetTrackID(it,ip,2);
 870
 871             if (idx0>=0 && dig>=zero) count[idx0]+=1;
 872             if (idx1>=0 && dig>=zero) count[idx1]+=1;
 873             if (idx2>=0 && dig>=zero) count[idx2]+=1;
 874           } while (digits->Next());
 875           for (Int_t j=0; j<np; j++)
 876             {
 877               TParticle *part = gAlice->Particle(j);
 878               if(part->Pt() < fMinGoodPt) continue;
 879               if(part->GetFirstMother() >= 0) continue;
 880               if (count[j]>1) //at least two digits at this padrow
 881                 {
 882                   crossed++;
 883                   count[j]=0;
 884                 }
 885             }
 886
 887           Int_t patch = AliL3Transform::GetPatch(sr);
 888           if(sp==kTRUE)
 889             patch=0;
 890           AliL3SpacePointData *points = fClusters[sl][patch];
 891           if(!points)
 892             continue;
 893           for(UInt_t k=0; k<fNcl[sl][patch]; k++)
 894             {
 895               if(points[k].fPadRow!=sr) continue;
 896               recs++;
 897             }
 898           ntuppel->Fill(sl,sr,crossed,recs);
 899         }
 900
 901       TD->Delete();
 902       delete[] count;
 903     }
 904   TFile *file = TFile::Open(outfile,"RECREATE");
 905   ntuppel->Write();
 906   file->Close();
 907
 908   rfile->Close();
 909   dfile->Close();
 910 #endif
 911 }
 912
 913 Float_t AliL3Evaluate::GetCrossingAngle(TParticle *part,Int_t slice,Int_t padrow,Float_t *xyz)
 914 {
 915   //Calculate the padrow crossing angle of the particle
 916
 917   Double_t kappa = AliL3Transform::GetBField()*AliL3Transform::GetBFact()/part->Pt();
 918
 919   Double_t radius = 1/fabs(kappa);
 920   if(part->GetPdgCode() > 0) kappa = -kappa;
 921
 922   Float_t angl[1] = {part->Phi()};
 923
 924   AliL3Transform::Global2LocalAngle(angl,slice);
 925
 926   Double_t charge = -1.*kappa;
 927
 928   Double_t trackPhi0 = angl[0] + charge*0.5*AliL3Transform::Pi()/fabs(charge);
 929
 930   Double_t x0=0;
 931   Double_t y0=0;
 932   Double_t xc = x0 - radius * cos(trackPhi0);
 933   Double_t yc = y0 - radius * sin(trackPhi0);
 934
 935   Double_t tangent[2];
 936   tangent[0] = -1.*(xyz[1] - yc)/radius;
 937   tangent[1] = (xyz[0] - xc)/radius;
 938
 939   Double_t perp_padrow[2] = {1,0}; //locally in slice
 940
 941   Double_t cos_beta = fabs(tangent[0]*perp_padrow[0] + tangent[1]*perp_padrow[1]);
 942   if(cos_beta > 1) cos_beta=1;
 943   return acos(cos_beta);
 944 }
 945
 946 Int_t AliL3Evaluate::FindPrimaries(Int_t nparticles)
 947 {
 948   // cuts:
 949   Double_t vertcut = 0.001;
 950   Double_t decacut = 3.;
 951   Double_t timecut = 0.;
 952   Int_t nprch1=0;
 953   TParticle * part = gAlice->Particle(0);
 954   Double_t xori = part->Vx();
 955   Double_t yori = part->Vy();
 956   Double_t zori = part->Vz();
 957   for(Int_t iprim = 0; iprim<nparticles; iprim++){   //loop on  tracks
 958
 959     part = gAlice->Particle(iprim);
 960     char *xxx=strstr(part->GetName(),"XXX");
 961     if(xxx)continue;
 962
 963     TParticlePDG *ppdg = part->GetPDG();
 964     if(TMath::Abs(ppdg->Charge())!=3)continue;  // only charged (no quarks)
 965
 966     Double_t dist=TMath::Sqrt((part->Vx()-xori)*(part->Vx()-xori)+(part->Vy()-yori)*(part->Vy()-yori)+(part->Vz()-zori)*(part->Vz()-zori));
 967     if(dist>vertcut)continue;  // cut on the vertex
 968
 969     if(part->T()>timecut)continue;
 970
 971     Double_t ptot=TMath::Sqrt(part->Px()*part->Px()+part->Py()*part->Py()+part->Pz()*part->Pz());
 972     if(ptot==(TMath::Abs(part->Pz())))continue; // no beam particles
 973
 974     Bool_t prmch = kTRUE;   // candidate primary track
 975     Int_t fidau=part->GetFirstDaughter();  // cut on daughters
 976     Int_t lasdau=0;
 977     Int_t ndau=0;
 978     if(fidau>=0){
 979       lasdau=part->GetLastDaughter();
 980       ndau=lasdau-fidau+1;
 981     }
 982     if(ndau>0){
 983       for(Int_t j=fidau;j<=lasdau;j++){
 984         TParticle *dau=gAlice->Particle(j);
 985         Double_t distd=TMath::Sqrt((dau->Vx()-xori)*(dau->Vx()-xori)+(dau->Vy()-yori)*(dau->Vy()-yori)+(dau->Vz()-zori)*(dau->Vz()-zori));
 986         if(distd<decacut)prmch=kFALSE;  // eliminate if the decay is near the vertex
 987       }
 988     }
 989
 990     if(prmch){
 991       nprch1++;
 992       part->SetBit(kPrimaryCharged);
 993     }
 994   }
 995
 996   return nprch1;
 997 }