initial commit

DrcPidFast.cxx

@@ -0,0 +1,75 @@
+#include "DrcPidFast.h"
+
+DrcPidFast::DrcPidFast(){
+
+  fMass[0]=0.000511;
+  fMass[1]=0.105658;
+  fMass[2]=0.139570;
+  fMass[3]=0.49368;
+  fMass[4]=0.938272;
+
+  // read Cherenkov track resolution map
+  TFile* file = TFile::Open("ctr_map.root");
+  fTrrMap = new TH2F();
+  file->GetObject("htrr", fTrrMap);  
+}
+
+DrcPidInfo DrcPidFast::GetInfo(int pdg,TVector3 mom, double err){
+
+  const int max = 5;
+  DrcPidInfo info;
+  int pid = get_pid(pdg);
+  double m = mom.Mag();
+  double theta = mom.Theta()*180/TMath::Pi();
+
+  // set default values
+  for(int i=0; i<max; i++){
+    info.probability[i]=0.25;
+    info.sigma[i]=100;
+  }
+
+  // check range
+  if(m>10) m = 10;
+  if(theta<25 || theta>140) return info;
+
+  int bin = fTrrMap->FindBin(theta,m);
+  double trr = fTrrMap->GetBinContent(bin);
+  double ctr = sqrt(trr*trr+err*err)*0.001;
+
+  // 1.46907 - fused silica
+  double true_cangle = acos(sqrt(m*m + fMass[pid]*fMass[pid])/m/1.46907);  
+  true_cangle += fRand.Gaus(0,ctr);
+
+  // return default values if momentum below Cherenkov threshold (true_cangle is NaN)
+  if(true_cangle != true_cangle) return info;
+
+  double cangle,sum=0,fsum=0;
+  double delta[max]={0}, probability[max]={0};
+
+  for(int i=0; i<max; i++){
+    cangle = acos(sqrt(m*m + fMass[i]*fMass[i])/m/1.46907);
+    if(cangle != cangle) continue;
+    delta[i] = fabs(cangle-true_cangle);
+    sum += delta[i];
+    info.sigma[i]=(cangle-true_cangle)/ctr;
+  }
+
+  // normalization
+  for(int i=0; i<max; i++){
+    if(delta[i]>0) info.probability[i] = sum/delta[i];
+    fsum += info.probability[i];
+  }
+  for(int i=0; i<max; i++) info.probability[i] /= fsum;
+
+  return info;
+}
+
+int DrcPidFast::get_pid(int pdg){
+  int pid=0;
+  if(pdg==11)   pid=0; //e
+  if(pdg==13)   pid=1; //mu
+  if(pdg==211)  pid=2; //pi
+  if(pdg==321)  pid=3; //K
+  if(pdg==2212) pid=4; //p
+  return pid;
+}

DrcPidFast.h

@@ -0,0 +1,39 @@
+// Fast reconstruction for EIC Barrel DIRC
+// original author: r.dzhygadlo at gsi.de
+
+#ifndef DrcPidFast_h
+#define DrcPidFast_h 1
+
+#include "TFile.h"
+#include "TH2F.h"
+#include "TVector3.h"
+#include "TRandom.h"
+#include <iostream>
+
+// prob - normalized to 1 probability for e,mu,pi,k,p
+// sigma - deviation of the determined Cherenkov angle from expected in terms of Cherenkov track resolution
+struct DrcPidInfo {
+    double probability[5];
+    double sigma[5];
+};
+
+class DrcPidFast{
+
+ public:
+  DrcPidFast();
+  ~DrcPidFast(){}
+
+  // pdg - Particle Data Group code of the particle
+  // mom - 3-momentum of the particle [GeV/c]
+  // err - error assosiated with track direction [mrad]
+  DrcPidInfo GetInfo(int pdg,TVector3 mom, double err=0);
+  TH2F *GetTrrMap(){ return fTrrMap; }
+
+ private:
+  int get_pid(int pdg);
+  TH2F *fTrrMap;
+  double fMass[5];
+  TRandom  fRand;
+};
+
+#endif

example.C

@@ -0,0 +1,25 @@
+#include "DrcPidFast.h"
+
+void example(int pdg=211,double mom=6){
+
+  DrcPidFast pid;
+
+  TVector3 momentum(0,0,mom);
+  momentum.RotateX(25/180.*TMath::Pi());
+
+  TH1F *hPi = new TH1F("hPi","",200,-10,10);
+  TH1F *hK = new TH1F("hK","",200,-10,10);
+
+  DrcPidInfo info;
+  for(int i=0; i<1000; i++){
+    info = pid.GetInfo(pdg,momentum,0);
+    hPi->Fill(info.sigma[2]);
+    hK->Fill(info.sigma[3]);
+  }
+
+  hPi->Draw();
+  hPi->SetLineColor(kBlue);
+  hK->SetLineColor(kRed);
+  hK->Draw("same");
+
+}

loadlib.C

@@ -0,0 +1,3 @@
+void loadlib(){
+  gROOT->ProcessLine(".L DrcPidFast.cxx+");
+}

plot_map.C

@@ -0,0 +1,81 @@
+#include "../prttools/prttools.C"
+#include "DrcPidFast.h"
+
+void plot_map(int pdg=211){
+
+  prt_setRootPalette(1);
+  TH2F *hsep = new TH2F("hsep",";polar angle [deg];momentum [GeV/c]",116,24.5,140.5,50,0.1,10.1);
+  TH1F *hPi = new TH1F("hPi","hPi",500,-100,100);
+  TH1F *hK = new TH1F("hK","hK",500,-100,100);
+  double m1=0,m2=0,s1=0,s2=0,sep=0;
+  DrcPidFast pid;    
+  DrcPidInfo info;
+  TF1 *ff;
+  //TCanvas *cc = new TCanvas("cc","cc",800,500);
+  for(double m=0.6; m<=10; m=m+0.2){
+    if(m>3){
+      delete hK;
+      delete hPi;
+      hK = new TH1F("hK","hK",200,-10,10);
+      hPi = new TH1F("hPi","hPi",200,-10,10);
+    }
+    if(m<1){
+      delete hK;
+      delete hPi;
+      hK = new TH1F("hK","hK",1000,-500,500);
+      hPi = new TH1F("hPi","hPi",1000,-500,500);
+    }
+    std::cout<<"m "<<m<<std::endl;
+
+    for(int t=25; t<=140; t++){
+      TVector3 mom(0,0,m);
+      mom.RotateX(t/180.*TMath::Pi());
+      for(int i=0; i<1000; i++){
+	info = pid.GetInfo(pdg,mom,0.5);
+	hPi->Fill(info.sigma[2]);
+	hK->Fill(info.sigma[3]);    
+      }
+
+      {
+	hK->Fit("gaus","Q0");
+	ff = hK->GetFunction("gaus");
+	if(ff) m1=ff->GetParameter(1);
+	if(ff) s1=ff->GetParameter(2);
+	hPi->Fit("gaus","Q0");
+	ff = hPi->GetFunction("gaus");
+	if(ff) m2=ff->GetParameter(1);
+	if(ff) s2=ff->GetParameter(2);
+
+	if(s1+s2>0) sep = (fabs(m2-m1))/(0.5*(s1+s2));
+      }
+
+      // hPi->Draw();
+      // hPi->SetLineColor(kBlue);
+      // hK->SetLineColor(kRed);
+      // hK->Draw("same");
+      // cc->Update();
+      // cc->WaitPrimitive();
+
+      hsep->Fill(t,m,sep);
+
+      hPi->Reset();
+      hK->Reset();
+    }
+  }
+
+  // prt_canvasAdd("diff",800,500);
+  // hPi->Draw();
+  // hPi->SetLineColor(kBlue);
+  // hK->SetLineColor(kRed);
+  // hK->Draw("same");
+
+  gStyle->SetOptStat(0);
+  prt_canvasAdd("map_sep",800,500);
+  hsep->GetYaxis()->SetRangeUser(0.6,10);
+  hsep->SetMaximum(10);
+  gPad->SetLogz();
+  hsep->Draw("colz");
+
+  prt_savepath ="data/plot_map";
+  prt_canvasSave();  
+}