123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905906907908909910911912913914915916917918919920921922923924925926927928929930931932933934935936937938939940941942943944945946947948949950951952953954955956957958959960961962963964965966967968969970971972973974975976977978979980981982983984985986987988989990991992993994995996997998999100010011002100310041005100610071008100910101011101210131014101510161017101810191020102110221023102410251026102710281029103010311032103310341035103610371038103910401041104210431044104510461047104810491050105110521053105410551056105710581059106010611062106310641065106610671068106910701071107210731074107510761077107810791080108110821083108410851086108710881089109010911092109310941095109610971098109911001101110211031104110511061107110811091110 |
- //
- // ================================================================================
- // Straw tube response simulation
- //
- //
- //
- // mandatory call at the beginning
- // once to set constants
- //
- // void TConst(Double_t Radius, Double_t pSTP, Double_t ArP, Double_t CO2P);
- //
- // Pysical (garfield) values:
- // Radius Press (pSTP) % Ar (ArP) % CO2 (CO2P)
- // 0.4 1 0.9 0.1
- // 0.5 1 0.9 0.1
- // 0.5 2 0.8 0.2
- //
- //
- // initialization at each track for MC applications
- //
- // define particle and its in-out hits
- // void TInit(Double_t Mass, Double_t Momentum, Double_t InOut[])
- //
- // void PutWireXYZ(Double_t w1, Double_t w2, Double_t w3,
- // Double_t w4, Double_t w5, Double_t w6); define the wire
- //
- // Double_t PartToTime(Double_t Mass, Double_t Momentum,
- // Double_t InOut[] ); straw time
- //
- // Double_t PartToADC(); ADC charge
- // signal corresponding to the energy loss of StrawCharge
- // for energy loss simulation
- //
- // where:
- //
- // pSTP = pressure in bar
- // Radius = tube radius
- // ArP and CO2P percentages of Argon and CO2
- // w1:w6 = coordinates of the extremes of the straw wire (straw axis)
- // Mass, Momentm = mass and momentum of the MC particle
- // InOut[0:5] = straw input and output coordinates traversed from the MC particle
- //
- // Routine for off line reconstruction
- //
- // Double_t TimnsToDiscm(Double_t time); from time (ns) to
- // radius in cm
- // where:
- // time: measured drift time in ns
- //
- //
- //
- // For other applications see the routines listed in TStraw.h
- //
- // Author: A. Rotondi (november 2005, revised and extended in August 2006)
- //
- // ================================================================================
- //
- #include "TRandom.h"
- #include "TMath.h"
- #include "TStraw.h"
- #include "TImage.h"
- #include "TVector3.h"
- #include "TVectorD.h"
- #include "TMatrixD.h"
- #include "TMatrixDEigen.h"
- ClassImp(TStraw);
- // ============================================================
- TStraw::TStraw() {
- CNumb=0;
- memset(CumClus,0,sizeof(CumClus));
- memset(PolyaCum,0,sizeof(PolyaCum));
- memset(Xs,0,sizeof(Xs));
- // class constructor
- // clear
- CNumb=0;
- CDist.clear();
- CDistC.clear();
- CNele.clear();
- CNeleT.clear();
- TeleTime.clear();
- AmplSig.clear();
- Pulse.clear();
- PulseT.clear();
- WDist.clear();
- }
- // ----------------------------------------------------------------------
- void TStraw::TConst(Double_t R1, Double_t P1, Double_t A1, Double_t C1) {
- // set constants for the simulation
- // Input
- // P1 = tube absolute pressure pSTP
- // R1 = tube radius Radius
- // A1 = argon percentage ArPerc
- // C1 = C=2 percentage C=2Perc
- //
- Radius=R1;
- pSTP=P1;
- // Input for the media (volume percentages)
- ArPerc = A1;
- CO2Perc = C1;
- // cluster dimensions in Ar and CO2 (experimantal values)
- // Double_t PClus[20] =
- // {0., .656, .150, .064, .035, .0225, .0155, .0105, .0081,
- // .0061, .0049, .0039, .0030, .0025, .0020, .0016, .0012,
- // .00095, .00075, .00063}; // Fischle fit
- Double_t PClus[20] =
- {0., .802, .0707, .020, .013, .008, .006, .005, .006,
- .008, .009, .007, .0050, .0040, .0033, .0029, .0025,
- .0023, .0022, .002}; // Lapique 1st calculation
- // Double_t PClus[20] =
- // {0., .841, .0340, .021, .013, .008, .006, .004, .003,
- // .008, .013, .008, .0050, .004, .0030, .0028, .0025,
- // .0023, .0022, .002}; // Lapique 2nd calculation
- // Double_t PClus[20] =
- // {0., .656, .150, .064, .035, .0225, .0155, .0105, .0081,
- // .0061, .0049, .0039, .0030, .0025, .0020, .0016, .0012,
- // .00080, .00059, .00045}; // Fischle empirical
- // PDouble_t Clus[20] =
- // {0., .656, .148, .0649, .0337, .0244, .0141, .0078, .0095,
- // .0063, .0062, .0042, .0028, .0018, .0023, .0017, .0014,
- // .00060, .00050, .00063}; // Fischle exp
- Double_t CO2Clus[20] =
- {0., .730, .162, .038, .020, .0110, .0147, .0060, .0084,
- .0052, .0020, .0042, .0021, .0025, .0038, .0021, .0009,
- .00013, .00064, .00048}; // Fischle exp
- // Double_t CH4Clus[20] =
- // {0., .786, .120, .032, .013, .0098, .0055, .0057, .0027,
- // .0029, .0020, .0016, .0013, .0010, .0012, .0006, .0005,
- // .00042, .00037, .00033}; // Fischle exp
-
- CH4Perc = 0.07;
- // -----------------------------------------------------
- // gain of the avalanche
- // Ar/CO2 90/10 1 bar (NTP) MAGY GARFIELD 20-7-2006
- GasGain=90000.;
- // argon ----------------------------------------------------
- AAr = 39.948; // Argon (39.948)
- ZAr= 18.0; // Argon (18)
- RhoAr = pSTP*1.78*1.e-03; // g/cm3 (1.78 mg/cm3)
- IAr = 188*1.e-09; // ionization potential (GeV) (188 eV)
- WiAr =26.7; // energy to create an ion pair (standard 26.7 eV)
- NclAr = 23.; // cluster/cm in Argon
- // CO2 -----------------------------------------------------
- ACO2 = 44; // CO2
- ZCO2 = 22.; // CO2
- RhoCO2 = pSTP*1.98*1.e-03; // g/cm3 CO2 (1.98 mg/cm3)
- ICO2 = 95.8*1.e-09; // ionization potential (GeV) (188 eV)
- WiCO2 = 33.0; // energy to create an ion pair (33 eV)
- NclCO2 = 35.5; // clusters/cm CO2 35.5
- // Methane CH4 ---------------------------------------------------------
- ACH4 = 16; // CO2 (39.948)
- ZCH4 = 10.; // CO2 (18)
- RhoCH4 = pSTP*0.71*1.e-03; // g/cm3 CO2 (0.71 mg/cm3)
- ICH4 = 40.6*1.e-09; // ionization potential (GeV) (188 eV)
- WiCH4 = 28.0; // energy to create an ion pair
- NclCH4 = 25.0;
- // Input for the media (weight percentages) ----------------------------
- ArWPerc = ArPerc *AAr /(ArPerc*AAr + CO2Perc*ACO2);
- CO2WPerc = CO2Perc*ACO2/(ArPerc*AAr + CO2Perc*ACO2);
-
- // mixture densiies ----------------------------------------------------
- RhoMixCO2 = 1./((ArWPerc/RhoAr) + (CO2WPerc/RhoCO2));
- RhoMixCH4 = 1./((ArWPerc/RhoAr) + (CH4WPerc/RhoCH4));
- //----------------------------------------------------------------------
- // particles (Gev, energy losses in Kev)
- PZeta = 1; // projectile charge
- piMass = 0.139; // particle mass (GeV)
- eMass = 0.511/1000.; // electron mass (GeV) (0.511 MeV)
- prMass = 0.93827; // proton mass (GeV)
- // ---------------------------------------------------------------------
- // thresholds for the straw tubes (default values) see TInit for current values
- Thresh1=10;
- Thresh2=30;
- // channels for the signal
- Nchann = 500;
- // ---------------------------------------------Emin------------------------
- NPolya= 100; // steps for the calculation of the Polya distributions
- Xmax=0.; // Polya istribution is calculated between o and Xmax (see Polya)
- bPolya = 0.5;
- Polya(bPolya); // cumulative of the Polya distribution
- // -----------------------------------------------------------------------
- // cumulative for the number of electron per cluster
- Double_t Wnorm = (ArPerc*NclAr + CO2Perc*NclCO2);
- CumClus[0]=(ArPerc*NclAr*PClus[0] + CO2Perc*NclCO2*CO2Clus[0])/Wnorm;
- //for(Int_t i=1; i<=20; i++) --->> out of bounds was changed
- for(Int_t i=1; i<20; i++) {
- CumClus[i]=(ArPerc*NclAr*PClus[i] + CO2Perc*NclCO2*CO2Clus[i])/Wnorm
- + CumClus[i-1];
- }
- CumClus[20]=1.;
- Double_t sum=0.;
- for(Int_t i=0; i<=19; i++) {
- sum += PClus[i];
- //cout<<" PClus["<<i<<"] = "<<PClus[i]<<endl;
- }
- //for(Int_t i=0;i<=20;i++) {cout<<"CumClus["<<i<<"] = "<<CumClus[i]<<endl;}
- //cout<<" Sum of Probabilities = "<<sum<<endl;
- }
- // ==============================================================
- void TStraw::PutWireXYZ(Double_t w1, Double_t w2, Double_t w3,
- Double_t w4, Double_t w5, Double_t w6){
- // get wire coordinates
- Wx1=w1;
- Wy1=w2;
- Wz1=w3;
-
- Wx2=w4;
- Wy2=w5;
- Wz2=w6;
- }
- // =============================================================
- void TStraw::TInit(Double_t xPMass, Double_t xPMom, Double_t InOut[]) {
- // initialization of the constants for each track and each straw
- // transfer of data
- // track geometrical quantities
- Xin=InOut[0]; Yin=InOut[1]; Zin=InOut[2];
- Xout=InOut[3]; Yout=InOut[4]; Zout=InOut[5];
- PMass = xPMass;
- PMom = xPMom;
-
- // path into the straw
- Dx = sqrt((Xout-Xin)*(Xout-Xin) +
- (Yout-Yin)*(Yout-Yin) +
- (Zout-Zin)*(Zout-Zin));
- // clear
- CNumb=0;
- PulseMax=0;
- PulseTime=0;
- CDist.clear();
- CDistC.clear();
- CNele.clear();
- CNeleT.clear();
- WDist.clear();
- TeleTime.clear();
- AmplSig.clear();
- Pulse.clear();
- PulseT.clear();
- // ---------------------------------------------------------------------
- // initialization of the constants
- // maximum energy transfer Emax (GeV) and related quantities
-
- PEn = sqrt(PMom*PMom + PMass*PMass); // particle energy GeV
- beta = PMom/PEn;
- gamma= 1/sqrt(1.-beta*beta);
- Double_t begam = beta*gamma; // local
- Double_t mratio= eMass/PMass; // local
- // calculation of the polarization Sternheimer factor for Argon
- Double_t ms = 2.80;
- Double_t Xst = log10(begam*sqrt(pSTP));
- Double_t Cs =-11.92;
- Double_t X1 = 4;
- Double_t Xo = 1.96;
- Double_t as = 0.389;
- Delta=0;
- if(Xo <= Xst && Xst <= X1) Delta = 4.6051702*Xst +Cs +as*pow(X1-Xst,ms);
- else if(X1< Xst) Delta = 4.6051702*Xst +Cs;
- if(Delta < 0) Delta=0;
- // calculation of other typical quantites
- Emax = (1.022*begam*begam/(1.+2.*gamma*mratio+mratio*mratio))/1000.; // GeV
- CsiAr = pSTP*0.5*0.3071*PZeta*PZeta*ZAr*RhoAr/(beta*beta*AAr)*1e-03; //GeV
- Double_t fact = 2.*eMass*beta*beta*gamma*gamma*Emax/(IAr*IAr);
- EmedAr = 2.*CsiAr*(0.5*TMath::Log(fact)-beta*beta - 0.5*Delta); // GeV/cm
- EminAr = pSTP*2.70*1.e-06; // GeV/cm (2.5 keV)
- // most prob energy (GeV)
- EmpAr = EmedAr + CsiAr*(0.422784 + beta*beta + log(CsiAr/Emax));
- CsiCO2 = pSTP*0.5*0.3071*PZeta*PZeta*ZCO2*RhoCO2/(beta*beta*ACO2)*1e-03; //GeV
- fact = 2.*eMass*beta*beta*gamma*gamma*Emax/(ICO2*ICO2);
- EmedCO2 = 2.*CsiCO2*(0.5*TMath::Log(fact)-beta*beta - 0.5*Delta); // GeV/cm
- EminCO2 = pSTP*3.60*1.e-06; // GeV/cm (2.5 keV)
- // most prob energy (GeV)
- EmpCO2 = EmedCO2 + CsiCO2*(0.422784 + beta*beta + log(CsiCO2/Emax));
- Csi = RhoMixCO2 * ((ArWPerc*CsiAr/RhoAr) + (CO2WPerc*CsiCO2/RhoCO2));;
- Emed = RhoMixCO2 * ((ArWPerc*EmedAr/RhoAr) + (CO2WPerc*EmedCO2/RhoCO2));
- Emin = RhoMixCO2 * ((ArWPerc*EminAr/RhoAr) + (CO2WPerc*EminCO2/RhoCO2));
- Emp = RhoMixCO2 * ((ArWPerc*EmpAr/RhoAr) + (CO2WPerc*EmpCO2/RhoCO2));
- // mean weighted interaction
- Wi = (ArPerc*NclAr*WiAr + CO2Perc*NclCO2*WiCO2)
- /(ArPerc*NclAr + CO2Perc*NclCO2);
-
- Ncl=pSTP*((ArPerc*NclAr) + (CO2Perc*NclCO2))*Emed/Emin; // <--- Cluster/cm
-
- Lcl=1./Ncl; // mean free path between clusters (cm)
- Ecl = 2.8; // mean number of electrons per clusters
- Cutoff = Ncl*Ecl*20; // limit to the number of primary electrons
- Ntote = Ecl*Ncl; // total mean number of electrons per cm
- // ---------------------------------------------------------------------
- // thresholds for the straw tubes (current values)
- // total number of electrons * scaling factor (max of signal x electron=1)
- Thresh1=Ncl*Ecl* 0.05;
- Thresh2=Ncl*Ecl* 0.15;
- // -----------------------------------------------------------------------
- TDirCos(); // director cosine of the track
-
- // control prints
- // cout<<" Dx "<<Dx<<" Csi "<<Csi*1.e+09<<" Emax MeV "<<
- // Emax*1.e+03<<" PMass "<<PMass<<" gamma "<<
- // gamma<<endl;
- // cout<<" RRise "<<RRise(gamma)<<endl;
- // cout<<" Thresh1 = "<<Thresh1<<" Thresh2 = "<<Thresh2<<endl;
- }
- // =============================================================
- // energy loss in GeV
- Double_t TStraw::STEloss() {
- return gRandom->Landau(Emp*Dx,Csi*Dx); // in GeV
- }
- // =============================================================
- Double_t TStraw::StrawCharge() {
- // sampling of the ionization energy loss, number of
- // primary electrons and their wire distances for each track
- // energy loss in GeV
-
- // clear
- Double_t TotEnEle = 0., ClusEner=0., Ecurr;
- CNeleT.clear();
- WDist.clear();
- // threshold for delta rays (>1.5 KeV from NIM A301(1991)202)
- Double_t thr=Csi*Dx*1.e+09/(Ntote*1500.);
- // total energy released by the electrons
- //Number of cluster in the length
- NNClus = Cluster(); // set vectors CDist CDistC e CNele
- // calculation of total energy released by the electrons
- for(Int_t j=1;j<=NNClus;j++){
- ClusEner = 0.;
- for(Int_t jj=1; jj<=(Int_t)(CNele.at(j-1)); jj++){
- Ecurr = Wi;
- TotEnEle += Ecurr;
- ClusEner += Ecurr;
- // delta rays simulation over 1.keV. Upper threshol 12 keV
- if(gRandom->Uniform() <= thr) {
- TotEnEle += int(Wi/(1.002- gRandom->Uniform()));
- ClusEner += int(Wi/(1.002- gRandom->Uniform()));
- }
- }
- // effective number of electrons per cluster
- CNeleT.push_back(int(ClusEner/Wi));
- }
- // record the distance from the wire of any electron
- // adding the diffusion effects
- Int_t owfl =0;
- for(Int_t k=1; k<=NNClus; k++) {
- for (Int_t kk=0; kk<(Int_t)CNeleT.at(k-1); kk++) {
- if(owfl > (Int_t) Cutoff) break;
- owfl++;
- WDist.push_back(WDistCalc(CDistC.at(k-1)));
- }
- }
-
- return TotEnEle*1.e-09; // in Gev
- }
- // =============================================================
- Int_t TStraw::Cluster() {
- // calculate the number of clusters CNumb
- // their distance CDist end
- // their number of electrons CNele
- CNumb=0;
- CDist.clear();
- CDistC.clear();
- CNele.clear();
- Double_t DisTot = 0;
- for(Int_t k=1; k<=1000; k++) {
- // distance
-
- Double_t path = -Lcl*log(1-gRandom->Uniform()+0.000001);
- DisTot+= path;
- if(DisTot>Dx) break;
- CNumb=k;
- CDist.push_back(path);
- CDistC.push_back(DisTot);
- CNele.push_back((Double_t)Eject() );
- }
-
- return CNumb;
- }
- // =============================================================
- Int_t TStraw::Eject() {
- // find the number of electrons in a cluster
- Int_t Inelect;
- Double_t nelect=0.;
-
- nelect= TMath::BinarySearch(21,CumClus,gRandom->Uniform());
- if(nelect<19) {
- return Inelect = (Int_t) (nelect) +1;
- }
- else {
- return Inelect = (Int_t)(20./(1.01-gRandom->Uniform()));
- }
- }
- // =============================================================
- Double_t TStraw::RRise(Double_t gamma_r) {
- // interpolate the relativisic rise of the
- // number of cluster per cm, starting from the one
- // measured at the ionization minimum
- Double_t Rise;
- Double_t lg = log10(gamma_r);
- if(1.<=gamma_r && gamma_r <= 2.2){
- Rise = -2.159*lg +1.7;
- }
- else if(2.2<=gamma_r && gamma_r <= 6.){
- Rise = 1.;
- }
- else if(6.<=gamma_r && gamma_r <= 200.){
- Rise = 0.302*lg + 0.765;
- }
- else if(200.<=gamma_r && gamma_r <= 1000.){
- Rise = 0.1431*lg + 1.131;
- }
- else{
- Rise= 1.54;
- }
- return Rise;
- }
- // ==========================================================================
- void TStraw::Polya(Double_t bpar) {
- // calculate the cumulative of the Polya distribution
- // for samplig the gain fluctuations
- Double_t eps = 0.0001;
- Double_t Dxx = 0.05, xx=0., x1,x2;
- Double_t PMax;
-
- Double_t k=1./bpar -1.;
- // find Xmax
- PMax = eps * pow(k,k)*exp(-k);
- Double_t value = 1.e+06;
- Xmax =2*k;
- while(value > PMax){
- Xmax +=Dxx;
- value = pow(Xmax,k)*exp(-Xmax);
- }
- Xmax += -0.5*Dxx;
- // calculate the cumulative
- Double_t dx= Xmax/NPolya;
- Xs[0]=0;
- for(Int_t i=1; i<NPolya; i++){
- x1 = xx;
- xx += dx;
- x2= xx;
- Xs[i]=x2;
- if(i>1)
- PolyaCum[i] = PolyaCum[i-1] +
- 0.5*dx*(pow(x1,k)*exp(-x1) + pow(x2,k)*exp(-x2))/TMath::Gamma(k+1);
- else
- PolyaCum[i] = PolyaCum[i-1] +
- 0.5*dx* pow(x2,k)*exp(-x2)/TMath::Gamma(k+1);
- }
- // adjust the normalization
- for(Int_t ii=0; ii<NPolya; ii++) PolyaCum[ii] /= PolyaCum[NPolya-1];
-
- }
- // =====================================================================
- Double_t TStraw::PolyaSamp() {
- // sampling a wire gain fluctuation in the gas
- Double_t xr = gRandom->Uniform();
- Int_t n = TMath::BinarySearch(NPolya,PolyaCum,xr);
- Double_t xsamp = Xmax;
- if(n<NPolya-1){
- xsamp= Xs[n] +
- ((xr-PolyaCum[n])/(PolyaCum[n+1] - PolyaCum[n])) *(Xs[n+1]-Xs[n]);
- }
- return xsamp;
- }
- // =====================================================================
- TVector3 TStraw::WDistCalc(Double_t DisTot) {
- // calculation of the distance of the cluster from the wire
- // DisTot is the cluster distance from the track entrance
- // diffusion effects are considered
- // wire director cosines
- Double_t Wlength = sqrt( (Wx2-Wx1)*(Wx2-Wx1) +
- (Wy2-Wy1)*(Wy2-Wy1) +
- (Wz2-Wz1)*(Wz2-Wz1) );
- Wp = (Wx2-Wx1)/Wlength;
- Wq = (Wy2-Wy1)/Wlength;
- Wr = (Wz2-Wz1)/Wlength;
- // current track coordinates
- Double_t xcor = Calpha * DisTot + Xin;
- Double_t ycor = Cbeta * DisTot + Yin;
- Double_t zcor = Cgamma * DisTot + Zin;
-
- // cross product VV1= (wire x track) for the distance
- Double_t XX1 = Wr*(ycor-Wy1) - Wq*(zcor-Wz1);
- Double_t YY1 = Wp*(zcor-Wz1) - Wr*(xcor-Wx1);
- Double_t ZZ1 = Wq*(xcor-Wx1) - Wp*(ycor-Wy1);
-
- // vector for the 3-D distance from the wire (from wire x VV1)
- Double_t XX = Wq*ZZ1 - Wr*YY1;
- Double_t YY = Wr*XX1 - Wp*ZZ1;
- Double_t ZZ = Wp*YY1 - Wq*XX1;
- //cout<<" XYZ "<<XX<<" "<<YY<<" "<<ZZ<<endl;
- Double_t DDistcm = sqrt(XX*XX + YY*YY +ZZ*ZZ);
- //director cosines of the distance vector
- Double_t cosx = XX/DDistcm;
- Double_t cosy = YY/DDistcm;
- Double_t cosz = ZZ/DDistcm;
-
- // sampling of the diffusion
- //Double_t DDist = DDistcm*10.; // distance in mm
- // longitudinal coefficient of diffusion (GARFIELD) cm --> micron
- // MAGY Ar/CO2 90/10 20-7-2006 GARFIELD
- Double_t SigL = DiffLong(DDistcm);
- // ... in cm
- SigL *= 1.e-04; // in cm
- // tranverse coefficient (GARFIELD) cm--> micron
- // MAGY Ar/CO2 90/10 20-7-2006 GARFIELD
- Double_t SigT = DiffTran(DDistcm);
- SigT *= 1.e-04; // in cm
- // sampling of Longitudinal and Transverse diffusion
- Double_t difL = gRandom->Gaus(0.,SigL);
- Double_t difT = gRandom->Gaus(0.,SigT);
-
- // vector addition to the distance
- // the transverse component has the same dir cos of the wire
- XX += difL*cosx + difT*Wp;
- YY += difL*cosy + difT*Wq;
- ZZ += difL*cosz + difT*Wr;
- //cout<<" XYZ+ dif "<<XX<<" "<<YY<<" "<<ZZ<<endl;
- //cout<<" --------------------------------------------------- "<<endl;
- TVector3 TDist(XX,YY,ZZ);
- return TDist;
- }
- // =====================================================================
- void TStraw::TDirCos() {
- // director cosines of the track (called for each track)
- //path into the straw
- Rpath = sqrt((Xout-Xin)*(Xout-Xin) +
- (Yout-Yin)*(Yout-Yin) +
- (Zout-Zin)*(Zout-Zin));
- //director cosines
- Calpha = (Xout-Xin)/Rpath;
- Cbeta = (Yout-Yin)/Rpath;
- Cgamma = (Zout-Zin)/Rpath;
- }
- // =====================================================================
- Int_t TStraw::TimeEle() {
- TeleTime.clear();
- // calculate the arrival times (ns) for each electron in TeleTime
- // from cm to ns <---------------------------------
- // return the number of electrons
- // Author:A. Rotondi 20-7-2006
- Int_t ido = WDist.size();
- Double_t etime;
- // cutoff the total charge to 20 times the average
- if(ido > (Int_t) Cutoff) ido = (Int_t) Cutoff;
- for(Int_t k=0; k<ido; k++) {
- Double_t dst = WDist.at(k).Mag(); // distance in cm
- // space to time calculation from GARFIELD
- if(pSTP <1.9){
- if(Radius<0.5){
- // V=1600V diameter 4 mm 90/10
- etime= -1.624e-05 + 0.1258*dst + 0.8079*pow(dst,2) -2.918*pow(dst,3)
- + 10.33*pow(dst,4) -10.84*pow(dst,5);
- }
- else{
- // V=1600V diameter 5 mm 90/10
- etime= -6.763e-05 + 0.1471*dst +0.3625*pow(dst,2) +0.3876*pow(dst,3)
- +1.04*pow(dst,4) -1.693*pow(dst,5);
- }
- }
- else{
- // 2 bar 2000V, 5 mm, 80/20
- etime= -0.0001014 + 0.1463*dst -0.1694*pow(dst,2) +2.4248*pow(dst,3)
- -1.793*pow(dst,4);
- }
- etime*=1000.; // nano seconds
- TeleTime.push_back(etime);
- }
- return TeleTime.size();
- }
- // =====================================================================
- Double_t TStraw::Signal(Double_t t, Double_t t0) {
- // electric signal at time t of a cluster arriving at t0
- Double_t elesig = 0;
- Double_t x= t - t0;
- if (x>0){
- Double_t A = 1.03e-03;
- Double_t B = 3.95;
- Double_t C = 0.228;
- Double_t D = -3.839e-02;
- Double_t E = 1.148e-03;
- elesig = A*exp(B*log(x)-C*x)*(1+D*x+E*x*x);
- }
- return elesig;
- }
- // =====================================================================
- Int_t TStraw::StrawSignal(Int_t nsteps) {
- // creation of nstep values of
- // the straw global Pulse (sum on all clusters)
- // return the number of primary electrons
- PulseMax=0;
- Pulse.clear();
- PulseT.clear();
- AmplSig.clear();
- Int_t neltot = TimeEle(); // creation and size of TeleTime (electron times)
- Double_t Tmax = 1.e-25;
- for(Int_t k=0; k< neltot; k++) if(Tmax<TeleTime.at(k)) Tmax=TeleTime.at(k);
- Tmax += 100.;
-
- Double_t Dt = Tmax/nsteps; // number of steps of the signal
- //AmplSig is the amplitude of each electron
- for(Int_t j=0; j< neltot; j++){
- AmplSig.push_back(bPolya*PolyaSamp());
- }
- // creation of the signal Pulse(PulseT) time in ns
- for(Int_t j=0; j< nsteps; j++){
- Double_t te = j*Dt;
- Double_t sumele=0.;
- for(Int_t jj=0; jj< neltot; jj++){
- Double_t te0 = TeleTime.at(jj);
- sumele += AmplSig.at(jj)*Signal(te,te0);
- }
-
- Pulse.push_back(sumele);
- PulseT.push_back(te);
- }
- // add a random noise (3% of maximum) plus
- // a 1% of 200 ns periodic signal
- Double_t Pmax = 1.e-25;
- for(Int_t k=0; k<(Int_t)Pulse.size(); k++)
- if(Pmax<Pulse.at(k)) Pmax=Pulse.at(k);
- for(Int_t k=0; k<(Int_t)Pulse.size(); k++){
- Pulse.at(k) += 0.03*Pmax*gRandom->Uniform()
- + 0.01*Pmax*(sin(6.28*PulseT.at(k)/120.));
- // if(Pulse[k]<0) Pulse[k] *= -1.;
- }
- PulseMax=Pmax;
- // set variable threshold for the signals
- Thresh1=0.05*Pmax;
- Thresh2=0.15*Pmax;
- return neltot;
- }
- // =====================================================================
- Int_t TStraw::StrawTime() {
- // simulate the discrimination of the straw signal
- // and give the time
- // discriminator technique: set 2 threshold, select the first one
- // (the low one) only if the second one is fired (FINUDA system)
- PulseTime=0.;
- Int_t ind=0;
- Int_t flag1=0;
- Int_t flag2=0;
- for(Int_t k=0; k<(Int_t) Pulse.size(); k++){
- if(flag1==0 && Pulse[k]>Thresh1) {
- flag1=1;
- ind=k;
- }
- if(Pulse[k]<Thresh1) {flag1=0; ind=0;} // reset if signal decreases
- if(flag1==1 && Pulse[k]>Thresh2) {
- flag2=1;
- break;
- }
- }
- if(flag1==1 && flag2==1) PulseTime=PulseT.at(ind);
-
- return ind;
- }
- // =====================================================================
- Double_t TStraw::TrueDist(Double_t Point[]) {
- // service routine that finds the distance in cm from the wire
- // by knowing the wire coordinates (class variables)
- // and the input-output points Point[6]
- Double_t truedist = 0;
- // wire director cosines
- Double_t Wlength = sqrt( (Wx2-Wx1)*(Wx2-Wx1) +
- (Wy2-Wy1)*(Wy2-Wy1) +
- (Wz2-Wz1)*(Wz2-Wz1) );
- Wp = (Wx2-Wx1)/Wlength;
- Wq = (Wy2-Wy1)/Wlength;
- Wr = (Wz2-Wz1)/Wlength;
- // director cosines of the given track
- Double_t Modu = sqrt( (Point[3]*Point[0])*(Point[3]*Point[0]) +
- (Point[4]*Point[1])*(Point[4]*Point[1]) +
- (Point[5]*Point[2])*(Point[5]*Point[2]) );
- Double_t dcx = (Point[3]-Point[0])/Modu;
- Double_t dcy = (Point[4]-Point[1])/Modu;
- Double_t dcz = (Point[5]-Point[2])/Modu;
-
- //distance formula
- Double_t p1 = (Point[0]-Wx1)*(dcy*Wr-dcz*Wq);
- Double_t p2 =-(Point[1]-Wy1)*(dcx*Wr-dcz*Wp);
- Double_t p3 = (Point[2]-Wz1)*(dcx*Wq-dcy*Wp);
- Double_t Det = p1+p2+p3;
- Double_t Disc = sqrt( (dcy*Wr-dcz*Wq)*(dcy*Wr-dcz*Wq) +
- (dcz*Wp-dcx*Wr)*(dcz*Wp-dcx*Wr) +
- (dcx*Wq-dcy*Wp)*(dcx*Wq-dcy*Wp) );
- if(Disc >0) truedist = TMath::Abs(Det/Disc);
-
- return truedist; // distance in cm
- }
- // =====================================================================
- Double_t TStraw::TimnsToDiscm(Double_t time) {
- // distance in cm from time in ns for pSTP =1 and pSTP=2 atm
- // utility routine for the track reconstruction
- //last update: A. Rotondi 3-3-2007
- Double_t drift, truer;
- // 1 absolute atm (NTP) 20-7-2006 GARFIELd Ar/CO2 90/10 MAGY
- // ns --> mm
-
- if(pSTP < 1.9) {
- if(Radius < 0.5){
- // 1600 V 4 mm variable theshold
- // drift = -5.68738e-02 +7.06505e-03*time
- // +3.70423e-03*pow(time,2) -1.39795e-04*pow(time,3)
- // +2.59863e-06*pow(time,4) -2.75729e-08*pow(time,5)
- // +1.69315e-10*pow(time,6) -5.59127e-13*pow(time,7)
- // +7.66371e-16*pow(time,8);
- // 1600 V 4 mm fixed threshold
- // drift = 7.50700e-02 -4.58303e-02*time
- // +6.01476e-03*pow(time,2) -1.78491e-04*pow(time,3)
- // +2.82377e-06*pow(time,4) -2.60709e-08*pow(time,5)
- // +1.40082e-10*pow(time,6) -4.04386e-13*pow(time,7)
- // +4.81504e-16*pow(time,8);
- drift = -0.070 -1.18253e-02
- -1.65315e-02 *time
- +6.18870e-03 *pow(time,2)
- -2.50647e-04 *pow(time,3)
- +5.15779e-06 *pow(time,4)
- -6.04826e-08 *pow(time,5)
- +4.07345e-10 *pow(time,6)
- -1.46437e-12 *pow(time,7)
- +2.17349e-15 *pow(time,8);
- }
- else{
- // 1600 V 5 mm
- drift = - 7.43091e-02
- + 6.34826e-03 *time
- + 2.70409e-03 *pow(time,2)
- - 7.96339e-05 *pow(time,3)
- + 1.12156e-06 *pow(time,4)
- - 8.76461e-09 *pow(time,5)
- + 3.87359e-11 *pow(time,6)
- - 9.05398e-14 *pow(time,7)
- + 8.69003e-17 *pow(time,8);
- }
- }
- else {
- // 2 absolute atm 5 mm 80/20 (1 bar overpressure)
- if(time < 100.) {
- drift = + 1.71849e-02
- - 3.99942e-02 *time
- + 8.96815e-03 *pow(time,2)
- - 3.23438e-04 *pow(time,3)
- + 4.88349e-06 *pow(time,4)
- - 2.27659e-08 *pow(time,5)
- - 1.00527e-10 *pow(time,6)
- - 6.18934e-13 *pow(time,7)
- + 2.69373e-14 *pow(time,8)
- - 1.26629e-16 *pow(time,9);
- }
- else{
- drift = - 2.65873e-01
- + 6.20514e-02 *time
- + 3.42564e-05 *pow(time,2)
- - 1.17871e-05 *pow(time,3)
- + 1.59485e-07 *pow(time,4)
- - 9.61794e-10 *pow(time,5)
- + 2.79989e-12 *pow(time,6)
- - 3.20470e-15 *pow(time,7);
- }
- // correction of the systematic shift for r=0.5 cm 2 atm
- // if( drift<=4.8){
- // truer= 2.36461e-01 -5.01598e-01*drift + 3.89468*pow(drift,2)
- // -4.17036*pow(drift,3) +2.19117*pow(drift,4)
- // -6.02374e-01*pow(drift,5) +8.32138e-02*pow(drift,6)
- // -4.55919e-03*pow(drift,7);
-
- // drift= TMath::Abs(truer);
- // }
- }
-
- drift = 0.1*drift;
- if(drift < 0.) drift=0.;
- return drift;
- }
- // --------------------------------------------------------------------------------
- Double_t TStraw::PartToTime(Double_t xPMass, Double_t xPMom, Double_t InOut[]) {
- // find the time of a particle of mass xPmass, momentum xPMom, with
- // input-output coordinate InOut[6]
- // Useful for MC pplication after a call to PutWireXYZ
- TInit(xPMass, xPMom, InOut); // start the event
- Out1 = StrawCharge(); // energy loss (GeV) to generate charge
- Out2 = StrawSignal(Nchann); // generate the straw signal
- Out3 = StrawTime(); // find the straw drift time PulseTime
-
- return PulseTime;
-
- }
- // --------------------------------------------------------------------------------
- Double_t TStraw::PartToADC() {
- // return the energy loss in the tube as charge signal
- // taking into account the Polya fluctuations
- Double_t ADCsignal=0.;
- for(Int_t j=1; j<= NNClus; j++){
-
- for(Int_t jc=1; jc<=(Int_t) CNeleT[j-1]; jc++)
- ADCsignal += bPolya * GasGain * PolyaSamp();
- }
- return ADCsignal;
- }
- // --------------------------------------------------------------------------------
- Double_t TStraw::DiffLong(Double_t Distcm) {
- // return the longitudinal diffusion
- // from cm to microns
- // MAGY GARFIELD Ar/CO2 90/10 1 bar (NTP) 20-7-2006
- //
- Double_t DiffMic;
- if(pSTP<1.9){
- if(Radius < 0.5){
- // 4 mm 1600 V 90/10
- DiffMic = 0.896 + 1387.*Distcm - 1.888e+04*pow(Distcm,2)
- + 1.799e+05*pow(Distcm,3) - 9.848e+05*pow(Distcm,4)
- + 3.009e+06*pow(Distcm,5) - 4.777e+06*pow(Distcm,6)
- + 3.074e+06*pow(Distcm,7);
- }
- else{
- // 5 mm 1600 V
- DiffMic = 1.537 + 1246.*Distcm - 1.357e+04*pow(Distcm,2)
- + 1.049e+05*pow(Distcm,3) - 4.755e+05*pow(Distcm,4)
- + 1.211e+06*pow(Distcm,5) - 1.6e+06*pow(Distcm,6)
- + 8.533e+05*pow(Distcm,7);
- }
- }
-
- else{
- // 2000V 5 mm 2 bar 80/20
- DiffMic = 2.135 +818.*Distcm - 1.044e+04*pow(Distcm,2)
- + 8.31e+04*pow(Distcm,3) - 3.492e+05*pow(Distcm,4)
- + 7.959e+05*pow(Distcm,5) - 9.378e+05*pow(Distcm,6)
- + 4.492e+05*pow(Distcm,7);
- }
- return DiffMic;
- }
- // --------------------------------------------------------------------------------
- Double_t TStraw::DiffTran(Double_t Distcm) {
- // return the transverse diffusion in microns
- // from cm to microns
- // MAGY GARFIELD Ar/CO2 90/10 1 bar (NTP) 20-7-2006
- //
- Double_t DiffMic;
- if(pSTP<1.9){
- if(Radius < 0.5){
- // 4 mm 1600 V 90/10
- // DiffMic = + 0.8513 + 1648.*Distcm - 1.085e+04*pow(Distcm,2)
- // + 7.38e+04*pow(Distcm,3) - 3.025e+05*pow(Distcm,4)
- // + 6.067e+05*pow(Distcm,5) - 4.643e+04*pow(Distcm,6);
- DiffMic = + 1.482 + 1529.*Distcm - 6755.*pow(Distcm,2)
- + 2.924e+04*pow(Distcm,3) - 0.9246e+05*pow(Distcm,4)
- + 1.548e+05*pow(Distcm,5) - 1.002e+05*pow(Distcm,6);
- }
- else{
- // 5 mm 1600 V 90/10
- DiffMic = + 1.482 + 1529.*Distcm - 6755.*pow(Distcm,2)
- + 2.924e+04*pow(Distcm,3) - 0.9246e+05*pow(Distcm,4)
- + 1.548e+05*pow(Distcm,5) - 1.002e+05*pow(Distcm,6);
- }
- }
- else{
- // 5 mm 2000 V 2 bar 80/20
- DiffMic = +2.094 + 1138.*Distcm - 7557.*pow(Distcm,2)
- + 2.968e+04*pow(Distcm,3) - 6.577e+04*pow(Distcm,4)
- + 7.581e+04*pow(Distcm,5) - 3.497e+04*pow(Distcm,6);
- }
- return DiffMic;
- }
- // --------------------------------------------------------
- Double_t TStraw::DistEle(Double_t time) {
- // dist in cm from time in ns for pSTP =1 and pSTP=2
- // utility routine for the >SINGLE< electron reconstruction
- //last update: A. Rotondi 20-7-2006
- Double_t drift;
- // 1 absolute atm (NTP) 20-7-2006 GARFIELd Ar/CO2 90/10 MAGY
- // ns --> cm
- time *= 0.001; // time in micro s
- if(pSTP < 2.) {
- if(Radius <0.5){
- // 1600 V 4 mm 90/10
- drift = 0.001629 + 6.194*time
- - 56.55*pow(time,2) + 355.8*pow(time,3)
- - 903.2*pow(time,4);
- }
- else{
- // 1600 V 5 mm 90/10
- drift = 0.003365 + 5.734*time
- - 41.88*pow(time,2) + 191.2*pow(time,3)
- - 333.4*pow(time,4) ;
- }
- }
- else {
-
- // 2 absolute atm 2000 V 5 mm 80/20 (1 bar overpressure)
- drift = 0.003365 + 7.056*time
- - 62.28*pow(time,2) + 306.1*pow(time,3)
- - 558.7*pow(time,4);
- }
- return drift; // in cm
- }
|