From 4de874d19c7239e7e86e8b8c82722faa42fcf63d Mon Sep 17 00:00:00 2001 From: vestbo Date: Thu, 22 Mar 2001 12:51:56 +0000 Subject: [PATCH 1/1] This commit was generated by cvs2svn to compensate for changes in r3174, which included commits to RCS files with non-trunk default branches. --- HLT/hough/AliL3Defs.h | 12 + HLT/hough/AliL3Histogram.cxx | 17 + HLT/hough/AliL3Histogram.h | 20 + HLT/hough/AliL3HoughEval.cxx | 205 ++++++++++ HLT/hough/AliL3HoughEval.h | 37 ++ HLT/hough/AliL3HoughLinkDef.h | 15 + HLT/hough/AliL3HoughMaxFinder.cxx | 113 ++++++ HLT/hough/AliL3HoughMaxFinder.h | 31 ++ HLT/hough/AliL3HoughMerge.cxx | 71 ++++ HLT/hough/AliL3HoughMerge.h | 31 ++ HLT/hough/AliL3HoughPixel.cxx | 45 +++ HLT/hough/AliL3HoughPixel.h | 45 +++ HLT/hough/AliL3HoughTransformer.cxx | 580 ++++++++++++++++++++++++++++ HLT/hough/AliL3HoughTransformer.h | 68 ++++ HLT/hough/AliL3Transform.cxx | 574 +++++++++++++++++++++++++++ HLT/hough/AliL3Transform.h | 65 ++++ HLT/hough/Makefile | 79 ++++ HLT/hough/hough.C | 143 +++++++ HLT/hough/hough_line.C | 126 ++++++ HLT/hough/hough_line_merge.C | 205 ++++++++++ HLT/hough/hough_merge.C | 148 +++++++ HLT/hough/hough_mergehistos.C | 63 +++ HLT/hough/rootlogon.C | 21 + 23 files changed, 2714 insertions(+) create mode 100644 HLT/hough/AliL3Defs.h create mode 100644 HLT/hough/AliL3Histogram.cxx create mode 100644 HLT/hough/AliL3Histogram.h create mode 100644 HLT/hough/AliL3HoughEval.cxx create mode 100644 HLT/hough/AliL3HoughEval.h create mode 100644 HLT/hough/AliL3HoughLinkDef.h create mode 100644 HLT/hough/AliL3HoughMaxFinder.cxx create mode 100644 HLT/hough/AliL3HoughMaxFinder.h create mode 100644 HLT/hough/AliL3HoughMerge.cxx create mode 100644 HLT/hough/AliL3HoughMerge.h create mode 100644 HLT/hough/AliL3HoughPixel.cxx create mode 100644 HLT/hough/AliL3HoughPixel.h create mode 100644 HLT/hough/AliL3HoughTransformer.cxx create mode 100644 HLT/hough/AliL3HoughTransformer.h create mode 100644 HLT/hough/AliL3Transform.cxx create mode 100644 HLT/hough/AliL3Transform.h create mode 100644 HLT/hough/Makefile create mode 100644 HLT/hough/hough.C create mode 100644 HLT/hough/hough_line.C create mode 100644 HLT/hough/hough_line_merge.C create mode 100644 HLT/hough/hough_merge.C create mode 100644 HLT/hough/hough_mergehistos.C create mode 100644 HLT/hough/rootlogon.C diff --git a/HLT/hough/AliL3Defs.h b/HLT/hough/AliL3Defs.h new file mode 100644 index 00000000000..c9bfdf052e5 --- /dev/null +++ b/HLT/hough/AliL3Defs.h @@ -0,0 +1,12 @@ +#ifndef _ALIL3DEFS_H_ +#define _ALIL3DEFS_H_ + +#include + +const Int_t NRowsSlice = 174; +const Int_t NRows[5][2] = {{ 0, 45},{46,77},{78,109},{110,141},{142,173}}; +const Double_t Pi = 3.14159265358979323846; +const Double_t ToRad = Pi/180.; + + +#endif /* _ALIL3DEFS_H_ */ diff --git a/HLT/hough/AliL3Histogram.cxx b/HLT/hough/AliL3Histogram.cxx new file mode 100644 index 00000000000..ccc5460968a --- /dev/null +++ b/HLT/hough/AliL3Histogram.cxx @@ -0,0 +1,17 @@ + +#include "AliL3Histogram.h" + +ClassImp(AliL3Histogram) + + +AliL3Histogram::AliL3Histogram() +{ + +} + + +AliL3Histogram::~AliL3Histogram() +{ + //Destructor + +} diff --git a/HLT/hough/AliL3Histogram.h b/HLT/hough/AliL3Histogram.h new file mode 100644 index 00000000000..b31e1c433a0 --- /dev/null +++ b/HLT/hough/AliL3Histogram.h @@ -0,0 +1,20 @@ +#ifndef ALIL3_HOUGHHISTOGRAM +#define ALIL3_HOUGHHISTOGRAM + +#include "AliL3RootTypes.h" + +class AliL3HoughTransformer : public TH2F { + + private: + + + + public: + AliL3Histogram(); + + + ClassDef(AliL3Histogram,1) + +}; + +#endif diff --git a/HLT/hough/AliL3HoughEval.cxx b/HLT/hough/AliL3HoughEval.cxx new file mode 100644 index 00000000000..9e0aaca1c96 --- /dev/null +++ b/HLT/hough/AliL3HoughEval.cxx @@ -0,0 +1,205 @@ +#include +#include +#include +#include +#include +#include + +#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,TH2F *hist,TH2F *fake) +{ + //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 = 2; + + for(Int_t i=0; iGetNTracks(); 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; + } + + if(fake) + fake->Fill(xyz[0],xyz[1],1); + 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(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); + } + + } + } + +} + + +void AliL3HoughEval::CompareMC(Char_t *rootfile,AliL3TrackArray *merged_tracks,Float_t *eta) +{ + + Int_t slice = 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; jUncheckedAt(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; tGetNTracks(); 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; + + +} diff --git a/HLT/hough/AliL3HoughEval.h b/HLT/hough/AliL3HoughEval.h new file mode 100644 index 00000000000..42ee1293ee5 --- /dev/null +++ b/HLT/hough/AliL3HoughEval.h @@ -0,0 +1,37 @@ +#ifndef ALIL3_HOUGH_Eval +#define ALIL3_HOUGH_Eval + +#include "AliL3RootTypes.h" + +class TClonesArray; +class TObjArray; +class TH2F; +class TH1F; +class AliL3Transform; +class AliL3TrackArray; +class AliL3HoughTransformer; + +class AliL3HoughEval : public TObject { + + private: + + Int_t fSlice; + AliL3HoughTransformer *fHoughTransformer; + AliL3Transform *fTransform; //! + + public: + AliL3HoughEval(); + AliL3HoughEval(AliL3HoughTransformer *transformer); + virtual ~AliL3HoughEval(); + + void LookInsideRoad(AliL3TrackArray *tracks,TH2F *hist=0,TH2F *fake=0); + void DisplaySlice(TH2F *hist); + void CompareMC(Char_t *rootfile,AliL3TrackArray *merged_tracks,Float_t *eta); + + void SetTransformer(AliL3HoughTransformer *t) {fHoughTransformer=t;} + + ClassDef(AliL3HoughEval,1) + +}; + +#endif diff --git a/HLT/hough/AliL3HoughLinkDef.h b/HLT/hough/AliL3HoughLinkDef.h new file mode 100644 index 00000000000..534862b547f --- /dev/null +++ b/HLT/hough/AliL3HoughLinkDef.h @@ -0,0 +1,15 @@ +#ifdef __CINT__ + +#pragma link off all globals; +#pragma link off all classes; +#pragma link off all functions; + +#pragma link C++ class AliL3HoughTransformer; +#pragma link C++ class AliL3HoughPixel; +//#pragma link C++ class AliL3HoughTrack; +#pragma link C++ class AliL3HoughMaxFinder; +#pragma link C++ class AliL3HoughEval; +#pragma link C++ class AliL3HoughMerge; + +#endif + diff --git a/HLT/hough/AliL3HoughMaxFinder.cxx b/HLT/hough/AliL3HoughMaxFinder.cxx new file mode 100644 index 00000000000..d36b28c6e36 --- /dev/null +++ b/HLT/hough/AliL3HoughMaxFinder.cxx @@ -0,0 +1,113 @@ +#include + +#include + +#include "AliL3TrackArray.h" +#include "AliL3HoughTrack.h" +#include "AliL3HoughMaxFinder.h" + +ClassImp(AliL3HoughMaxFinder) + + +AliL3HoughMaxFinder::AliL3HoughMaxFinder() +{ + //Default constructor + fThreshold = 0; + //fTracks = 0; + fHistoType=0; +} + + +AliL3HoughMaxFinder::AliL3HoughMaxFinder(Char_t *histotype) +{ + //Constructor + + //fTracks = new AliL3TrackArray("AliL3HoughTrack"); + if(strcmp(histotype,"KappaPhi")==0) fHistoType='c'; + if(strcmp(histotype,"DPsi")==0) fHistoType='l'; + +} + + +AliL3HoughMaxFinder::~AliL3HoughMaxFinder() +{ + //Destructor + +} + + +AliL3TrackArray *AliL3HoughMaxFinder::FindMaxima(TH2F *hist,Int_t *rowrange,Int_t ref_row) +{ + //Locate all the maxima in input histogram. + //Maxima is defined as bins with more entries than the + //immediately neighbouring bins. + + Int_t xmin = hist->GetXaxis()->GetFirst(); + Int_t xmax = hist->GetXaxis()->GetLast(); + Int_t ymin = hist->GetYaxis()->GetFirst(); + Int_t ymax = hist->GetYaxis()->GetLast(); + Int_t bin[9],track_counter=0; + Stat_t value[9]; + + AliL3TrackArray *tracks = new AliL3TrackArray("AliL3HoughTrack"); + AliL3HoughTrack *track; + + for(Int_t xbin=xmin+1; xbinGetBin(xbin-1,ybin-1); + bin[1] = hist->GetBin(xbin,ybin-1); + bin[2] = hist->GetBin(xbin+1,ybin-1); + bin[3] = hist->GetBin(xbin-1,ybin); + bin[4] = hist->GetBin(xbin,ybin); + bin[5] = hist->GetBin(xbin+1,ybin); + bin[6] = hist->GetBin(xbin-1,ybin+1); + bin[7] = hist->GetBin(xbin,ybin+1); + bin[8] = hist->GetBin(xbin+1,ybin+1); + value[0] = hist->GetBinContent(bin[0]); + value[1] = hist->GetBinContent(bin[1]); + value[2] = hist->GetBinContent(bin[2]); + value[3] = hist->GetBinContent(bin[3]); + value[4] = hist->GetBinContent(bin[4]); + value[5] = hist->GetBinContent(bin[5]); + value[6] = hist->GetBinContent(bin[6]); + value[7] = hist->GetBinContent(bin[7]); + value[8] = hist->GetBinContent(bin[8]); + + if(value[4] <= fThreshold) continue;//central bin below threshold + + if(value[4]>value[0] && value[4]>value[1] && value[4]>value[2] + && value[4]>value[3] && value[4]>value[5] && value[4]>value[6] + && value[4]>value[7] && value[4]>value[8]) + { + //Found a local maxima + Float_t max_x = hist->GetXaxis()->GetBinCenter(xbin); + Float_t max_y = hist->GetYaxis()->GetBinCenter(ybin); + + track = (AliL3HoughTrack*)tracks->NextTrack(); + + + switch(fHistoType) + { + case 'c': + track->SetTrackParameters(max_x,max_y,(Int_t)value[4]); + break; + case 'l': + track->SetLineParameters(max_x,max_y,(Int_t)value[4],rowrange,ref_row); + break; + default: + printf("AliL3HoughMaxFinder: Error in tracktype\n"); + } + + track_counter++; + + + } + else + continue; //not a maxima + } + } + tracks->QSort(); + return tracks; +} diff --git a/HLT/hough/AliL3HoughMaxFinder.h b/HLT/hough/AliL3HoughMaxFinder.h new file mode 100644 index 00000000000..e937dda1cc6 --- /dev/null +++ b/HLT/hough/AliL3HoughMaxFinder.h @@ -0,0 +1,31 @@ +#ifndef ALIL3_HOUGH_MaxFinder +#define ALIL3_HOUGH_MaxFinder + +#include "AliL3RootTypes.h" + +class AliL3TrackArray; + +class AliL3HoughMaxFinder : public TObject { + + private: + + Int_t fThreshold; + //AliL3TrackArray *fTracks; //! + + Char_t fHistoType; + + public: + AliL3HoughMaxFinder(); + AliL3HoughMaxFinder(Char_t *histotype); + virtual ~AliL3HoughMaxFinder(); + + AliL3TrackArray *FindMaxima(TH2F *hist,Int_t *rowrange=0,Int_t ref_row=0); + void SetThreshold(Int_t f) {fThreshold = f;} + + //AliL3TrackArray *GetTracks() {return fTracks;} + + ClassDef(AliL3HoughMaxFinder,1) + +}; + +#endif diff --git a/HLT/hough/AliL3HoughMerge.cxx b/HLT/hough/AliL3HoughMerge.cxx new file mode 100644 index 00000000000..534e329c8c3 --- /dev/null +++ b/HLT/hough/AliL3HoughMerge.cxx @@ -0,0 +1,71 @@ +#include + +#include "AliL3TrackArray.h" +#include "AliL3HoughTrack.h" +#include "AliL3HoughMerge.h" +#include "AliL3HoughTransformer.h" + +ClassImp(AliL3HoughMerge) + + +AliL3HoughMerge::AliL3HoughMerge() +{ + //Default constructor + + fInTracks = NULL; + fOutTracks = NULL; + fNPatches = 5; +} + + +AliL3HoughMerge::AliL3HoughMerge(Int_t slice,Int_t row_patches) +{ + //Constructor + + //fInTracks = (AliL3TrackArray**)new Byte_t[row_patches*sizeof(AliL3TrackArray*)]; + fInTracks = new AliL3TrackArray*[row_patches]; + fNPatches = row_patches; + for(Int_t i=0; iAddTracks(tracks); //copies tracks to new trackarray. Does not delete the track objects. +} + +void AliL3HoughMerge::FillHisto(TH2F *merge_hist) +{ + + for(Int_t pat=0; pat < fNPatches; pat++) + { + for(Int_t t=0; tGetNTracks(); t++) + { + AliL3HoughTrack *tr = (AliL3HoughTrack*)fInTracks[pat]->GetCheckedTrack(t); + if(!tr) {printf("AliL3HoughMerge NO TRACK\n"); continue;} + merge_hist->Fill(tr->GetKappa(),tr->GetPhi0(),tr->GetNHits()); + } + } +} + +/* +void AliL3HoughMerge::MergeLines() +{ + + + +} +*/ diff --git a/HLT/hough/AliL3HoughMerge.h b/HLT/hough/AliL3HoughMerge.h new file mode 100644 index 00000000000..cf3653b064e --- /dev/null +++ b/HLT/hough/AliL3HoughMerge.h @@ -0,0 +1,31 @@ +#ifndef ALIL3_HOUGH_Merge +#define ALIL3_HOUGH_Merge + +#include "AliL3RootTypes.h" + +class TClonesArray; +class TObjArray; +class TH2F; +class AliL3TrackArray; + +class AliL3HoughMerge : public TObject { + + private: + + AliL3TrackArray **fInTracks; //! + AliL3TrackArray *fOutTracks; //! + + Int_t fNPatches; + public: + AliL3HoughMerge(); + AliL3HoughMerge(Int_t slice,Int_t row_patches=5); + virtual ~AliL3HoughMerge(); + + void FillTracks(AliL3TrackArray *tracks,Int_t patch); + void FillHisto(TH2F *merge_hist); + + ClassDef(AliL3HoughMerge,1) + +}; + +#endif diff --git a/HLT/hough/AliL3HoughPixel.cxx b/HLT/hough/AliL3HoughPixel.cxx new file mode 100644 index 00000000000..f7c89a88510 --- /dev/null +++ b/HLT/hough/AliL3HoughPixel.cxx @@ -0,0 +1,45 @@ + +#include "AliL3Transform.h" +#include "AliL3HoughPixel.h" + +ClassImp(AliL3HoughPixel) + +AliL3HoughPixel::AliL3HoughPixel() +{ + +} + +AliL3HoughPixel::AliL3HoughPixel(Int_t pad,Int_t time,Short_t signal) + // Int_t slice,Int_t padrow) +{ + //Constructor + fPad = pad; + fTime = time; + fSignal = signal; + /*fSlice = slice; + fPadrow = padrow; + */ +} + +AliL3HoughPixel::~AliL3HoughPixel() +{ + +} +/* +void AliL3HoughPixel::SetCoordinates() +{ + + AliL3Transform *transform = new AliL3Transform(); + + Float_t xyz[3]; + Int_t sector,row; + transform->Slice2Sector(fSlice,fPadrow,sector,row); + transform->Raw2Global(xyz,sector,row,fPad,fTime); + fX = xyz[0]; + fY = xyz[1]; + fZ = xyz[2]; + fPhi = transform->GetPhi(xyz); + fEta = transform->GetEta(xyz); + +} +*/ diff --git a/HLT/hough/AliL3HoughPixel.h b/HLT/hough/AliL3HoughPixel.h new file mode 100644 index 00000000000..734ac771929 --- /dev/null +++ b/HLT/hough/AliL3HoughPixel.h @@ -0,0 +1,45 @@ +#ifndef ALI_PIXEL +#define ALI_PIXEL + +#include + +class AliL3HoughPixel : public TObject { + + private: + Short_t fSignal; + Int_t fPad; + Int_t fTime; + /*Int_t fSlice; + Int_t fPadrow; + + Float_t fX; + Float_t fY; + Float_t fZ; + Double_t fPhi; + Double_t fEta; + */ + public: + + AliL3HoughPixel(); + //AliL3HoughPixel(Int_t pad,Int_t time,Short_t signal,Int_t slice,Int_t padrow); + AliL3HoughPixel(Int_t pad,Int_t time,Short_t signal); + virtual ~AliL3HoughPixel(); + + /*Int_t GetPad() {return fPad;} + Int_t GetTime() {return fTime;} + Short_t GetSignal() {return fSignal;} + Double_t GetPhi() {return fPhi;} + Double_t GetEta() {return fEta;} + Float_t GetX() {return fX;} + Float_t GetY() {return fY;} + Float_t GetZ() {return fZ;} + + void SetCoordinates(); + */ + AliL3HoughPixel *nextRowPixel; + + ClassDef(AliL3HoughPixel,1) +}; + + +#endif diff --git a/HLT/hough/AliL3HoughTransformer.cxx b/HLT/hough/AliL3HoughTransformer.cxx new file mode 100644 index 00000000000..5dfeccdea15 --- /dev/null +++ b/HLT/hough/AliL3HoughTransformer.cxx @@ -0,0 +1,580 @@ +#include + +#include +#include +#include + +#include "AliTPCParam.h" +#include "AliSimDigits.h" + +#include "AliL3Defs.h" +#include "AliL3Transform.h" +#include "AliL3HoughPixel.h" +#include "AliL3HoughTransformer.h" +#include "AliL3HoughTrack.h" +#include "AliL3TrackArray.h" + +ClassImp(AliL3HoughTransformer) + +AliL3HoughTransformer::AliL3HoughTransformer() +{ + //Default constructor + fTransform = 0; + fEtaMin = 0; + fEtaMax = 0; + fSlice = 0; + fPatch = 0; + fRowContainer = 0; + fPhiRowContainer = 0; + fNumOfPadRows=0; + fContainerBounds=0; + fNDigits=0; + fIndex = 0; +} + + +AliL3HoughTransformer::AliL3HoughTransformer(Int_t slice,Int_t patch,Float_t *etarange,Int_t phi_segments) +{ + //Constructor + + fTransform = new AliL3Transform(); + + + fEtaMin = etarange[0]; + fEtaMax = etarange[1]; + fSlice = slice; + fPatch = patch; + fNPhiSegments = phi_segments; + fNumOfPadRows=NRowsSlice; +} + + +AliL3HoughTransformer::~AliL3HoughTransformer() +{ + //Destructor + if(fRowContainer) + delete [] fRowContainer; + if(fTransform) + delete fTransform; + if(fPhiRowContainer) + delete [] fPhiRowContainer; + if(fIndex) + { + for(Int_t i=0; iGetYaxis()->GetFirst(); + Int_t ymax = hist->GetYaxis()->GetLast(); + Int_t nbinsy = hist->GetNbinsY(); + + fIndex = new Int_t*[fNDigits]; + for(Int_t i=0; inextRowPixel) + { + Float_t xyz[3]; + fTransform->Slice2Sector(fSlice,padrow,sector,row); + fTransform->Raw2Local(xyz,sector,row,pixel->fPad,pixel->fTime); + + Double_t r_pix = sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]); + + Double_t phi_pix = fTransform->GetPhi(xyz); + Short_t signal = pixel->fCharge; + Int_t index = pixel->fIndex; + if(index >= fNDigits) + printf("AliL3HoughTransformer::InitTemplates : Index error! %d\n",index); + for(Int_t p=ymin; p<=ymax; p++) + { + + Double_t phi0 = hist->GetYaxis()->GetBinCenter(p); + Double_t kappa = 2*sin(phi_pix-phi0)/r_pix; + + Int_t bin = hist->FindBin(kappa,phi0); + if(fIndex[index][p]!=0) + printf("AliL3HoughTransformer::InitTemplates : Overlapping indexes\n"); + fIndex[index][p] = bin; + } + } + } + +} + + +void AliL3HoughTransformer::CountBins() +{ + + Int_t middle_row = 87; //middle of the slice + + Double_t r_in_bundle = fTransform->Row2X(middle_row); + // Double_t phi_min = (fSlice*20 - 10)*ToRad; + //Double_t phi_max = (fSlice*20 + 10)*ToRad; + Double_t phi_min = -15*ToRad; + Double_t phi_max = 15*ToRad; + + Double_t phi_slice = (phi_max - phi_min)/fNPhiSegments; + Double_t min_phi0 = 10000; + Double_t max_phi0 = 0; + Double_t min_kappa = 100000; + Double_t max_kappa = 0; + + Int_t xbin = 60; + Int_t ybin = 60; + Float_t xrange[2] = {-0.006 , 0.006}; //Pt 0.2-> + Float_t yrange[2] = {-0.26 , 0.26}; //slice 2 0.55->0.88 + + TH2F *histo = new TH2F("histo","Parameter space",xbin,xrange[0],xrange[1],ybin,yrange[0],yrange[1]); + + for(Int_t padrow=NRows[fPatch][0]; padrow <= NRows[fPatch][1]; padrow++) + { + for(Int_t pad=0; pad < fTransform->GetNPads(padrow); pad++) + { + for(Int_t time = 0; time < fTransform->GetNTimeBins(); time++) + { + Float_t xyz[3]; + Int_t sector,row; + fTransform->Slice2Sector(fSlice,padrow,sector,row); + fTransform->Raw2Global(xyz,sector,row,pad,time); + Double_t eta = fTransform->GetEta(xyz); + if(eta < fEtaMin || eta > fEtaMax) continue; + fTransform->Global2Local(xyz,sector); + Double_t r_pix = sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]); + Double_t phi_pix = fTransform->GetPhi(xyz); + + for(Int_t p=0; p<=fNPhiSegments; p++) + { + Double_t phi_in_bundle = phi_min + p*phi_slice; + + Double_t tanPhi0 = (r_pix*sin(phi_in_bundle)-r_in_bundle*sin(phi_pix))/(r_pix*cos(phi_in_bundle)-r_in_bundle*cos(phi_pix)); + + Double_t phi0 = atan(tanPhi0); + // if(phi0 < 0.55 || phi0 > 0.88) continue; + + //if(phi0 < 0) phi0 = phi0 +2*Pi; + //Double_t kappa = sin(phi_in_bundle - phi0)*2/r_in_bundle; + + Double_t angle = phi_pix - phi0; + Double_t kappa = 2*sin(angle)/r_pix; + histo->Fill(kappa,phi0,1); + if(phi0 < min_phi0) + min_phi0 = phi0; + if(phi0 > max_phi0) + max_phi0 = phi0; + if(kappa < min_kappa) + min_kappa = kappa; + if(kappa > max_kappa) + max_kappa = kappa; + + } + + } + + } + + } + Int_t count=0,bi=0; + + Int_t xmin = histo->GetXaxis()->GetFirst(); + Int_t xmax = histo->GetXaxis()->GetLast(); + Int_t ymin = histo->GetYaxis()->GetFirst(); + Int_t ymax = histo->GetYaxis()->GetLast(); + + + for(Int_t xbin=xmin+1; xbinGetBin(xbin,ybin); + if(histo->GetBinContent(bin)>0) + count++; + } + } + + + printf("Number of possible tracks in this region %d, bins %d\n",count,bi); + printf("Phi, min %f max %f\n",min_phi0,max_phi0); + printf("Kappa, min %f max %f\n",min_kappa,max_kappa); + histo->Draw("box"); +} + + +void AliL3HoughTransformer::Transform2Circle(TH2F *hist,Int_t middle_row) +{ + //Transformation is done with respect to local coordinates in slice. + + + AliL3Digits *pix1; + Int_t sector,row; + + //Define a common point + + /* + Double_t rowdist1 = fTransform->Row2X(middle_row-1); + Double_t rowdist2 = fTransform->Row2X(middle_row); + Double_t r_in_bundle = rowdist1 + (rowdist1-rowdist2)/2; + */ + + Double_t r_in_bundle = fTransform->Row2X(middle_row); + //Make overlap between slices + Double_t phi_min = -15*ToRad; + Double_t phi_max = 15*ToRad; + + Double_t phi_slice = (phi_max - phi_min)/fNPhiSegments; + + for(Int_t p=0; p <= fNPhiSegments; p++) + { + Double_t phi_in_bundle = phi_min + p*phi_slice; + //printf("phi %f in slice %d patch %d middle row %f\n",phi_in_bundle/ToRad,fSlice,fPatch,r_in_bundle); + + for(Int_t padrow = NRows[fPatch][0]; padrow <= NRows[fPatch][1]; padrow++) + //for(Int_t padrow = middle_row; padrow <= 173; padrow++) + //for(Int_t padrow = 0; padrow <= middle_row; padrow++) + { + + for(pix1=(AliL3Digits*)fRowContainer[padrow].first; pix1!=0; pix1=(AliL3Digits*)pix1->nextRowPixel) + { + + Float_t xyz[3]; + fTransform->Slice2Sector(fSlice,padrow,sector,row); + fTransform->Raw2Local(xyz,sector,row,pix1->fPad,pix1->fTime); + //fTransform->Raw2Global(xyz,sector,row,pix1->fPad,pix1->fTime); + + Double_t r_pix = sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]); + + Double_t phi_pix = fTransform->GetPhi(xyz); + + //Double_t tanPhi0 = (r_pix*sin(phi_in_bundle)-r_in_bundle*sin(phi_pix))/(r_pix*cos(phi_in_bundle)-r_in_bundle*cos(phi_pix)); + Double_t tanPhi0 = (r_in_bundle*sin(phi_pix)-r_pix*sin(phi_in_bundle))/(r_in_bundle*cos(phi_pix)-r_pix*cos(phi_in_bundle)); + + Double_t phi0 = atan(tanPhi0); + //if(padrow > middle_row) + //phi0 = -phi0; + //if(phi0 < 0.55 || phi0 > 0.88) continue; + + Double_t angle = phi_pix - phi0; + Double_t kappa = 2*sin(angle)/r_pix; + + //Double_t angle = phi_in_bundle - phi0; + //Double_t kappa = 2*sin(angle)/r_in_bundle; + + //if(kappa < -0.006 || kappa > 0.006) continue; + + Short_t signal = pix1->fCharge; + + hist->Fill(kappa,phi0,signal); + + + } + + } + } + +} + + +void AliL3HoughTransformer::Transform2Circle(TH2F *hist) +{ + //Transformation is done with respect to local coordinates in slice. + //Transform every pixel into whole phirange, using parametrisation: + //kappa = 2*sin(phi-phi0)/R + + printf("Transforming 1 pixel only\n"); + + AliL3Digits *pix1; + Int_t sector,row; + + Int_t ymin = hist->GetYaxis()->GetFirst(); + Int_t ymax = hist->GetYaxis()->GetLast(); + Int_t nbinsy = hist->GetNbinsY(); + + for(Int_t padrow = NRows[fPatch][0]; padrow <= NRows[fPatch][1]; padrow++) + { + + for(pix1=(AliL3Digits*)fRowContainer[padrow].first; pix1!=0; pix1=(AliL3Digits*)pix1->nextRowPixel) + { + + Short_t signal = pix1->fCharge; + Int_t index = pix1->fIndex; + + for(Int_t p=0; p <= nbinsy; p++) + hist->AddBinContent(fIndex[index][p],signal); + + } + + } + + +} + +/* + void AliL3HoughTransformer::Transform2Circle(TH2F *hist) + { + //Transformation is done with respect to local coordinates in slice. + //Transform every pixel into whole phirange, using parametrisation: + //kappa = 2*sin(phi-phi0)/R + + printf("Transforming 1 pixel only\n"); + + AliL3Digits *pix1; + Int_t sector,row; + + Int_t ymin = hist->GetYaxis()->GetFirst(); + Int_t ymax = hist->GetYaxis()->GetLast(); + + for(Int_t padrow = NRows[fPatch][0]; padrow <= NRows[fPatch][1]; padrow++) + { + + for(pix1=(AliL3Digits*)fRowContainer[padrow].first; pix1!=0; pix1=(AliL3Digits*)pix1->nextRowPixel) + { + + Float_t xyz[3]; + fTransform->Slice2Sector(fSlice,padrow,sector,row); + fTransform->Raw2Local(xyz,sector,row,pix1->fPad,pix1->fTime); + + Double_t r_pix = sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]); + + Double_t phi_pix = fTransform->GetPhi(xyz); + Short_t signal = pix1->fCharge; + + for(Int_t p=ymin+1; p<=ymax; p++) + { + Double_t phi0 = hist->GetYaxis()->GetBinCenter(p); + Double_t kappa = 2*sin(phi_pix-phi0)/r_pix; + hist->Fill(kappa,phi0,signal); + } + + } + + } + + + } +*/ + +void AliL3HoughTransformer::TransformLines2Circle(TH2F *hist,AliL3TrackArray *tracks) +{ + + for(Int_t i=0; iGetNTracks(); i++) + { + AliL3HoughTrack *track = (AliL3HoughTrack*)tracks->GetCheckedTrack(i); + if(!track) {printf("AliL3HoughTransformer::TransformLines2Circle : NO TRACK IN ARRAY\n"); continue;} + + Double_t xmin = fTransform->Row2X(track->GetFirstRow()); + Double_t xmax = fTransform->Row2X(track->GetLastRow()); + + Double_t a = -1./tan(track->GetPsiLine()); + Double_t b = track->GetDLine()/sin(track->GetPsiLine()); + + Double_t ymin = a*xmin + b; + Double_t ymax = a*xmax + b; + + Double_t middle_x = xmin + (xmax-xmin)/2; + Double_t middle_y = ymin + (ymax-ymin)/2; + + Double_t r_middle = sqrt(middle_x*middle_x + middle_y*middle_y); + Double_t phi = atan2(middle_y,middle_x); + Double_t phi0 = 2*phi - track->GetPsiLine(); + Double_t kappa = 2*sin(phi-phi0)/r_middle; + hist->Fill(kappa,phi0,track->GetWeight()); + + } + + +} + +void AliL3HoughTransformer::Transform2Line(TH2F *hist,Int_t ref_row,Int_t *rowrange,Double_t *phirange,TH2F *raw) +{ + //Transform every pixel into histogram, using parametrisation: + //D = x*cos(psi) + y*sin(psi) + + // printf("In Transform; rowrange %d %d ref_row %d phirange %f %f\n",rowrange[0],rowrange[1],ref_row,phirange[0],phirange[1]); + + AliL3Digits *pix1; + + Int_t xmin = hist->GetXaxis()->GetFirst(); + Int_t xmax = hist->GetXaxis()->GetLast(); + + Double_t x0 = fTransform->Row2X(ref_row); + Double_t y0 = 0; + + Int_t sector,row; + //for(Int_t padrow = NRows[fPatch][0]; padrow <= NRows[fPatch][1]; padrow++) + + Double_t phi_min = -10*ToRad; + Double_t phi_max = 10*ToRad; + Double_t delta_phi = (phi_max-phi_min)/fNPhiSegments; + + + Int_t phi_min_index = (Int_t)((phirange[0]*ToRad-phi_min)/delta_phi); + Int_t phi_max_index = (Int_t)((phirange[1]*ToRad-phi_min)/delta_phi); + + + Int_t index; + + for(Int_t phi=phi_min_index; phi <= phi_max_index; phi++) + { + for(Int_t padrow = rowrange[0]; padrow <= rowrange[1]; padrow++) + { + + index = phi*fNumOfPadRows + padrow; + //printf("Looping index %d\n",index); + if(index > fContainerBounds || index < 0) + { + printf("AliL3HoughTransformer::Transform2Line : index %d out of range \n",index); + return; + } + //for(pix1=(AliL3Digits*)fRowContainer[padrow].first; pix1!=0; pix1=(AliL3Digits*)pix1->nextRowPixel) + for(pix1=(AliL3Digits*)fPhiRowContainer[index].first; pix1!=0; pix1=(AliL3Digits*)pix1->nextPhiRowPixel) + { + //printf("Transforming pixel in index %d pad %d time %d padrow %d\n",index,pix1->fPad,pix1->fTime,padrow); + Float_t xyz[3]; + fTransform->Slice2Sector(fSlice,padrow,sector,row); + fTransform->Raw2Local(xyz,sector,row,pix1->fPad,pix1->fTime); + + if(raw) + raw->Fill(xyz[0],xyz[1],pix1->fCharge); + + xyz[0] = xyz[0]-x0; + xyz[1] = xyz[1]-y0; + + //printf("calculating..."); + for(Int_t d=xmin+1; d<=xmax; d++) + { + Double_t psi = hist->GetXaxis()->GetBinCenter(d); + Double_t D = xyz[0]*cos(psi) + xyz[1]*sin(psi); + + Short_t signal = pix1->fCharge; + hist->Fill(psi,D,signal); + //printf("Filling histo, psi %f D %f\n",psi,D); + } + //printf("done\n"); + } + //printf(" done\n"); + } + + } + +} + +void AliL3HoughTransformer::GetPixels(Char_t *rootfile,TH2F *hist) +{ + + TFile *file = new TFile(rootfile); + file->cd(); + + AliTPCParam *param=(AliTPCParam *)file->Get("75x40_100x60"); + TTree *t=(TTree*)file->Get("TreeD_75x40_100x60"); + + AliSimDigits da, *digarr=&da; + t->GetBranch("Segment")->SetAddress(&digarr); + Stat_t num_of_entries=t->GetEntries(); + + Int_t digit_counter=0; + Float_t xyz[3]; + Double_t eta,phi; + + Int_t nrows = NRows[fPatch][1] - NRows[fPatch][0] + 1; + printf("nrows %d slice %d patch %d\n",nrows,fSlice,fPatch); + + if(fRowContainer) + delete [] fRowContainer; + fRowContainer = new AliL3HoughContainer[fNumOfPadRows]; + memset(fRowContainer,0,fNumOfPadRows*sizeof(AliL3HoughContainer)); + + + fContainerBounds = (fNPhiSegments+1)*(fNumOfPadRows+1); + printf("Allocating %d bytes to container of size %d\n",fContainerBounds*sizeof(AliL3HoughContainer),fContainerBounds); + if(fPhiRowContainer) + delete [] fPhiRowContainer; + fPhiRowContainer = new AliL3HoughContainer[fContainerBounds]; + memset(fPhiRowContainer,0,fContainerBounds*sizeof(AliL3HoughContainer)); + + Double_t phi_min = -10*ToRad; + Double_t phi_max = 10*ToRad; + Double_t delta_phi = (phi_max-phi_min)/fNPhiSegments; + Int_t index; + digit_counter=0; + + for (Int_t i=0; iGetEvent(i); + Int_t sector; + Int_t row; + param->AdjustSectorRow(digarr->GetID(),sector,row); + Int_t slice,padrow; + fTransform->Sector2Slice(slice,padrow,sector,row); + if(slice != fSlice) continue; + if(padrow < NRows[fPatch][0]) continue; + if(padrow > NRows[fPatch][1]) break; + digarr->First(); + do { + Int_t time=digarr->CurrentRow(); + Int_t pad=digarr->CurrentColumn(); + Short_t signal=digarr->CurrentDigit(); + if(time < param->GetMaxTBin()-1 && time > 0) + if(digarr->GetDigit(time-1,pad) <= param->GetZeroSup() + && digarr->GetDigit(time+1,pad) <= param->GetZeroSup()) + continue; + + + fTransform->Raw2Global(xyz,sector,row,pad,time); + eta = fTransform->GetEta(xyz); + if(eta < fEtaMin || eta > fEtaMax) continue; + fTransform->Global2Local(xyz,sector); + + + phi = fTransform->GetPhi(xyz); + if(hist) + hist->Fill(xyz[0],xyz[1],signal); + + AliL3Digits *dig = new AliL3Digits; + dig->fIndex = digit_counter; + digit_counter++; + dig->fCharge = signal; + dig->fPad = pad; + dig->fTime = time; + + if(fRowContainer[padrow].first == NULL) + fRowContainer[padrow].first = (void*)dig; + else + ((AliL3Digits*)(fRowContainer[padrow].last))->nextRowPixel=dig; + fRowContainer[padrow].last = (void*)dig; + + Int_t phi_index = (Int_t)((phi-phi_min)/delta_phi); + index = phi_index*fNumOfPadRows + padrow; + if(phi_index > fContainerBounds || phi_index < 0) + { + printf("AliL3HoughTransform::GetPixels : index out of range %d\n",phi_index); + continue; + } + + if(fPhiRowContainer[index].first == NULL) + fPhiRowContainer[index].first = (void*)dig; + else + ((AliL3Digits*)(fPhiRowContainer[index].last))->nextPhiRowPixel = dig; + fPhiRowContainer[index].last=(void*)dig; + + }while (digarr->Next()); + + } + + fNDigits = digit_counter; + printf("digitcounter %d\n",digit_counter); + printf("Allocated %d bytes to pixels\n",digit_counter*sizeof(AliL3Digits)); + file->Close(); + delete file; + +} diff --git a/HLT/hough/AliL3HoughTransformer.h b/HLT/hough/AliL3HoughTransformer.h new file mode 100644 index 00000000000..2833e6a4023 --- /dev/null +++ b/HLT/hough/AliL3HoughTransformer.h @@ -0,0 +1,68 @@ +#ifndef ALIL3_HOUGHTRANSFORMER +#define ALIL3_HOUGHTRANSFORMER + +#include "AliL3RootTypes.h" + +struct AliL3Digits +{ + UShort_t fCharge; + UChar_t fPad; + UShort_t fTime; + Int_t fIndex; + AliL3Digits *nextPhiRowPixel; + AliL3Digits *nextRowPixel; +}; + +struct AliL3HoughContainer +{ + void *first; + void *last; +}; + +class TH2F; +class AliL3Transform; +class AliL3HoughEvaluate; +class AliL3TrackArray; + +class AliL3HoughTransformer : public TObject { + + private: + + friend class AliL3HoughEval; + + AliL3Transform *fTransform; //! + + Int_t fNPhiSegments; //Number of patches in phi. + Float_t fEtaMin; + Float_t fEtaMax; + Int_t fNumEtaSegments; + Int_t fNumOfPadRows; + + AliL3HoughContainer *fRowContainer; //! + AliL3HoughContainer *fPhiRowContainer; //! + Int_t fContainerBounds; + Int_t fNDigits; + Int_t **fIndex; //! + + + Int_t fSlice; + Int_t fPatch; + + public: + AliL3HoughTransformer(); + AliL3HoughTransformer(Int_t slice,Int_t patch,Float_t *etarange,Int_t phi_segments); + virtual ~AliL3HoughTransformer(); + + void Transform2Circle(TH2F *hist,Int_t middle_row); + void Transform2Circle(TH2F *hist); + void Transform2Line(TH2F *hist,Int_t ref_row,Int_t *rowrange,Double_t *phirange,TH2F *raw=0); + void TransformLines2Circle(TH2F *hist,AliL3TrackArray *tracks); + void GetPixels(Char_t *rootfile,TH2F *hist=0); + void InitTemplates(TH2F *hist); + void CountBins(); + + ClassDef(AliL3HoughTransformer,1) + +}; + +#endif diff --git a/HLT/hough/AliL3Transform.cxx b/HLT/hough/AliL3Transform.cxx new file mode 100644 index 00000000000..cdcd6ba8235 --- /dev/null +++ b/HLT/hough/AliL3Transform.cxx @@ -0,0 +1,574 @@ + +//Author: Anders Strand Vestbo +//Author: Uli Frankenfeld +//Last Modified: 05.12.2000 + +#include "AliL3Logging.h" +#include "AliL3Transform.h" +//#include +#include +//___________________________ +// AliL3Transform +// +// Transformation class for ALICE TPC. + +ClassImp(AliL3Transform); + + +AliL3Transform::AliL3Transform(){ + //constructor + Init(); +} + + +AliL3Transform::~AliL3Transform(){ +} + +void AliL3Transform::Init(){ + //sector: + fNTimeBins = 512; //?uli + fNRowLow = 55; + fNRowUp = 119; + fNSectorLow = 36; + fNSectorUp = 36; + fNSector = 72; + fPadPitchWidthLow = 0.400000; + fPadPitchWidthUp = 0.600000; + fZWidth = 0.56599998474121093750; + fZSigma = 0.22880849748219134199; + + //slices: + fNSlice = 36; + fNRow = 174; + fPi = 3.14159265358979323846; + for(Int_t i=0;i<36;i++){ + fCos[i] = cos(fPi/9*i); + fSin[i] = sin(fPi/9*i); + } + + fX[0] = 88.850006103515625; + fX[1] = 89.600006103515625; + fX[2] = 90.350006103515625; + fX[3] = 91.100006103515625; + fX[4] = 91.850006103515625; + fX[5] = 92.600006103515625; + fX[6] = 93.350006103515625; + fX[7] = 94.100006103515625; + fX[8] = 94.850006103515625; + fX[9] = 95.600006103515625; + fX[10] = 96.350006103515625; + fX[11] = 97.100006103515625; + fX[12] = 97.850006103515625; + fX[13] = 98.600006103515625; + fX[14] = 99.350006103515625; + fX[15] = 100.100006103515625; + fX[16] = 100.850006103515625; + fX[17] = 101.600006103515625; + fX[18] = 102.350006103515625; + fX[19] = 103.100006103515625; + fX[20] = 103.850006103515625; + fX[21] = 104.600006103515625; + fX[22] = 105.350006103515625; + fX[23] = 106.100006103515625; + fX[24] = 106.850006103515625; + fX[25] = 107.600006103515625; + fX[26] = 108.350006103515625; + fX[27] = 109.100006103515625; + fX[28] = 109.850006103515625; + fX[29] = 110.600006103515625; + fX[30] = 111.350006103515625; + fX[31] = 112.100006103515625; + fX[32] = 112.850006103515625; + fX[33] = 113.600006103515625; + fX[34] = 114.350006103515625; + fX[35] = 115.100006103515625; + fX[36] = 115.850006103515625; + fX[37] = 116.600006103515625; + fX[38] = 117.350006103515625; + fX[39] = 118.100006103515625; + fX[40] = 118.850006103515625; + fX[41] = 119.600006103515625; + fX[42] = 120.350006103515625; + fX[43] = 121.100006103515625; + fX[44] = 121.850006103515625; + fX[45] = 122.600006103515625; + fX[46] = 123.350006103515625; + fX[47] = 124.100006103515625; + fX[48] = 124.850006103515625; + fX[49] = 125.600006103515625; + fX[50] = 126.350006103515625; + fX[51] = 127.100006103515625; + fX[52] = 127.850006103515625; + fX[53] = 128.600006103515625; + fX[54] = 129.350006103515625; + fX[55] = 132.574996948242188; + fX[56] = 133.574996948242188; + fX[57] = 134.574996948242188; + fX[58] = 135.574996948242188; + fX[59] = 136.574996948242188; + fX[60] = 137.574996948242188; + fX[61] = 138.574996948242188; + fX[62] = 139.574996948242188; + fX[63] = 140.574996948242188; + fX[64] = 141.574996948242188; + fX[65] = 142.574996948242188; + fX[66] = 143.574996948242188; + fX[67] = 144.574996948242188; + fX[68] = 145.574996948242188; + fX[69] = 146.574996948242188; + fX[70] = 147.574996948242188; + fX[71] = 148.574996948242188; + fX[72] = 149.574996948242188; + fX[73] = 150.574996948242188; + fX[74] = 151.574996948242188; + fX[75] = 152.574996948242188; + fX[76] = 153.574996948242188; + fX[77] = 154.574996948242188; + fX[78] = 155.574996948242188; + fX[79] = 156.574996948242188; + fX[80] = 157.574996948242188; + fX[81] = 158.574996948242188; + fX[82] = 159.574996948242188; + fX[83] = 160.574996948242188; + fX[84] = 161.574996948242188; + fX[85] = 162.574996948242188; + fX[86] = 163.574996948242188; + fX[87] = 164.574996948242188; + fX[88] = 165.574996948242188; + fX[89] = 166.574996948242188; + fX[90] = 167.574996948242188; + fX[91] = 168.574996948242188; + fX[92] = 169.574996948242188; + fX[93] = 170.574996948242188; + fX[94] = 171.574996948242188; + fX[95] = 172.574996948242188; + fX[96] = 173.574996948242188; + fX[97] = 174.574996948242188; + fX[98] = 175.574996948242188; + fX[99] = 176.574996948242188; + fX[100] = 177.574996948242188; + fX[101] = 178.574996948242188; + fX[102] = 179.574996948242188; + fX[103] = 180.574996948242188; + fX[104] = 181.574996948242188; + fX[105] = 182.574996948242188; + fX[106] = 183.574996948242188; + fX[107] = 184.574996948242188; + fX[108] = 185.574996948242188; + fX[109] = 186.574996948242188; + fX[110] = 187.574996948242188; + fX[111] = 188.574996948242188; + fX[112] = 189.574996948242188; + fX[113] = 190.574996948242188; + fX[114] = 191.574996948242188; + fX[115] = 192.574996948242188; + fX[116] = 193.574996948242188; + fX[117] = 194.574996948242188; + fX[118] = 195.574996948242188; + fX[119] = 196.574996948242188; + fX[120] = 197.574996948242188; + fX[121] = 198.574996948242188; + fX[122] = 199.574996948242188; + fX[123] = 200.574996948242188; + fX[124] = 201.574996948242188; + fX[125] = 202.574996948242188; + fX[126] = 203.574996948242188; + fX[127] = 204.574996948242188; + fX[128] = 205.574996948242188; + fX[129] = 206.574996948242188; + fX[130] = 207.574996948242188; + fX[131] = 208.574996948242188; + fX[132] = 209.574996948242188; + fX[133] = 210.574996948242188; + fX[134] = 211.574996948242188; + fX[135] = 212.574996948242188; + fX[136] = 213.574996948242188; + fX[137] = 214.574996948242188; + fX[138] = 215.574996948242188; + fX[139] = 216.574996948242188; + fX[140] = 217.574996948242188; + fX[141] = 218.574996948242188; + fX[142] = 219.574996948242188; + fX[143] = 220.574996948242188; + fX[144] = 221.574996948242188; + fX[145] = 222.574996948242188; + fX[146] = 223.574996948242188; + fX[147] = 224.574996948242188; + fX[148] = 225.574996948242188; + fX[149] = 226.574996948242188; + fX[150] = 227.574996948242188; + fX[151] = 228.574996948242188; + fX[152] = 229.574996948242188; + fX[153] = 230.574996948242188; + fX[154] = 231.574996948242188; + fX[155] = 232.574996948242188; + fX[156] = 233.574996948242188; + fX[157] = 234.574996948242188; + fX[158] = 235.574996948242188; + fX[159] = 236.574996948242188; + fX[160] = 237.574996948242188; + fX[161] = 238.574996948242188; + fX[162] = 239.574996948242188; + fX[163] = 240.574996948242188; + fX[164] = 241.574996948242188; + fX[165] = 242.574996948242188; + fX[166] = 243.574996948242188; + fX[167] = 244.574996948242188; + fX[168] = 245.574996948242188; + fX[169] = 246.574996948242188; + fX[170] = 247.574996948242188; + fX[171] = 248.574996948242188; + fX[172] = 249.574996948242188; + fX[173] = 250.574996948242188; + fNPads[0] = 71; + fNPads[1] = 71; + fNPads[2] = 71; + fNPads[3] = 73; + fNPads[4] = 73; + fNPads[5] = 73; + fNPads[6] = 75; + fNPads[7] = 75; + fNPads[8] = 75; + fNPads[9] = 77; + fNPads[10] = 77; + fNPads[11] = 77; + fNPads[12] = 79; + fNPads[13] = 79; + fNPads[14] = 79; + fNPads[15] = 81; + fNPads[16] = 81; + fNPads[17] = 81; + fNPads[18] = 83; + fNPads[19] = 83; + fNPads[20] = 83; + fNPads[21] = 85; + fNPads[22] = 85; + fNPads[23] = 85; + fNPads[24] = 87; + fNPads[25] = 87; + fNPads[26] = 87; + fNPads[27] = 89; + fNPads[28] = 89; + fNPads[29] = 89; + fNPads[30] = 89; + fNPads[31] = 91; + fNPads[32] = 91; + fNPads[33] = 91; + fNPads[34] = 93; + fNPads[35] = 93; + fNPads[36] = 93; + fNPads[37] = 95; + fNPads[38] = 95; + fNPads[39] = 95; + fNPads[40] = 97; + fNPads[41] = 97; + fNPads[42] = 97; + fNPads[43] = 99; + fNPads[44] = 99; + fNPads[45] = 99; + fNPads[46] = 101; + fNPads[47] = 101; + fNPads[48] = 101; + fNPads[49] = 103; + fNPads[50] = 103; + fNPads[51] = 103; + fNPads[52] = 105; + fNPads[53] = 105; + fNPads[54] = 105; + fNPads[55] = 73; + fNPads[56] = 73; + fNPads[57] = 73; + fNPads[58] = 75; + fNPads[59] = 75; + fNPads[60] = 75; + fNPads[61] = 75; + fNPads[62] = 77; + fNPads[63] = 77; + fNPads[64] = 77; + fNPads[65] = 79; + fNPads[66] = 79; + fNPads[67] = 79; + fNPads[68] = 81; + fNPads[69] = 81; + fNPads[70] = 81; + fNPads[71] = 81; + fNPads[72] = 83; + fNPads[73] = 83; + fNPads[74] = 83; + fNPads[75] = 85; + fNPads[76] = 85; + fNPads[77] = 85; + fNPads[78] = 85; + fNPads[79] = 87; + fNPads[80] = 87; + fNPads[81] = 87; + fNPads[82] = 89; + fNPads[83] = 89; + fNPads[84] = 89; + fNPads[85] = 91; + fNPads[86] = 91; + fNPads[87] = 91; + fNPads[88] = 91; + fNPads[89] = 93; + fNPads[90] = 93; + fNPads[91] = 93; + fNPads[92] = 95; + fNPads[93] = 95; + fNPads[94] = 95; + fNPads[95] = 95; + fNPads[96] = 97; + fNPads[97] = 97; + fNPads[98] = 97; + fNPads[99] = 99; + fNPads[100] = 99; + fNPads[101] = 99; + fNPads[102] = 101; + fNPads[103] = 101; + fNPads[104] = 101; + fNPads[105] = 101; + fNPads[106] = 103; + fNPads[107] = 103; + fNPads[108] = 103; + fNPads[109] = 105; + fNPads[110] = 105; + fNPads[111] = 105; + fNPads[112] = 105; + fNPads[113] = 107; + fNPads[114] = 107; + fNPads[115] = 107; + fNPads[116] = 109; + fNPads[117] = 109; + fNPads[118] = 109; + fNPads[119] = 111; + fNPads[120] = 111; + fNPads[121] = 111; + fNPads[122] = 111; + fNPads[123] = 113; + fNPads[124] = 113; + fNPads[125] = 113; + fNPads[126] = 115; + fNPads[127] = 115; + fNPads[128] = 115; + fNPads[129] = 115; + fNPads[130] = 117; + fNPads[131] = 117; + fNPads[132] = 117; + fNPads[133] = 119; + fNPads[134] = 119; + fNPads[135] = 119; + fNPads[136] = 121; + fNPads[137] = 121; + fNPads[138] = 121; + fNPads[139] = 121; + fNPads[140] = 123; + fNPads[141] = 123; + fNPads[142] = 123; + fNPads[143] = 125; + fNPads[144] = 125; + fNPads[145] = 125; + fNPads[146] = 125; + fNPads[147] = 127; + fNPads[148] = 127; + fNPads[149] = 127; + fNPads[150] = 129; + fNPads[151] = 129; + fNPads[152] = 129; + fNPads[153] = 129; + fNPads[154] = 131; + fNPads[155] = 131; + fNPads[156] = 131; + fNPads[157] = 133; + fNPads[158] = 133; + fNPads[159] = 133; + fNPads[160] = 135; + fNPads[161] = 135; + fNPads[162] = 135; + fNPads[163] = 135; + fNPads[164] = 137; + fNPads[165] = 137; + fNPads[166] = 137; + fNPads[167] = 139; + fNPads[168] = 139; + fNPads[169] = 139; + fNPads[170] = 139; + fNPads[171] = 141; + fNPads[172] = 141; + fNPads[173] = 141; +} + + +Double_t AliL3Transform::GetEta(Float_t *xyz) +{ + Double_t r3 = sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]+xyz[2]*xyz[2]); + Double_t eta = 0.5 * log((r3+xyz[2])/(r3-xyz[2])); + return eta; +} + +Double_t AliL3Transform::GetPhi(Float_t *xyz) +{ + + Double_t phi = atan2(xyz[1],xyz[0]); + //if(phi<0) phi=phi+2*TMath::Pi(); + return phi; +} + + +Bool_t AliL3Transform::Slice2Sector(Int_t slice, Int_t slicerow, Int_t & sector, Int_t &row) const{ + if(slicerow<0&&slicerow>=fNRow) return kFALSE; + if(slice<0||slice>=fNSlice) return kFALSE; + + if(slicerow=fNSector) return kFALSE; + if(sector=fNSector||row<0) return kFALSE; + if(sector=fNRowLow) return kFALSE; + slice = sector; + slicerow = row; + } + else{ + if(row>=fNRowUp) return kFALSE; + slice = sector - fNSectorLow; + slicerow = row + fNRowLow; + } + return kTRUE; +} + +Double_t AliL3Transform::Row2X(Int_t slicerow){ + if(slicerow<0||slicerow>=fNRow) return 0; + return fX[slicerow]; +} + +void AliL3Transform::Local2Global(Float_t *xyz,Int_t slice) +{ + //Transformation to global coordinate system + Float_t x0 = xyz[0]; + Float_t y0 = xyz[1]; + Float_t cs,sn; + cs = fCos[slice]; + sn = fSin[slice]; + xyz[0]=x0*cs-y0*sn; + xyz[1]=x0*sn+y0*cs; + xyz[2]=xyz[2];//global z=local z +} + +void AliL3Transform::Local2GlobalAngle(Float_t *angle,Int_t slice){ + angle[0] = fmod(angle[0]+slice*fPi/9,2*fPi); +} + +void AliL3Transform::Global2LocalAngle(Float_t *angle,Int_t slice){ + angle[0] = angle[0]-slice*fPi/9; + if(angle[0]<0) angle[0]+=2*fPi; +} + +void AliL3Transform::Raw2Local(Float_t *xyz,Int_t sector,Int_t row,Float_t pad,Float_t time) +{ + //Transformation from rawdata to local coordinate system + + Int_t slice,slicerow; + Sector2Slice(slice, slicerow, sector, row); + + xyz[0]=Row2X(slicerow); + Int_t npads= fNPads[slicerow]; + if(sectorSetOptStat(0); + + Int_t xbin = 60; + Int_t ybin = 60; + Float_t xrange[2] = {-0.006 , 0.006}; //Pt 0.2-> + //Float_t yrange[2] = {-0.17 , 0.17}; //slice 2 0.55->0.88 + Float_t yrange[2] = {-0.26 , 0.26};// -15 - 15 + //Float_t yrange[2] = {0.55 , 0.88}; //slice 2 0.55->0.88 + TH2F *hist = new TH2F("hist","Parameter space",xbin,xrange[0],xrange[1],ybin,yrange[0],yrange[1]); + hist->GetXaxis()->SetTitle("#kappa [cm^{-1}]"); + hist->GetYaxis()->SetTitle("#Phi [rad]"); + SetTH1Options(hist); + + Int_t xr[2] = {0,250}; + Int_t yr[2] = {-125,125}; + TH1F *ntracks = new TH1F("ntracks","",100,0,200); + TH2F *raw = new TH2F("raw","",250,xr[0],xr[1],250,yr[0],yr[1]); + TH2F *road = new TH2F("road","",250,0,250,250,yr[0],yr[1]); + TH2F *fake = new TH2F("fake","",250,0,250,250,-125,125); + TH2F *peaks = new TH2F("peaks","",xbin,xrange[0],xrange[1],ybin,yrange[0],yrange[1]); + peaks->SetMarkerStyle(3); + peaks->SetMarkerColor(2); + + int slice = 2,n_phi_segments=30; + //float eta[2] = {0.3,0.4}; + float eta[2] = {0.03,0.04}; + + AliL3HoughTransformer *a = new AliL3HoughTransformer(slice,patch,eta,n_phi_segments); + a->GetPixels(rootfile,raw); + a->InitTemplates(hist); + + //a->Transform2Circle(hist,78); + TStopwatch sw; + a->Transform2Circle(hist); + sw.Stop(); + printf("Transformation done in %f ms\n",sw.CpuTime()*1000); + + AliL3HoughMaxFinder *b = new AliL3HoughMaxFinder("KappaPhi"); + //b->SetThreshold(10000); + AliL3TrackArray *tracks = b->FindMaxima(hist); + + AliL3HoughEval *c = new AliL3HoughEval(a); + printf("number of peaks before %d\n",tracks->GetNTracks()); + + /* + for(int i=0; iGetNTracks(); i++) + { + AliL3HoughTrack *track = (AliL3HoughTrack*)tracks->GetCheckedTrack(i); + if(!track) printf("no track\n"); + printf("Pt %f Phi %f kappa %f weigth %d charge %d\n",track->GetPt(),track->GetPhi0(),track->GetKappa(),track->GetNHits(),track->GetCharge()); + } + */ + //c->LookInsideRoad(a,tracks,road,fake,ntracks); + c->LookInsideRoad(tracks); + printf("Number of trackcandidates %d\n",tracks->GetNTracks()); + + AliL3Transform *transform = new AliL3Transform(); + Int_t NRows[5][2] = {{ 0, 45},{46,77},{78,109},{110,141},{142,173}}; + for(Int_t i=0; iGetNTracks(); i++) + { + AliL3HoughTrack *track = (AliL3HoughTrack*)tracks->GetCheckedTrack(i); + if(!track) {printf("NO TRACK\n"); continue;} + peaks->Fill(track->GetKappa(),track->GetPhi0(),1); + //for(int row=0; row<174; row++) + for(int row=NRows[patch][0]; row<=NRows[patch][1]; row++) + { + Float_t xyz[3]; + track->GetCrossingPoint(row,xyz); + //track->GetCrossingPoint(slice,row,xyz); + //printf("Track does not cross line at row %d\n",row); + //transform->Local2Global(xyz,slice); + road->Fill(xyz[0],xyz[1],1); + } + //float an[2] = {track->GetPhi0(),0}; + //transform->Local2GlobalAngle(an,slice); + printf("pt %f phi0 %f weight %d\n",track->GetPt(),track->GetPhi0(),track->GetWeight()); + } + + + TCanvas *c1 = new TCanvas("c1","",900,900); + SetCanvasOptions(c1); + c1->Divide(2,2); + c1->cd(1); + hist->Draw("lego1"); + // peaks->Draw("same"); + + // TCanvas *c2 = new TCanvas("c2","",2); + c1->cd(3); + raw->Draw(""); + + //TCanvas *c3 = new TCanvas("c3","",2); + //road->SetMarkerStyle(5); + //road->SetMarkerColor(3); + c1->cd(4); + road->SetMarkerStyle(9); + road->Draw(); + /* + TCanvas *c4 = new TCanvas("c4","",2); + fake->Draw("cont1"); + + TCanvas *c5 = new TCanvas("c5","",2); + ntracks->Draw(); + */ + delete tracks; + delete a; + delete b; + delete c; + + //---------------------------------------------------------------------------- + if(MC) + { + TFile *file = new TFile(rootfile); + file->cd(); + + 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; jUncheckedAt(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; + + printf("Generated particle %d, with pt %f phi %f\n",p->GetPdgCode(),ptg,phi_part); + } + file->Close(); + delete file; + } + +} diff --git a/HLT/hough/hough_line.C b/HLT/hough/hough_line.C new file mode 100644 index 00000000000..89e8f41a768 --- /dev/null +++ b/HLT/hough/hough_line.C @@ -0,0 +1,126 @@ +void hough_line(char *rootfile,int patch,Bool_t MC=false) +{ + gStyle->SetOptStat(0); + + + Double_t pi = TMath::Pi(); + Int_t xbin = 60; + Int_t ybin = 60; + Float_t xrange[2] = {-pi/2,pi/2}; + Float_t yrange[2] = {-40,40}; + + TH2F *hist = new TH2F("hist","Parameter space",xbin,xrange[0],xrange[1],ybin,yrange[0],yrange[1]); + hist->GetXaxis()->SetTitle("#Psi"); + hist->GetYaxis()->SetTitle("D"); + SetTH1Options(hist); + + Int_t xr[2] = {0,250}; + Int_t yr[2] = {-125,125}; + TH1F *ntracks = new TH1F("ntracks","",100,0,200); + TH2F *raw = new TH2F("raw","",250,xr[0],xr[1],250,yr[0],yr[1]); + TH2F *road = new TH2F("road","",250,0,250,250,yr[0],yr[1]); + TH2F *container = new TH2F("container","",250,0,250,250,yr[0],yr[1]); + TH2F *peaks = new TH2F("peaks","",xbin,xrange[0],xrange[1],ybin,yrange[0],yrange[1]); + peaks->SetMarkerStyle(3); + peaks->SetMarkerColor(2); + + int slice = 2,n_phi_segments=20; + float eta[2] = {0.3,0.4}; + //float eta[2] = {0.05,0.06}; + + AliL3HoughTransformer *a = new AliL3HoughTransformer(slice,patch,eta,n_phi_segments); + + const Int_t NRows[5][2] = {{ 0, 45},{46,77},{78,109},{110,141},{142,173}}; + const int ref_row[5]={22,60,93,125,157}; + int rowrange[2] = {6,12}; + double phirange[2] = {-10.,10.}; + a->GetPixels(rootfile,raw); + a->Transform2Line(hist,9,rowrange,phirange); + + AliL3HoughMaxFinder *b = new AliL3HoughMaxFinder("DPsi"); + AliL3TrackArray *tracks = b->FindMaxima(hist,rowrange,5); + + TH2F *cross = new TH2F("cross","",250,xr[0],xr[1],250,yr[0],yr[1]); + for(int i=0; i<2; i++) + //for(int i=0; iGetNTracks(); i++) + { + AliL3HoughTrack *track = (AliL3HoughTrack*)tracks->GetCheckedTrack(i); + if(!track) continue; + Double_t xy[2]; + for(int j=0; j<174; j++) + { + track->GetLineCrossingPoint(j,xy); + cross->Fill(xy[0],xy[1],1); + } + printf("Segment %f %f weight %d\n",track->GetPsiLine(),track->GetDLine(),track->GetWeight()); + } + cross->SetMarkerStyle(6); + cross->SetMarkerColor(2); + + + TCanvas *c1 = new TCanvas("c1","",900,900); + SetCanvasOptions(c1); + c1->Divide(2,2); + c1->cd(1); + hist->Draw("box"); + peaks->SetMarkerStyle(3); + peaks->SetMarkerColor(2); + //peaks->Draw("same"); + + c1->cd(2); + container->Draw(); + + c1->cd(3); + raw->Draw(""); + cross->Draw("same"); + return; + //TCanvas *c3 = new TCanvas("c3","",2); + //road->SetMarkerStyle(5); + //road->SetMarkerColor(3); + c1->cd(4); + road->SetMarkerStyle(9); + road->Draw(); + /* + TCanvas *c4 = new TCanvas("c4","",2); + fake->Draw("cont1"); + + TCanvas *c5 = new TCanvas("c5","",2); + ntracks->Draw(); + */ + delete a; + delete b; + delete c; + + //---------------------------------------------------------------------------- + if(MC) + { + TFile *file = new TFile(rootfile); + file->cd(); + + 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; jUncheckedAt(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; + + printf("Generated particle %d, with pt %f phi %f\n",p->GetPdgCode(),ptg,phi_part); + } + file->Close(); + delete file; + } + +} diff --git a/HLT/hough/hough_line_merge.C b/HLT/hough/hough_line_merge.C new file mode 100644 index 00000000000..f24f0694e77 --- /dev/null +++ b/HLT/hough/hough_line_merge.C @@ -0,0 +1,205 @@ +void hough_line_merge(char *rootfile,int patch,Bool_t MC=false) +{ + gStyle->SetOptStat(0); + + + Double_t pi = TMath::Pi(); + Int_t xbin = 60; + Int_t ybin = 60; + Float_t xrange[2] = {-pi/2,pi/2}; + Float_t yrange[2] = {-40,40}; + + TH2F *hist; + + Int_t xr[2] = {0,250}; + Int_t yr[2] = {-125,125}; + TH1F *ntracks = new TH1F("ntracks","",100,0,200); + TH2F *raw = new TH2F("raw","",250,xr[0],xr[1],250,yr[0],yr[1]); + TH2F *road = new TH2F("road","",250,0,250,250,yr[0],yr[1]); + TH2F *fake = new TH2F("fake","",250,0,250,250,-125,125); + TH2F *peaks = new TH2F("peaks","",xbin,xrange[0],xrange[1],ybin,yrange[0],yrange[1]); + peaks->SetMarkerStyle(3); + peaks->SetMarkerColor(2); + + int slice = 2,n_phi_segments=20; + float eta[2] = {0.3,0.4}; + //float eta[2] = {0.05,0.06}; + + AliL3HoughTransformer *a = new AliL3HoughTransformer(slice,patch,eta,n_phi_segments); + a->GetPixels(rootfile,raw); + + AliL3HoughMaxFinder *b = new AliL3HoughMaxFinder("DPsi"); + + const Int_t NRows[5][2] = {{ 0, 45},{46,77},{78,109},{110,141},{142,173}}; + //const int ref_row[4]={5,15,25,35}; + //int rowrange[4][2]={{0,10},{10,20},{20,30},{30,45}}; + int ref_row[7] = {3,9,15,21,27,33,39}; + int rowrange[7][2] = {{0,6},{6,12},{12,18},{18,24},{24,30},{30,36},{36,42}}; + + double phirange[7][2] = {{-10.,-5.},{-7.5,-2.5},{-5.,0.},{-2.5,2.5},{0.,5.},{2.5,7.5},{5.,10.}}; + int npatches = 7; + int count=0; + + //AliL3TrackArray **track_pieces = new AliL3TrackArray*[npatches]; + AliL3TrackArray *track_pieces = new AliL3TrackArray("AliL3HoughTrack"); + + double pran[2] = {-10,10}; + + // int pa=3; + for(int row_pat=0; row_pat < 7; row_pat++) + { + + for(int pat=0; pat<7; pat++) + { + + TH2F *histo = new TH2F("histo","Parameter space",xbin,xrange[0],xrange[1],ybin,yrange[0],yrange[1]); + + printf("transforming ref_row %d rowrange %d %d phirange %f %f\n",ref_row[row_pat],rowrange[row_pat][0],rowrange[row_pat][1],phirange[pat][0],phirange[pat][1]); + + a->Transform2Line(histo,ref_row[row_pat],rowrange[row_pat],phirange[pat]); + //a->Transform2Line(histo,ref_row[row_pat],rowrange[row_pat],pran); + + AliL3TrackArray *tmp = (AliL3TrackArray*)b->FindMaxima(histo,rowrange[row_pat],ref_row[row_pat]); + + track_pieces->AddTracks(tmp); + + //AliL3HoughTrack *tra = (AliL3HoughTrack*)b->GetTracks()->GetCheckedTrack(0); + //printf("psi %f rowrange %d %d\n",tra->GetPsiLine(),tra->GetFirstRow(),tra->GetLastRow()); + + delete histo; + delete tmp; + } + + } + + + printf("Transformation complete, genereted %d tracks\n",track_pieces->GetNTracks()); + + + //track_pieces->Sort(); + TH2F *cross = new TH2F("cross","",250,xr[0],xr[1],250,yr[0],yr[1]); + Double_t xy[2]; + //for(int i=0; iGetNTracks()==0) continue; + AliL3HoughTrack *track = (AliL3HoughTrack*)track_pieces->GetCheckedTrack(i); + if(!track) {printf("NO TRACK\n"); continue;} + //for(int j=rowrange[pa][0]; j<=rowrange[pa][1]; j++) + /* + for(int j=0; j<174; j++) + { + track->GetLineCrossingPoint(j,xy); + cross->Fill(xy[0],xy[1],1); + } + */ + //printf("weight %d psi %f D %f rowrange %d %d\n",track->GetWeight(),track->GetPsiLine(),track->GetDLine(),track->GetFirstRow(),track->GetLastRow()); + + } + cross->SetMarkerStyle(6); + cross->SetMarkerColor(2); + //cross->Draw(); + + + + AliL3HoughTransformer *a2 = new AliL3HoughTransformer(slice,pat,eta,n_phi_segments); + + xbin = 60 ; + ybin = 60; + xrange[0] = -0.006; + xrange[1] = 0.006; + yrange[0] = -0.26; + yrange[1] = 0.26; + hist = new TH2F("hist","Parameter space",xbin,xrange[0],xrange[1],ybin,yrange[0],yrange[1]); + + TH2F *cross_helix = new TH2F("cross_helix","",250,xr[0],xr[1],250,yr[0],yr[1]); + + a2->TransformLines2Circle(hist,(AliL3TrackArray*)track_pieces); + + AliL3HoughMaxFinder *b2 = new AliL3HoughMaxFinder("KappaPhi"); + AliL3TrackArray *final_tracks = b2->FindMaxima(hist); + + printf("Number of track candidates before checking %d\n",final_tracks->GetNTracks()); + AliL3HoughEval *eval = new AliL3HoughEval(slice); + eval->SetTransformer(a); + eval->LookInsideRoad(final_tracks); + + printf("Total number of tracks after checking %d\n",final_tracks->GetNTracks()); + for(int i=0; i<3; i++)//final_tracks->GetNTracks(); i++) + { + + AliL3HoughTrack *track = (AliL3HoughTrack*)final_tracks->GetCheckedTrack(i); + if(!track) continue; + + for(int r=0; r<174; r++) + { + Float_t xyz_helix[3]; + track->GetCrossingPoint(r,xyz_helix); + cross_helix->Fill(xyz_helix[0],xyz_helix[1],1); + } + + peaks->Fill(track->GetKappa(),track->GetPhi0(),1); + printf("Found track, pt %f phi0 %f\n",track->GetPt(),track->GetPhi0()); + } + cross_helix->SetMarkerStyle(6); + cross_helix->SetMarkerColor(2); + TCanvas *c1 = new TCanvas("c1","",900,900); + SetCanvasOptions(c1); + c1->Divide(2,2); + c1->cd(1); + hist->Draw("box"); + + peaks->SetMarkerStyle(3); + peaks->SetMarkerColor(2); + //peaks->Draw("same"); + + c1->cd(2); + cross_helix->Draw(); + + c1->cd(3); + raw->Draw(""); + cross_helix->DrawCopy("same"); + //cross->Draw("same"); + + c1->cd(4); + //hist->Draw("box"); + raw->DrawCopy(); + + delete track_pieces; + delete a; + delete b; + delete a2; + delete eval; + //---------------------------------------------------------------------------- + if(MC) + { + TFile *file = new TFile(rootfile); + file->cd(); + + 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; jUncheckedAt(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; + + printf("Generated particle %d, with pt %f phi %f\n",p->GetPdgCode(),ptg,phi_part); + } + file->Close(); + delete file; + } + +} diff --git a/HLT/hough/hough_merge.C b/HLT/hough/hough_merge.C new file mode 100644 index 00000000000..9eaffcb9c7a --- /dev/null +++ b/HLT/hough/hough_merge.C @@ -0,0 +1,148 @@ +void hough_merge(char *rootfile,Bool_t MC=false) +{ + gStyle->SetOptStat(0); + + Int_t xbin = 60; + Int_t ybin = 60; + Float_t xrange[2] = {-0.006 , 0.006}; //Pt 0.1->4GeV 0.006-0.006 + //Float_t yrange[2] = {-0.17 , 0.17}; //slice 2 0.55->0.88 + Float_t yrange[2] = {-0.26 , 0.26}; //slice 2 0.55->0.88 + //Float_t yrange[2] = {0.55 , 0.88}; //slice 2 0.55->0.88 + + Int_t xr[2] = {0,250}; + Int_t yr[2] = {-125,125}; + TH1F *ntracks = new TH1F("ntracks","",100,0,200); + TH2F *raw = new TH2F("raw","",250,xr[0],xr[1],250,yr[0],yr[1]); + TH2F *road = new TH2F("road","",250,0,250,250,yr[0],yr[1]); + TH2F *fake = new TH2F("fake","",300,0,300,250,-125,125); + TH2F *peaks = new TH2F("peaks","",xbin,xrange[0],xrange[1],ybin,yrange[0],yrange[1]); + peaks->SetMarkerStyle(3); + peaks->SetMarkerColor(2); + + int slice = 2,patch=0,n_phi_segments=40; + float eta[2] = {0.3,0.4}; + // float eta[2] = {0.04,0.05}; + + AliL3HoughTransformer *a; + AliL3HoughMaxFinder *b; + AliL3HoughEval *c = new AliL3HoughEval(slice); + AliL3HoughMerge *d = new AliL3HoughMerge(slice); + TH2F *hist; + for(int pat=0; pat<5; pat++) + { + a = new AliL3HoughTransformer(slice,pat,eta,n_phi_segments); + hist = new TH2F("hist","Parameter space",xbin,xrange[0],xrange[1],ybin,yrange[0],yrange[1]); + a->GetPixels(rootfile,raw); + a->Transform2Circle(hist,87); + + b = new AliL3HoughMaxFinder("KappaPsi"); + AliL3TrackArray *tracks = b->FindMaxima(hist); + c->SetTransformer(a); + c->LookInsideRoad(tracks); + d->FillTracks(tracks,pat); + delete a; + delete b; + } + return; + + + TH2F *merge_hist = new TH2F("merge_hist","Parameter space",xbin,xrange[0],xrange[1],ybin,yrange[0],yrange[1]); + + d->FillHisto(merge_hist); + + + b = new AliL3HoughMaxFinder(slice,pat); + b->FindMaxima(merge_hist); + + AliL3TrackArray *mtrs = b->GetTracks(); + + c->CompareMC(rootfile,mtrs,eta); + + AliL3Transform *transform = new AliL3Transform(); + for(int tr=0; trGetNTracks(); tr++) + { + AliL3HoughTrack *t = (AliL3HoughTrack*)mtrs->GetCheckedTrack(tr); + if(!t) {printf("NO TRACK11\n"); continue;} + if(t->GetMCid() < 0) {printf("Track %d was fake, weigth %d\n",tr,t->GetNHits()); continue;} + printf("Found pt %f weigth %d\n",t->GetPt(),t->GetNHits()); + //printf("charge %f\n",t->GetCharge()); + for(Int_t j=0; j<174; j++) + { + Float_t xyz[3]; + if(!t->GetCrossingPoint(j,xyz)) + { + printf("Track does not cross line\n"); + continue; + } + road->Fill(xyz[0],xyz[1],1); + + } + + + } + + + TCanvas *c1 = new TCanvas("c1","",800,800); + SetCanvasOptions(c1); + c1->Divide(2,2); + c1->cd(1); + merge_hist->Draw("box"); + + //TCanvas *c2 = new TCanvas("c2","",2); + c1->cd(3); + raw->Draw(); + + //TCanvas *c3 = new TCanvas("c3","",2); + c1->cd(4); + road->SetMarkerStyle(7); + road->Draw(); + + delete b; + delete d; + delete c; + //---------------------------------------------------------------------------- + if(MC) + { + TFile *file = new TFile(rootfile); + file->cd(); + + 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; jUncheckedAt(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; + int found = 0; + for(int m=0; mGetNTracks(); m++) + { + AliL3HoughTrack *tra = (AliL3HoughTrack*)mtrs->GetCheckedTrack(m); + if(tra->GetMCid() != j) continue; + found = 1; + double dPt = fabs(ptg-tra->GetPt()); + float phi0[2] = {tra->GetPhi0(),0}; + transform->Local2GlobalAngle(phi0,slice); + double dPhi0 = fabs(phi_part-phi0[0]); + printf("Difference Pt %f Phi0 %f\n",dPt,dPhi0); + + } + if(!found) printf("Track %d was not found\n",j); + printf("Generated particle %d, with pt %f phi %f\n",p->GetPdgCode(),ptg,phi_part); + } + file->Close(); + delete file; + } + +} diff --git a/HLT/hough/hough_mergehistos.C b/HLT/hough/hough_mergehistos.C new file mode 100644 index 00000000000..854b50b179f --- /dev/null +++ b/HLT/hough/hough_mergehistos.C @@ -0,0 +1,63 @@ +void hough_mergehistos(char *rootfile) +{ + gStyle->SetOptStat(0); + + Int_t xbin = 60; + Int_t ybin = 60; + Float_t xrange[2] = {-0.006 , 0.006}; //Pt 0.1->4GeV 0.006-0.006 + //Float_t yrange[2] = {-0.17 , 0.17}; //slice 2 0.55->0.88 + Float_t yrange[2] = {-0.26 , 0.26}; //slice 2 0.55->0.88 + //Float_t yrange[2] = {0.55 , 0.88}; //slice 2 0.55->0.88 + + Int_t xr[2] = {0,250}; + Int_t yr[2] = {-125,125}; + TH1F *ntracks = new TH1F("ntracks","",100,0,200); + TH2F *raw = new TH2F("raw","",250,xr[0],xr[1],250,yr[0],yr[1]); + TH2F *road = new TH2F("road","",250,0,250,250,yr[0],yr[1]); + TH2F *fake = new TH2F("fake","",300,0,300,250,-125,125); + TH2F *peaks = new TH2F("peaks","",xbin,xrange[0],xrange[1],ybin,yrange[0],yrange[1]); + peaks->SetMarkerStyle(3); + peaks->SetMarkerColor(2); + + int slice = 2,patch=0,n_phi_segments=40; + float eta[2] = {0.3,0.4}; + // float eta[2] = {0.04,0.05}; + + AliL3HoughTransformer *a; + AliL3HoughMaxFinder *b = new AliL3HoughMaxFinder("KappaPsi"); + + TH2F **hist = new TH2F*[5]; + + for(int pat=0; pat<5; pat++) + { + a = new AliL3HoughTransformer(slice,pat,eta,n_phi_segments); + hist[pat] = new TH2F("hist","Parameter space",xbin,xrange[0],xrange[1],ybin,yrange[0],yrange[1]); + a->GetPixels(rootfile,raw); + a->InitTemplates(hist[pat]); + a->Transform2Circle(hist[pat]); + + //AliL3TrackArray *tracks = b->FindMaxima(hist[pat]); + + delete a; + + } + + + TCanvas *c1 = new TCanvas("c1","",800,800); + c1->Divide(2,2); + c1->cd(1); + hist[0]->DrawCopy("box"); + c1->cd(2); + hist[4]->DrawCopy("box"); + c1->cd(3); + hist[2]->DrawCopy("box"); + + c1->cd(4); + hist[4]->Add(hist[0],hist[2]); + hist[4]->Draw("box"); + + delete b; + + + delete [] hist; +} diff --git a/HLT/hough/rootlogon.C b/HLT/hough/rootlogon.C new file mode 100644 index 00000000000..82a2c4d6f04 --- /dev/null +++ b/HLT/hough/rootlogon.C @@ -0,0 +1,21 @@ + +{ + printf("\nWELCOME to the magic world of Level3\n\n"); + + //load library neccesary for TPC alone + gSystem->Load("$(ROOTSYS)/lib/libPhysics"); + gSystem->Load("$(ROOTSYS)/lib/libEG"); + gSystem->Load("$(ALICE_ROOT)/lib/tgt_$(ALICE_TARGET)/libSTEER"); + gSystem->Load("$(ALICE_ROOT)/lib/tgt_$(ALICE_TARGET)/libCONTAINERS"); + gSystem->Load("$(ALICE_ROOT)/lib/tgt_$(ALICE_TARGET)/libTPC"); + gSystem->Load("/usr/local/franken/lib/MLUC/lib/linux-i386/libMLUC.so"); + + gSystem->Load("$(ALICE)/mylibs/libAliL3"); + gSystem->Load("$(ALICE)/mylibs/libAliL3Hough"); + + gROOT->LoadMacro("$(ALICE)/dev/XFunct.C"); + gStyle->SetStatBorderSize(1); + gStyle->SetTitleBorderSize(0); + printf("TPC libraries loaded\n"); + printf("Level3 libraries loaded\n"); +} -- 2.43.0