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