[u/mrichter/AliRoot.git] / HLT / hough / AliL3HoughEval.cxx

//Author:        Anders Strand Vestbo
//Last Modified: 28.6.01

#include <TClonesArray.h>
#include <TH2.h>
#include <TH1.h>
#include <TFile.h>
#include <AliRun.h>
#include <TParticle.h>
#include <TTree.h>

#include "AliTPCParam.h"
#include "AliSimDigits.h"
#include "AliL3TrackArray.h"
#include "AliL3Transform.h"
#include "AliL3HoughTransformer.h"
#include "AliL3Defs.h"
#include "AliL3HoughTrack.h"
#include "AliL3HoughEval.h"

ClassImp(AliL3HoughEval)


AliL3HoughEval::AliL3HoughEval()
{
  //Default constructor
  fTransform = new AliL3Transform();
  fHoughTransformer = 0;
  fNumOfRowsToMiss = 1;
  fNumOfPadsToLook = 1;
}


AliL3HoughEval::AliL3HoughEval(AliL3HoughTransformer *transformer)
{
  //Constructor
  fHoughTransformer = transformer;
  fTransform = new AliL3Transform();
  fNumOfRowsToMiss = 1;
  fNumOfPadsToLook = 1;
}


AliL3HoughEval::~AliL3HoughEval()
{
  //Destructor
  if(fTransform)
    delete fTransform;
  if(fMcTrackTable)
    delete [] fMcTrackTable;
}


Bool_t AliL3HoughEval::LookInsideRoad(AliL3HoughTrack *track,Int_t eta_index,Bool_t remove)
{
  //Look at rawdata along the road specified by the track candidates.
  
  if(!fHoughTransformer)
    {
      printf("\nAliL3HoughEval: No transformer object\n");
      return kFALSE;
    }
  
  Int_t patch = fHoughTransformer->fPatch;
  Int_t slice = fHoughTransformer->fSlice;

  
  Int_t sector,row;
  Int_t lut_index;
  Char_t **track_lut = fHoughTransformer->fTrackTable;
  Int_t nrow=0,npixs=0;
  Float_t xyz[3];
    
  for(Int_t padrow = NRows[patch][0]; padrow <= NRows[patch][1]; padrow++)
    {
      Int_t prow = padrow - NRows[patch][0];
      if(!track->GetCrossingPoint(padrow,xyz))  
	{
	  printf("AliL3HoughEval::LookInsideRoad : Track does not cross line!!\n");
	  continue;
	}
      
      fTransform->Slice2Sector(slice,padrow,sector,row);
      fTransform->Local2Raw(xyz,sector,row);
      npixs=0;
      //Look at both sides of the crossing point:
      for(Int_t p=(Int_t)xyz[1]-fNumOfPadsToLook; p<=(Int_t)xyz[1]+fNumOfPadsToLook; p++)
	{
	  if(p<0 || p>fTransform->GetNPads(padrow)) continue;
	  lut_index = (prow<<8) + p;
	  if(track_lut[eta_index][lut_index]>0) //There was a signal here
	    npixs++;
	}
      
      if(npixs > 0)
	{
	  nrow++;
	}	  
    }
  if(nrow >= NRows[patch][1]-NRows[patch][0]-fNumOfRowsToMiss)//this was a good track
    {
      track->SetEtaIndex(eta_index);
      if(remove)
	RemoveTrackFromImage(track,eta_index);
      return kTRUE;
    }
  else
    return kFALSE;
}

void AliL3HoughEval::LookInsideRawRoad(AliL3TrackArray *tracks,Int_t eta_index,Bool_t remove)
{
  //Evalaute the track candidates by looking along the trajectory.
  //If remove is on, the pixels along the track will be removed. 


  if(!fHoughTransformer)
    {
      printf("AliL3HoughEval: No transformer object\n");
      return;
    }
  
  Int_t patch = fHoughTransformer->fPatch;
  Int_t slice = fHoughTransformer->fSlice;

  
  Int_t sector,row;
  Int_t lut_index;
  Char_t **track_lut = fHoughTransformer->fTrackTable;
  Int_t nrow,npixs;
  Float_t xyz[3];
  
  for(Int_t i=0; i<tracks->GetNTracks(); i++)
    {
      AliL3HoughTrack *track = (AliL3HoughTrack*)tracks->GetCheckedTrack(i);
      if(!track) {printf("No track\n"); return;}
      
      nrow = 0;
      
      for(Int_t padrow = NRows[patch][0]; padrow <= NRows[patch][1]; padrow++)
	{
	  Int_t prow = padrow - NRows[patch][0];
	  if(!track->GetCrossingPoint(padrow,xyz))  
	    {
	      printf("AliL3HoughEval::LookInsideRoad : Track does not cross line!!\n");
	      continue;
	    }
	  	  
	  fTransform->Slice2Sector(slice,padrow,sector,row);
	  fTransform->Local2Raw(xyz,sector,row);
	  npixs=0;
	  //Look at both sides of the crossing point:
	  for(Int_t p=(Int_t)xyz[1]-fNumOfPadsToLook; p<=(Int_t)xyz[1]+fNumOfPadsToLook; p++)
	    {
	      if(p<0 || p>fTransform->GetNPads(padrow)) continue;
	      lut_index = (prow<<8) + p;
	      if(track_lut[eta_index][lut_index]>0) //There was a signal here
		npixs++;
	    }
	  
	  if(npixs > 0)
	    {
	      nrow++;
	    }	  
	}
      if(nrow < NRows[patch][1]-NRows[patch][0]-fNumOfRowsToMiss)
	tracks->Remove(i); //this was not a good enough track
      else if(remove)
	RemoveTrackFromImage(track,eta_index); //this was a good track, so remove it from the image.
    }
  
  tracks->Compress();
}

void AliL3HoughEval::RemoveTrackFromImage(AliL3HoughTrack *track,Int_t eta_index)
{
  //Remove the pixels along the track in the image. 
  
  Int_t patch = fHoughTransformer->fPatch;
  Int_t slice = fHoughTransformer->fSlice;

  Int_t maxnpads = fTransform->GetNPads(NRows[patch][1]);
  
  Int_t lut_index;
  Char_t **track_lut = fHoughTransformer->fTrackTable;
  Int_t sector,row;
  Float_t xyz[3];
  for(Int_t padrow = NRows[patch][0]; padrow<=NRows[patch][1]; padrow++)
    {
      Int_t prow = padrow - NRows[patch][0];
      if(!track->GetCrossingPoint(padrow,xyz))  
	{
	  printf("AliL3HoughEval::LookInsideRoad : Track does not cross line!!\n");
	  continue;
	}
      
      fTransform->Slice2Sector(slice,padrow,sector,row);
      fTransform->Local2Raw(xyz,sector,row);
      
      for(Int_t p=(Int_t)xyz[1]-fNumOfPadsToLook; p<=(Int_t)xyz[1]+fNumOfPadsToLook; p++)
	{
	  if(p<0 || p>fTransform->GetNPads(padrow)) continue;
	  lut_index = (prow<<8) + p;
	  track_lut[eta_index][lut_index] = -1; //remove it!!
	}
    }
  
  
}

void AliL3HoughEval::DisplaySlice(TH2F *hist)
{

  AliL3Digits *pixel;
  
  for(Int_t padrow = 0; padrow < 174; padrow++)
    {
      
      for(pixel=(AliL3Digits*)fHoughTransformer->fRowContainer[padrow].first; pixel!=0; pixel=(AliL3Digits*)pixel->nextRowPixel)
	{
	  
	  Int_t sector,row;
	  Float_t xyz_pix[3];
	  fTransform->Slice2Sector(fHoughTransformer->fSlice,padrow,sector,row);
	  fTransform->Raw2Local(xyz_pix,sector,row,pixel->fPad,pixel->fTime); //y alone pad
	  
	  if(hist)
	    {
	      hist->Fill(xyz_pix[0],xyz_pix[1],pixel->fCharge);
	    }    
	  
	}
    }
  
}

void AliL3HoughEval::DefineGoodParticles(Char_t *rootfile,Double_t pet)
{
  //define the particles that produce good enough signals to be recognized in the transform

  Int_t num_eta_segments = fHoughTransformer->fNumEtaSegments;
  Double_t eta_max = fHoughTransformer->fEtaMax;
  fMcTrackTable = new Int_t[num_eta_segments];
  for(Int_t i=0; i<num_eta_segments; i++)
    fMcTrackTable[i]=0;
  Double_t etaslice = eta_max/num_eta_segments;

  Int_t patch = fHoughTransformer->fPatch;
  Int_t slice = fHoughTransformer->fSlice;
  
  TFile *file = new TFile(rootfile);
  file->cd();

  AliRun *gAlice = (AliRun*)file->Get("gAlice");
  gAlice->GetEvent(0);
  
  TClonesArray *particles=gAlice->Particles();
  Int_t np = particles->GetEntriesFast();
  
  Int_t row_to_miss = fNumOfRowsToMiss;
  AliTPCParam *param = (AliTPCParam*)file->Get("75x40_100x60");
  Int_t zero = param->GetZeroSup();
  TTree *TD=(TTree*)gDirectory->Get("TreeD_75x40_100x60");
  AliSimDigits da, *digits=&da;
  TD->GetBranch("Segment")->SetAddress(&digits); //Return pointer to branch segment.
  Int_t *good=new Int_t[np];
  Int_t *count = new Int_t[np]; //np number of particles.
  Int_t good_number = NRows[patch][1]-NRows[patch][0]-row_to_miss;
  Int_t i;
  for (i=0; i<np; i++) count[i]=0;
  Int_t sectors_by_rows=(Int_t)TD->GetEntries();
  for (i=0; i<sectors_by_rows; i++) 
    {
      if (!TD->GetEvent(i)) continue;
      Int_t sec,row;
      param->AdjustSectorRow(digits->GetID(),sec,row);
      Int_t ss,sr;
      fTransform->Sector2Slice(ss,sr,sec,row);
      if(ss!=slice) continue;
      if(sr < NRows[patch][0]) continue;
      if(sr > NRows[patch][1]) break;
      digits->First();
      while (digits->Next()) {
	Int_t it=digits->CurrentRow(), ip=digits->CurrentColumn();
	Short_t dig = digits->GetDigit(it,ip);
	Int_t idx0=digits->GetTrackID(it,ip,0); 
	Int_t idx1=digits->GetTrackID(it,ip,1);
	Int_t idx2=digits->GetTrackID(it,ip,2);
	if (idx0>=0 && dig>=zero) count[idx0]+=1;
	if (idx1>=0 && dig>=zero) count[idx1]+=1;
	if (idx2>=0 && dig>=zero) count[idx2]+=1;
      }
      for (Int_t j=0; j<np; j++) 
	{
	  if (count[j]>1) {//at least two digits at this padrow
	    good[j]++;
	  }
	  count[j]=0;
	}
    }
  delete[] count;
  
  Int_t good_one=0;
  //TObjArray *part = new TObjArray(0,0);
  for(i=0; i<np; i++)
   {
     TParticle *p = (TParticle*)particles->UncheckedAt(i);
     if(p->GetFirstMother()>0) continue; //secondary particle
     if(good[i] < good_number) continue; //too few padrows
     Double_t ptg=p->Pt();
     if(ptg<pet) continue;
     Int_t eta_index = (Int_t)(p->Eta()/etaslice);
     if(eta_index < 0 || eta_index > num_eta_segments)
       continue;
     //fMcTrackTable[eta_index]++;
     //part->AddLast(p);
     good_one++;
   }
  
  printf("nparticles %d\n",good_one);
  file->Close();
  delete [] good;
  delete file;
  //return part;
}


void AliL3HoughEval::CompareMC(Char_t *rootfile,AliL3TrackArray *merged_tracks,Float_t *eta)
{
  
  Int_t slice = fHoughTransformer->fSlice;
  
  TFile *file = new TFile(rootfile);
  file->cd();
  
  AliRun *gAlice = (AliRun*)file->Get("gAlice");
  gAlice->GetEvent(0);
  
  TClonesArray *particles=gAlice->Particles();  
  Int_t n=particles->GetEntriesFast();
  Float_t torad=TMath::Pi()/180;
  Float_t phi_min = slice*20 - 10;
  Float_t phi_max = slice*20 + 10;
  
  for (Int_t j=0; j<n; j++) {
    TParticle *p=(TParticle*)particles->UncheckedAt(j);
    if (p->GetFirstMother()>=0) continue;  //secondary particle
    if(p->Eta() < eta[0] || p->Eta() > eta[1]) continue;
    Double_t ptg=p->Pt(),pxg=p->Px(),pyg=p->Py();//,pzg=p->Pz();
    Double_t phi_part = TMath::ATan2(pyg,pxg);
    if (phi_part < 0) phi_part += 2*TMath::Pi();
    
    if(phi_part < phi_min*torad || phi_part > phi_max*torad) {continue;}
    if(ptg<0.100) continue;
    
    AliL3HoughTrack *sel_track=0;
    
    Double_t min_dist = 10000;
    
    for(Int_t t=0; t<merged_tracks->GetNTracks(); t++)
      {
	AliL3HoughTrack *track = (AliL3HoughTrack*)merged_tracks->GetCheckedTrack(t);
	if(!track) {printf("AliL3HoughEval: NO TRACK\n"); continue;}
	Float_t phi0[2] = {track->GetPhi0(),0};
	fTransform->Local2GlobalAngle(phi0,slice);
	Double_t dPt = ptg - track->GetPt();
	Double_t dPhi0 = phi_part - phi0[0];
	Double_t distance = pow(dPt,2) + pow(dPhi0,2);
	if(distance < min_dist)
	  {
	    min_dist = distance;
	    sel_track = track;
	  }      
	
      }
    if(sel_track)
      sel_track->SetBestMCid(j,min_dist);
    
    Double_t dpt = fabs(ptg-sel_track->GetPt());
    Double_t dphi0 = fabs(phi_part-sel_track->GetPhi0());
    //printf("Found match, min_dist %f dPt %f dPhi0 %f\n",min_dist,dpt,dphi0);
  }
  file->Close();
  delete file;
  
}
Commit	Line	Data
f80b98cb	1	//Author: Anders Strand Vestbo
	2	//Last Modified: 28.6.01
	3
4de874d1	4	#include <TClonesArray.h>
	5	#include <TH2.h>
	6	#include <TH1.h>
	7	#include <TFile.h>
	8	#include <AliRun.h>
	9	#include <TParticle.h>
f80b98cb	10	#include <TTree.h>
4de874d1	11
f80b98cb	12	#include "AliTPCParam.h"
f80b98cb	13	#include "AliSimDigits.h"
4de874d1	14	#include "AliL3TrackArray.h"
	15	#include "AliL3Transform.h"
	16	#include "AliL3HoughTransformer.h"
	17	#include "AliL3Defs.h"
	18	#include "AliL3HoughTrack.h"
	19	#include "AliL3HoughEval.h"
	20
	21	ClassImp(AliL3HoughEval)
	22
	23
	24	AliL3HoughEval::AliL3HoughEval()
	25	{
	26	//Default constructor
	27	fTransform = new AliL3Transform();
f80b98cb	28	fHoughTransformer = 0;
	29	fNumOfRowsToMiss = 1;
	30	fNumOfPadsToLook = 1;
4de874d1	31	}
	32
	33
	34	AliL3HoughEval::AliL3HoughEval(AliL3HoughTransformer *transformer)
	35	{
	36	//Constructor
	37	fHoughTransformer = transformer;
	38	fTransform = new AliL3Transform();
f80b98cb	39	fNumOfRowsToMiss = 1;
f80b98cb	40	fNumOfPadsToLook = 1;
4de874d1	41	}
	42
	43
	44	AliL3HoughEval::~AliL3HoughEval()
	45	{
	46	//Destructor
	47	if(fTransform)
	48	delete fTransform;
f80b98cb	49	if(fMcTrackTable)
f80b98cb	50	delete [] fMcTrackTable;
4de874d1	51	}
	52
	53
f80b98cb	54	Bool_t AliL3HoughEval::LookInsideRoad(AliL3HoughTrack *track,Int_t eta_index,Bool_t remove)
4de874d1	55	{
4de874d1	56	//Look at rawdata along the road specified by the track candidates.
f80b98cb	57
	58	if(!fHoughTransformer)
	59	{
	60	printf("\nAliL3HoughEval: No transformer object\n");
	61	return kFALSE;
	62	}
	63
	64	Int_t patch = fHoughTransformer->fPatch;
	65	Int_t slice = fHoughTransformer->fSlice;
	66
	67
	68	Int_t sector,row;
	69	Int_t lut_index;
	70	Char_t **track_lut = fHoughTransformer->fTrackTable;
	71	Int_t nrow=0,npixs=0;
	72	Float_t xyz[3];
	73
	74	for(Int_t padrow = NRows[patch][0]; padrow <= NRows[patch][1]; padrow++)
	75	{
	76	Int_t prow = padrow - NRows[patch][0];
	77	if(!track->GetCrossingPoint(padrow,xyz))
	78	{
	79	printf("AliL3HoughEval::LookInsideRoad : Track does not cross line!!\n");
	80	continue;
	81	}
	82
	83	fTransform->Slice2Sector(slice,padrow,sector,row);
	84	fTransform->Local2Raw(xyz,sector,row);
	85	npixs=0;
	86	//Look at both sides of the crossing point:
	87	for(Int_t p=(Int_t)xyz[1]-fNumOfPadsToLook; p<=(Int_t)xyz[1]+fNumOfPadsToLook; p++)
	88	{
	89	if(p<0 \|\| p>fTransform->GetNPads(padrow)) continue;
	90	lut_index = (prow<<8) + p;
	91	if(track_lut[eta_index][lut_index]>0) //There was a signal here
	92	npixs++;
	93	}
	94
	95	if(npixs > 0)
	96	{
	97	nrow++;
	98	}
	99	}
	100	if(nrow >= NRows[patch][1]-NRows[patch][0]-fNumOfRowsToMiss)//this was a good track
	101	{
	102	track->SetEtaIndex(eta_index);
	103	if(remove)
	104	RemoveTrackFromImage(track,eta_index);
	105	return kTRUE;
	106	}
	107	else
	108	return kFALSE;
	109	}
	110
	111	void AliL3HoughEval::LookInsideRawRoad(AliL3TrackArray *tracks,Int_t eta_index,Bool_t remove)
	112	{
	113	//Evalaute the track candidates by looking along the trajectory.
	114	//If remove is on, the pixels along the track will be removed.
	115
4de874d1	116
	117	if(!fHoughTransformer)
	118	{
	119	printf("AliL3HoughEval: No transformer object\n");
	120	return;
	121	}
	122
4de874d1	123	Int_t patch = fHoughTransformer->fPatch;
4de874d1	124	Int_t slice = fHoughTransformer->fSlice;
4de874d1	125
f80b98cb	126
	127	Int_t sector,row;
	128	Int_t lut_index;
	129	Char_t **track_lut = fHoughTransformer->fTrackTable;
	130	Int_t nrow,npixs;
	131	Float_t xyz[3];
	132
4de874d1	133	for(Int_t i=0; i<tracks->GetNTracks(); i++)
	134	{
	135	AliL3HoughTrack track = (AliL3HoughTrack)tracks->GetCheckedTrack(i);
	136	if(!track) {printf("No track\n"); return;}
	137
f80b98cb	138	nrow = 0;
4de874d1	139
	140	for(Int_t padrow = NRows[patch][0]; padrow <= NRows[patch][1]; padrow++)
	141	{
f80b98cb	142	Int_t prow = padrow - NRows[patch][0];
4de874d1	143	if(!track->GetCrossingPoint(padrow,xyz))
	144	{
	145	printf("AliL3HoughEval::LookInsideRoad : Track does not cross line!!\n");
	146	continue;
	147	}
f80b98cb	148
	149	fTransform->Slice2Sector(slice,padrow,sector,row);
	150	fTransform->Local2Raw(xyz,sector,row);
	151	npixs=0;
	152	//Look at both sides of the crossing point:
	153	for(Int_t p=(Int_t)xyz[1]-fNumOfPadsToLook; p<=(Int_t)xyz[1]+fNumOfPadsToLook; p++)
4de874d1	154	{
f80b98cb	155	if(p<0 \|\| p>fTransform->GetNPads(padrow)) continue;
	156	lut_index = (prow<<8) + p;
	157	if(track_lut[eta_index][lut_index]>0) //There was a signal here
	158	npixs++;
4de874d1	159	}
f80b98cb	160
4de874d1	161	if(npixs > 0)
	162	{
	163	nrow++;
	164	}
	165	}
f80b98cb	166	if(nrow < NRows[patch][1]-NRows[patch][0]-fNumOfRowsToMiss)
	167	tracks->Remove(i); //this was not a good enough track
	168	else if(remove)
	169	RemoveTrackFromImage(track,eta_index); //this was a good track, so remove it from the image.
4de874d1	170	}
	171
	172	tracks->Compress();
f80b98cb	173	}
	174
	175	void AliL3HoughEval::RemoveTrackFromImage(AliL3HoughTrack *track,Int_t eta_index)
	176	{
	177	//Remove the pixels along the track in the image.
	178
	179	Int_t patch = fHoughTransformer->fPatch;
	180	Int_t slice = fHoughTransformer->fSlice;
4de874d1	181
f80b98cb	182	Int_t maxnpads = fTransform->GetNPads(NRows[patch][1]);
	183
	184	Int_t lut_index;
	185	Char_t **track_lut = fHoughTransformer->fTrackTable;
	186	Int_t sector,row;
	187	Float_t xyz[3];
	188	for(Int_t padrow = NRows[patch][0]; padrow<=NRows[patch][1]; padrow++)
	189	{
	190	Int_t prow = padrow - NRows[patch][0];
	191	if(!track->GetCrossingPoint(padrow,xyz))
	192	{
	193	printf("AliL3HoughEval::LookInsideRoad : Track does not cross line!!\n");
	194	continue;
	195	}
	196
	197	fTransform->Slice2Sector(slice,padrow,sector,row);
	198	fTransform->Local2Raw(xyz,sector,row);
	199
	200	for(Int_t p=(Int_t)xyz[1]-fNumOfPadsToLook; p<=(Int_t)xyz[1]+fNumOfPadsToLook; p++)
	201	{
	202	if(p<0 \|\| p>fTransform->GetNPads(padrow)) continue;
	203	lut_index = (prow<<8) + p;
	204	track_lut[eta_index][lut_index] = -1; //remove it!!
	205	}
	206	}
	207
4de874d1	208
	209	}
	210
	211	void AliL3HoughEval::DisplaySlice(TH2F *hist)
	212	{
	213
	214	AliL3Digits *pixel;
	215
	216	for(Int_t padrow = 0; padrow < 174; padrow++)
	217	{
	218
	219	for(pixel=(AliL3Digits)fHoughTransformer->fRowContainer[padrow].first; pixel!=0; pixel=(AliL3Digits)pixel->nextRowPixel)
	220	{
	221
	222	Int_t sector,row;
	223	Float_t xyz_pix[3];
	224	fTransform->Slice2Sector(fHoughTransformer->fSlice,padrow,sector,row);
	225	fTransform->Raw2Local(xyz_pix,sector,row,pixel->fPad,pixel->fTime); //y alone pad
	226
	227	if(hist)
	228	{
	229	hist->Fill(xyz_pix[0],xyz_pix[1],pixel->fCharge);
	230	}
	231
	232	}
	233	}
	234
	235	}
	236
f80b98cb	237	void AliL3HoughEval::DefineGoodParticles(Char_t *rootfile,Double_t pet)
4de874d1	238	{
f80b98cb	239	//define the particles that produce good enough signals to be recognized in the transform
	240
	241	Int_t num_eta_segments = fHoughTransformer->fNumEtaSegments;
	242	Double_t eta_max = fHoughTransformer->fEtaMax;
	243	fMcTrackTable = new Int_t[num_eta_segments];
	244	for(Int_t i=0; i<num_eta_segments; i++)
	245	fMcTrackTable[i]=0;
	246	Double_t etaslice = eta_max/num_eta_segments;
	247
	248	Int_t patch = fHoughTransformer->fPatch;
	249	Int_t slice = fHoughTransformer->fSlice;
4de874d1	250
4de874d1	251	TFile *file = new TFile(rootfile);
	252	file->cd();
	253
	254	AliRun gAlice = (AliRun)file->Get("gAlice");
	255	gAlice->GetEvent(0);
	256
	257	TClonesArray *particles=gAlice->Particles();
f80b98cb	258	Int_t np = particles->GetEntriesFast();
	259
	260	Int_t row_to_miss = fNumOfRowsToMiss;
	261	AliTPCParam param = (AliTPCParam)file->Get("75x40_100x60");
	262	Int_t zero = param->GetZeroSup();
	263	TTree TD=(TTree)gDirectory->Get("TreeD_75x40_100x60");
	264	AliSimDigits da, *digits=&da;
	265	TD->GetBranch("Segment")->SetAddress(&digits); //Return pointer to branch segment.
	266	Int_t *good=new Int_t[np];
	267	Int_t *count = new Int_t[np]; //np number of particles.
	268	Int_t good_number = NRows[patch][1]-NRows[patch][0]-row_to_miss;
	269	Int_t i;
	270	for (i=0; i<np; i++) count[i]=0;
	271	Int_t sectors_by_rows=(Int_t)TD->GetEntries();
	272	for (i=0; i<sectors_by_rows; i++)
	273	{
	274	if (!TD->GetEvent(i)) continue;
	275	Int_t sec,row;
	276	param->AdjustSectorRow(digits->GetID(),sec,row);
	277	Int_t ss,sr;
	278	fTransform->Sector2Slice(ss,sr,sec,row);
	279	if(ss!=slice) continue;
	280	if(sr < NRows[patch][0]) continue;
	281	if(sr > NRows[patch][1]) break;
	282	digits->First();
	283	while (digits->Next()) {
	284	Int_t it=digits->CurrentRow(), ip=digits->CurrentColumn();
	285	Short_t dig = digits->GetDigit(it,ip);
	286	Int_t idx0=digits->GetTrackID(it,ip,0);
	287	Int_t idx1=digits->GetTrackID(it,ip,1);
	288	Int_t idx2=digits->GetTrackID(it,ip,2);
	289	if (idx0>=0 && dig>=zero) count[idx0]+=1;
	290	if (idx1>=0 && dig>=zero) count[idx1]+=1;
	291	if (idx2>=0 && dig>=zero) count[idx2]+=1;
	292	}
	293	for (Int_t j=0; j<np; j++)
	294	{
	295	if (count[j]>1) {//at least two digits at this padrow
	296	good[j]++;
	297	}
	298	count[j]=0;
	299	}
	300	}
	301	delete[] count;
	302
	303	Int_t good_one=0;
	304	//TObjArray *part = new TObjArray(0,0);
	305	for(i=0; i<np; i++)
	306	{
	307	TParticle p = (TParticle)particles->UncheckedAt(i);
	308	if(p->GetFirstMother()>0) continue; //secondary particle
	309	if(good[i] < good_number) continue; //too few padrows
	310	Double_t ptg=p->Pt();
	311	if(ptg<pet) continue;
	312	Int_t eta_index = (Int_t)(p->Eta()/etaslice);
	313	if(eta_index < 0 \|\| eta_index > num_eta_segments)
	314	continue;
	315	//fMcTrackTable[eta_index]++;
	316	//part->AddLast(p);
	317	good_one++;
	318	}
	319
	320	printf("nparticles %d\n",good_one);
9ee44abe	321	file->Close();
f80b98cb	322	delete [] good;
9ee44abe	323	delete file;
f80b98cb	324	//return part;
9ee44abe	325	}
	326
	327
	328	void AliL3HoughEval::CompareMC(Char_t rootfile,AliL3TrackArray merged_tracks,Float_t *eta)
	329	{
	330
f80b98cb	331	Int_t slice = fHoughTransformer->fSlice;
	332
	333	TFile *file = new TFile(rootfile);
	334	file->cd();
	335
	336	AliRun gAlice = (AliRun)file->Get("gAlice");
	337	gAlice->GetEvent(0);
9ee44abe	338
f80b98cb	339	TClonesArray *particles=gAlice->Particles();
4de874d1	340	Int_t n=particles->GetEntriesFast();
	341	Float_t torad=TMath::Pi()/180;
	342	Float_t phi_min = slice*20 - 10;
	343	Float_t phi_max = slice*20 + 10;
	344
	345	for (Int_t j=0; j<n; j++) {
	346	TParticle p=(TParticle)particles->UncheckedAt(j);
	347	if (p->GetFirstMother()>=0) continue; //secondary particle
	348	if(p->Eta() < eta[0] \|\| p->Eta() > eta[1]) continue;
	349	Double_t ptg=p->Pt(),pxg=p->Px(),pyg=p->Py();//,pzg=p->Pz();
	350	Double_t phi_part = TMath::ATan2(pyg,pxg);
	351	if (phi_part < 0) phi_part += 2*TMath::Pi();
	352
	353	if(phi_part < phi_mintorad \|\| phi_part > phi_maxtorad) {continue;}
	354	if(ptg<0.100) continue;
	355
	356	AliL3HoughTrack *sel_track=0;
	357
	358	Double_t min_dist = 10000;
	359
	360	for(Int_t t=0; t<merged_tracks->GetNTracks(); t++)
	361	{
	362	AliL3HoughTrack track = (AliL3HoughTrack)merged_tracks->GetCheckedTrack(t);
	363	if(!track) {printf("AliL3HoughEval: NO TRACK\n"); continue;}
	364	Float_t phi0[2] = {track->GetPhi0(),0};
	365	fTransform->Local2GlobalAngle(phi0,slice);
	366	Double_t dPt = ptg - track->GetPt();
	367	Double_t dPhi0 = phi_part - phi0[0];
	368	Double_t distance = pow(dPt,2) + pow(dPhi0,2);
	369	if(distance < min_dist)
	370	{
	371	min_dist = distance;
	372	sel_track = track;
	373	}
	374
	375	}
	376	if(sel_track)
	377	sel_track->SetBestMCid(j,min_dist);
	378
	379	Double_t dpt = fabs(ptg-sel_track->GetPt());
	380	Double_t dphi0 = fabs(phi_part-sel_track->GetPhi0());
	381	//printf("Found match, min_dist %f dPt %f dPhi0 %f\n",min_dist,dpt,dphi0);
	382	}
f80b98cb	383	file->Close();
	384	delete file;
	385
4de874d1	386	}