HLT/src/AliL3Evaluate.cxx

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