-//Author: Anders Strand Vestbo
-//Last Modified: 28.6.01
+//$Id$
-#include <TClonesArray.h>
-#include <TH2.h>
+// Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
+//*-- Copyright © ASV
+
+#include "AliL3StandardIncludes.h"
+
+#ifdef use_root
#include <TH1.h>
#include <TFile.h>
-#include <AliRun.h>
-#include <TParticle.h>
-#include <TTree.h>
+#endif
-#include "AliTPCParam.h"
-#include "AliSimDigits.h"
+#include "AliL3Logging.h"
+#include "AliL3HoughEval.h"
+#include "AliL3MemHandler.h"
#include "AliL3TrackArray.h"
-#include "AliL3Transform.h"
-#include "AliL3HoughTransformer.h"
-#include "AliL3Defs.h"
+#include "AliL3HoughBaseTransformer.h"
+#include "AliL3DigitData.h"
#include "AliL3HoughTrack.h"
-#include "AliL3HoughEval.h"
+#include "AliL3Transform.h"
+#include "AliL3Histogram.h"
+#include "AliL3Histogram1D.h"
-ClassImp(AliL3HoughEval)
+#if GCCVERSION == 3
+using namespace std;
+#endif
+
+//_____________________________________________________________
+// AliL3HoughEval
+//
+// Evaluation class for tracklets produced by the Hough transform.
+ClassImp(AliL3HoughEval)
AliL3HoughEval::AliL3HoughEval()
{
- //Default constructor
- fTransform = new AliL3Transform();
- fHoughTransformer = 0;
- fNumOfRowsToMiss = 1;
+
+ fRemoveFoundTracks = kFALSE;
fNumOfPadsToLook = 1;
-}
-
-
-AliL3HoughEval::AliL3HoughEval(AliL3HoughTransformer *transformer)
-{
- //Constructor
- fHoughTransformer = transformer;
- fTransform = new AliL3Transform();
fNumOfRowsToMiss = 1;
- fNumOfPadsToLook = 1;
+ fEtaHistos=0;
+ fRowPointers = 0;
}
AliL3HoughEval::~AliL3HoughEval()
{
- //Destructor
- if(fTransform)
- delete fTransform;
- if(fMcTrackTable)
- delete [] fMcTrackTable;
+ fHoughTransformer = 0;
+ if(fRowPointers)
+ {
+ for(Int_t i=0; i<fNrows; i++)
+ fRowPointers[i] = 0;
+ delete [] fRowPointers;
+ }
}
+void AliL3HoughEval::InitTransformer(AliL3HoughBaseTransformer *transformer)
+{
+ fHoughTransformer = transformer;
+ fSlice = fHoughTransformer->GetSlice();
+ fPatch = fHoughTransformer->GetPatch();
+ fNrows = AliL3Transform::GetLastRow(fPatch) - AliL3Transform::GetFirstRow(fPatch) + 1;
+ fNEtaSegments = fHoughTransformer->GetNEtaSegments();
+ fEtaMin = fHoughTransformer->GetEtaMin();
+ fEtaMax = fHoughTransformer->GetEtaMax();
+ GenerateLUT();
+}
-Bool_t AliL3HoughEval::LookInsideRoad(AliL3HoughTrack *track,Int_t eta_index,Bool_t remove)
+void AliL3HoughEval::GenerateLUT()
{
- //Look at rawdata along the road specified by the track candidates.
+ //Generate a Look-up table, to limit the access to raw data
+
+ if(!fRowPointers)
+ fRowPointers = new AliL3DigitRowData*[fNrows];
+
+ AliL3DigitRowData *tempPt = (AliL3DigitRowData*)fHoughTransformer->GetDataPointer();
+ if(!tempPt)
+ printf("\nAliL3HoughEval::GenerateLUT : Zero data pointer\n");
- if(!fHoughTransformer)
+ for(Int_t i=AliL3Transform::GetFirstRow(fPatch); i<=AliL3Transform::GetLastRow(fPatch); i++)
{
- printf("\nAliL3HoughEval: No transformer object\n");
- return kFALSE;
+ Int_t prow = i - AliL3Transform::GetFirstRow(fPatch);
+ fRowPointers[prow] = tempPt;
+ AliL3MemHandler::UpdateRowPointer(tempPt);
}
- Int_t patch = fHoughTransformer->fPatch;
- Int_t slice = fHoughTransformer->fSlice;
+}
+Bool_t AliL3HoughEval::LookInsideRoad(AliL3HoughTrack *track,Int_t &nrows_crossed,Bool_t remove)
+{
+ //Look at rawdata along the road specified by the track candidates.
+ //If track is good, return true, if not return false.
Int_t sector,row;
- Int_t lut_index;
- Char_t **track_lut = fHoughTransformer->fTrackTable;
- Int_t nrow=0,npixs=0;
+
+ Int_t nrow=0,npixs=0;//,rows_crossed=0;
Float_t xyz[3];
+
+ Int_t total_charge=0;//total charge along the road
- for(Int_t padrow = NRows[patch][0]; padrow <= NRows[patch][1]; padrow++)
+ for(Int_t padrow = AliL3Transform::GetFirstRow(fPatch); padrow <= AliL3Transform::GetLastRow(fPatch); padrow++)
{
- Int_t prow = padrow - NRows[patch][0];
+ Int_t prow = padrow - AliL3Transform::GetFirstRow(fPatch);
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);
+ AliL3Transform::Slice2Sector(fSlice,padrow,sector,row);
+ AliL3Transform::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++)
+
+ //Get the timebins for this pad
+ AliL3DigitRowData *tempPt = fRowPointers[prow];
+ if(!tempPt)
{
- 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++;
+ printf("AliL3HoughEval::LookInsideRoad : Zero data pointer\n");
+ continue;
}
- if(npixs > 0)
+ //Look at both sides of the pad:
+ for(Int_t p=(Int_t)rint(xyz[1])-fNumOfPadsToLook; p<=(Int_t)rint(xyz[1])+fNumOfPadsToLook; p++)
+ {
+ AliL3DigitData *digPt = tempPt->fDigitData;
+ for(UInt_t j=0; j<tempPt->fNDigit; j++)
+ {
+ UChar_t pad = digPt[j].fPad;
+ Int_t charge = digPt[j].fCharge;
+ if(charge <= fHoughTransformer->GetLowerThreshold()) continue;
+ if(pad < p) continue;
+ if(pad > p) break;
+ UShort_t time = digPt[j].fTime;
+ Double_t eta = AliL3Transform::GetEta(fSlice,padrow,pad,time);
+ Int_t pixel_index = fHoughTransformer->GetEtaIndex(eta);
+ if(pixel_index != track->GetEtaIndex()) continue;
+ total_charge += digPt[j].fCharge;
+ if(remove)
+ digPt[j].fCharge = 0; //Erase the track from image
+ npixs++;
+ }
+ }
+
+ if(npixs > 1)//At least 2 digits on this padrow
{
nrow++;
}
}
- if(nrow >= NRows[patch][1]-NRows[patch][0]-fNumOfRowsToMiss)//this was a good track
+ if(remove)
+ return kTRUE;
+
+ nrows_crossed += nrow; //Update the number of rows crossed.
+
+ if(nrow >= AliL3Transform::GetNRows(fPatch) - fNumOfRowsToMiss)//this was a good track
{
- track->SetEtaIndex(eta_index);
- if(remove)
- RemoveTrackFromImage(track,eta_index);
+ if(fRemoveFoundTracks)
+ {
+ Int_t dummy=0;
+ LookInsideRoad(track,dummy,kTRUE);
+ }
return kTRUE;
}
else
return kFALSE;
}
-void AliL3HoughEval::LookInsideRawRoad(AliL3TrackArray *tracks,Int_t eta_index,Bool_t remove)
+void AliL3HoughEval::FindEta(AliL3TrackArray *tracks)
{
- //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++)
+ Int_t ntracks = tracks->GetNTracks();
+ fEtaHistos = new AliL3Histogram1D*[ntracks];
+
+ Char_t hname[100];
+ for(Int_t i=0; i<ntracks; 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++)
+ sprintf(hname,"etahist_%d",i);
+ fEtaHistos[i] = new AliL3Histogram1D(hname,hname,100,0,1);
+ }
+ Double_t etaslice = (fEtaMax - fEtaMin)/fNEtaSegments;
+
+ for(Int_t ntr=0; ntr<ntracks; ntr++)
+ {
+ AliL3HoughTrack *track = (AliL3HoughTrack*)tracks->GetCheckedTrack(ntr);
+ if(!track) continue;
+ for(Int_t padrow = AliL3Transform::GetFirstRow(fPatch); padrow <= AliL3Transform::GetLastRow(fPatch); padrow++)
{
- Int_t prow = padrow - NRows[patch][0];
+ Int_t prow = padrow - AliL3Transform::GetFirstRow(fPatch);
+
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++)
+
+ AliL3Transform::Slice2Sector(fSlice,padrow,sector,row);
+ AliL3Transform::Local2Raw(xyz,sector,row);
+
+ //Get the timebins for this pad
+ AliL3DigitRowData *tempPt = fRowPointers[prow];
+ if(!tempPt)
{
- 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++;
+ printf("AliL3HoughEval::LookInsideRoad : Zero data pointer\n");
+ continue;
}
- if(npixs > 0)
+ //Look at both sides of the pad:
+ for(Int_t p=(Int_t)rint(xyz[1])-fNumOfPadsToLook; p<=(Int_t)rint(xyz[1])+fNumOfPadsToLook; p++)
{
- nrow++;
- }
+ AliL3DigitData *digPt = tempPt->fDigitData;
+ for(UInt_t j=0; j<tempPt->fNDigit; j++)
+ {
+ UChar_t pad = digPt[j].fPad;
+ Int_t charge = digPt[j].fCharge;
+ if(charge <= fHoughTransformer->GetLowerThreshold()) continue;
+ if(pad < p) continue;
+ if(pad > p) break;
+ UShort_t time = digPt[j].fTime;
+ Double_t eta = AliL3Transform::GetEta(fSlice,padrow,pad,time);
+ Int_t pixel_index = (Int_t)(eta/etaslice);
+ if(pixel_index > track->GetEtaIndex()+1) continue;
+ if(pixel_index < track->GetEtaIndex()-1) break;
+ fEtaHistos[ntr]->Fill(eta,digPt[j].fCharge);
+ }
+ }
}
- 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.
}
+ for(Int_t i=0; i<ntracks; i++)
+ {
+ AliL3Histogram1D *hist = fEtaHistos[i];
+ Int_t max_bin = hist->GetMaximumBin();
+ Double_t max_value = hist->GetBinContent(max_bin);
+ AliL3HoughTrack *track = (AliL3HoughTrack*)tracks->GetCheckedTrack(i);
+ if(!track) continue;
+ if(hist->GetBinContent(max_bin-1)<max_value && hist->GetBinContent(max_bin+1)<max_value)
+ {
+ track->SetWeight((Int_t)max_value,kTRUE);
+ track->SetEta(hist->GetBinCenter(max_bin));
+ track->SetNHits(track->GetWeight());
+ }
+ else
+ {
+ track->SetWeight(0);
+ tracks->Remove(i); //remove this track, because it was not a peak
+ }
+ }
tracks->Compress();
+
+ //for(Int_t i=0; i<ntracks; i++)
+ //delete fEtaHistos[i];
+ //delete [] fEtaHistos;
}
-void AliL3HoughEval::RemoveTrackFromImage(AliL3HoughTrack *track,Int_t eta_index)
+void AliL3HoughEval::DisplayEtaSlice(Int_t eta_index,AliL3Histogram *hist)
{
- //Remove the pixels along the track in the image.
-
- Int_t patch = fHoughTransformer->fPatch;
- Int_t slice = fHoughTransformer->fSlice;
+ //Display the current raw data inside the (slice,patch)
- Int_t maxnpads = fTransform->GetNPads(NRows[patch][1]);
+ if(!hist)
+ {
+ printf("AliL3HoughEval::DisplayEtaSlice : No input histogram!\n");
+ return;
+ }
- 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++)
+ for(Int_t padrow = AliL3Transform::GetFirstRow(fPatch); padrow <= AliL3Transform::GetLastRow(fPatch); padrow++)
{
- Int_t prow = padrow - NRows[patch][0];
- if(!track->GetCrossingPoint(padrow,xyz))
+ Int_t prow = padrow - AliL3Transform::GetFirstRow(fPatch);
+
+ AliL3DigitRowData *tempPt = fRowPointers[prow];
+ if(!tempPt)
{
- printf("AliL3HoughEval::LookInsideRoad : Track does not cross line!!\n");
+ printf("AliL3HoughEval::DisplayEtaSlice : Zero data pointer\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++)
+ AliL3DigitData *digPt = tempPt->fDigitData;
+ if((Int_t)tempPt->fRow != padrow)
+ {
+ printf("\nAliL3HoughEval::DisplayEtaSlice : Mismatching padrows!!!\n");
+ return;
+ }
+ for(UInt_t j=0; j<tempPt->fNDigit; j++)
{
- if(p<0 || p>fTransform->GetNPads(padrow)) continue;
- lut_index = (prow<<8) + p;
- track_lut[eta_index][lut_index] = -1; //remove it!!
+ UChar_t pad = digPt[j].fPad;
+ UChar_t charge = digPt[j].fCharge;
+ UShort_t time = digPt[j].fTime;
+ if((Int_t)charge <= fHoughTransformer->GetLowerThreshold() || (Int_t)charge >= fHoughTransformer->GetUpperThreshold()) continue;
+ Float_t xyz[3];
+ Int_t sector,row;
+ AliL3Transform::Slice2Sector(fSlice,padrow,sector,row);
+ AliL3Transform::Raw2Local(xyz,sector,row,pad,time);
+ Double_t eta = AliL3Transform::GetEta(xyz);
+ Int_t pixel_index = fHoughTransformer->GetEtaIndex(eta);//(Int_t)(eta/etaslice);
+ if(pixel_index != eta_index) continue;
+ hist->Fill(xyz[0],xyz[1],charge);
}
}
-
}
-void AliL3HoughEval::DisplaySlice(TH2F *hist)
+#ifdef use_root
+void AliL3HoughEval::CompareMC(AliL3TrackArray *tracks,Char_t *trackfile,Int_t threshold)
{
-
- AliL3Digits *pixel;
-
- for(Int_t padrow = 0; padrow < 174; padrow++)
+ /*
+ struct GoodTrack goodtracks[15000];
+ Int_t nt=0;
+ ifstream in(trackfile);
+ if(in)
{
-
- for(pixel=(AliL3Digits*)fHoughTransformer->fRowContainer[padrow].first; pixel!=0; pixel=(AliL3Digits*)pixel->nextRowPixel)
+ printf("Reading good tracks from file %s\n",trackfile);
+ while (in>>goodtracks[nt].label>>goodtracks[nt].code>>
+ goodtracks[nt].px>>goodtracks[nt].py>>goodtracks[nt].pz>>
+ goodtracks[nt].pt>>goodtracks[nt].eta>>goodtracks[nt].nhits)
{
-
- 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)
+ nt++;
+ if (nt==15000)
{
- hist->Fill(xyz_pix[0],xyz_pix[1],pixel->fCharge);
- }
-
+ cerr<<"Too many good tracks"<<endl;
+ break;
+ }
+ }
+ if (!in.eof())
+ {
+ LOG(AliL3Log::kError,"AliL3HoughEval::CompareMC","Input file")
+ <<"Error in file reading"<<ENDLOG;
+ return;
}
}
+ else
+ {
+ LOG(AliL3Log::kError,"AliL3HoughEval::CompareMC","Input")
+ <<"No input trackfile "<<trackfile<<ENDLOG;
+ }
-}
-
-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);
+ Int_t *particles = new Int_t[fNEtaSegments];
+ Int_t *ftracks = new Int_t[fNEtaSegments];
+ for(Int_t i=0; i<fNEtaSegments; i++)
+ {
+ particles[i]=0;
+ ftracks[i]=0;
+ }
- TClonesArray *particles=gAlice->Particles();
- Int_t np = particles->GetEntriesFast();
+ TH1F *ptgood = new TH1F("ptgood","ptgood",5,0,2);
+ TH1F *ptfound = new TH1F("ptfound","ptgood",5,0,2);
+ TH1F *pteff = new TH1F("pteff","pteff",5,0,2);
+ TH1F *etafound = new TH1F("etafound","etafound",5,0,1);
+ TH1F *etagood = new TH1F("etagood","etagood",5,0,1);
+ TH1F *etaeff = new TH1F("etaeff","etaeff",5,0,1);
- 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++)
+ Double_t etaslice = (fEtaMax - fEtaMin)/fNEtaSegments;
+ for(Int_t i=0; i<tracks->GetNTracks(); 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;
- }
+ AliL3HoughTrack *tr = (AliL3HoughTrack*)tracks->GetCheckedTrack(i);
+ if(!tr) continue;
+ if(tr->GetWeight()<threshold) continue;
+ Int_t trackindex = tr->GetEtaIndex();
+ if(trackindex <0 || trackindex >= fNEtaSegments) continue;
+ ftracks[trackindex]++;
+ ptfound->Fill(tr->GetPt());
+ etafound->Fill(tr->GetEta());
+ }
+ for(Int_t i=0; i<nt; i++)
+ {
+ if(goodtracks[i].nhits < 174) continue;
+ if(goodtracks[i].pt < 0.2) continue;
+ Int_t particleindex = (Int_t)(goodtracks[i].eta/etaslice);
+ if(particleindex < 0 || particleindex >= fNEtaSegments) continue;
+ particles[particleindex]++;
+ ptgood->Fill(goodtracks[i].pt);
+ etagood->Fill(goodtracks[i].eta);
}
- 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;
-}
-
+ Double_t found=0;
+ Double_t good =0;
+ for(Int_t i=0; i<fNEtaSegments; i++)
+ {
+ //printf("Slice %d : Found tracks %d, good tracks %d\n",i,ftracks[i],particles[i]);
+ found += ftracks[i];
+ good += particles[i];
+ }
+ printf("And the total efficiency was: %f\n",found/good);
-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);
- }
+ ptgood->Sumw2(); ptfound->Sumw2();
+ etagood->Sumw2(); etafound->Sumw2();
+ pteff->Divide(ptfound,ptgood,1,1,"b");
+ etaeff->Divide(etafound,etagood,1,1,"b");
+ TFile *file = TFile::Open("eff.root","RECREATE");
+ ptgood->Write();
+ ptfound->Write();
+ pteff->Write();
+ etafound->Write();
+ etagood->Write();
+ etaeff->Write();
file->Close();
- delete file;
+ delete [] particles;
+ delete [] ftracks;
+ */
}
+
+#endif
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff