| 1 | #include <iostream>
|
|---|
| 2 | #include <TMath.h>
|
|---|
| 3 | #include <TVectorD.h>
|
|---|
| 4 | #include <TVector3.h>
|
|---|
| 5 | #include <TH2.h>
|
|---|
| 6 | #include <TF1.h>
|
|---|
| 7 | #include <TGraph.h>
|
|---|
| 8 | #include <TLine.h>
|
|---|
| 9 | #include <TRandom.h>
|
|---|
| 10 | #include <TCanvas.h>
|
|---|
| 11 | #include <TStyle.h>
|
|---|
| 12 | #include <TLegend.h>
|
|---|
| 13 | #include "external/TrackCovariance/TrkUtil.h"
|
|---|
| 14 |
|
|---|
| 15 |
|
|---|
| 16 | void ExampleClCount(Int_t Nev = 1000, Double_t ang = 90., Int_t Opt=0)
|
|---|
| 17 | {
|
|---|
| 18 | //
|
|---|
| 19 | // Constants
|
|---|
| 20 | Double_t angRad = ang*TMath::Pi() / 180.; // To radians
|
|---|
| 21 | Double_t KpMass = 0.493677; // Charged Kaon mass
|
|---|
| 22 | Double_t PiMass = 0.13957018; // Charged Pion mass
|
|---|
| 23 | Double_t Bz = 2.0; // Field
|
|---|
| 24 | TVector3 x(0.0, 0.0, 0.0); // Track origin
|
|---|
| 25 | //
|
|---|
| 26 | Double_t pmin = 0.1; // Momentum range
|
|---|
| 27 | Double_t pmax = 50.;
|
|---|
| 28 | //
|
|---|
| 29 | // Setup chamber
|
|---|
| 30 | //
|
|---|
| 31 | TrkUtil *TU = new TrkUtil(Bz);
|
|---|
| 32 | //cout << "TU done" << endl;
|
|---|
| 33 | Double_t Rmin = 0.35;
|
|---|
| 34 | Double_t Rmax = 2.0;
|
|---|
| 35 | Double_t Zmin = -2.0;
|
|---|
| 36 | Double_t Zmax = 2.0;
|
|---|
| 37 | TU->SetDchBoundaries(Rmin, Rmax, Zmin, Zmax);
|
|---|
| 38 | TU->SetGasMix(Opt); // gas selection
|
|---|
| 39 | char title[50];
|
|---|
| 40 | Double_t y_min = 0.0;
|
|---|
| 41 | Double_t y_max = 0.0;
|
|---|
| 42 | switch (Opt)
|
|---|
| 43 | {
|
|---|
| 44 | case 0:
|
|---|
| 45 | {
|
|---|
| 46 | sprintf(title, "Helium 90 - Isobutane 10"); // He-Isobutane
|
|---|
| 47 | y_min = 1500.; y_max = 3500;
|
|---|
| 48 | break;
|
|---|
| 49 | }
|
|---|
| 50 | case 1:
|
|---|
| 51 | {
|
|---|
| 52 | sprintf(title, "Helium 100"); // pure He
|
|---|
| 53 | y_min = 400.; y_max = 1000;
|
|---|
| 54 | break;
|
|---|
| 55 | }
|
|---|
| 56 | case 2:
|
|---|
| 57 | {
|
|---|
| 58 | sprintf(title, "Argon 50 - Ethane 50"); // Argon - Ethane
|
|---|
| 59 | y_min = 6000.; y_max = 8500.;
|
|---|
| 60 | break;
|
|---|
| 61 | }
|
|---|
| 62 | case 3:
|
|---|
| 63 | {
|
|---|
| 64 | sprintf(title, "Argon 100"); // pure Argon
|
|---|
| 65 | y_min = 4000.; y_max = 6000.;
|
|---|
| 66 | break;
|
|---|
| 67 | }
|
|---|
| 68 | }
|
|---|
| 69 | //cout << "DchBoundaries set" << endl;
|
|---|
| 70 | //
|
|---|
| 71 | // Histograms
|
|---|
| 72 | //
|
|---|
| 73 | TH2F *h_kaon = new TH2F("h_kaon", "Nclusters distribution", 200, 0., 50., 50, y_min, y_max);
|
|---|
| 74 | TH2F *h_pion = new TH2F("h_pion", "Nclusters distribution", 200, 0., 50., 50, y_min, y_max);
|
|---|
| 75 | //cout << "histograms booked" << endl;
|
|---|
| 76 | //
|
|---|
| 77 | // Fill plots
|
|---|
| 78 | //
|
|---|
| 79 | for (Int_t n = 0; n < Nev; n++)
|
|---|
| 80 | {
|
|---|
| 81 | Double_t R = gRandom->Rndm();
|
|---|
| 82 | Double_t pval = pmin + R*(pmax - pmin);
|
|---|
| 83 | Double_t px = pval*TMath::Sin(angRad);
|
|---|
| 84 | Double_t pz = pval*TMath::Cos(angRad);
|
|---|
| 85 | TVector3 p(px, 0.0, pz);
|
|---|
| 86 | Double_t Q = 1.0;
|
|---|
| 87 | TVectorD Par = TrkUtil::XPtoPar(x, p, Q, Bz);
|
|---|
| 88 | //cout << "Parameters done. pval = "<< pval << endl;
|
|---|
| 89 | //
|
|---|
| 90 | Double_t Ncl = 0.0;
|
|---|
| 91 | if (TU->IonClusters(Ncl, KpMass, Par))
|
|---|
| 92 | {
|
|---|
| 93 | //cout << "Kaon extracted. pval = "<< pval << endl;
|
|---|
| 94 | h_kaon->Fill(pval, Ncl);
|
|---|
| 95 | //cout << "Kaon filled. Ncl = "<<Ncl << endl;
|
|---|
| 96 | }
|
|---|
| 97 | if (TU->IonClusters(Ncl, PiMass, Par))
|
|---|
| 98 | {
|
|---|
| 99 | //cout << "Pion extracted. pval = "<<pval << endl;
|
|---|
| 100 | h_pion->Fill(pval, Ncl);
|
|---|
| 101 | //cout << "Pion filled. Ncl = "<<Ncl << endl;
|
|---|
| 102 | }
|
|---|
| 103 | }
|
|---|
| 104 | //
|
|---|
| 105 | // graph K/pi separation
|
|---|
| 106 | //
|
|---|
| 107 | const Int_t Nint = 501;
|
|---|
| 108 | Double_t pmom[Nint];
|
|---|
| 109 | Double_t pmn = 0.2; Double_t pmx = 100.;
|
|---|
| 110 | Double_t stp = (TMath::Log(pmx) - TMath::Log(pmn)) / (Double_t)(Nint - 1);
|
|---|
| 111 | for (Int_t i = 0; i < Nint; i++)pmom[i] = TMath::Exp(i * stp + TMath::Log(pmn));
|
|---|
| 112 | Double_t SigDiff[Nint];
|
|---|
| 113 | for (Int_t i = 0; i < Nint; i++)
|
|---|
| 114 | {
|
|---|
| 115 | Double_t bgK = pmom[i] / KpMass;
|
|---|
| 116 | Double_t bgP = pmom[i] / PiMass;
|
|---|
| 117 | Double_t px = pmom[i]*TMath::Sin(angRad);
|
|---|
| 118 | Double_t pz = pmom[i]*TMath::Cos(angRad);
|
|---|
| 119 | TVector3 p(px, 0.0, pz);
|
|---|
| 120 | Double_t Q = 1.0;
|
|---|
| 121 | TVectorD Par = TrkUtil::XPtoPar(x, p, Q, Bz);
|
|---|
| 122 | Double_t tLen = TU->TrkLen(Par);
|
|---|
| 123 | Double_t CluKp = TrkUtil::Nclusters(bgK,Opt)*tLen;
|
|---|
| 124 | Double_t CluKpErr = TMath::Sqrt(CluKp);
|
|---|
| 125 | Double_t CluPi = TrkUtil::Nclusters(bgP,Opt)*tLen;
|
|---|
| 126 | Double_t CluPiErr = TMath::Sqrt(CluPi);
|
|---|
| 127 | //cout << "Momentum = " << pmom[i] << ", length= " << tLen
|
|---|
| 128 | // << ", CluKp = " << CluKp << ", CluPi = " << CluPi << endl;
|
|---|
| 129 | //if (CluKp + CluPi > 0.0)SigDiff[i] = TMath::Abs(CluPi - CluKp) / TMath::Sqrt(CluKp + CluPi);
|
|---|
| 130 | if (CluKp + CluPi > 0.0)SigDiff[i] = TMath::Abs(CluPi - CluKp) / (0.5*(CluKpErr + CluPiErr));
|
|---|
| 131 | }
|
|---|
| 132 | //
|
|---|
| 133 | // Plots
|
|---|
| 134 | //
|
|---|
| 135 | TCanvas *cnv = new TCanvas("cnv", title, 50, 50, 800, 500);
|
|---|
| 136 | cnv->Divide(2, 1);
|
|---|
| 137 | cnv->cd(1);
|
|---|
| 138 | gStyle->SetOptStat(0);
|
|---|
| 139 | h_kaon->SetMarkerColor(kRed);
|
|---|
| 140 | h_kaon->SetLineColor(kRed);
|
|---|
| 141 | h_kaon->GetXaxis()->SetTitle("Momentum (GeV)");
|
|---|
| 142 | h_kaon->Draw();
|
|---|
| 143 | h_pion->SetMarkerColor(kBlack);
|
|---|
| 144 | h_pion->SetLineColor(kBlack);
|
|---|
| 145 | h_pion->Draw("SAME");
|
|---|
| 146 | //
|
|---|
| 147 | TLegend* lg = new TLegend(0.1, 0.9, 0.3, 0.70);
|
|---|
| 148 | TString LgTitle = "Particle type:";
|
|---|
| 149 | lg->SetHeader(LgTitle);
|
|---|
| 150 | lg->AddEntry(h_pion, "#pi", "L");
|
|---|
| 151 | lg->AddEntry(h_kaon, "#color[2]{K}", "L");
|
|---|
| 152 | lg->Draw();
|
|---|
| 153 | cnv->cd(2);
|
|---|
| 154 | gPad->SetLogx();
|
|---|
| 155 | TGraph *gr = new TGraph(Nint, pmom, SigDiff);
|
|---|
| 156 | gr->SetMinimum(0.0);
|
|---|
| 157 | gr->SetMaximum(8.0);
|
|---|
| 158 | gr->GetXaxis()->SetLimits(0.,100.);
|
|---|
| 159 | gr->SetTitle("K/#pi separation in nr. of #sigma");
|
|---|
| 160 | gr->GetXaxis()->SetTitle("Momentum (GeV)");
|
|---|
| 161 | gr->SetLineColor(kBlue);
|
|---|
| 162 | gr->Draw("APC");
|
|---|
| 163 | TLine *line = new TLine(0.0, 3.0, 100., 3.0);
|
|---|
| 164 | line->Draw("SAME");
|
|---|
| 165 | //
|
|---|
| 166 | TCanvas* cnvf = new TCanvas("cnvf", "Interpolating function", 100, 100, 800, 500);
|
|---|
| 167 | cnvf->Divide(1, 1);
|
|---|
| 168 | cnvf->cd(1);
|
|---|
| 169 | gPad->SetLogx();
|
|---|
| 170 | TF1* f_ncl = new TF1("f_ncl", TU, &TrkUtil::funcNcl, 0.5, 1000., 1);
|
|---|
| 171 | f_ncl->SetNpx(500);
|
|---|
| 172 | f_ncl->Draw();
|
|---|
| 173 | }
|
|---|
