** Warning **
Issuing rollback() due to DESTROY without explicit disconnect() of DBD::mysql::db handle dbname=lxr at /lxr/lib/LXR/Common.pm line 1113.
/CMSSW_15_1_X_2025-06-20-2300/SimMuon/DTDigitizer/test/plotDriftTime.r
Warning, /SimMuon/DTDigitizer/test/plotDriftTime.r is written in an unsupported language. File is not indexed.
0001
0002 /**********************************************************************
0003 * Various tests of the parametrization funtions.
0004 * Execute with: .x plotDriftTime.r
0005 *
0006 * Author: G. Bevilacqua, N. Amapane
0007 *
0008 **********************************************************************/
0009
0010 #include <iostream>
0011 #include "TGraph.h"
0012 #include "TString.h"
0013 #include "TPostScript.h"
0014 #include "TCanvas.h"
0015
0016 // in cm!
0017 #define DT_Cell_HalfWidth 2.1
0018
0019 typedef struct {
0020 double v_drift, t_drift, delta_t, t_width_m, t_width_p ;
0021 } drift_time ;
0022
0023 typedef struct {
0024 double v_drift, x_drift, delta_x, x_width_m, x_width_p ;
0025 } drift_distance;
0026
0027
0028 // Create graphs for parametrized drift time as a function one of the parameters.
0029 TGraphAsymmErrors * doGraph(int funct, // which of the pars varies (0->3)
0030 Double_t pars[3],// The other pars
0031 int Nx, Double_t * xv, // x array
0032 bool old, // use the old parametrization
0033 short interpol // use interpolation
0034 ){
0035 drift_time * DT;
0036
0037 static const int maxsize = 200;
0038 Double_t tv[maxsize];
0039 Double_t stm[maxsize];
0040 Double_t stp[maxsize];
0041
0042 if (Nx>maxsize) return 0;
0043
0044 Double_t par[4]; // x,alpha,bWire,bNorm,interpol
0045
0046 // Assign fixed parameters
0047 int idum = 0;
0048 for (int i=0; i<4; i++) {
0049 if (i != funct) {
0050 par[i] = pars[idum];
0051 idum++;
0052 }
0053 }
0054
0055 for(int i=0; i<Nx; i++) {
0056 par[funct] = xv[i];
0057 if(old){
0058 tv[i] = oldParametrization(par[0], par[1], par[2], par[3]);
0059 } else {
0060 DT = driftTime( par[0], par[1], par[2], par[3], interpol);
0061 tv[i] = DT->t_drift;
0062 stm[i] = DT->t_width_m;
0063 stp[i] = DT->t_width_p;
0064 }
0065 }
0066
0067 TGraphAsymmErrors *grx;
0068 if(old){
0069 grx = new TGraphAsymmErrors(Nx,xv,tv,0,0,0,0);
0070 } else {
0071 grx = new TGraphAsymmErrors(Nx,xv,tv,0,0,stm,stp);
0072 }
0073 return grx;
0074 }
0075
0076 // Plot graphs of the drift time (created with doGraph)
0077 void plotGraph(TString &title,
0078 TGraphAsymmErrors *grx1,
0079 TGraphAsymmErrors *grx2,
0080 TGraphAsymmErrors *grx3)
0081 {
0082 grx1->GetYaxis()->SetTitle("time (ns)");
0083 grx1->SetTitle(title.Data());
0084 grx1->SetLineColor(1);
0085 grx1->SetMarkerColor(1);
0086 grx1->SetMarkerStyle(2);
0087 grx1->SetMarkerSize(.4);
0088 grx1->SetMinimum(0);
0089 grx1->SetMaximum(500);
0090 grx1->Draw("AP");
0091
0092 grx2->SetLineColor(4);
0093 grx2->SetMarkerColor(4);
0094 grx2->SetMarkerStyle(22);
0095 grx2->SetMarkerSize(.7);
0096 grx2->Draw("P");
0097
0098 grx3->SetLineColor(2);
0099 grx3->SetMarkerColor(2);
0100 grx3->SetMarkerStyle(21);
0101 grx3->SetMarkerSize(.4);
0102 grx3->Draw("P");
0103 }
0104
0105
0106 // Option 1: Plot drift time as a function of x
0107 void DtVSx(Double_t alpha, Double_t bWire, Double_t bNorm, Int_t ifl, Int_t flagCanvas, Int_t flagFile)
0108 {
0109 static const Int_t N_x = 84; //passo x = 0.5 mm
0110 Double_t xarray[N_x];
0111
0112 float max = DT_Cell_HalfWidth+0.1;
0113
0114 float step = max*2./N_x;
0115
0116 for(Int_t i=0; i<N_x; i++) {
0117 xarray[i] = i*step - max;
0118 }
0119
0120 Double_t par[3];
0121 par[0] = alpha;
0122 par[1] = bWire;
0123 par[2] = bNorm;
0124
0125 TGraphAsymmErrors *gOld = doGraph(0,par,N_x,xarray,true,0);
0126 TGraphAsymmErrors *gNew = doGraph(0,par,N_x,xarray,false,0);
0127 TGraphAsymmErrors *gNewI = doGraph(0,par,N_x,xarray,false,1);
0128
0129 TString title = "Drift time for #alpha=" + dToTString(alpha) +
0130 ", Bnorm =" + dToTString(bNorm) +
0131 ", Bwire =" + dToTString(bWire);
0132 TCanvas *cx=0;
0133 if(flagCanvas==1) {
0134 cx = new TCanvas(title,"Drift time as a function of x", 600,400);
0135 }
0136
0137 plotGraph(title, gOld,gNew,gNewI);
0138
0139 if(flagFile==1) {
0140 TString path ="Plots/PlotDtVSx/"+title+".ps";
0141 path.ReplaceAll(" ","_");
0142 path.ReplaceAll("#","");
0143 cx->Print(path.Data());
0144 }
0145 }
0146
0147 // Option 2: Plot drift time as a function of alpha
0148 void DtVSalpha(Double_t x, Double_t bWire, Double_t bNorm, Int_t ifl, Int_t flagCanvas, Int_t flagFile)
0149 {
0150
0151 static const Int_t N_alpha = 120;
0152 float max = 60.;
0153
0154 //---
0155 float step = max*2./N_alpha;
0156 Double_t xarray[N_alpha];
0157
0158 for(Int_t i=0; i<N_alpha; i++) {
0159 xarray[i] = i*step - max;
0160 }
0161
0162 Double_t par[3];
0163 par[0] = x;
0164 par[1] = bWire;
0165 par[2] = bNorm;
0166
0167 TGraphAsymmErrors *gNew = doGraph(1,par,N_alpha,xarray,false,0);
0168 TGraphAsymmErrors *gNewI = doGraph(1,par,N_alpha,xarray,false,1);
0169 TGraphAsymmErrors *gOld = doGraph(1,par,N_alpha,xarray,true,0);
0170
0171
0172 TString title = "Drift time for x=" + dToTString(x) +
0173 ", Bnorm =" + dToTString(bNorm) +
0174 ", Bwire =" + dToTString(bWire);
0175 TCanvas *ca=0;
0176 if(flagCanvas==1) {
0177 ca = new TCanvas(title,"Drift time as a function of alpha", 600,400);
0178 }
0179
0180 plotGraph(title, gOld,gNew,gNewI);
0181
0182 if(flagFile==1) {
0183 TString path ="Plots/PlotDtVSalpha/"+title+".ps";
0184 path.ReplaceAll(" ","_");
0185 path.ReplaceAll("#","");
0186 ca->Print(path.Data());
0187 }
0188 }
0189
0190 // Option 3: Plot drift time as a function of Bwire
0191 void DtVSbWire(Double_t x, Double_t alpha, Double_t bNorm, Int_t ifl, Int_t flagCanvas, Int_t flagFile)
0192 {
0193 static const Int_t Nx = 50;
0194 float max = 0.5;
0195
0196 //---
0197 float step = max*2./Nx;
0198 Double_t xarray[Nx];
0199
0200 for(Int_t i=0; i<Nx; i++) {
0201 xarray[i] = i*step - max;
0202 }
0203
0204 Double_t par[3];
0205 par[0] = x;
0206 par[1] = alpha;
0207 par[2] = bNorm;
0208
0209 TGraphAsymmErrors *gOld = doGraph(2,par,Nx,xarray,true,0);
0210 TGraphAsymmErrors *gNew = doGraph(2,par,Nx,xarray,false,0);
0211 TGraphAsymmErrors *gNewI = doGraph(2,par,Nx,xarray,false,1);
0212
0213 TString title = "Drift time for x=" + dToTString(x) +
0214 ", #alpha=" + dToTString(alpha) +
0215 ", Bnorm =" + dToTString(bNorm);
0216 TCanvas *cbWire=0;
0217 if (flagCanvas==1) {
0218 cbWire = new TCanvas(title,"Drift time as a function of bWire", 600,400);
0219 }
0220
0221 plotGraph(title, gOld,gNew,gNewI);
0222
0223 if(flagFile==1)
0224 {
0225 TString path ="Plots/PlotDtVSbWire/"+title+".ps";
0226 path.ReplaceAll(" ","_");
0227 path.ReplaceAll("#","");
0228 cbWire->Print(path.Data());
0229 }
0230
0231 }
0232
0233 // Option 4: Plot drift time as a function of Bnorm
0234 void DtVSbNorm(Double_t x, Double_t alpha, Double_t bWire, Int_t ifl, Int_t flagCanvas, Int_t flagFile)
0235 {
0236
0237 static const Int_t Nx = 85;
0238 float max = 0.85;
0239
0240 //---
0241 float step = max/Nx;
0242 Double_t xarray[Nx];
0243
0244 for(Int_t i=0; i<Nx; i++) {
0245 xarray[i] = i*step;
0246 }
0247
0248 Double_t par[3];
0249 par[0] = x;
0250 par[1] = alpha;
0251 par[2] = bWire;
0252
0253 TGraphAsymmErrors *gOld = doGraph(3,par,Nx,xarray,true,0);
0254 TGraphAsymmErrors *gNew = doGraph(3,par,Nx,xarray,false,0);
0255 TGraphAsymmErrors *gNewI = doGraph(3,par,Nx,xarray,false,1);
0256
0257 TString title = "Drift time for x=" + dToTString(x) +
0258 ", #alpha=" + dToTString(alpha) +
0259 ", Bwire =" + dToTString(bWire);
0260 TCanvas *cbNorm=0;
0261 if(flagCanvas==1) {
0262 cbNorm = new TCanvas(title,"Drift time as a function of bNorm", 600,400);
0263 }
0264
0265 plotGraph(title, gOld,gNew,gNewI);
0266
0267 if(flagFile==1)
0268 {
0269 TString path ="Plots/PlotDtVSbNorm/"+title+".ps";
0270 path.ReplaceAll(" ","_");
0271 path.ReplaceAll("#","");
0272 cbNorm->Print(path.Data());
0273 }
0274
0275 }
0276
0277
0278
0279 // Option 6: Simulate a time box
0280 void timeBox(double alpha, double bWire, double bNorm, bool old){
0281 short interpol = 1;
0282
0283 int smear = 1;
0284 if (!old) {
0285 cout << "Use smearing of times? [0/1]" <<endl;
0286 cin >> smear;
0287 }
0288
0289 int ntracks = 50000;
0290
0291 cout << "Simulation of " << ntracks << " tracks with " << (old?"old":"new")
0292 << " parametrization " << endl
0293 << " Smearing of times " << ((smear&&(!old))?"ON":"OFF") << endl;
0294
0295
0296 TH1D *hTBox = new TH1D("hTBox","Time box",275,-50,500);
0297 TRandom rnd;
0298 drift_time * DT;
0299 for (int i = 0; i<ntracks; i++) {
0300 float x = rnd.Uniform(-2.1,2.1);
0301 if (old) {
0302 hTBox->Fill(oldParametrization(x,alpha,bWire,bNorm));
0303 } else {
0304 float dt;
0305 if (smear) {
0306 dt = smearedTime(x,alpha,bWire,bNorm, interpol);
0307 } else {
0308 DT = driftTime(x,alpha,bWire,bNorm, interpol);
0309 dt = DT->t_drift;
0310 }
0311 hTBox->Fill(dt);
0312 }
0313 }
0314 TCanvas *cx=new TCanvas("c_hTBox","Time box");
0315 hTBox->Draw();
0316 }
0317
0318 // Option 5: Print value computed by the parametrization
0319 void printDt(double x, double alpha, double bWire, double bNorm, int ifl)
0320 {
0321 short interpol = 0;
0322
0323 drift_time * DT;
0324 DT = driftTime(x,alpha,bWire,bNorm, interpol);
0325
0326 cout<<endl<<"driftTime(x="<<x<<",alpha="<<alpha<<",bWire="<<bWire<<",bNorm="<<bNorm<<",ifl="<<ifl<<") = "<<DT->t_drift<<" ns"<<endl;
0327
0328 cout<<"sigma_r = "<<DT->t_width_p<<" ns"<<endl;
0329 cout<<"sigma_l = "<<DT->t_width_m<<" ns"<<endl;
0330
0331 }
0332
0333
0334 // Used by option 9,
0335 float plotNoSmearing(float alpha, float bWire, float bNorm, short interpol) {
0336 drift_time * DT;
0337 drift_distance * DX;
0338
0339 TMarker * m = new TMarker();
0340 m->SetMarkerColor(kBlue);
0341 m->SetMarkerStyle(2);
0342 for (int i =0; i < 210; i++){
0343 float x = i/100.;
0344 DT = driftTime(x, alpha, bWire, bNorm, interpol);
0345 float dt = DT->t_drift;
0346 DX = trackDistance(dt,alpha,bWire,bNorm, interpol);
0347 m->DrawMarker(x,DX->x_drift/10.-x);
0348 }
0349 }
0350
0351 // Used by option 9,
0352 void resPullPlots(float alpha=0, float bWire=0, float bNorm=0, short interpol=1) {
0353 gROOT->LoadMacro("macros.C");
0354
0355 TStyle * style = getStyle();
0356 style->SetOptStat("OURMEN");
0357 // style->SetOptStat("RME");
0358 style->SetOptFit(101);
0359 style->cd();
0360
0361 int form = 1;
0362
0363 // TFile * f = new TFile("MBHitAnalysis_NewDigis.root");
0364 // TFile * f = new TFile("MBHitAnalysis_NewDigis_onlyMu.root");
0365
0366 TCanvas * c1;
0367
0368 // goto resvspos2;
0369
0370 TCanvas *cx=new TCanvas("c_hTBox","Time box");
0371 hTBox->Draw();
0372
0373 c1 = newCanvas("c_hPos",2,1,form);
0374 c1->cd(1);
0375 hPos->Draw();
0376 c1->cd(2);
0377 hPosC->Draw();
0378
0379 c1 = newCanvas("c_hRes",2,1,form);
0380 c1->cd(1);
0381 drawGFit(&hRes, -0.2,0.2,-0.1,0.1);
0382 c1->cd(2);
0383 drawGFit(&hResC, -0.2,0.2,-0.1,0.1);
0384
0385
0386 c1 = newCanvas("c_hPull",2,1,form);
0387 c1->cd(1);
0388 drawGFit(&hPull, -5,5);
0389 c1->cd(2);
0390 drawGFit(&hPullC, -5,5);
0391
0392 resvspos:
0393
0394 c1 = newCanvas("c_hResVsPos");
0395 plotAndProfileX(&hResVsPos,-0.1,0.1,true);
0396 plotNoSmearing(alpha, bWire, bNorm, interpol);
0397
0398 c1 = newCanvas("c_hResVsPosC");
0399 plotAndProfileX(&hResVsPosC,-0.1,0.1,true);
0400 plotNoSmearing(alpha, bWire, bNorm, interpol);
0401
0402
0403 pullvspos:
0404 c1 = newCanvas("c_hPullVsPos");
0405 plotAndProfileX(&hPullVsPos,-5,8, true);
0406 c1 = newCanvas("c_hPullVsPosC");
0407 plotAndProfileX(&hPullVsPosC,-5,8, true);
0408
0409 end:
0410 }
0411
0412 // Option 9: simulate resolution and pulls
0413 void resPull(double alpha, double bWire, double bNorm){
0414
0415 int ntracks = 50000;
0416 bool smear=true;
0417 //smear = false;
0418 short interpol = 1;
0419 bool secondLevelCorr = true;
0420 secondLevelCorr = false;
0421
0422 TH1F * hTBox = new TH1F ("hTBox","Time box",275,-50,500);
0423 TH1F * hPos = new TH1F ("hPos", "RHit position", 100, 0,2.5);
0424 TH1F * hRes = new TH1F ("hRes", "RHit residual", 1290, -4.3,4.3);
0425 TH1F * hPull = new TH1F ("hPull", "RHit pull", 100, -5, 5);
0426 TH2F * hResVsPos = new TH2F("hResVsPos", "RHit residual vs position",
0427 100, 0,2.5, 1290, -4.3,4.3);
0428 TH2F * hPullVsPos = new TH2F("hPullVsPos", "RHit pull vs position",
0429 100, 0,2.5, 130, -5.,8.);
0430
0431 TH2F * hTPosVsPos = new TH2F("hTPosVsPos", "True pos vs rec pos",
0432 220, 0,2.2, 220, 0,2.2);
0433 TH2F * hratio = new TH2F("hratio", "RHit residual vs position",
0434 220, 0,2.2, 200, 0., 1.5);
0435
0436 // TH2F * hTPosVsPos1 = new TH2F("hTPosVsPos1", "True pos vs rec pos",
0437 // 220, 0,2.2, 220, 0,2.2);
0438 // TH2F * hResVsPos1 = new TH2F("hResVsPos1", "RHit residual vs position",
0439 // 100, 0,2.5, 1290, -4.3,4.3);
0440
0441
0442
0443 TH1F * hPosC = new TH1F ("hPosC", "RHit position", 100, 0,2.5);
0444 TH1F * hResC = new TH1F ("hResC", "RHit residual", 1290, -4.3,4.3);
0445 TH1F * hPullC = new TH1F ("hPullC", "RHit pull", 100, -5, 5);
0446 TH1F * hPullC_fixed = new TH1F ("hPullC_fixed", "RHit pull", 100, -5, 5);
0447 TH2F * hResVsPosC = new TH2F("hResVsPosC", "RHit residual vs position",
0448 100, 0,2.5, 860, -4.3,4.3);
0449 TH2F * hPullVsPosC = new TH2F("hPullVsPosC", "RHit pull vs position",
0450 100, 0,2.5, 130, -5.,8.);
0451
0452
0453 TRandom rnd;
0454 drift_time * DT;
0455 drift_distance * DX;
0456 for (int i = 0; i<ntracks; i++) {
0457 float x = rnd.Uniform(-2.1,2.1);
0458 float dt;
0459 if(smear) {
0460 dt = smearedTime(x,alpha,bWire,bNorm, interpol);
0461 } else {
0462 DT = driftTime(x, alpha, bWire, bNorm, interpol);
0463 dt = DT->t_drift;
0464 }
0465 DX = trackDistance(dt,alpha,bWire,bNorm, interpol);
0466 float rX = DX->x_drift/10.;
0467 float res = rX - fabs(x);
0468 float pull = res*20./(DX->x_width_m + DX->x_width_p);
0469 hTBox->Fill(dt);
0470 hPos->Fill(rX);
0471 hRes->Fill(res);
0472 hPull->Fill(pull);
0473 hResVsPos->Fill(fabs(x),res);
0474 hPullVsPos->Fill(fabs(x),pull);
0475
0476
0477 hTPosVsPos->Fill(fabs(x),rX);
0478 hratio->Fill(rX,fabs(x)/rX);
0479
0480 // TF1 * fun = new TF1("fun", "pol2", 0, 0.2);
0481 // fun->SetParameter(0,-3.80066e-01);
0482 // fun->SetParameter(1, 1.36233e+01);
0483 // fun->SetParameter(2,-3.49056e+01);
0484
0485 // float corX = rX;
0486 // if (corX<0.2) corX *= fun->Eval(rX);
0487
0488 // hTPosVsPos1->Fill(fabs(x),corX);
0489 // hResVsPos1->Fill(fabs(x),corX-fabs(x));
0490
0491
0492 // correct for difference mean - mean value
0493 float dX = 0.;
0494 dX = (DX->x_width_p - DX->x_width_m)*sqrt(2/TMath::Pi())/10.;
0495 float rXC = TMath::Max(0,rX-dX);
0496 float sigma = (DX->x_width_m + DX->x_width_p)/20.;
0497 if (secondLevelCorr && fabs(x)<0.4) { //FIXME!!!
0498 DT = driftTime(rXC, alpha, bWire, bNorm, interpol); //SIGN
0499 float dX1 = (DT->t_width_p - DT->t_width_m)*sqrt(2/TMath::Pi())*DT->v_drift/10.;
0500 rXC = TMath::Max(0,rX-dX1);
0501 // sigma = (DT->t_width_p * DT->t_width_m)*DT->v_drift/20.;
0502 }
0503
0504 res = rXC - fabs(x);
0505 pull = res/sigma;
0506 hPosC->Fill(rXC);
0507 hResC->Fill(res);
0508 hPullC->Fill(pull);
0509 hPullC_fixed->Fill(res/0.01423);
0510 hResVsPosC->Fill(fabs(x),res);
0511 hPullVsPosC->Fill(fabs(x),pull);
0512 }
0513
0514 resPullPlots(alpha, bWire, bNorm, interpol);
0515
0516 cout << endl
0517 << "number of tracks: " << ntracks << endl
0518 << "time smearing: " << (smear?"ON":"OFF") << endl
0519 << "second-level correction: " << (secondLevelCorr?"ON":"OFF") << endl
0520 << "interpolation: " << (interpol?"ON":"OFF") << endl
0521 << endl
0522 << "Plots with a name ending by C include peak-mean correction." << endl
0523 << "red markers: profile of scatter plots" << endl
0524 << "blue markers: behaviour without smearing" << endl
0525 << endl;
0526 }
0527
0528
0529 // Option 10: Verify the peak-mean correction for the asymmetric arrival time distribution
0530
0531 void plotAsymDist(double mean, double sigma1, double sigma2) {
0532 int ntracks = 500000;
0533
0534 TH1D *hParam = new TH1D("hParam","Double half-gaussian",500,0,500);
0535
0536 for (int i = 0; i<ntracks; i++) {
0537 double x = asymGausSample(mean, sigma1, sigma2);
0538 hParam->Fill(x);
0539 }
0540
0541 TCanvas *c1 = new TCanvas("c_Param","Double half-gaussian sampling");
0542 hParam->Draw();
0543
0544 cout << endl << "The mean calculated from the parametrized t, sigma_L, sigma_R is:" << endl
0545 << "mean+(sigma_R-sigma_L)*sqrt(2./pi) = "
0546 << mean+(sigma2-sigma1)*sqrt(2./TMath::Pi()) <<endl;
0547
0548 }
0549
0550 // Option 11: Simulate the mean timer distribution
0551 void meanTimer(double alpha, double bWire, double bNorm, bool old){
0552
0553 int smear = 1;
0554 if (!old) {
0555 cout << "Use smearing of times? [0/1]" <<endl;
0556 cin >> smear;
0557 }
0558
0559 int ntracks = 50000;
0560
0561 cout << "Simulation of " << ntracks << " tracks with " << (old?"old":"new")
0562 << " parametrization " << endl
0563 << " Smearing of times " << ((smear&&(!old))?"ON":"OFF") << endl;
0564
0565
0566 TH1D *hMeanTimer = new TH1D("hMeanTimer","Mean timer",275,340,420);
0567 TRandom rnd;
0568
0569 for (int i = 0; i<ntracks; i++) {
0570 float x = rnd.Uniform(0,2.1);
0571 float t1=DriftTimeChoosing(x,alpha,bWire,bNorm,old,smear);
0572 float t2=DriftTimeChoosing(2.1-x,alpha,bWire,bNorm,old,smear);
0573 float t3=DriftTimeChoosing(x,alpha,bWire,bNorm,old,smear);
0574
0575 hMeanTimer->Fill( (t1+t3)/2 + t2 );
0576 }
0577 TCanvas *cx=new TCanvas("c_hMeanTimer","Mean Timer");
0578 hMeanTimer->Draw();
0579 }
0580
0581 double DriftTimeChoosing(double x,double alpha, double bWire, double bNorm, bool old, bool smear){
0582 short interpol = 1;
0583 drift_time * DT;
0584 if (old)
0585 return oldParametrization(x,alpha,bWire,bNorm);
0586 else {
0587 float dt;
0588 if (smear)
0589 dt = smearedTime(x,alpha,bWire,bNorm, interpol);
0590 else {
0591 DT = driftTime(x,alpha,bWire,bNorm, interpol);
0592 dt = DT->t_drift;
0593 }
0594 return dt;
0595 }
0596 }
0597
0598
0599
0600 int plot()
0601 {
0602 cout<<endl<<"Drift time in a muon chamber"<<endl;
0603 cout<<endl<<"[x] = cm , -2.1 cm < x < 2.1 cm"<<endl;
0604 cout<<"[dt] = ns"<<endl;
0605 cout<<"[alpha] = degrees , -30 < alpha < 30"<<endl;
0606 cout<<"[bWire] = T, |bWire| < 0.4"<<endl;
0607 cout<<"[bNorm] = T, |bNorm| < 0.75 (symmetric)"<<endl;
0608
0609 Double_t dt, x, alpha, bWire, bNorm;
0610 Int_t ifl, flc, flf, fli;
0611 ifl = 0;
0612 flc = 1; //flagCanvas (for DtVS<...>): 0 = no creating canvas
0613 flf = 0; //flagFile (for DtVS<...>): 0 = no saving; 1 = save to file
0614 fli = 0; //flagIsto (for ShowGaussian): 0 = curve; 1 = curve + 500ev sampling
0615
0616 cout << endl
0617 << "Choose: " << endl
0618 << " 1 Plot drift time as a function of x"<<endl
0619 << " 2 Plot drift time as a function of alpha"<<endl
0620 << " 3 Plot drift time as a function of bWire"<<endl
0621 << " 4 Plot drift time as a function of bNorm"<<endl
0622 << " 5 Predefined collection of the above plots"<<endl
0623 << " 6 Compute drift time for the given parameters"<<endl
0624 << " 7 Simulated time box (new parametrization)"<<endl
0625 << " 8 Simulated time box (old parametrization)"<<endl
0626 << " 9 Simulated resolution and pulls"<<endl
0627 << " 10 Verify the peak-mean correction for the arrival time distribution"<<endl
0628 << " 11 Mean Timer distribution"<<endl;
0629
0630 Int_t n1 = 0;
0631 cin>>n1;
0632
0633 if(n1==1) {
0634 cout<<"Insert alpha value [deg]: ";
0635 cin>>alpha;
0636 cout<<"Insert bWire value [T]: ";
0637 cin>>bWire;
0638 cout<<"Insert bNorm value [T]: ";
0639 cin>>bNorm;
0640
0641 cout<<endl<<" alpha = "<<alpha<<endl;
0642 cout<<" bWire = "<<bWire<<endl;
0643 cout<<" bNorm = "<<bNorm<<endl;
0644 cout<<" ifl = "<<ifl<<endl;
0645
0646 DtVSx(alpha,bWire,bNorm,ifl,flc,flf);
0647
0648 helpDt();
0649
0650 }
0651
0652 else if(n1==2) {
0653 cout<<"Insert x value [cm from wire]: ";
0654 cin>>x;
0655 cout<<"Insert bWire value: [T]";
0656 cin>>bWire;
0657 cout<<"Insert bNorm value: [T]";
0658 cin>>bNorm;
0659
0660 cout<<endl<<" x = "<<x<<endl;
0661 cout<<" bWire = "<<bWire<<endl;
0662 cout<<" bNorm = "<<bNorm<<endl;
0663 cout<<" ifl = "<<ifl<<endl;
0664
0665 DtVSalpha(x,bWire,bNorm,ifl,flc,flf);
0666 helpDt();
0667
0668 }
0669
0670 else if(n1==3) {
0671 cout<<"Insert x value [cm from wire]: ";
0672 cin>>x;
0673 cout<<"Insert alpha value [deg]: ";
0674 cin>>alpha;
0675 cout<<"Insert bNorm value [T]: ";
0676 cin>>bNorm;
0677
0678 cout<<endl<<" x = "<<x<<endl;
0679 cout<<" alpha = "<<alpha<<endl;
0680 cout<<" bNorm = "<<bNorm<<endl;
0681 cout<<" ifl = "<<ifl<<endl;
0682
0683 DtVSbWire(x,alpha,bNorm,ifl,flc,flf);
0684 helpDt();
0685
0686 }
0687
0688 else if(n1==4) {
0689 cout<<"Insert x value [cm from wire]: ";
0690 cin>>x;
0691 cout<<"Insert alpha value [deg]: ";
0692 cin>>alpha;
0693 cout<<"Insert bWire value [T]: ";
0694 cin>>bWire;
0695
0696 cout<<endl<<" x = "<<x<<endl;
0697 cout<<" alpha = "<<alpha<<endl;
0698 cout<<" bWire = "<<bWire<<endl;
0699 cout<<" ifl = "<<ifl<<endl;
0700
0701 DtVSbNorm(x,alpha,bWire,ifl,flc,flf);
0702 helpDt();
0703
0704 }
0705
0706 else if(n1==5) {
0707 Double_t x, alpha, bWire, bNorm;
0708
0709
0710 cX = new TCanvas("DTvsX","Drift time as a function of x", 600,400);
0711 int i=1;
0712 cX->Divide(4,3);
0713 cX->cd(i++);
0714 alpha=0; bWire=0; bNorm=0;
0715 DtVSx(alpha,bWire,bNorm,0,0,0);
0716 cX->cd(i++);
0717 alpha=0; bWire=0.4; bNorm=0.7;
0718 DtVSx(alpha,bWire,bNorm,0,0,0);
0719 cX->cd(i++);
0720 alpha=0; bWire=0; bNorm=0;
0721 DtVSx(alpha,bWire,bNorm,0,0,0);
0722 cX->cd(i++);
0723 alpha=0; bWire=0.4; bNorm=0.7;
0724 DtVSx(alpha,bWire,bNorm,0,0,0);
0725
0726 cX->cd(i++);
0727 alpha=30; bWire=0; bNorm=0;
0728 DtVSx(alpha,bWire,bNorm,0,0,0);
0729 cX->cd(i++);
0730 alpha=30; bWire=0.4; bNorm=0.7;
0731 DtVSx(alpha,bWire,bNorm,0,0,0);
0732 cX->cd(i++);
0733 alpha=30; bWire=0; bNorm=0;
0734 DtVSx(alpha,bWire,bNorm,0,0,0);
0735 cX->cd(i++);
0736 alpha=30; bWire=0.4; bNorm=0.7;
0737 DtVSx(alpha,bWire,bNorm,0,0,0);
0738
0739 cX->cd(i++);
0740 alpha=55; bWire=0; bNorm=0;
0741 DtVSx(alpha,bWire,bNorm,0,0,0);
0742 cX->cd(i++);
0743 alpha=55; bWire=0.4; bNorm=0.7;
0744 DtVSx(alpha,bWire,bNorm,0,0,0);
0745 cX->cd(i++);
0746 alpha=55; bWire=0; bNorm=0;
0747 DtVSx(alpha,bWire,bNorm,0,0,0);
0748 cX->cd(i++);
0749 alpha=55; bWire=0.4; bNorm=0.7;
0750 DtVSx(alpha,bWire,bNorm,0,0,0);
0751
0752 cX = new TCanvas("DTvsalpha","Drift time as a function of alpha", 600,400);
0753 i=1;
0754 cX->Divide(4,3);
0755 cX->cd(i++);
0756 x=0.5; bWire=0; bNorm=0;
0757 DtVSalpha(x,bWire,bNorm,0,0,0);
0758 cX->cd(i++);
0759 x=0.5; bWire=0.4; bNorm=0;
0760 DtVSalpha(x,bWire,bNorm,0,0,0);
0761 cX->cd(i++);
0762 x=0.5; bWire=0; bNorm=0.75;
0763 DtVSalpha(x,bWire,bNorm,0,0,0);
0764 cX->cd(i++);
0765 x=0.5; bWire=0.4; bNorm=0.75;
0766 DtVSalpha(x,bWire,bNorm,0,0,0);
0767
0768 cX->cd(i++);
0769 x=1; bWire=0; bNorm=0;
0770 DtVSalpha(x,bWire,bNorm,0,0,0);
0771 cX->cd(i++);
0772 x=1; bWire=0.4; bNorm=0;
0773 DtVSalpha(x,bWire,bNorm,0,0,0);
0774 cX->cd(i++);
0775 x=1; bWire=0; bNorm=0.75;
0776 DtVSalpha(x,bWire,bNorm,0,0,0);
0777 cX->cd(i++);
0778 x=1; bWire=0.4; bNorm=0.75;
0779 DtVSalpha(x,bWire,bNorm,0,0,0);
0780
0781 cX->cd(i++);
0782 x=1.5; bWire=0; bNorm=0;
0783 DtVSalpha(x,bWire,bNorm,0,0,0);
0784 cX->cd(i++);
0785 x=1.5; bWire=0.4; bNorm=0;
0786 DtVSalpha(x,bWire,bNorm,0,0,0);
0787 cX->cd(i++);
0788 x=1.5; bWire=0; bNorm=0.75;
0789 DtVSalpha(x,bWire,bNorm,0,0,0);
0790 cX->cd(i++);
0791 x=1.5; bWire=0.4; bNorm=0.75;
0792 DtVSalpha(x,bWire,bNorm,0,0,0);
0793
0794
0795 cX = new TCanvas("DTvsbWire","Drift time as a function of bWire", 600,400);
0796 i=1;
0797 cX->Divide(4,3);
0798 cX->cd(i++);
0799 x=0.5; alpha=0; bNorm=0;
0800 DtVSbWire(x,alpha,bNorm,0,0,0);
0801 cX->cd(i++);
0802 x=0.5; alpha=30; bNorm=0;
0803 DtVSbWire(x,alpha,bNorm,0,0,0);
0804 cX->cd(i++);
0805 x=0.5; alpha=0; bNorm=0.4;
0806 DtVSbWire(x,alpha,bNorm,0,0,0);
0807 cX->cd(i++);
0808 x=0.5; alpha=30; bNorm=0.4;
0809 DtVSbWire(x,alpha,bNorm,0,0,0);
0810
0811 cX->cd(i++);
0812 x=1; alpha=0; bNorm=0;
0813 DtVSbWire(x,alpha,bNorm,0,0,0);
0814 cX->cd(i++);
0815 x=1; alpha=30; bNorm=0;
0816 DtVSbWire(x,alpha,bNorm,0,0,0);
0817 cX->cd(i++);
0818 x=1; alpha=0; bNorm=0.4;
0819 DtVSbWire(x,alpha,bNorm,0,0,0);
0820 cX->cd(i++);
0821 x=1; alpha=30; bNorm=0.4;
0822 DtVSbWire(x,alpha,bNorm,0,0,0);
0823
0824 cX->cd(i++);
0825 x=1.5; alpha=0; bNorm=0;
0826 DtVSbWire(x,alpha,bNorm,0,0,0);
0827 cX->cd(i++);
0828 x=1.5; alpha=30; bNorm=0;
0829 DtVSbWire(x,alpha,bNorm,0,0,0);
0830 cX->cd(i++);
0831 x=1.5; alpha=0; bNorm=0.4;
0832 DtVSbWire(x,alpha,bNorm,0,0,0);
0833 cX->cd(i++);
0834 x=1.5; alpha=30; bNorm=0.4;
0835 DtVSbWire(x,alpha,bNorm,0,0,0);
0836
0837
0838 cX = new TCanvas("DTvsbNorm","Drift time as a function of bNorm", 600,400);
0839 i=1;
0840 cX->Divide(4,3);
0841 cX->cd(i++);
0842 x=0.5; alpha=0; bWire=0;
0843 DtVSbNorm(x,alpha,bWire,0,0,0);
0844 cX->cd(i++);
0845 alpha=30; bWire=0;
0846 DtVSbNorm(x,alpha,bWire,0,0,0);
0847 cX->cd(i++);
0848 alpha=0; bWire=0.4;
0849 DtVSbNorm(x,alpha,bWire,0,0,0);
0850 cX->cd(i++);
0851 alpha=30; bWire=0.4;
0852 DtVSbNorm(x,alpha,bWire,0,0,0);
0853
0854 cX->cd(i++);
0855 x=1; alpha=0; bWire=0;
0856 DtVSbNorm(x,alpha,bWire,0,0,0);
0857 cX->cd(i++);
0858 alpha=30; bWire=0;
0859 DtVSbNorm(x,alpha,bWire,0,0,0);
0860 cX->cd(i++);
0861 alpha=0; bWire=0.4;
0862 DtVSbNorm(x,alpha,bWire,0,0,0);
0863 cX->cd(i++);
0864 alpha=30; bWire=0.4;
0865 DtVSbNorm(x,alpha,bWire,0,0,0);
0866
0867 cX->cd(i++);
0868 x=1.5; alpha=0; bWire=0;
0869 DtVSbNorm(x,alpha,bWire,0,0,0);
0870 cX->cd(i++);
0871 alpha=30; bWire=0;
0872 DtVSbNorm(x,alpha,bWire,0,0,0);
0873 cX->cd(i++);
0874 alpha=0; bWire=0.4;
0875 DtVSbNorm(x,alpha,bWire,0,0,0);
0876 cX->cd(i++);
0877 alpha=30; bWire=0.4;
0878 DtVSbNorm(x,alpha,bWire,0,0,0);
0879 }
0880
0881 else if(n1==6) {
0882 cout<<"Insert x value [cm from wire]: ";
0883 cin>>x;
0884 cout<<"Insert alpha value [deg]: ";
0885 cin>>alpha;
0886 cout<<"Insert bWire value [T]: ";
0887 cin>>bWire;
0888 cout<<"Insert bNorm value [T]: ";
0889 cin>>bNorm;
0890
0891 cout<<endl<<" x = "<<x<<endl;
0892 cout<<" alpha = "<<alpha<<endl;
0893 cout<<" bWire = "<<bWire<<endl;
0894 cout<<" bNorm = "<<bNorm<<endl;
0895 cout<<" ifl = "<<ifl<<endl;
0896
0897 // ... To be implemented: draw the reconstructed arrival time distribution.
0898
0899 printDt(x,alpha,bWire,bNorm,ifl);
0900
0901 }
0902
0903 else if(n1==7||n1==8) {
0904 cout<<"Insert alpha value [deg]: ";
0905 cin>>alpha;
0906 cout<<"Insert bWire value [T]: ";
0907 cin>>bWire;
0908 cout<<"Insert bNorm value [T]: ";
0909 cin>>bNorm;
0910
0911 cout<<endl<<" alpha = "<<alpha<<endl;
0912 cout<<" bWire = "<<bWire<<endl;
0913 cout<<" bNorm = "<<bNorm<<endl;
0914
0915 bool old = false;
0916 if (n1==8) old = true;
0917 timeBox(alpha,bWire,bNorm,old);
0918
0919 }
0920
0921 else if(n1==9) {
0922 cout<<"Insert alpha value [deg]: ";
0923 cin>>alpha;
0924 cout<<"Insert bWire value [T]: ";
0925 cin>>bWire;
0926 cout<<"Insert bNorm value [T]: ";
0927 cin>>bNorm;
0928
0929 cout<<endl<<" alpha = "<<alpha<<endl;
0930 cout<<" bWire = "<<bWire<<endl;
0931 cout<<" bNorm = "<<bNorm<<endl;
0932
0933 resPull(alpha, bWire, bNorm);
0934 }
0935
0936 else if(n1==10) {
0937
0938 double mean = 200;
0939 double sigma1 = 10;
0940 double sigma2 = 100;
0941
0942 plotAsymDist(mean, sigma1, sigma2);
0943 }
0944
0945 else if(n1==11) {
0946 int old= 0;
0947 cout<<"Insert alpha value [deg]: ";
0948 cin>>alpha;
0949 cout<<"Insert bWire value [T]: ";
0950 cin>>bWire;
0951 cout<<"Insert bNorm value [T]: ";
0952 cin>>bNorm;
0953 cout<<"Old[1] or New[0] parametrization? ";
0954 cin >> old;
0955 cout<<endl;
0956 cout<<" alpha = "<<alpha<<endl;
0957 cout<<" bWire = "<<bWire<<endl;
0958 cout<<" bNorm = "<<bNorm<<endl;
0959 cout<<" ifl = "<<ifl<<endl;
0960
0961 meanTimer(alpha,bWire,bNorm,old);
0962
0963 }
0964
0965
0966
0967
0968
0969
0970 else {
0971 cout<<"There's something wrong. n1= " << n1 << endl;
0972 n1 = 0;
0973 return 1;
0974 }
0975
0976
0977 cout<<endl<<"Execution successful"<<endl;
0978 return 0;
0979 }
0980
0981
0982 // Print a simple legend for plots...
0983 void helpDt() {
0984 cout << endl << "Legend:" <<endl
0985 << "Red: new parametrization (Puerta, Garcia-Abia), with interpolation" << endl
0986 << "Blue: new parametrization (Puerta, Garcia-Abia), no interpolation" << endl
0987 << "Black: old parametrization (Gresele, Rovelli)" <<endl;
0988 }
0989
0990
0991 // The entry point.
0992 // just call .x plotDriftTime.r
0993 void plotDriftTime() {
0994 // gROOT->GetList()->Delete();
0995 // gROOT->GetListOfCanvases()->Delete();
0996
0997 gSystem->Load("libCLHEP");
0998 gSystem->Load("libtestSimMuonDTDigitizer.so");
0999 plot();
1000 }
1001
1002
1003
1004
