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 = 70;
   6   Int_t ybin = 70;
   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   ntracks = new TH1F("ntracks","",100,0,200);
  15   raw = new TH2F("raw","",250,xr[0],xr[1],250,yr[0],yr[1]);
  16   road = new TH2F("road","",250,0,250,250,yr[0],yr[1]);
  17   fake = new TH2F("fake","",300,0,300,250,-125,125);
  18   peaks = new TH2F("peaks","",xbin,xrange[0],xrange[1],ybin,yrange[0],yrange[1]);
  19   peaks->SetMarkerStyle(3);
  20   peaks->SetMarkerColor(2);
  21   road->SetMarkerStyle(6);
  22
  23   int slice = 2;
  24   //float eta[2] = {0.3,0.4};
  25   float eta[2] = {0.04,0.05};
  26
  27   b = new AliL3HoughMaxFinder("KappaPhi");
  28   //hist = new TH2F("hist","Parameter space",xbin,xrange[0],xrange[1],ybin,yrange[0],yrange[1]);
  29   TH2F **histos = new TH2F*[5];
  30   for(int pat=0; pat<5; pat++)
  31     {
  32       a = new AliL3HoughTransformer(slice,pat,eta);
  33       histos[pat] = new TH2F("hist","Parameter space",xbin,xrange[0],xrange[1],ybin,yrange[0],yrange[1]);
  34       a->GetPixels(rootfile,raw);
  35       a->InitTemplates(hist);
  36       a->Transform2Circle(hist);
  37       delete a;
  38     }
  39
  40   hist = new TH2F("hist","Parameter space",xbin,xrange[0],xrange[1],ybin,yrange[0],yrange[1]);
  41   for(Int_t i=0; i<5; i++)
  42     {
  43       hist->Add(histos[i],1);
  44     }
  45   tracks = (AliL3TrackArray*)b->FindPeak(hist,4,0.95,5);
  46   track = (AliL3HoughTrack*)tracks->GetCheckedTrack(0);
  47   peaks->Fill(track->GetKappa(),track->GetPhi0(),1);
  48
  49
  50   for(Int_t padrow=0; padrow<174; padrow++)
  51     {
  52       Float_t xyz[3];
  53       track->GetCrossingPoint(padrow,xyz);
  54       road->Fill(xyz[0],xyz[1],1);
  55     }
  56
  57   c1 = new TCanvas("c1","",1000,1000);
  58   c1->Divide(2,2);
  59   c1->cd(1);
  60   hist->Draw("box");
  61   peaks->Draw("same");
  62   c1->cd(2);
  63   raw->Draw();
  64   c1->cd(3);
  65   road->Draw();
  66   printf("Pt %f phi0 %f\n",track->GetPt(),track->GetPhi0());
  67   return;
  68
  69
  70   delete b;
  71   delete d;
  72   delete c;
  73   //----------------------------------------------------------------------------
  74   if(MC)
  75     {
  76       TFile *file = new TFile(rootfile);
  77       file->cd();
  78
  79       gAlice = (AliRun*)file->Get("gAlice");
  80       gAlice->GetEvent(0);
  81
  82       TClonesArray *particles=gAlice->Particles();
  83       Int_t n=particles->GetEntriesFast();
  84       Float_t torad=TMath::Pi()/180;
  85       Float_t phi_min = slice*20 - 10;
  86       Float_t phi_max = slice*20 + 10;
  87
  88       for (Int_t j=0; j<n; j++) {
  89         TParticle *p=(TParticle*)particles->UncheckedAt(j);
  90         if (p->GetFirstMother()>=0) continue;  //secondary particle
  91         if(p->Eta() < eta[0] || p->Eta() > eta[1]) continue;
  92         Double_t ptg=p->Pt(),pxg=p->Px(),pyg=p->Py(),pzg=p->Pz();
  93         Double_t phi_part = TMath::ATan2(pyg,pxg);
  94         if (phi_part < 0) phi_part += 2*TMath::Pi();
  95
  96         if(phi_part < phi_min*torad || phi_part > phi_max*torad) {continue;}
  97         if(ptg<0.100) continue;
  98         int found = 0;
  99         for(int m=0; m<mtrs->GetNTracks(); m++)
 100           {
 101             AliL3HoughTrack *tra = (AliL3HoughTrack*)mtrs->GetCheckedTrack(m);
 102             if(tra->GetMCid() != j) continue;
 103             found = 1;
 104             double dPt = fabs(ptg-tra->GetPt());
 105             float  phi0[2] = {tra->GetPhi0(),0};
 106             transform->Local2GlobalAngle(phi0,slice);
 107             double dPhi0 = fabs(phi_part-phi0[0]);
 108             printf("Difference Pt %f Phi0 %f\n",dPt,dPhi0);
 109
 110           }
 111         if(!found) printf("Track %d was not found\n",j);
 112         printf("Generated particle %d, with pt %f phi %f\n",p->GetPdgCode(),ptg,phi_part);
 113       }
 114       file->Close();
 115       delete file;
 116     }
 117
 118 }