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