HLT/hough/hough_line_merge.C

   1 void hough_line_merge(char *rootfile,int patch,Bool_t MC=false)
   2 {
   3   gStyle->SetOptStat(0);
   4
   5
   6   Double_t pi = TMath::Pi();
   7   Int_t xbin = 60;
   8   Int_t ybin = 60;
   9   Float_t xrange[2] = {-pi/2,pi/2};
  10   Float_t yrange[2] = {-40,40};
  11
  12   TH2F *hist;
  13
  14   Int_t xr[2] = {0,250};
  15   Int_t yr[2] = {-125,125};
  16   TH1F *ntracks = new TH1F("ntracks","",100,0,200);
  17   TH2F *raw = new TH2F("raw","",250,xr[0],xr[1],250,yr[0],yr[1]);
  18   TH2F *road = new TH2F("road","",250,0,250,250,yr[0],yr[1]);
  19   TH2F *fake = new TH2F("fake","",250,0,250,250,-125,125);
  20   TH2F *peaks = new TH2F("peaks","",xbin,xrange[0],xrange[1],ybin,yrange[0],yrange[1]);
  21   peaks->SetMarkerStyle(3);
  22   peaks->SetMarkerColor(2);
  23
  24   int slice = 2,n_phi_segments=20;
  25   float eta[2] = {0.3,0.4};
  26   //float eta[2] = {0.05,0.06};
  27
  28   AliL3HoughTransformer *a = new AliL3HoughTransformer(slice,patch,eta,n_phi_segments);
  29   a->GetPixels(rootfile,raw);
  30
  31   AliL3HoughMaxFinder *b = new AliL3HoughMaxFinder("DPsi");
  32
  33   const Int_t NRows[5][2] = {{ 0, 45},{46,77},{78,109},{110,141},{142,173}};
  34   //const int ref_row[4]={5,15,25,35};
  35   //int rowrange[4][2]={{0,10},{10,20},{20,30},{30,45}};
  36   int ref_row[7] = {3,9,15,21,27,33,39};
  37   int rowrange[7][2] = {{0,6},{6,12},{12,18},{18,24},{24,30},{30,36},{36,42}};
  38
  39   double phirange[7][2] = {{-10.,-5.},{-7.5,-2.5},{-5.,0.},{-2.5,2.5},{0.,5.},{2.5,7.5},{5.,10.}};
  40   int npatches = 7;
  41   int count=0;
  42
  43   //AliL3TrackArray **track_pieces = new AliL3TrackArray*[npatches];
  44   AliL3TrackArray *track_pieces = new AliL3TrackArray("AliL3HoughTrack");
  45
  46   double pran[2] = {-10,10};
  47
  48   // int pa=3;
  49   for(int row_pat=0; row_pat < 7; row_pat++)
  50     {
  51
  52       for(int pat=0; pat<7; pat++)
  53         {
  54
  55           TH2F *histo = new TH2F("histo","Parameter space",xbin,xrange[0],xrange[1],ybin,yrange[0],yrange[1]);
  56
  57           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]);
  58
  59           a->Transform2Line(histo,ref_row[row_pat],rowrange[row_pat],phirange[pat]);
  60           //a->Transform2Line(histo,ref_row[row_pat],rowrange[row_pat],pran);
  61
  62           AliL3TrackArray *tmp = (AliL3TrackArray*)b->FindMaxima(histo,rowrange[row_pat],ref_row[row_pat]);
  63
  64           track_pieces->AddTracks(tmp);
  65
  66           //AliL3HoughTrack *tra = (AliL3HoughTrack*)b->GetTracks()->GetCheckedTrack(0);
  67           //printf("psi %f rowrange %d %d\n",tra->GetPsiLine(),tra->GetFirstRow(),tra->GetLastRow());
  68
  69           delete histo;
  70           delete tmp;
  71         }
  72
  73     }
  74
  75
  76   printf("Transformation complete, genereted %d tracks\n",track_pieces->GetNTracks());
  77
  78
  79   //track_pieces->Sort();
  80   TH2F *cross = new TH2F("cross","",250,xr[0],xr[1],250,yr[0],yr[1]);
  81   Double_t xy[2];
  82   //for(int i=0; i<npatches; i++)
  83   for(int i=0; i<4; i++)
  84     {
  85       if(track_pieces->GetNTracks()==0) continue;
  86       AliL3HoughTrack *track = (AliL3HoughTrack*)track_pieces->GetCheckedTrack(i);
  87       if(!track) {printf("NO TRACK\n"); continue;}
  88       //for(int j=rowrange[pa][0]; j<=rowrange[pa][1]; j++)
  89       /*
  90       for(int j=0; j<174; j++)
  91         {
  92           track->GetLineCrossingPoint(j,xy);
  93           cross->Fill(xy[0],xy[1],1);
  94         }
  95       */
  96       //printf("weight %d psi %f D %f rowrange %d %d\n",track->GetWeight(),track->GetPsiLine(),track->GetDLine(),track->GetFirstRow(),track->GetLastRow());
  97
  98     }
  99   cross->SetMarkerStyle(6);
 100   cross->SetMarkerColor(2);
 101   //cross->Draw();
 102
 103
 104
 105   AliL3HoughTransformer *a2 = new AliL3HoughTransformer(slice,pat,eta,n_phi_segments);
 106
 107   xbin = 60 ;
 108   ybin = 60;
 109   xrange[0] = -0.006;
 110   xrange[1] = 0.006;
 111   yrange[0] = -0.26;
 112   yrange[1] = 0.26;
 113   hist = new TH2F("hist","Parameter space",xbin,xrange[0],xrange[1],ybin,yrange[0],yrange[1]);
 114
 115   TH2F *cross_helix = new TH2F("cross_helix","",250,xr[0],xr[1],250,yr[0],yr[1]);
 116
 117   a2->TransformLines2Circle(hist,(AliL3TrackArray*)track_pieces);
 118
 119   AliL3HoughMaxFinder *b2 = new AliL3HoughMaxFinder("KappaPhi");
 120   AliL3TrackArray *final_tracks = b2->FindMaxima(hist);
 121
 122   printf("Number of track candidates before checking %d\n",final_tracks->GetNTracks());
 123   AliL3HoughEval *eval = new AliL3HoughEval(slice);
 124   eval->SetTransformer(a);
 125   eval->LookInsideRoad(final_tracks);
 126
 127   printf("Total number of tracks after checking %d\n",final_tracks->GetNTracks());
 128   for(int i=0; i<3; i++)//final_tracks->GetNTracks(); i++)
 129     {
 130
 131       AliL3HoughTrack *track = (AliL3HoughTrack*)final_tracks->GetCheckedTrack(i);
 132       if(!track) continue;
 133
 134       for(int r=0; r<174; r++)
 135         {
 136           Float_t xyz_helix[3];
 137           track->GetCrossingPoint(r,xyz_helix);
 138           cross_helix->Fill(xyz_helix[0],xyz_helix[1],1);
 139         }
 140
 141       peaks->Fill(track->GetKappa(),track->GetPhi0(),1);
 142       printf("Found track, pt %f phi0 %f\n",track->GetPt(),track->GetPhi0());
 143     }
 144   cross_helix->SetMarkerStyle(6);
 145   cross_helix->SetMarkerColor(2);
 146   TCanvas *c1 = new TCanvas("c1","",900,900);
 147   SetCanvasOptions(c1);
 148   c1->Divide(2,2);
 149   c1->cd(1);
 150   hist->Draw("box");
 151
 152   peaks->SetMarkerStyle(3);
 153   peaks->SetMarkerColor(2);
 154   //peaks->Draw("same");
 155
 156   c1->cd(2);
 157   cross_helix->Draw();
 158
 159   c1->cd(3);
 160   raw->Draw("");
 161   cross_helix->DrawCopy("same");
 162   //cross->Draw("same");
 163
 164   c1->cd(4);
 165   //hist->Draw("box");
 166   raw->DrawCopy();
 167
 168   delete track_pieces;
 169   delete a;
 170   delete b;
 171   delete a2;
 172   delete eval;
 173   //----------------------------------------------------------------------------
 174   if(MC)
 175     {
 176       TFile *file = new TFile(rootfile);
 177       file->cd();
 178
 179       gAlice = (AliRun*)file->Get("gAlice");
 180       gAlice->GetEvent(0);
 181
 182       TClonesArray *particles=gAlice->Particles();
 183       Int_t n=particles->GetEntriesFast();
 184       Float_t torad=TMath::Pi()/180;
 185       Float_t phi_min = slice*20 - 10;
 186       Float_t phi_max = slice*20 + 10;
 187
 188       for (Int_t j=0; j<n; j++) {
 189         TParticle *p=(TParticle*)particles->UncheckedAt(j);
 190         if (p->GetFirstMother()>=0) continue;  //secondary particle
 191         if(p->Eta() < eta[0] || p->Eta() > eta[1]) continue;
 192         Double_t ptg=p->Pt(),pxg=p->Px(),pyg=p->Py(),pzg=p->Pz();
 193         Double_t phi_part = TMath::ATan2(pyg,pxg);
 194         if (phi_part < 0) phi_part += 2*TMath::Pi();
 195
 196         if(phi_part < phi_min*torad || phi_part > phi_max*torad) {continue;}
 197         if(ptg<0.100) continue;
 198
 199         printf("Generated particle %d, with pt %f phi %f\n",p->GetPdgCode(),ptg,phi_part);
 200       }
 201       file->Close();
 202       delete file;
 203     }
 204
 205 }