HLT/hough/AliL3HoughEval.cxx

   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
  43 void AliL3HoughEval::LookInsideRoad(AliL3TrackArray *tracks,Bool_t remove,TH2F *hist)
  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;
  57   Int_t rows_to_miss = 1;
  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
  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
  94               if(remove) //Remove this pixel from image
  95                 pixel->fIndex = -1;
  96
  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
 146 TClonesArray *AliL3HoughEval::GetParticles(Char_t *rootfile)
 147 {
 148
 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();
 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);
 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_min*torad || phi_part > phi_max*torad) {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   }
 212
 213 }