HLT/hough/hough_merge.C

   1 void hough_merge(char *rootfile,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.1->4GeV 0.006-0.006
   8   //Float_t yrange[2] = {-0.17 , 0.17}; //slice 2 0.55->0.88
   9   Float_t yrange[2] = {-0.26 , 0.26}; //slice 2 0.55->0.88
  10   //Float_t yrange[2] = {0.55 , 0.88}; //slice 2 0.55->0.88
  11
  12   Int_t xr[2] = {0,250};
  13   Int_t yr[2] = {-125,125};
  14   TH1F *ntracks = new TH1F("ntracks","",100,0,200);
  15   TH2F *raw = new TH2F("raw","",250,xr[0],xr[1],250,yr[0],yr[1]);
  16   TH2F *road = new TH2F("road","",250,0,250,250,yr[0],yr[1]);
  17   TH2F *fake = new TH2F("fake","",300,0,300,250,-125,125);
  18   TH2F *peaks = new TH2F("peaks","",xbin,xrange[0],xrange[1],ybin,yrange[0],yrange[1]);
  19   peaks->SetMarkerStyle(3);
  20   peaks->SetMarkerColor(2);
  21
  22   int slice = 2,patch=0,n_phi_segments=40;
  23   float eta[2] = {0.3,0.4};
  24   // float eta[2] = {0.04,0.05};
  25
  26   AliL3HoughTransformer *a;
  27   AliL3HoughMaxFinder *b;
  28   AliL3HoughEval *c = new AliL3HoughEval(slice);
  29   AliL3HoughMerge *d = new AliL3HoughMerge(slice);
  30   TH2F *hist;
  31   for(int pat=0; pat<5; pat++)
  32     {
  33       a = new AliL3HoughTransformer(slice,pat,eta,n_phi_segments);
  34       hist = new TH2F("hist","Parameter space",xbin,xrange[0],xrange[1],ybin,yrange[0],yrange[1]);
  35       a->GetPixels(rootfile,raw);
  36       a->Transform2Circle(hist,87);
  37
  38       b = new AliL3HoughMaxFinder("KappaPsi");
  39       AliL3TrackArray *tracks = b->FindMaxima(hist);
  40       c->SetTransformer(a);
  41       c->LookInsideRoad(tracks);
  42       d->FillTracks(tracks,pat);
  43       delete a;
  44       delete b;
  45     }
  46   return;
  47
  48
  49   TH2F *merge_hist = new TH2F("merge_hist","Parameter space",xbin,xrange[0],xrange[1],ybin,yrange[0],yrange[1]);
  50
  51   d->FillHisto(merge_hist);
  52
  53
  54   b = new AliL3HoughMaxFinder(slice,pat);
  55   b->FindMaxima(merge_hist);
  56
  57   AliL3TrackArray *mtrs = b->GetTracks();
  58
  59   c->CompareMC(rootfile,mtrs,eta);
  60
  61   AliL3Transform *transform = new AliL3Transform();
  62   for(int tr=0; tr<mtrs->GetNTracks(); tr++)
  63     {
  64       AliL3HoughTrack *t = (AliL3HoughTrack*)mtrs->GetCheckedTrack(tr);
  65       if(!t) {printf("NO TRACK11\n"); continue;}
  66       if(t->GetMCid() < 0) {printf("Track %d was fake, weigth %d\n",tr,t->GetNHits()); continue;}
  67       printf("Found pt %f weigth %d\n",t->GetPt(),t->GetNHits());
  68       //printf("charge %f\n",t->GetCharge());
  69       for(Int_t j=0; j<174; j++)
  70         {
  71           Float_t xyz[3];
  72           if(!t->GetCrossingPoint(j,xyz))
  73             {
  74               printf("Track does not cross line\n");
  75               continue;
  76             }
  77           road->Fill(xyz[0],xyz[1],1);
  78
  79         }
  80
  81
  82     }
  83
  84
  85   TCanvas *c1 = new TCanvas("c1","",800,800);
  86   SetCanvasOptions(c1);
  87   c1->Divide(2,2);
  88   c1->cd(1);
  89   merge_hist->Draw("box");
  90
  91   //TCanvas *c2 = new TCanvas("c2","",2);
  92   c1->cd(3);
  93   raw->Draw();
  94
  95   //TCanvas *c3 = new TCanvas("c3","",2);
  96   c1->cd(4);
  97   road->SetMarkerStyle(7);
  98   road->Draw();
  99
 100   delete b;
 101   delete d;
 102   delete c;
 103   //----------------------------------------------------------------------------
 104   if(MC)
 105     {
 106       TFile *file = new TFile(rootfile);
 107       file->cd();
 108
 109       gAlice = (AliRun*)file->Get("gAlice");
 110       gAlice->GetEvent(0);
 111
 112       TClonesArray *particles=gAlice->Particles();
 113       Int_t n=particles->GetEntriesFast();
 114       Float_t torad=TMath::Pi()/180;
 115       Float_t phi_min = slice*20 - 10;
 116       Float_t phi_max = slice*20 + 10;
 117
 118       for (Int_t j=0; j<n; j++) {
 119         TParticle *p=(TParticle*)particles->UncheckedAt(j);
 120         if (p->GetFirstMother()>=0) continue;  //secondary particle
 121         if(p->Eta() < eta[0] || p->Eta() > eta[1]) continue;
 122         Double_t ptg=p->Pt(),pxg=p->Px(),pyg=p->Py(),pzg=p->Pz();
 123         Double_t phi_part = TMath::ATan2(pyg,pxg);
 124         if (phi_part < 0) phi_part += 2*TMath::Pi();
 125
 126         if(phi_part < phi_min*torad || phi_part > phi_max*torad) {continue;}
 127         if(ptg<0.100) continue;
 128         int found = 0;
 129         for(int m=0; m<mtrs->GetNTracks(); m++)
 130           {
 131             AliL3HoughTrack *tra = (AliL3HoughTrack*)mtrs->GetCheckedTrack(m);
 132             if(tra->GetMCid() != j) continue;
 133             found = 1;
 134             double dPt = fabs(ptg-tra->GetPt());
 135             float  phi0[2] = {tra->GetPhi0(),0};
 136             transform->Local2GlobalAngle(phi0,slice);
 137             double dPhi0 = fabs(phi_part-phi0[0]);
 138             printf("Difference Pt %f Phi0 %f\n",dPt,dPhi0);
 139
 140           }
 141         if(!found) printf("Track %d was not found\n",j);
 142         printf("Generated particle %d, with pt %f phi %f\n",p->GetPdgCode(),ptg,phi_part);
 143       }
 144       file->Close();
 145       delete file;
 146     }
 147
 148 }