TStraw.cxx 31 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905906907908909910911912913914915916917918919920921922923924925926927928929930931932933934935936937938939940941942943944945946947948949950951952953954955956957958959960961962963964965966967968969970971972973974975976977978979980981982983984985986987988989990991992993994995996997998999100010011002100310041005100610071008100910101011101210131014101510161017101810191020102110221023102410251026102710281029103010311032103310341035103610371038103910401041104210431044104510461047104810491050105110521053105410551056105710581059106010611062106310641065106610671068106910701071107210731074107510761077107810791080108110821083108410851086108710881089109010911092109310941095109610971098109911001101110211031104110511061107110811091110
  1. //
  2. // ================================================================================
  3. // Straw tube response simulation
  4. //
  5. //
  6. //
  7. // mandatory call at the beginning
  8. // once to set constants
  9. //
  10. // void TConst(Double_t Radius, Double_t pSTP, Double_t ArP, Double_t CO2P);
  11. //
  12. // Pysical (garfield) values:
  13. // Radius Press (pSTP) % Ar (ArP) % CO2 (CO2P)
  14. // 0.4 1 0.9 0.1
  15. // 0.5 1 0.9 0.1
  16. // 0.5 2 0.8 0.2
  17. //
  18. //
  19. // initialization at each track for MC applications
  20. //
  21. // define particle and its in-out hits
  22. // void TInit(Double_t Mass, Double_t Momentum, Double_t InOut[])
  23. //
  24. // void PutWireXYZ(Double_t w1, Double_t w2, Double_t w3,
  25. // Double_t w4, Double_t w5, Double_t w6); define the wire
  26. //
  27. // Double_t PartToTime(Double_t Mass, Double_t Momentum,
  28. // Double_t InOut[] ); straw time
  29. //
  30. // Double_t PartToADC(); ADC charge
  31. // signal corresponding to the energy loss of StrawCharge
  32. // for energy loss simulation
  33. //
  34. // where:
  35. //
  36. // pSTP = pressure in bar
  37. // Radius = tube radius
  38. // ArP and CO2P percentages of Argon and CO2
  39. // w1:w6 = coordinates of the extremes of the straw wire (straw axis)
  40. // Mass, Momentm = mass and momentum of the MC particle
  41. // InOut[0:5] = straw input and output coordinates traversed from the MC particle
  42. //
  43. // Routine for off line reconstruction
  44. //
  45. // Double_t TimnsToDiscm(Double_t time); from time (ns) to
  46. // radius in cm
  47. // where:
  48. // time: measured drift time in ns
  49. //
  50. //
  51. //
  52. // For other applications see the routines listed in TStraw.h
  53. //
  54. // Author: A. Rotondi (november 2005, revised and extended in August 2006)
  55. //
  56. // ================================================================================
  57. //
  58. #include "TRandom.h"
  59. #include "TMath.h"
  60. #include "TStraw.h"
  61. #include "TImage.h"
  62. #include "TVector3.h"
  63. #include "TVectorD.h"
  64. #include "TMatrixD.h"
  65. #include "TMatrixDEigen.h"
  66. ClassImp(TStraw);
  67. // ============================================================
  68. TStraw::TStraw() {
  69. CNumb=0;
  70. memset(CumClus,0,sizeof(CumClus));
  71. memset(PolyaCum,0,sizeof(PolyaCum));
  72. memset(Xs,0,sizeof(Xs));
  73. // class constructor
  74. // clear
  75. CNumb=0;
  76. CDist.clear();
  77. CDistC.clear();
  78. CNele.clear();
  79. CNeleT.clear();
  80. TeleTime.clear();
  81. AmplSig.clear();
  82. Pulse.clear();
  83. PulseT.clear();
  84. WDist.clear();
  85. }
  86. // ----------------------------------------------------------------------
  87. void TStraw::TConst(Double_t R1, Double_t P1, Double_t A1, Double_t C1) {
  88. // set constants for the simulation
  89. // Input
  90. // P1 = tube absolute pressure pSTP
  91. // R1 = tube radius Radius
  92. // A1 = argon percentage ArPerc
  93. // C1 = C=2 percentage C=2Perc
  94. //
  95. Radius=R1;
  96. pSTP=P1;
  97. // Input for the media (volume percentages)
  98. ArPerc = A1;
  99. CO2Perc = C1;
  100. // cluster dimensions in Ar and CO2 (experimantal values)
  101. // Double_t PClus[20] =
  102. // {0., .656, .150, .064, .035, .0225, .0155, .0105, .0081,
  103. // .0061, .0049, .0039, .0030, .0025, .0020, .0016, .0012,
  104. // .00095, .00075, .00063}; // Fischle fit
  105. Double_t PClus[20] =
  106. {0., .802, .0707, .020, .013, .008, .006, .005, .006,
  107. .008, .009, .007, .0050, .0040, .0033, .0029, .0025,
  108. .0023, .0022, .002}; // Lapique 1st calculation
  109. // Double_t PClus[20] =
  110. // {0., .841, .0340, .021, .013, .008, .006, .004, .003,
  111. // .008, .013, .008, .0050, .004, .0030, .0028, .0025,
  112. // .0023, .0022, .002}; // Lapique 2nd calculation
  113. // Double_t PClus[20] =
  114. // {0., .656, .150, .064, .035, .0225, .0155, .0105, .0081,
  115. // .0061, .0049, .0039, .0030, .0025, .0020, .0016, .0012,
  116. // .00080, .00059, .00045}; // Fischle empirical
  117. // PDouble_t Clus[20] =
  118. // {0., .656, .148, .0649, .0337, .0244, .0141, .0078, .0095,
  119. // .0063, .0062, .0042, .0028, .0018, .0023, .0017, .0014,
  120. // .00060, .00050, .00063}; // Fischle exp
  121. Double_t CO2Clus[20] =
  122. {0., .730, .162, .038, .020, .0110, .0147, .0060, .0084,
  123. .0052, .0020, .0042, .0021, .0025, .0038, .0021, .0009,
  124. .00013, .00064, .00048}; // Fischle exp
  125. // Double_t CH4Clus[20] =
  126. // {0., .786, .120, .032, .013, .0098, .0055, .0057, .0027,
  127. // .0029, .0020, .0016, .0013, .0010, .0012, .0006, .0005,
  128. // .00042, .00037, .00033}; // Fischle exp
  129. CH4Perc = 0.07;
  130. // -----------------------------------------------------
  131. // gain of the avalanche
  132. // Ar/CO2 90/10 1 bar (NTP) MAGY GARFIELD 20-7-2006
  133. GasGain=90000.;
  134. // argon ----------------------------------------------------
  135. AAr = 39.948; // Argon (39.948)
  136. ZAr= 18.0; // Argon (18)
  137. RhoAr = pSTP*1.78*1.e-03; // g/cm3 (1.78 mg/cm3)
  138. IAr = 188*1.e-09; // ionization potential (GeV) (188 eV)
  139. WiAr =26.7; // energy to create an ion pair (standard 26.7 eV)
  140. NclAr = 23.; // cluster/cm in Argon
  141. // CO2 -----------------------------------------------------
  142. ACO2 = 44; // CO2
  143. ZCO2 = 22.; // CO2
  144. RhoCO2 = pSTP*1.98*1.e-03; // g/cm3 CO2 (1.98 mg/cm3)
  145. ICO2 = 95.8*1.e-09; // ionization potential (GeV) (188 eV)
  146. WiCO2 = 33.0; // energy to create an ion pair (33 eV)
  147. NclCO2 = 35.5; // clusters/cm CO2 35.5
  148. // Methane CH4 ---------------------------------------------------------
  149. ACH4 = 16; // CO2 (39.948)
  150. ZCH4 = 10.; // CO2 (18)
  151. RhoCH4 = pSTP*0.71*1.e-03; // g/cm3 CO2 (0.71 mg/cm3)
  152. ICH4 = 40.6*1.e-09; // ionization potential (GeV) (188 eV)
  153. WiCH4 = 28.0; // energy to create an ion pair
  154. NclCH4 = 25.0;
  155. // Input for the media (weight percentages) ----------------------------
  156. ArWPerc = ArPerc *AAr /(ArPerc*AAr + CO2Perc*ACO2);
  157. CO2WPerc = CO2Perc*ACO2/(ArPerc*AAr + CO2Perc*ACO2);
  158. // mixture densiies ----------------------------------------------------
  159. RhoMixCO2 = 1./((ArWPerc/RhoAr) + (CO2WPerc/RhoCO2));
  160. RhoMixCH4 = 1./((ArWPerc/RhoAr) + (CH4WPerc/RhoCH4));
  161. //----------------------------------------------------------------------
  162. // particles (Gev, energy losses in Kev)
  163. PZeta = 1; // projectile charge
  164. piMass = 0.139; // particle mass (GeV)
  165. eMass = 0.511/1000.; // electron mass (GeV) (0.511 MeV)
  166. prMass = 0.93827; // proton mass (GeV)
  167. // ---------------------------------------------------------------------
  168. // thresholds for the straw tubes (default values) see TInit for current values
  169. Thresh1=10;
  170. Thresh2=30;
  171. // channels for the signal
  172. Nchann = 500;
  173. // ---------------------------------------------Emin------------------------
  174. NPolya= 100; // steps for the calculation of the Polya distributions
  175. Xmax=0.; // Polya istribution is calculated between o and Xmax (see Polya)
  176. bPolya = 0.5;
  177. Polya(bPolya); // cumulative of the Polya distribution
  178. // -----------------------------------------------------------------------
  179. // cumulative for the number of electron per cluster
  180. Double_t Wnorm = (ArPerc*NclAr + CO2Perc*NclCO2);
  181. CumClus[0]=(ArPerc*NclAr*PClus[0] + CO2Perc*NclCO2*CO2Clus[0])/Wnorm;
  182. //for(Int_t i=1; i<=20; i++) --->> out of bounds was changed
  183. for(Int_t i=1; i<20; i++) {
  184. CumClus[i]=(ArPerc*NclAr*PClus[i] + CO2Perc*NclCO2*CO2Clus[i])/Wnorm
  185. + CumClus[i-1];
  186. }
  187. CumClus[20]=1.;
  188. Double_t sum=0.;
  189. for(Int_t i=0; i<=19; i++) {
  190. sum += PClus[i];
  191. //cout<<" PClus["<<i<<"] = "<<PClus[i]<<endl;
  192. }
  193. //for(Int_t i=0;i<=20;i++) {cout<<"CumClus["<<i<<"] = "<<CumClus[i]<<endl;}
  194. //cout<<" Sum of Probabilities = "<<sum<<endl;
  195. }
  196. // ==============================================================
  197. void TStraw::PutWireXYZ(Double_t w1, Double_t w2, Double_t w3,
  198. Double_t w4, Double_t w5, Double_t w6){
  199. // get wire coordinates
  200. Wx1=w1;
  201. Wy1=w2;
  202. Wz1=w3;
  203. Wx2=w4;
  204. Wy2=w5;
  205. Wz2=w6;
  206. }
  207. // =============================================================
  208. void TStraw::TInit(Double_t xPMass, Double_t xPMom, Double_t InOut[]) {
  209. // initialization of the constants for each track and each straw
  210. // transfer of data
  211. // track geometrical quantities
  212. Xin=InOut[0]; Yin=InOut[1]; Zin=InOut[2];
  213. Xout=InOut[3]; Yout=InOut[4]; Zout=InOut[5];
  214. PMass = xPMass;
  215. PMom = xPMom;
  216. // path into the straw
  217. Dx = sqrt((Xout-Xin)*(Xout-Xin) +
  218. (Yout-Yin)*(Yout-Yin) +
  219. (Zout-Zin)*(Zout-Zin));
  220. // clear
  221. CNumb=0;
  222. PulseMax=0;
  223. PulseTime=0;
  224. CDist.clear();
  225. CDistC.clear();
  226. CNele.clear();
  227. CNeleT.clear();
  228. WDist.clear();
  229. TeleTime.clear();
  230. AmplSig.clear();
  231. Pulse.clear();
  232. PulseT.clear();
  233. // ---------------------------------------------------------------------
  234. // initialization of the constants
  235. // maximum energy transfer Emax (GeV) and related quantities
  236. PEn = sqrt(PMom*PMom + PMass*PMass); // particle energy GeV
  237. beta = PMom/PEn;
  238. gamma= 1/sqrt(1.-beta*beta);
  239. Double_t begam = beta*gamma; // local
  240. Double_t mratio= eMass/PMass; // local
  241. // calculation of the polarization Sternheimer factor for Argon
  242. Double_t ms = 2.80;
  243. Double_t Xst = log10(begam*sqrt(pSTP));
  244. Double_t Cs =-11.92;
  245. Double_t X1 = 4;
  246. Double_t Xo = 1.96;
  247. Double_t as = 0.389;
  248. Delta=0;
  249. if(Xo <= Xst && Xst <= X1) Delta = 4.6051702*Xst +Cs +as*pow(X1-Xst,ms);
  250. else if(X1< Xst) Delta = 4.6051702*Xst +Cs;
  251. if(Delta < 0) Delta=0;
  252. // calculation of other typical quantites
  253. Emax = (1.022*begam*begam/(1.+2.*gamma*mratio+mratio*mratio))/1000.; // GeV
  254. CsiAr = pSTP*0.5*0.3071*PZeta*PZeta*ZAr*RhoAr/(beta*beta*AAr)*1e-03; //GeV
  255. Double_t fact = 2.*eMass*beta*beta*gamma*gamma*Emax/(IAr*IAr);
  256. EmedAr = 2.*CsiAr*(0.5*TMath::Log(fact)-beta*beta - 0.5*Delta); // GeV/cm
  257. EminAr = pSTP*2.70*1.e-06; // GeV/cm (2.5 keV)
  258. // most prob energy (GeV)
  259. EmpAr = EmedAr + CsiAr*(0.422784 + beta*beta + log(CsiAr/Emax));
  260. CsiCO2 = pSTP*0.5*0.3071*PZeta*PZeta*ZCO2*RhoCO2/(beta*beta*ACO2)*1e-03; //GeV
  261. fact = 2.*eMass*beta*beta*gamma*gamma*Emax/(ICO2*ICO2);
  262. EmedCO2 = 2.*CsiCO2*(0.5*TMath::Log(fact)-beta*beta - 0.5*Delta); // GeV/cm
  263. EminCO2 = pSTP*3.60*1.e-06; // GeV/cm (2.5 keV)
  264. // most prob energy (GeV)
  265. EmpCO2 = EmedCO2 + CsiCO2*(0.422784 + beta*beta + log(CsiCO2/Emax));
  266. Csi = RhoMixCO2 * ((ArWPerc*CsiAr/RhoAr) + (CO2WPerc*CsiCO2/RhoCO2));;
  267. Emed = RhoMixCO2 * ((ArWPerc*EmedAr/RhoAr) + (CO2WPerc*EmedCO2/RhoCO2));
  268. Emin = RhoMixCO2 * ((ArWPerc*EminAr/RhoAr) + (CO2WPerc*EminCO2/RhoCO2));
  269. Emp = RhoMixCO2 * ((ArWPerc*EmpAr/RhoAr) + (CO2WPerc*EmpCO2/RhoCO2));
  270. // mean weighted interaction
  271. Wi = (ArPerc*NclAr*WiAr + CO2Perc*NclCO2*WiCO2)
  272. /(ArPerc*NclAr + CO2Perc*NclCO2);
  273. Ncl=pSTP*((ArPerc*NclAr) + (CO2Perc*NclCO2))*Emed/Emin; // <--- Cluster/cm
  274. Lcl=1./Ncl; // mean free path between clusters (cm)
  275. Ecl = 2.8; // mean number of electrons per clusters
  276. Cutoff = Ncl*Ecl*20; // limit to the number of primary electrons
  277. Ntote = Ecl*Ncl; // total mean number of electrons per cm
  278. // ---------------------------------------------------------------------
  279. // thresholds for the straw tubes (current values)
  280. // total number of electrons * scaling factor (max of signal x electron=1)
  281. Thresh1=Ncl*Ecl* 0.05;
  282. Thresh2=Ncl*Ecl* 0.15;
  283. // -----------------------------------------------------------------------
  284. TDirCos(); // director cosine of the track
  285. // control prints
  286. // cout<<" Dx "<<Dx<<" Csi "<<Csi*1.e+09<<" Emax MeV "<<
  287. // Emax*1.e+03<<" PMass "<<PMass<<" gamma "<<
  288. // gamma<<endl;
  289. // cout<<" RRise "<<RRise(gamma)<<endl;
  290. // cout<<" Thresh1 = "<<Thresh1<<" Thresh2 = "<<Thresh2<<endl;
  291. }
  292. // =============================================================
  293. // energy loss in GeV
  294. Double_t TStraw::STEloss() {
  295. return gRandom->Landau(Emp*Dx,Csi*Dx); // in GeV
  296. }
  297. // =============================================================
  298. Double_t TStraw::StrawCharge() {
  299. // sampling of the ionization energy loss, number of
  300. // primary electrons and their wire distances for each track
  301. // energy loss in GeV
  302. // clear
  303. Double_t TotEnEle = 0., ClusEner=0., Ecurr;
  304. CNeleT.clear();
  305. WDist.clear();
  306. // threshold for delta rays (>1.5 KeV from NIM A301(1991)202)
  307. Double_t thr=Csi*Dx*1.e+09/(Ntote*1500.);
  308. // total energy released by the electrons
  309. //Number of cluster in the length
  310. NNClus = Cluster(); // set vectors CDist CDistC e CNele
  311. // calculation of total energy released by the electrons
  312. for(Int_t j=1;j<=NNClus;j++){
  313. ClusEner = 0.;
  314. for(Int_t jj=1; jj<=(Int_t)(CNele.at(j-1)); jj++){
  315. Ecurr = Wi;
  316. TotEnEle += Ecurr;
  317. ClusEner += Ecurr;
  318. // delta rays simulation over 1.keV. Upper threshol 12 keV
  319. if(gRandom->Uniform() <= thr) {
  320. TotEnEle += int(Wi/(1.002- gRandom->Uniform()));
  321. ClusEner += int(Wi/(1.002- gRandom->Uniform()));
  322. }
  323. }
  324. // effective number of electrons per cluster
  325. CNeleT.push_back(int(ClusEner/Wi));
  326. }
  327. // record the distance from the wire of any electron
  328. // adding the diffusion effects
  329. Int_t owfl =0;
  330. for(Int_t k=1; k<=NNClus; k++) {
  331. for (Int_t kk=0; kk<(Int_t)CNeleT.at(k-1); kk++) {
  332. if(owfl > (Int_t) Cutoff) break;
  333. owfl++;
  334. WDist.push_back(WDistCalc(CDistC.at(k-1)));
  335. }
  336. }
  337. return TotEnEle*1.e-09; // in Gev
  338. }
  339. // =============================================================
  340. Int_t TStraw::Cluster() {
  341. // calculate the number of clusters CNumb
  342. // their distance CDist end
  343. // their number of electrons CNele
  344. CNumb=0;
  345. CDist.clear();
  346. CDistC.clear();
  347. CNele.clear();
  348. Double_t DisTot = 0;
  349. for(Int_t k=1; k<=1000; k++) {
  350. // distance
  351. Double_t path = -Lcl*log(1-gRandom->Uniform()+0.000001);
  352. DisTot+= path;
  353. if(DisTot>Dx) break;
  354. CNumb=k;
  355. CDist.push_back(path);
  356. CDistC.push_back(DisTot);
  357. CNele.push_back((Double_t)Eject() );
  358. }
  359. return CNumb;
  360. }
  361. // =============================================================
  362. Int_t TStraw::Eject() {
  363. // find the number of electrons in a cluster
  364. Int_t Inelect;
  365. Double_t nelect=0.;
  366. nelect= TMath::BinarySearch(21,CumClus,gRandom->Uniform());
  367. if(nelect<19) {
  368. return Inelect = (Int_t) (nelect) +1;
  369. }
  370. else {
  371. return Inelect = (Int_t)(20./(1.01-gRandom->Uniform()));
  372. }
  373. }
  374. // =============================================================
  375. Double_t TStraw::RRise(Double_t gamma_r) {
  376. // interpolate the relativisic rise of the
  377. // number of cluster per cm, starting from the one
  378. // measured at the ionization minimum
  379. Double_t Rise;
  380. Double_t lg = log10(gamma_r);
  381. if(1.<=gamma_r && gamma_r <= 2.2){
  382. Rise = -2.159*lg +1.7;
  383. }
  384. else if(2.2<=gamma_r && gamma_r <= 6.){
  385. Rise = 1.;
  386. }
  387. else if(6.<=gamma_r && gamma_r <= 200.){
  388. Rise = 0.302*lg + 0.765;
  389. }
  390. else if(200.<=gamma_r && gamma_r <= 1000.){
  391. Rise = 0.1431*lg + 1.131;
  392. }
  393. else{
  394. Rise= 1.54;
  395. }
  396. return Rise;
  397. }
  398. // ==========================================================================
  399. void TStraw::Polya(Double_t bpar) {
  400. // calculate the cumulative of the Polya distribution
  401. // for samplig the gain fluctuations
  402. Double_t eps = 0.0001;
  403. Double_t Dxx = 0.05, xx=0., x1,x2;
  404. Double_t PMax;
  405. Double_t k=1./bpar -1.;
  406. // find Xmax
  407. PMax = eps * pow(k,k)*exp(-k);
  408. Double_t value = 1.e+06;
  409. Xmax =2*k;
  410. while(value > PMax){
  411. Xmax +=Dxx;
  412. value = pow(Xmax,k)*exp(-Xmax);
  413. }
  414. Xmax += -0.5*Dxx;
  415. // calculate the cumulative
  416. Double_t dx= Xmax/NPolya;
  417. Xs[0]=0;
  418. for(Int_t i=1; i<NPolya; i++){
  419. x1 = xx;
  420. xx += dx;
  421. x2= xx;
  422. Xs[i]=x2;
  423. if(i>1)
  424. PolyaCum[i] = PolyaCum[i-1] +
  425. 0.5*dx*(pow(x1,k)*exp(-x1) + pow(x2,k)*exp(-x2))/TMath::Gamma(k+1);
  426. else
  427. PolyaCum[i] = PolyaCum[i-1] +
  428. 0.5*dx* pow(x2,k)*exp(-x2)/TMath::Gamma(k+1);
  429. }
  430. // adjust the normalization
  431. for(Int_t ii=0; ii<NPolya; ii++) PolyaCum[ii] /= PolyaCum[NPolya-1];
  432. }
  433. // =====================================================================
  434. Double_t TStraw::PolyaSamp() {
  435. // sampling a wire gain fluctuation in the gas
  436. Double_t xr = gRandom->Uniform();
  437. Int_t n = TMath::BinarySearch(NPolya,PolyaCum,xr);
  438. Double_t xsamp = Xmax;
  439. if(n<NPolya-1){
  440. xsamp= Xs[n] +
  441. ((xr-PolyaCum[n])/(PolyaCum[n+1] - PolyaCum[n])) *(Xs[n+1]-Xs[n]);
  442. }
  443. return xsamp;
  444. }
  445. // =====================================================================
  446. TVector3 TStraw::WDistCalc(Double_t DisTot) {
  447. // calculation of the distance of the cluster from the wire
  448. // DisTot is the cluster distance from the track entrance
  449. // diffusion effects are considered
  450. // wire director cosines
  451. Double_t Wlength = sqrt( (Wx2-Wx1)*(Wx2-Wx1) +
  452. (Wy2-Wy1)*(Wy2-Wy1) +
  453. (Wz2-Wz1)*(Wz2-Wz1) );
  454. Wp = (Wx2-Wx1)/Wlength;
  455. Wq = (Wy2-Wy1)/Wlength;
  456. Wr = (Wz2-Wz1)/Wlength;
  457. // current track coordinates
  458. Double_t xcor = Calpha * DisTot + Xin;
  459. Double_t ycor = Cbeta * DisTot + Yin;
  460. Double_t zcor = Cgamma * DisTot + Zin;
  461. // cross product VV1= (wire x track) for the distance
  462. Double_t XX1 = Wr*(ycor-Wy1) - Wq*(zcor-Wz1);
  463. Double_t YY1 = Wp*(zcor-Wz1) - Wr*(xcor-Wx1);
  464. Double_t ZZ1 = Wq*(xcor-Wx1) - Wp*(ycor-Wy1);
  465. // vector for the 3-D distance from the wire (from wire x VV1)
  466. Double_t XX = Wq*ZZ1 - Wr*YY1;
  467. Double_t YY = Wr*XX1 - Wp*ZZ1;
  468. Double_t ZZ = Wp*YY1 - Wq*XX1;
  469. //cout<<" XYZ "<<XX<<" "<<YY<<" "<<ZZ<<endl;
  470. Double_t DDistcm = sqrt(XX*XX + YY*YY +ZZ*ZZ);
  471. //director cosines of the distance vector
  472. Double_t cosx = XX/DDistcm;
  473. Double_t cosy = YY/DDistcm;
  474. Double_t cosz = ZZ/DDistcm;
  475. // sampling of the diffusion
  476. //Double_t DDist = DDistcm*10.; // distance in mm
  477. // longitudinal coefficient of diffusion (GARFIELD) cm --> micron
  478. // MAGY Ar/CO2 90/10 20-7-2006 GARFIELD
  479. Double_t SigL = DiffLong(DDistcm);
  480. // ... in cm
  481. SigL *= 1.e-04; // in cm
  482. // tranverse coefficient (GARFIELD) cm--> micron
  483. // MAGY Ar/CO2 90/10 20-7-2006 GARFIELD
  484. Double_t SigT = DiffTran(DDistcm);
  485. SigT *= 1.e-04; // in cm
  486. // sampling of Longitudinal and Transverse diffusion
  487. Double_t difL = gRandom->Gaus(0.,SigL);
  488. Double_t difT = gRandom->Gaus(0.,SigT);
  489. // vector addition to the distance
  490. // the transverse component has the same dir cos of the wire
  491. XX += difL*cosx + difT*Wp;
  492. YY += difL*cosy + difT*Wq;
  493. ZZ += difL*cosz + difT*Wr;
  494. //cout<<" XYZ+ dif "<<XX<<" "<<YY<<" "<<ZZ<<endl;
  495. //cout<<" --------------------------------------------------- "<<endl;
  496. TVector3 TDist(XX,YY,ZZ);
  497. return TDist;
  498. }
  499. // =====================================================================
  500. void TStraw::TDirCos() {
  501. // director cosines of the track (called for each track)
  502. //path into the straw
  503. Rpath = sqrt((Xout-Xin)*(Xout-Xin) +
  504. (Yout-Yin)*(Yout-Yin) +
  505. (Zout-Zin)*(Zout-Zin));
  506. //director cosines
  507. Calpha = (Xout-Xin)/Rpath;
  508. Cbeta = (Yout-Yin)/Rpath;
  509. Cgamma = (Zout-Zin)/Rpath;
  510. }
  511. // =====================================================================
  512. Int_t TStraw::TimeEle() {
  513. TeleTime.clear();
  514. // calculate the arrival times (ns) for each electron in TeleTime
  515. // from cm to ns <---------------------------------
  516. // return the number of electrons
  517. // Author:A. Rotondi 20-7-2006
  518. Int_t ido = WDist.size();
  519. Double_t etime;
  520. // cutoff the total charge to 20 times the average
  521. if(ido > (Int_t) Cutoff) ido = (Int_t) Cutoff;
  522. for(Int_t k=0; k<ido; k++) {
  523. Double_t dst = WDist.at(k).Mag(); // distance in cm
  524. // space to time calculation from GARFIELD
  525. if(pSTP <1.9){
  526. if(Radius<0.5){
  527. // V=1600V diameter 4 mm 90/10
  528. etime= -1.624e-05 + 0.1258*dst + 0.8079*pow(dst,2) -2.918*pow(dst,3)
  529. + 10.33*pow(dst,4) -10.84*pow(dst,5);
  530. }
  531. else{
  532. // V=1600V diameter 5 mm 90/10
  533. etime= -6.763e-05 + 0.1471*dst +0.3625*pow(dst,2) +0.3876*pow(dst,3)
  534. +1.04*pow(dst,4) -1.693*pow(dst,5);
  535. }
  536. }
  537. else{
  538. // 2 bar 2000V, 5 mm, 80/20
  539. etime= -0.0001014 + 0.1463*dst -0.1694*pow(dst,2) +2.4248*pow(dst,3)
  540. -1.793*pow(dst,4);
  541. }
  542. etime*=1000.; // nano seconds
  543. TeleTime.push_back(etime);
  544. }
  545. return TeleTime.size();
  546. }
  547. // =====================================================================
  548. Double_t TStraw::Signal(Double_t t, Double_t t0) {
  549. // electric signal at time t of a cluster arriving at t0
  550. Double_t elesig = 0;
  551. Double_t x= t - t0;
  552. if (x>0){
  553. Double_t A = 1.03e-03;
  554. Double_t B = 3.95;
  555. Double_t C = 0.228;
  556. Double_t D = -3.839e-02;
  557. Double_t E = 1.148e-03;
  558. elesig = A*exp(B*log(x)-C*x)*(1+D*x+E*x*x);
  559. }
  560. return elesig;
  561. }
  562. // =====================================================================
  563. Int_t TStraw::StrawSignal(Int_t nsteps) {
  564. // creation of nstep values of
  565. // the straw global Pulse (sum on all clusters)
  566. // return the number of primary electrons
  567. PulseMax=0;
  568. Pulse.clear();
  569. PulseT.clear();
  570. AmplSig.clear();
  571. Int_t neltot = TimeEle(); // creation and size of TeleTime (electron times)
  572. Double_t Tmax = 1.e-25;
  573. for(Int_t k=0; k< neltot; k++) if(Tmax<TeleTime.at(k)) Tmax=TeleTime.at(k);
  574. Tmax += 100.;
  575. Double_t Dt = Tmax/nsteps; // number of steps of the signal
  576. //AmplSig is the amplitude of each electron
  577. for(Int_t j=0; j< neltot; j++){
  578. AmplSig.push_back(bPolya*PolyaSamp());
  579. }
  580. // creation of the signal Pulse(PulseT) time in ns
  581. for(Int_t j=0; j< nsteps; j++){
  582. Double_t te = j*Dt;
  583. Double_t sumele=0.;
  584. for(Int_t jj=0; jj< neltot; jj++){
  585. Double_t te0 = TeleTime.at(jj);
  586. sumele += AmplSig.at(jj)*Signal(te,te0);
  587. }
  588. Pulse.push_back(sumele);
  589. PulseT.push_back(te);
  590. }
  591. // add a random noise (3% of maximum) plus
  592. // a 1% of 200 ns periodic signal
  593. Double_t Pmax = 1.e-25;
  594. for(Int_t k=0; k<(Int_t)Pulse.size(); k++)
  595. if(Pmax<Pulse.at(k)) Pmax=Pulse.at(k);
  596. for(Int_t k=0; k<(Int_t)Pulse.size(); k++){
  597. Pulse.at(k) += 0.03*Pmax*gRandom->Uniform()
  598. + 0.01*Pmax*(sin(6.28*PulseT.at(k)/120.));
  599. // if(Pulse[k]<0) Pulse[k] *= -1.;
  600. }
  601. PulseMax=Pmax;
  602. // set variable threshold for the signals
  603. Thresh1=0.05*Pmax;
  604. Thresh2=0.15*Pmax;
  605. return neltot;
  606. }
  607. // =====================================================================
  608. Int_t TStraw::StrawTime() {
  609. // simulate the discrimination of the straw signal
  610. // and give the time
  611. // discriminator technique: set 2 threshold, select the first one
  612. // (the low one) only if the second one is fired (FINUDA system)
  613. PulseTime=0.;
  614. Int_t ind=0;
  615. Int_t flag1=0;
  616. Int_t flag2=0;
  617. for(Int_t k=0; k<(Int_t) Pulse.size(); k++){
  618. if(flag1==0 && Pulse[k]>Thresh1) {
  619. flag1=1;
  620. ind=k;
  621. }
  622. if(Pulse[k]<Thresh1) {flag1=0; ind=0;} // reset if signal decreases
  623. if(flag1==1 && Pulse[k]>Thresh2) {
  624. flag2=1;
  625. break;
  626. }
  627. }
  628. if(flag1==1 && flag2==1) PulseTime=PulseT.at(ind);
  629. return ind;
  630. }
  631. // =====================================================================
  632. Double_t TStraw::TrueDist(Double_t Point[]) {
  633. // service routine that finds the distance in cm from the wire
  634. // by knowing the wire coordinates (class variables)
  635. // and the input-output points Point[6]
  636. Double_t truedist = 0;
  637. // wire director cosines
  638. Double_t Wlength = sqrt( (Wx2-Wx1)*(Wx2-Wx1) +
  639. (Wy2-Wy1)*(Wy2-Wy1) +
  640. (Wz2-Wz1)*(Wz2-Wz1) );
  641. Wp = (Wx2-Wx1)/Wlength;
  642. Wq = (Wy2-Wy1)/Wlength;
  643. Wr = (Wz2-Wz1)/Wlength;
  644. // director cosines of the given track
  645. Double_t Modu = sqrt( (Point[3]*Point[0])*(Point[3]*Point[0]) +
  646. (Point[4]*Point[1])*(Point[4]*Point[1]) +
  647. (Point[5]*Point[2])*(Point[5]*Point[2]) );
  648. Double_t dcx = (Point[3]-Point[0])/Modu;
  649. Double_t dcy = (Point[4]-Point[1])/Modu;
  650. Double_t dcz = (Point[5]-Point[2])/Modu;
  651. //distance formula
  652. Double_t p1 = (Point[0]-Wx1)*(dcy*Wr-dcz*Wq);
  653. Double_t p2 =-(Point[1]-Wy1)*(dcx*Wr-dcz*Wp);
  654. Double_t p3 = (Point[2]-Wz1)*(dcx*Wq-dcy*Wp);
  655. Double_t Det = p1+p2+p3;
  656. Double_t Disc = sqrt( (dcy*Wr-dcz*Wq)*(dcy*Wr-dcz*Wq) +
  657. (dcz*Wp-dcx*Wr)*(dcz*Wp-dcx*Wr) +
  658. (dcx*Wq-dcy*Wp)*(dcx*Wq-dcy*Wp) );
  659. if(Disc >0) truedist = TMath::Abs(Det/Disc);
  660. return truedist; // distance in cm
  661. }
  662. // =====================================================================
  663. Double_t TStraw::TimnsToDiscm(Double_t time) {
  664. // distance in cm from time in ns for pSTP =1 and pSTP=2 atm
  665. // utility routine for the track reconstruction
  666. //last update: A. Rotondi 3-3-2007
  667. Double_t drift, truer;
  668. // 1 absolute atm (NTP) 20-7-2006 GARFIELd Ar/CO2 90/10 MAGY
  669. // ns --> mm
  670. if(pSTP < 1.9) {
  671. if(Radius < 0.5){
  672. // 1600 V 4 mm variable theshold
  673. // drift = -5.68738e-02 +7.06505e-03*time
  674. // +3.70423e-03*pow(time,2) -1.39795e-04*pow(time,3)
  675. // +2.59863e-06*pow(time,4) -2.75729e-08*pow(time,5)
  676. // +1.69315e-10*pow(time,6) -5.59127e-13*pow(time,7)
  677. // +7.66371e-16*pow(time,8);
  678. // 1600 V 4 mm fixed threshold
  679. // drift = 7.50700e-02 -4.58303e-02*time
  680. // +6.01476e-03*pow(time,2) -1.78491e-04*pow(time,3)
  681. // +2.82377e-06*pow(time,4) -2.60709e-08*pow(time,5)
  682. // +1.40082e-10*pow(time,6) -4.04386e-13*pow(time,7)
  683. // +4.81504e-16*pow(time,8);
  684. drift = -0.070 -1.18253e-02
  685. -1.65315e-02 *time
  686. +6.18870e-03 *pow(time,2)
  687. -2.50647e-04 *pow(time,3)
  688. +5.15779e-06 *pow(time,4)
  689. -6.04826e-08 *pow(time,5)
  690. +4.07345e-10 *pow(time,6)
  691. -1.46437e-12 *pow(time,7)
  692. +2.17349e-15 *pow(time,8);
  693. }
  694. else{
  695. // 1600 V 5 mm
  696. drift = - 7.43091e-02
  697. + 6.34826e-03 *time
  698. + 2.70409e-03 *pow(time,2)
  699. - 7.96339e-05 *pow(time,3)
  700. + 1.12156e-06 *pow(time,4)
  701. - 8.76461e-09 *pow(time,5)
  702. + 3.87359e-11 *pow(time,6)
  703. - 9.05398e-14 *pow(time,7)
  704. + 8.69003e-17 *pow(time,8);
  705. }
  706. }
  707. else {
  708. // 2 absolute atm 5 mm 80/20 (1 bar overpressure)
  709. if(time < 100.) {
  710. drift = + 1.71849e-02
  711. - 3.99942e-02 *time
  712. + 8.96815e-03 *pow(time,2)
  713. - 3.23438e-04 *pow(time,3)
  714. + 4.88349e-06 *pow(time,4)
  715. - 2.27659e-08 *pow(time,5)
  716. - 1.00527e-10 *pow(time,6)
  717. - 6.18934e-13 *pow(time,7)
  718. + 2.69373e-14 *pow(time,8)
  719. - 1.26629e-16 *pow(time,9);
  720. }
  721. else{
  722. drift = - 2.65873e-01
  723. + 6.20514e-02 *time
  724. + 3.42564e-05 *pow(time,2)
  725. - 1.17871e-05 *pow(time,3)
  726. + 1.59485e-07 *pow(time,4)
  727. - 9.61794e-10 *pow(time,5)
  728. + 2.79989e-12 *pow(time,6)
  729. - 3.20470e-15 *pow(time,7);
  730. }
  731. // correction of the systematic shift for r=0.5 cm 2 atm
  732. // if( drift<=4.8){
  733. // truer= 2.36461e-01 -5.01598e-01*drift + 3.89468*pow(drift,2)
  734. // -4.17036*pow(drift,3) +2.19117*pow(drift,4)
  735. // -6.02374e-01*pow(drift,5) +8.32138e-02*pow(drift,6)
  736. // -4.55919e-03*pow(drift,7);
  737. // drift= TMath::Abs(truer);
  738. // }
  739. }
  740. drift = 0.1*drift;
  741. if(drift < 0.) drift=0.;
  742. return drift;
  743. }
  744. // --------------------------------------------------------------------------------
  745. Double_t TStraw::PartToTime(Double_t xPMass, Double_t xPMom, Double_t InOut[]) {
  746. // find the time of a particle of mass xPmass, momentum xPMom, with
  747. // input-output coordinate InOut[6]
  748. // Useful for MC pplication after a call to PutWireXYZ
  749. TInit(xPMass, xPMom, InOut); // start the event
  750. Out1 = StrawCharge(); // energy loss (GeV) to generate charge
  751. Out2 = StrawSignal(Nchann); // generate the straw signal
  752. Out3 = StrawTime(); // find the straw drift time PulseTime
  753. return PulseTime;
  754. }
  755. // --------------------------------------------------------------------------------
  756. Double_t TStraw::PartToADC() {
  757. // return the energy loss in the tube as charge signal
  758. // taking into account the Polya fluctuations
  759. Double_t ADCsignal=0.;
  760. for(Int_t j=1; j<= NNClus; j++){
  761. for(Int_t jc=1; jc<=(Int_t) CNeleT[j-1]; jc++)
  762. ADCsignal += bPolya * GasGain * PolyaSamp();
  763. }
  764. return ADCsignal;
  765. }
  766. // --------------------------------------------------------------------------------
  767. Double_t TStraw::DiffLong(Double_t Distcm) {
  768. // return the longitudinal diffusion
  769. // from cm to microns
  770. // MAGY GARFIELD Ar/CO2 90/10 1 bar (NTP) 20-7-2006
  771. //
  772. Double_t DiffMic;
  773. if(pSTP<1.9){
  774. if(Radius < 0.5){
  775. // 4 mm 1600 V 90/10
  776. DiffMic = 0.896 + 1387.*Distcm - 1.888e+04*pow(Distcm,2)
  777. + 1.799e+05*pow(Distcm,3) - 9.848e+05*pow(Distcm,4)
  778. + 3.009e+06*pow(Distcm,5) - 4.777e+06*pow(Distcm,6)
  779. + 3.074e+06*pow(Distcm,7);
  780. }
  781. else{
  782. // 5 mm 1600 V
  783. DiffMic = 1.537 + 1246.*Distcm - 1.357e+04*pow(Distcm,2)
  784. + 1.049e+05*pow(Distcm,3) - 4.755e+05*pow(Distcm,4)
  785. + 1.211e+06*pow(Distcm,5) - 1.6e+06*pow(Distcm,6)
  786. + 8.533e+05*pow(Distcm,7);
  787. }
  788. }
  789. else{
  790. // 2000V 5 mm 2 bar 80/20
  791. DiffMic = 2.135 +818.*Distcm - 1.044e+04*pow(Distcm,2)
  792. + 8.31e+04*pow(Distcm,3) - 3.492e+05*pow(Distcm,4)
  793. + 7.959e+05*pow(Distcm,5) - 9.378e+05*pow(Distcm,6)
  794. + 4.492e+05*pow(Distcm,7);
  795. }
  796. return DiffMic;
  797. }
  798. // --------------------------------------------------------------------------------
  799. Double_t TStraw::DiffTran(Double_t Distcm) {
  800. // return the transverse diffusion in microns
  801. // from cm to microns
  802. // MAGY GARFIELD Ar/CO2 90/10 1 bar (NTP) 20-7-2006
  803. //
  804. Double_t DiffMic;
  805. if(pSTP<1.9){
  806. if(Radius < 0.5){
  807. // 4 mm 1600 V 90/10
  808. // DiffMic = + 0.8513 + 1648.*Distcm - 1.085e+04*pow(Distcm,2)
  809. // + 7.38e+04*pow(Distcm,3) - 3.025e+05*pow(Distcm,4)
  810. // + 6.067e+05*pow(Distcm,5) - 4.643e+04*pow(Distcm,6);
  811. DiffMic = + 1.482 + 1529.*Distcm - 6755.*pow(Distcm,2)
  812. + 2.924e+04*pow(Distcm,3) - 0.9246e+05*pow(Distcm,4)
  813. + 1.548e+05*pow(Distcm,5) - 1.002e+05*pow(Distcm,6);
  814. }
  815. else{
  816. // 5 mm 1600 V 90/10
  817. DiffMic = + 1.482 + 1529.*Distcm - 6755.*pow(Distcm,2)
  818. + 2.924e+04*pow(Distcm,3) - 0.9246e+05*pow(Distcm,4)
  819. + 1.548e+05*pow(Distcm,5) - 1.002e+05*pow(Distcm,6);
  820. }
  821. }
  822. else{
  823. // 5 mm 2000 V 2 bar 80/20
  824. DiffMic = +2.094 + 1138.*Distcm - 7557.*pow(Distcm,2)
  825. + 2.968e+04*pow(Distcm,3) - 6.577e+04*pow(Distcm,4)
  826. + 7.581e+04*pow(Distcm,5) - 3.497e+04*pow(Distcm,6);
  827. }
  828. return DiffMic;
  829. }
  830. // --------------------------------------------------------
  831. Double_t TStraw::DistEle(Double_t time) {
  832. // dist in cm from time in ns for pSTP =1 and pSTP=2
  833. // utility routine for the >SINGLE< electron reconstruction
  834. //last update: A. Rotondi 20-7-2006
  835. Double_t drift;
  836. // 1 absolute atm (NTP) 20-7-2006 GARFIELd Ar/CO2 90/10 MAGY
  837. // ns --> cm
  838. time *= 0.001; // time in micro s
  839. if(pSTP < 2.) {
  840. if(Radius <0.5){
  841. // 1600 V 4 mm 90/10
  842. drift = 0.001629 + 6.194*time
  843. - 56.55*pow(time,2) + 355.8*pow(time,3)
  844. - 903.2*pow(time,4);
  845. }
  846. else{
  847. // 1600 V 5 mm 90/10
  848. drift = 0.003365 + 5.734*time
  849. - 41.88*pow(time,2) + 191.2*pow(time,3)
  850. - 333.4*pow(time,4) ;
  851. }
  852. }
  853. else {
  854. // 2 absolute atm 2000 V 5 mm 80/20 (1 bar overpressure)
  855. drift = 0.003365 + 7.056*time
  856. - 62.28*pow(time,2) + 306.1*pow(time,3)
  857. - 558.7*pow(time,4);
  858. }
  859. return drift; // in cm
  860. }