HLT/hough/hough_line.C

   1 void hough_line(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 = new TH2F("hist","Parameter space",xbin,xrange[0],xrange[1],ybin,yrange[0],yrange[1]);
  13   hist->GetXaxis()->SetTitle("#Psi");
  14   hist->GetYaxis()->SetTitle("D");
  15   SetTH1Options(hist);
  16
  17   Int_t xr[2] = {0,250};
  18   Int_t yr[2] = {-125,125};
  19   TH1F *ntracks = new TH1F("ntracks","",100,0,200);
  20   TH2F *raw = new TH2F("raw","",250,xr[0],xr[1],250,yr[0],yr[1]);
  21   TH2F *road = new TH2F("road","",250,0,250,250,yr[0],yr[1]);
  22   TH2F *container = new TH2F("container","",250,0,250,250,yr[0],yr[1]);
  23   TH2F *peaks = new TH2F("peaks","",xbin,xrange[0],xrange[1],ybin,yrange[0],yrange[1]);
  24   peaks->SetMarkerStyle(3);
  25   peaks->SetMarkerColor(2);
  26
  27   int slice = 2,n_phi_segments=20;
  28   float eta[2] = {0.3,0.4};
  29   //float eta[2] = {0.05,0.06};
  30
  31   AliL3HoughTransformer *a = new AliL3HoughTransformer(slice,patch,eta,n_phi_segments);
  32
  33   const Int_t NRows[5][2] = {{ 0, 45},{46,77},{78,109},{110,141},{142,173}};
  34   const int ref_row[5]={22,60,93,125,157};
  35   int rowrange[2] = {6,12};
  36   double phirange[2] = {-10.,10.};
  37   a->GetPixels(rootfile,raw);
  38   a->Transform2Line(hist,9,rowrange,phirange);
  39
  40   AliL3HoughMaxFinder *b = new AliL3HoughMaxFinder("DPsi");
  41   AliL3TrackArray *tracks = b->FindMaxima(hist,rowrange,5);
  42
  43   TH2F *cross = new TH2F("cross","",250,xr[0],xr[1],250,yr[0],yr[1]);
  44   for(int i=0; i<2; i++)
  45   //for(int i=0; i<tracks->GetNTracks(); i++)
  46     {
  47       AliL3HoughTrack *track = (AliL3HoughTrack*)tracks->GetCheckedTrack(i);
  48       if(!track) continue;
  49       Double_t xy[2];
  50       for(int j=0; j<174; j++)
  51         {
  52           track->GetLineCrossingPoint(j,xy);
  53           cross->Fill(xy[0],xy[1],1);
  54         }
  55       printf("Segment %f %f weight %d\n",track->GetPsiLine(),track->GetDLine(),track->GetWeight());
  56     }
  57   cross->SetMarkerStyle(6);
  58   cross->SetMarkerColor(2);
  59
  60
  61   TCanvas *c1 = new TCanvas("c1","",900,900);
  62   SetCanvasOptions(c1);
  63   c1->Divide(2,2);
  64   c1->cd(1);
  65   hist->Draw("box");
  66   peaks->SetMarkerStyle(3);
  67   peaks->SetMarkerColor(2);
  68   //peaks->Draw("same");
  69
  70   c1->cd(2);
  71   container->Draw();
  72
  73   c1->cd(3);
  74   raw->Draw("");
  75   cross->Draw("same");
  76   return;
  77   //TCanvas *c3 = new TCanvas("c3","",2);
  78   //road->SetMarkerStyle(5);
  79   //road->SetMarkerColor(3);
  80   c1->cd(4);
  81   road->SetMarkerStyle(9);
  82   road->Draw();
  83   /*
  84   TCanvas *c4 = new TCanvas("c4","",2);
  85   fake->Draw("cont1");
  86
  87   TCanvas *c5 = new TCanvas("c5","",2);
  88   ntracks->Draw();
  89   */
  90   delete a;
  91   delete b;
  92   delete c;
  93
  94   //----------------------------------------------------------------------------
  95   if(MC)
  96     {
  97       TFile *file = new TFile(rootfile);
  98       file->cd();
  99
 100       gAlice = (AliRun*)file->Get("gAlice");
 101       gAlice->GetEvent(0);
 102
 103       TClonesArray *particles=gAlice->Particles();
 104       Int_t n=particles->GetEntriesFast();
 105       Float_t torad=TMath::Pi()/180;
 106       Float_t phi_min = slice*20 - 10;
 107       Float_t phi_max = slice*20 + 10;
 108
 109       for (Int_t j=0; j<n; j++) {
 110         TParticle *p=(TParticle*)particles->UncheckedAt(j);
 111         if (p->GetFirstMother()>=0) continue;  //secondary particle
 112         if(p->Eta() < eta[0] || p->Eta() > eta[1]) continue;
 113         Double_t ptg=p->Pt(),pxg=p->Px(),pyg=p->Py(),pzg=p->Pz();
 114         Double_t phi_part = TMath::ATan2(pyg,pxg);
 115         if (phi_part < 0) phi_part += 2*TMath::Pi();
 116
 117         if(phi_part < phi_min*torad || phi_part > phi_max*torad) {continue;}
 118         if(ptg<0.100) continue;
 119
 120         printf("Generated particle %d, with pt %f phi %f\n",p->GetPdgCode(),ptg,phi_part);
 121       }
 122       file->Close();
 123       delete file;
 124     }
 125
 126 }