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

#include <TClonesArray.h>
#include <TH2.h>
#include <TH1.h>
#include <TFile.h>
#include <AliRun.h>
#include <TParticle.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 = NULL;
}


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


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


void AliL3HoughEval::LookInsideRoad(AliL3TrackArray *tracks,Bool_t remove,TH2F *hist)
{
  //Look at rawdata along the road specified by the track candidates.

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

  for(Int_t i=0; i<tracks->GetNTracks(); i++)
    {
      AliL3HoughTrack *track = (AliL3HoughTrack*)tracks->GetCheckedTrack(i);
      if(!track) {printf("No track\n"); return;}
      
      Int_t nrow = 0;
      
      for(Int_t padrow = NRows[patch][0]; padrow <= NRows[patch][1]; padrow++)
	{
	  
	  Float_t xyz[3];
	  //Get the crossing point of the track and current padrow:
	  //if(!track->GetCrossingPoint(slice,padrow,xyz))
	  if(!track->GetCrossingPoint(padrow,xyz))  
	    {
	      printf("AliL3HoughEval::LookInsideRoad : Track does not cross line!!\n");
	      continue;
	    }
	  
	  Int_t npixs = 0;
	  
	  //Get the pixels along the track candidate
	  for(pixel=(AliL3Digits*)fHoughTransformer->fRowContainer[padrow].first; pixel!=0; pixel=(AliL3Digits*)pixel->nextRowPixel)
	    {
	      
	      Int_t sector,row;
	      Float_t xyz_pix[3];
	      fTransform->Slice2Sector(slice,padrow,sector,row);
	      fTransform->Raw2Local(xyz_pix,sector,row,pixel->fPad,pixel->fTime); //y alone pad
	      
	      
	      //check if we are inside road
	      if(fabs(xyz_pix[1] - xyz[1]) > num_of_pads_to_look*fTransform->GetPadPitchWidthLow()) continue; 
	      npixs++;
	      
	      if(remove) //Remove this pixel from image
		pixel->fIndex = -1;

	      if(hist)
		{
		  //fTransform->Local2Global(xyz_pix,slice);
		  hist->Fill(xyz_pix[0],xyz_pix[1],pixel->fCharge);
		}    
	      
	    }
	  if(npixs > 0)
	    {
	      nrow++;
	    }	  
	}
      
      if(nrow < NRows[patch][1]-NRows[patch][0]-rows_to_miss)
	tracks->Remove(i);
	      
    }
  
  tracks->Compress();

  
}

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);
	    }    
	  
	}
    }
  
}

TClonesArray *AliL3HoughEval::GetParticles(Char_t *rootfile)
{
  
  TFile *file = new TFile(rootfile);
  file->cd();

  AliRun *gAlice = (AliRun*)file->Get("gAlice");
  gAlice->GetEvent(0);
  
  TClonesArray *particles=gAlice->Particles();
  return particles;

  file->Close();
  delete file;
}


void AliL3HoughEval::CompareMC(Char_t *rootfile,AliL3TrackArray *merged_tracks,Float_t *eta)
{
  
  Int_t slice = fSlice;
  
  TClonesArray *particles = GetParticles(rootfile);
  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);
  }
    
}
Commit	Line	Data
4de874d1	1	#include <TClonesArray.h>
	2	#include <TH2.h>
	3	#include <TH1.h>
	4	#include <TFile.h>
	5	#include <AliRun.h>
	6	#include <TParticle.h>
	7
	8	#include "AliL3TrackArray.h"
	9	#include "AliL3Transform.h"
	10	#include "AliL3HoughTransformer.h"
	11	#include "AliL3Defs.h"
	12	#include "AliL3HoughTrack.h"
	13	#include "AliL3HoughEval.h"
	14
	15	ClassImp(AliL3HoughEval)
	16
	17
	18	AliL3HoughEval::AliL3HoughEval()
	19	{
	20	//Default constructor
	21	fTransform = new AliL3Transform();
	22	fHoughTransformer = NULL;
	23	}
	24
	25
	26	AliL3HoughEval::AliL3HoughEval(AliL3HoughTransformer *transformer)
	27	{
	28	//Constructor
	29	fHoughTransformer = transformer;
	30	fTransform = new AliL3Transform();
	31	}
	32
	33
	34	AliL3HoughEval::~AliL3HoughEval()
	35	{
	36	//Destructor
	37	if(fTransform)
	38	delete fTransform;
	39
	40	}
	41
	42
52a2a604	43	void AliL3HoughEval::LookInsideRoad(AliL3TrackArray tracks,Bool_t remove,TH2F hist)
4de874d1	44	{
	45	//Look at rawdata along the road specified by the track candidates.
	46
	47	if(!fHoughTransformer)
	48	{
	49	printf("AliL3HoughEval: No transformer object\n");
	50	return;
	51	}
	52
	53	AliL3Digits *pixel;
	54	Int_t patch = fHoughTransformer->fPatch;
	55	Int_t slice = fHoughTransformer->fSlice;
	56	Int_t num_of_pads_to_look = 1;
9ee44abe	57	Int_t rows_to_miss = 1;
4de874d1	58
	59	for(Int_t i=0; i<tracks->GetNTracks(); i++)
	60	{
	61	AliL3HoughTrack track = (AliL3HoughTrack)tracks->GetCheckedTrack(i);
	62	if(!track) {printf("No track\n"); return;}
	63
	64	Int_t nrow = 0;
	65
	66	for(Int_t padrow = NRows[patch][0]; padrow <= NRows[patch][1]; padrow++)
	67	{
	68
	69	Float_t xyz[3];
	70	//Get the crossing point of the track and current padrow:
	71	//if(!track->GetCrossingPoint(slice,padrow,xyz))
	72	if(!track->GetCrossingPoint(padrow,xyz))
	73	{
	74	printf("AliL3HoughEval::LookInsideRoad : Track does not cross line!!\n");
	75	continue;
	76	}
	77
4de874d1	78	Int_t npixs = 0;
	79
	80	//Get the pixels along the track candidate
	81	for(pixel=(AliL3Digits)fHoughTransformer->fRowContainer[padrow].first; pixel!=0; pixel=(AliL3Digits)pixel->nextRowPixel)
	82	{
	83
	84	Int_t sector,row;
	85	Float_t xyz_pix[3];
	86	fTransform->Slice2Sector(slice,padrow,sector,row);
	87	fTransform->Raw2Local(xyz_pix,sector,row,pixel->fPad,pixel->fTime); //y alone pad
	88
	89
	90	//check if we are inside road
	91	if(fabs(xyz_pix[1] - xyz[1]) > num_of_pads_to_look*fTransform->GetPadPitchWidthLow()) continue;
	92	npixs++;
	93
52a2a604	94	if(remove) //Remove this pixel from image
	95	pixel->fIndex = -1;
	96
4de874d1	97	if(hist)
	98	{
	99	//fTransform->Local2Global(xyz_pix,slice);
	100	hist->Fill(xyz_pix[0],xyz_pix[1],pixel->fCharge);
	101	}
	102
	103	}
	104	if(npixs > 0)
	105	{
	106	nrow++;
	107	}
	108	}
	109
	110	if(nrow < NRows[patch][1]-NRows[patch][0]-rows_to_miss)
	111	tracks->Remove(i);
	112
	113	}
	114
	115	tracks->Compress();
	116
	117
	118	}
	119
	120	void AliL3HoughEval::DisplaySlice(TH2F *hist)
	121	{
	122
	123	AliL3Digits *pixel;
	124
	125	for(Int_t padrow = 0; padrow < 174; padrow++)
	126	{
	127
	128	for(pixel=(AliL3Digits)fHoughTransformer->fRowContainer[padrow].first; pixel!=0; pixel=(AliL3Digits)pixel->nextRowPixel)
	129	{
	130
	131	Int_t sector,row;
	132	Float_t xyz_pix[3];
	133	fTransform->Slice2Sector(fHoughTransformer->fSlice,padrow,sector,row);
	134	fTransform->Raw2Local(xyz_pix,sector,row,pixel->fPad,pixel->fTime); //y alone pad
	135
	136	if(hist)
	137	{
	138	hist->Fill(xyz_pix[0],xyz_pix[1],pixel->fCharge);
	139	}
	140
	141	}
	142	}
	143
	144	}
	145
9ee44abe	146	TClonesArray AliL3HoughEval::GetParticles(Char_t rootfile)
4de874d1	147	{
4de874d1	148
4de874d1	149	TFile *file = new TFile(rootfile);
	150	file->cd();
	151
	152	AliRun gAlice = (AliRun)file->Get("gAlice");
	153	gAlice->GetEvent(0);
	154
	155	TClonesArray *particles=gAlice->Particles();
9ee44abe	156	return particles;
	157
	158	file->Close();
	159	delete file;
	160	}
	161
	162
	163	void AliL3HoughEval::CompareMC(Char_t rootfile,AliL3TrackArray merged_tracks,Float_t *eta)
	164	{
	165
	166	Int_t slice = fSlice;
	167
	168	TClonesArray *particles = GetParticles(rootfile);
4de874d1	169	Int_t n=particles->GetEntriesFast();
	170	Float_t torad=TMath::Pi()/180;
	171	Float_t phi_min = slice*20 - 10;
	172	Float_t phi_max = slice*20 + 10;
	173
	174	for (Int_t j=0; j<n; j++) {
	175	TParticle p=(TParticle)particles->UncheckedAt(j);
	176	if (p->GetFirstMother()>=0) continue; //secondary particle
	177	if(p->Eta() < eta[0] \|\| p->Eta() > eta[1]) continue;
	178	Double_t ptg=p->Pt(),pxg=p->Px(),pyg=p->Py();//,pzg=p->Pz();
	179	Double_t phi_part = TMath::ATan2(pyg,pxg);
	180	if (phi_part < 0) phi_part += 2*TMath::Pi();
	181
	182	if(phi_part < phi_mintorad \|\| phi_part > phi_maxtorad) {continue;}
	183	if(ptg<0.100) continue;
	184
	185	AliL3HoughTrack *sel_track=0;
	186
	187	Double_t min_dist = 10000;
	188
	189	for(Int_t t=0; t<merged_tracks->GetNTracks(); t++)
	190	{
	191	AliL3HoughTrack track = (AliL3HoughTrack)merged_tracks->GetCheckedTrack(t);
	192	if(!track) {printf("AliL3HoughEval: NO TRACK\n"); continue;}
	193	Float_t phi0[2] = {track->GetPhi0(),0};
	194	fTransform->Local2GlobalAngle(phi0,slice);
	195	Double_t dPt = ptg - track->GetPt();
	196	Double_t dPhi0 = phi_part - phi0[0];
	197	Double_t distance = pow(dPt,2) + pow(dPhi0,2);
	198	if(distance < min_dist)
	199	{
	200	min_dist = distance;
	201	sel_track = track;
	202	}
	203
	204	}
	205	if(sel_track)
	206	sel_track->SetBestMCid(j,min_dist);
	207
	208	Double_t dpt = fabs(ptg-sel_track->GetPt());
	209	Double_t dphi0 = fabs(phi_part-sel_track->GetPhi0());
	210	//printf("Found match, min_dist %f dPt %f dPhi0 %f\n",min_dist,dpt,dphi0);
	211	}
9ee44abe	212
4de874d1	213	}