HLT/src/AliL3Evaluate.cxx

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