diff --git a/Construct.dat b/Construct.dat new file mode 100644 index 0000000..aa23497 --- /dev/null +++ b/Construct.dat @@ -0,0 +1,4 @@ +1 3 3 1 + +hw, hh, columns, rows +hw and hh can be float values diff --git a/ConvoluteLoop.C b/ConvoluteLoop.C index 93f0ce6..f21b06e 100644 --- a/ConvoluteLoop.C +++ b/ConvoluteLoop.C @@ -1,7 +1,9 @@ #include #include #include +#include #include +#include #include #include "Math/Integrator.h" #include "RLAfuncs.cxx" @@ -27,7 +29,7 @@ double Func_1(double *t, double *par) // The detector-response function double Func_2(double *t, double *par) // The charge dispersion function { double result; - result = Qfuncs::PadRect(t[0], par[0], par[1], par[2], par[3]); + result = Qfuncs::PadRect(t[0], par[0], par[1], par[2], par[3], par[4]); //result = Qfuncs::TestTri(t[0]); return result; } @@ -47,10 +49,14 @@ Other than that, this is exactly the convolution code that Klaus has posted on the skipper physics link. */ -void ConvoluteLoop(double * maspar, double * relativepos) +void ConvoluteLoop(double * maspar, double * relativepos, char name[], char filename[], char title[], bool newFile) { //unpacking the master parameter array. These are set in the parent class. - double w = 0.25; + TStopwatch clock; + clock.Start(); + double R = maspar[5]; //75 micrometers in PowerPoint + if(newFile==true){TFile f(filename,"RECREATE");} + TFile f(filename,"UPDATE"); int n = int(maspar[0]+.1); double tR = maspar[1]; double sL = maspar[2]; @@ -66,14 +72,12 @@ void ConvoluteLoop(double * maspar, double * relativepos) fQ = new TF1("fQ", Func_2, tmin, tmax, 5);// 5 params: xl, xh, yl, yh, w // To be adjusted fRLA->SetParameters(tR, sL, tr, tf);//(double t_R, sigmaL, t_r, t_f)=(50.,13.,40.,2000.) all in ns - fQ->SetParameters(xl, xh, yl, yh, w);//(xl, xh, yl, yh)=(-1, 1, -3, 3) all in mm + fQ->SetParameters(xl, xh, yl, yh, R);//(xl, xh, yl, yh)=(-1, 1, -3, 3) all in mm fRLA->SetNpx(n); fQ->SetNpx(n); double tcplus = 0; double dt = (tmax-tmin)/double(n); - char name[40]; - sprintf(name,"(%g,%g) to (%g,%g))",xl,yl,xh,yh); - TH1D *cumConv = new TH1D("Signal",name, n, tmin, tmax); + TH1D *cumConv = new TH1D(name,title, n, tmin, tmax); cumConv->Reset(); TH1D *cumInt = new TH1D("cumInt", "", n, tmin, tmax); cumInt->Reset(); @@ -101,10 +105,15 @@ void ConvoluteLoop(double * maspar, double * relativepos) } } - cumConv->SetBinContent(iT, -maxInt);//cumInt->GetBinContent(n)); + cumConv->SetBinContent(iT, maxInt);//cumInt->GetBinContent(n)); } - cumConv->SetBinContent(n, 1.0E-12);//cumInt->GetBinContent(n)); cumConv->GetXaxis()->SetTitle("Time in ns"); - cumConv->Draw(); - cumConv->GetYaxis()->SetRangeUser(-3500,0); + cumConv->Write(); + delete cumConv; + delete cumInt; + cumConv = NULL; + cumInt = NULL; + f.Close(); + clock.Stop(); + cout << " Overall time: "; clock.Print(); cout << endl; } diff --git a/Input.dat b/Input.dat new file mode 100644 index 0000000..62fdd22 --- /dev/null +++ b/Input.dat @@ -0,0 +1,5 @@ +7 3 .5 1 11 4 1000 -200 4 + +X-start, Y-start, x-stepsize, y-stepsize, steps-in-x, steps-in-y, R-start, R-step, steps-in-R + +Number of steps include the 0 step! diff --git a/Qfuncs.cxx b/Qfuncs.cxx index 4fc56a9..ad85a24 100644 --- a/Qfuncs.cxx +++ b/Qfuncs.cxx @@ -17,15 +17,16 @@ #include "Qfuncs.h" #include "TMath.h" +#include "iostream" using namespace TMath; //______________________________________________________________________________ -double Qfuncs::PadRect(double tt, double xl, double xh, double yl, double yh, double w) +double Qfuncs::PadRect(double tt, double xl, double xh, double yl, double yh, double R) { // Calculate dispersion of charge across 2-D resistive anode - double N = 8000;// in q_e, typical COMPASS GEM value - double R = 2*530000;// Ohms per square + double N = 10000;// in q_e, 8000 is typical COMPASS GEM value + double w = .250;//530000;// Ohms per square double C = 0.21E-12;// 0.21pF/mm^2 //double w = 250E-3;// in mm /* @@ -35,14 +36,15 @@ double Qfuncs::PadRect(double tt, double xl, double xh, double yl, double yh, do double yl = par[2];//1;// in mm double yh = par[3];//7;// in mm */ - double Qpad;// = 1/Sqrt(4*t/(R*C)+2*Power(w,2)); + double ttn=-999; if(tt>0) { - double S = Sqrt(4*(tt*Power(double(10),-9))/(R*C)+2*Power(w,2)); + ttn=tt/double(1000000000); + double S = Sqrt(4*ttn/(R*C)+2*Power(w,2)); double diffxErf = Erf(xh/S)-Erf(xl/S); double diffyErf = Erf(yh/S)-Erf(yl/S); - Qpad = N/4 * ( diffxErf*diffyErf); + Qpad = N/double(4) * ( diffxErf*diffyErf ); } else { diff --git a/Qfuncs.h b/Qfuncs.h index 3fdb144..88849d3 100644 --- a/Qfuncs.h +++ b/Qfuncs.h @@ -33,7 +33,7 @@ */ namespace Qfuncs { - double PadRect(double tt=0, double xl=1, double xh=1, double yl=1, double yh=1, double w = 1); + double PadRect(double tt=0, double xl=1, double xh=1, double yl=1, double yh=1, double R = 1); double TestBox(double tt=0); double TestTri(double tt=0); } diff --git a/aragorn.C b/aragorn.C index e46e34b..8ec2abd 100644 --- a/aragorn.C +++ b/aragorn.C @@ -1,41 +1,45 @@ #include "aragorn.h" #include "quink.C" #include -#include "TCanvas.h" +#include +#include +#include +#include +#include using namespace std; -double matpar [5] = {30, 60, 13, 40, 2000}; //Parameters of the materiasl and the detector +double matpar [6] = {40, 60, 13, 40, 2000, 0}; //Parameters of the materials and the detector +TDatime now; + +float aragorn::hw; +float aragorn::hh; +int aragorn::rows; +int aragorn::columns; aragorn *aragorn::__instance=0; aragorn::aragorn() {cout << "Their pace is quickening." << endl;} //Just to confirm it runs correctly -/* NONE OF THIS IS CORRECT ANYMORE, BUT I WILL FIX IT SOON -The following is an example of how the quink array would look, with the quinks numbered as if they have already -correctly been associated with each other for x and y equal to 4: - On top is electrode 0 - ------------------------- - | 42 43 44 45 46 47 48 | - | 35 36 37 38 39 40 41 | - | 28 29 30 31 32 33 34 | On the right is electrode 1 - On the left is electrode 3 | 21 22 23 24 25 26 27 | - | 14 15 16 17 18 19 20 | - | 7 8 9 10 11 12 13 | - | 0 1 2 3 4 5 6 | - ------------------------- - On bottom is electrode 2 -*/ void aragorn::Construct() { + ifstream in; + in.open("Construct.dat"); + int c,r; + double a,b; + in>>a>>b>>c>>r; + Sethw(a); + Sethh(b); + Setcol(c); + Setrow(r); if(TheQuinks.size()>0){cout<<"This setup has already been constructed.";} else { - for(int j=hh;j<2*hh*rows;j=j+2*hh) + for(int j=1;j<=2*rows;j+=2) { //Creating all quinks, one row at at time and then moving down a row. - for(int i=hw;i<2*hw*columns;i=i+2*hw) + for(int i=1;i<=2*columns;i+=2) { - quink * newest =(new quink(double(i-hw),double(i+hw),double(j-hh),double(j+hh))); + quink * newest =(new quink(double(i*hw-hw),double(i*hw+hw),double(j*hh-hh),double(j*hh+hh))); TheQuinks.push_back(newest); } } @@ -51,30 +55,170 @@ void aragorn::Report() void aragorn::SetChargeCloud(double a,double b) { - //eventually this function will hopefully + // eventually this function will hopefully // take some x and y position and turn it // into a circle of charge with a dense // center and diffuse radial charge. // for today however, it does nothing. } -void aragorn::Event() +// This is just a function to make sure that the setup was constructed correctly +void aragorn::CheckPad(int n) +{ + if(n>TheQuinks.size()-1){cout<<"This is not a valid pad number.";} + else if (n<0){cout<<"This is not a valid pad number";} + else{TheQuinks[n]->Report();} +} + +void aragorn::Observe() { - double x1,y1; //x and y coordinate of charge cloud - - cout<<"Where should the charge cloud hit?"<>x1>>y1; - TCanvas* c1 = new TCanvas("c1", "Readout", 800, 600); - c1->Divide(columns,rows); - for(int i =1;i<=TheQuinks.size();i++) + int n = 7; + vector < vector > values(n+1); + int step = 40/(n); //tmax divided by the number of time samples + TH1F * h1 = new TH1F(); + char buffer [50]; + TCanvas * c1 = new TCanvas("c1","Progression of Charge"); + //TFile f("C:\\Users\\Overlord\\Desktop\\Programs\\Research\\ChargeDiv\\PRF\\PRFData3.root"); + TFile f("C:\\root_v5.34.25\\macros\\SmallPads.root"); + for (int i = 0; i <= n; ++i) + { + for (int k = 1; k < TheQuinks.size()+1; ++k) + { + sprintf(buffer,"Pad_li_%i",k); + h1 = (TH1F*)f.Get(buffer); + values[i].push_back(h1->GetBinContent(i*step+1)); //Value at a time i(index)*n(step size) + delete(h1); + } + } + c1->Divide(4,2); //The product of these need to be the n value +1; + double max; + double myval; + for (int i = 0; i <= n; ++i) { - c1->cd(i); - TheQuinks[i-1]->Display(matpar, x1, y1); - cout<cd(i+1); + gPad->Divide(columns,rows,0,0); + for (int j = 0; j < values[i].size(); ++j) + { + gPad->cd(j+1); + myval = abs(values[i][j]); + if(myval>.95*max){gPad->SetFillColor(kRed);} + else if((myval<.95*max)&&(myval>=.5*max)){gPad->SetFillColor(kRed+1);} + else if((myval<.5*max)&&(myval>=.3*max)){gPad->SetFillColor(kRed+2);} + else if((myval<.3*max)&&(myval>=.1*max)){gPad->SetFillColor(kRed+3);} + else if((myval<.1*max)&&(myval>=.01*max)){gPad->SetFillColor(kRed+4);} + else {gPad->SetFillColor(kBlack);} + gPad->Draw(); + c1->cd(i+1); + } } } -// This is just a function to make sure that the setup was constructed correctly -void aragorn::CheckPad(int n) + +void aragorn::SingleEvent() +{ + ifstream in; + in.open("Input.dat"); + in.clear(); + double x1,y1; + in>>x1>>y1; + in.close(); + char name[40]; + bool newFile = true; + char filename[40]; + int timenow=now.Convert(); + sprintf(filename,"Single_%i.root",timenow); + for (int j = 0; jToFile(matpar, x1, y1, name, filename, newFile); + cout<<"Pad "<>x1>>y1>>stepx>>stepy>>xev>>yev; + in.close(); + char name[40]; + char filename[40]; + bool newFile = true; + int timenow=now.Convert(); + sprintf(filename,"PRF_%i.root",timenow); + TStopwatch loopclock; + loopclock.Start(); + for(int j=0;jToFile(matpar, xtemp, ytemp, name, filename, newFile); + } + cout<Report(); + ifstream in; + in.open("Input.dat"); + double x1,y1,stepx,stepy,xtemp,ytemp,r1,stepr; + int xev,yev,rev; + in>>x1>>y1>>stepx>>stepy>>xev>>yev>>r1>>stepr>>rev; + in.close(); + cout<<"R start "<ToFile(matpar, xtemp, ytemp, name, filename, newFile); + } + cout< > x) +{ + for(int i=0;i TheQuinks; - - static const int hw=1; //half of the width of each pad - static const int hh=3; //half of the height of each pad - static const int columns = 5; //How many pads in a row - static const int rows = 3; // How many pads in each column - + void Sethw(double a){hw=a;} + void Sethh(double b){hh=b;} + void Setrow(int r){rows=r;} + void Setcol(int c){columns=c;} void Report(); void Construct(); // - void Event(); + void SingleEvent(); void SetChargeCloud(double,double); + void Observe(); void CheckPad(int); + void MakePRF(); + void VaryResistance(); + void Examine(vector< vector >); aragorn(); static aragorn *__instance; diff --git a/quink.C b/quink.C index 98b4981..cf984fd 100644 --- a/quink.C +++ b/quink.C @@ -1,9 +1,7 @@ #include "quink.h" #include -#include #include "TCanvas.h" #include "ConvoluteLoop.C" - using namespace std; quink::quink(double Xl, double Xh,double Yl, double Yh) @@ -14,9 +12,16 @@ quink::quink(double Xl, double Xh,double Yl, double Yh) yh = Yh; } -void quink::Display(double t[5], double &x, double &y) //Start time, end time, sample interval, and position of charge cloud. +// This is just a function used to make sure it was constructed correctly. +void quink::Report() +{ + char name[40]; + sprintf(name,"(%2.2f,%2.2f) to (%2.2f,%2.2f)",xl,yl,xh,yh); + cout< Readings; double relativepos [4]; // creating the relative positions, taking the positions of the pad // and subtracting the position of the charge cloud hit @@ -24,14 +29,9 @@ void quink::Display(double t[5], double &x, double &y) //Start time, end time, s relativepos[1] = (xh-x); relativepos[2] = (yl-y); relativepos[3] = (yh-y); + char title[80]; + sprintf(title,"(%2.2f,%2.2f) to (%2.2f,%2.2f) with chg-ctr at (%2.2f,%2.2f)",relativepos[0],relativepos[1],relativepos[2],relativepos[3],x,y); // Pass the array with all of the material constants // and the relative position to the convolution - ConvoluteLoop(t,relativepos); -} -// This is just a function used to make sure it was constructed correctly. -void quink::Report() -{ - char name[40]; - sprintf(name,"(%g,%g) to (%g,%g)",xl,yl,xh,yh); - cout<