HLT/hough/hough.C

   1 void hough(char *rootfile,int patch,Bool_t MC=false)
   2 {
   3   gStyle->SetOptStat(0);
   4
   5   Int_t xbin = 60;
   6   Int_t ybin = 60;
   7   Float_t xrange[2] = {-0.006 , 0.006}; //Pt 0.2->
   8   //Float_t yrange[2] = {-0.17 , 0.17}; //slice 2 0.55->0.88
   9   Float_t yrange[2] = {-0.26 , 0.26};// -15 - 15
  10   //Float_t yrange[2] = {0.55 , 0.88}; //slice 2 0.55->0.88
  11   TH2F *hist = new TH2F("hist","Parameter space",xbin,xrange[0],xrange[1],ybin,yrange[0],yrange[1]);
  12   hist->GetXaxis()->SetTitle("#kappa [cm^{-1}]");
  13   hist->GetYaxis()->SetTitle("#Phi [rad]");
  14   SetTH1Options(hist);
  15
  16   Int_t xr[2] = {0,250};
  17   Int_t yr[2] = {-125,125};
  18   TH1F *ntracks = new TH1F("ntracks","",100,0,200);
  19   TH2F *raw = new TH2F("raw","",250,xr[0],xr[1],250,yr[0],yr[1]);
  20   TH2F *road = new TH2F("road","",250,0,250,250,yr[0],yr[1]);
  21   TH2F *fake = new TH2F("fake","",250,0,250,250,-125,125);
  22   TH2F *peaks = new TH2F("peaks","",xbin,xrange[0],xrange[1],ybin,yrange[0],yrange[1]);
  23   peaks->SetMarkerStyle(3);
  24   peaks->SetMarkerColor(2);
  25
  26   int slice = 2,n_phi_segments=30;
  27   float eta[2] = {0.3,0.4};
  28   //float eta[2] = {0.03,0.04};
  29
  30   AliL3HoughTransformer *a = new AliL3HoughTransformer(slice,patch,eta,n_phi_segments);
  31   a->GetPixels(rootfile,raw);
  32   a->InitTemplates(hist);
  33
  34   //a->Transform2Circle(hist,78);
  35   TStopwatch sw;
  36   a->Transform2Circle(hist);
  37   sw.Stop();
  38   printf("Transformation done in %f ms\n",sw.CpuTime()*1000);
  39
  40   AliL3HoughMaxFinder *b = new AliL3HoughMaxFinder("KappaPhi");
  41   //b->SetThreshold(10000);
  42   //AliL3TrackArray *tracks = b->FindMaxima(hist);
  43
  44   AliL3TrackArray *tracks = b->FindPeak(hist,2,0.95,2);
  45
  46   AliL3HoughEval *c = new AliL3HoughEval(a);
  47   printf("number of peaks before %d\n",tracks->GetNTracks());
  48
  49   /*
  50     for(int i=0; i<tracks->GetNTracks(); i++)
  51     {
  52       AliL3HoughTrack *track = (AliL3HoughTrack*)tracks->GetCheckedTrack(i);
  53       if(!track) printf("no track\n");
  54       printf("Pt %f Phi %f kappa %f weigth %d charge %d\n",track->GetPt(),track->GetPhi0(),track->GetKappa(),track->GetNHits(),track->GetCharge());
  55     }
  56   */
  57   //c->LookInsideRoad(a,tracks,road,fake,ntracks);
  58   c->LookInsideRoad(tracks);
  59   printf("Number of trackcandidates %d\n",tracks->GetNTracks());
  60
  61   AliL3Transform *transform = new AliL3Transform();
  62   Int_t NRows[5][2] = {{ 0, 45},{46,77},{78,109},{110,141},{142,173}};
  63   for(Int_t i=0; i<tracks->GetNTracks(); i++)
  64     {
  65       AliL3HoughTrack *track = (AliL3HoughTrack*)tracks->GetCheckedTrack(i);
  66       if(!track) {printf("NO TRACK\n"); continue;}
  67       peaks->Fill(track->GetKappa(),track->GetPhi0(),1);
  68       //for(int row=0; row<174; row++)
  69       for(int row=NRows[patch][0]; row<=NRows[patch][1]; row++)
  70         {
  71           Float_t xyz[3];
  72           track->GetCrossingPoint(row,xyz);
  73           //track->GetCrossingPoint(slice,row,xyz);
  74             //printf("Track does not cross line at row %d\n",row);
  75           //transform->Local2Global(xyz,slice);
  76           road->Fill(xyz[0],xyz[1],1);
  77         }
  78       //float an[2] = {track->GetPhi0(),0};
  79       //transform->Local2GlobalAngle(an,slice);
  80       printf("pt %f phi0 %f weight %d\n",track->GetPt(),track->GetPhi0(),track->GetWeight());
  81     }
  82
  83
  84   TCanvas *c1 = new TCanvas("c1","",900,900);
  85   SetCanvasOptions(c1);
  86   c1->Divide(2,2);
  87   c1->cd(1);
  88   hist->Draw("lego1");
  89   //  peaks->Draw("same");
  90
  91   //  TCanvas *c2 = new TCanvas("c2","",2);
  92   c1->cd(3);
  93   raw->Draw("");
  94
  95   //TCanvas *c3 = new TCanvas("c3","",2);
  96   //road->SetMarkerStyle(5);
  97   //road->SetMarkerColor(3);
  98   c1->cd(4);
  99   road->SetMarkerStyle(9);
 100   road->Draw();
 101   /*
 102   TCanvas *c4 = new TCanvas("c4","",2);
 103   fake->Draw("cont1");
 104
 105   TCanvas *c5 = new TCanvas("c5","",2);
 106   ntracks->Draw();
 107   */
 108   delete tracks;
 109   delete a;
 110   delete b;
 111   delete c;
 112
 113   //----------------------------------------------------------------------------
 114   if(MC)
 115     {
 116       TFile *file = new TFile(rootfile);
 117       file->cd();
 118
 119       gAlice = (AliRun*)file->Get("gAlice");
 120       gAlice->GetEvent(0);
 121
 122       TClonesArray *particles=gAlice->Particles();
 123       Int_t n=particles->GetEntriesFast();
 124       Float_t torad=TMath::Pi()/180;
 125       Float_t phi_min = slice*20 - 10;
 126       Float_t phi_max = slice*20 + 10;
 127
 128       for (Int_t j=0; j<n; j++) {
 129         TParticle *p=(TParticle*)particles->UncheckedAt(j);
 130         if (p->GetFirstMother()>=0) continue;  //secondary particle
 131         if(p->Eta() < eta[0] || p->Eta() > eta[1]) continue;
 132         Double_t ptg=p->Pt(),pxg=p->Px(),pyg=p->Py(),pzg=p->Pz();
 133         Double_t phi_part = TMath::ATan2(pyg,pxg);
 134         if (phi_part < 0) phi_part += 2*TMath::Pi();
 135
 136         if(phi_part < phi_min*torad || phi_part > phi_max*torad) {continue;}
 137         if(ptg<0.100) continue;
 138
 139         printf("Generated particle %d, with pt %f phi %f\n",p->GetPdgCode(),ptg,phi_part);
 140       }
 141       file->Close();
 142       delete file;
 143     }
 144
 145 }