StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StTpcRSMaker.cxx
1 // doxygen info here
4 /*
5  The maker's algorithms and formulae based on
6  http://www.inst.bnl.gov/programs/gasnobledet/publications/Mathieson's_Book.pdf,
7  and Photo Absorption Model
8  "Ionization energy loss in very thin absorbers.", V.M. Grishin, V.K. Ermilova, S.K. Kotelnikov Nucl.Instrum.Meth.A309:476-484,1991
9  "A method to improve tracking and particle identification in TPC's and silicon detectors.", Hans Bichsel, Nucl.Instrum.Meth.A562:154-197,2006
10  HEED: "Modeling of ionization produced by fast charged particles in gases", I.B. Smirnov, Nucl.Instrum.Meth.A55(2005) 747-493.
11 
12 */
13 #include <assert.h>
14 #include "StTpcRSMaker.h"
15 #include "Stiostream.h"
16 // SCL
17 #include "StGlobals.hh"
18 #include "StThreeVectorD.hh"
19 #include "StPhysicalHelixD.hh"
20 // ROOT
21 #include "TClassTable.h"
22 #include "TDataSetIter.h"
23 #include "TTableSorter.h"
24 #include "TRandom.h"
25 #include "TSystem.h"
26 #include "TFile.h"
27 #include "TBenchmark.h"
28 #include "TProfile2D.h"
29 #include "TH3.h"
30 #include "TVirtualMC.h"
31 #include "TInterpreter.h"
32 #include "Math/SpecFuncMathMore.h"
33 #include "StDbUtilities/StCoordinates.hh"
34 #include "StDbUtilities/StTpcCoordinateTransform.hh"
35 // Dave's Header file
36 #include "StDbUtilities/StMagUtilities.h"
37 //#include "StDaqLib/TPC/trans_table.hh"
38 #include "StDetectorDbMaker/St_tpcAltroParamsC.h"
39 #include "StDetectorDbMaker/St_asic_thresholdsC.h"
40 #include "StDetectorDbMaker/St_tss_tssparC.h"
41 #include "StDetectorDbMaker/St_tpcPadGainT0BC.h"
42 #include "StDetectorDbMaker/St_TpcResponseSimulatorC.h"
43 #include "StDetectorDbMaker/St_tpcAnodeHVavgC.h"
44 #include "StDetectorDbMaker/StDetectorDbTpcRDOMasks.h"
45 #include "StDetectorDbMaker/St_tpcPadConfigC.h"
46 #include "StDetectorDbMaker/St_tpcGainCorrectionC.h"
47 #include "StDetectorDbMaker/St_TpcAvgCurrentC.h"
48 #include "StDetectorDbMaker/St_TpcAvgPowerSupplyC.h"
49 #include "StDetectorDbMaker/St_trigDetSumsC.h"
50 #include "StDetectorDbMaker/St_tpcBXT0CorrEPDC.h"
51 #include "StDetectorDbMaker/St_TpcPadPedRMSC.h"
52 #include "StEventUtilities/StEbyET0.h"
53 #include "StDetectorDbMaker/St_beamInfoC.h"
54 #include "StDetectorDbMaker/St_GatingGridBC.h"
55 #include "StDetectorDbMaker/St_starTriggerDelayC.h"
56 #if 0
57 #include "StParticleTable.hh"
58 #include "StParticleDefinition.hh"
59 #endif
60 #include "Altro.h"
61 #include "TRVector.h"
62 #include "StBichsel/Bichsel.h"
63 #include "StBichsel/StdEdxModel.h"
64 #include "StdEdxY2Maker/StTpcdEdxCorrection.h"
65 // g2t tables
66 #include "tables/St_g2t_track_Table.h"
67 #include "tables/St_g2t_vertex_Table.h"
68 #include "tables/St_g2t_tpc_hit_Table.h"
69 #include "StEventTypes.h"
70 struct HitPoint_t {
71  Int_t indx;
72  Int_t TrackId;
73  Double_t s; // track length to current point
74  Double_t sMin, sMax;
75  g2t_tpc_hit_st *tpc_hitC;
76  StGlobalCoordinate xyzG;
77  StTpcLocalSectorCoordinate coorLS;
78  StTpcLocalSectorDirection dirLS, BLS;
79  StTpcPadCoordinate Pad;
80 };
81 //#define ElectronHack
82 //#define Old_dNdx_Table
83 //#define __ELECTRONS_TUPLE__
84 #ifdef __TFG__VERSION__
85 #define __DEBUG__
86 #define __CHECK_RDOMAP_AND_VOLTAGE__
87 #endif /* __TFG__VERSION__ */
88 #if defined(__DEBUG__)
89 #define PrPP(A,B) if (Debug()%10 > 2) {LOG_INFO << "StTpcRSMaker::" << (#A) << "\t" << (#B) << " = \t" << (B) << endm;}
90 #else
91 #define PrPP(A,B)
92 #endif
93 static const char rcsid[] = "$Id: StTpcRSMaker.cxx,v 1.92 2020/05/22 20:49:19 fisyak Exp $";
94 #ifndef __DEBUG__
95 static Bool_t ClusterProfile = kFALSE;
96 #else
97 static Bool_t ClusterProfile = kTRUE;
98 #endif
99 #define Laserino 170
100 #define Chasrino 171
101 // Inner Outer
102 static Double_t t0IO[2] = {1.20868e-9, 1.43615e-9}; // recalculated in InducedCharge
103 static const Double_t tauC[2] = {999.655e-9, 919.183e-9};
104 TF1* StTpcRSMaker::fgTimeShape3[2] = {0, 0};
105 TF1* StTpcRSMaker::fgTimeShape0[2] = {0, 0};
106 static Double_t fgTriggerT0 = 0;
107 static Double_t timeBinMin = -0.5;
108 static Double_t timeBinMax = 44.5;
109 //________________________________________________________________________________
110 static const Int_t nx[4] = {200,500, 145, 401};
111 static const Double_t xmin[4] = {-10., -6, -0.5, -0.5};
112 static const Double_t xmax[4] = { 10., 44,144.5, 400.5};
113 static const Int_t nz = 42;
114 static const Double_t zmin = -210;
115 static const Double_t zmax = -zmin;
116 // io pt
117 static TProfile2D *hist[4][2] = {0};
118 static TProfile2D *PixelHist[2] = {0};
119 static TProfile2D *GainHist[2] = {0};
120 static const Int_t nChecks = 22;
121 static TH1 *checkList[2][22] = {0};
122 static TString TpcMedium("TPCE_SENSITIVE_GAS");
123 Short_t StTpcRSMaker::mADCs[__MaxNumberOfTimeBins__];
124 //________________________________________________________________________________
125 ClassImp(StTpcRSMaker);
126 //________________________________________________________________________________
127 StTpcRSMaker::StTpcRSMaker(const char *name):
128  StMaker(name),
129  mLaserScale(1),
130  minSignal(1e-4),
131  ElectronRange(0.0055), // Electron Range(.055mm)
132  ElectronRangeEnergy(3000), // eV
133  ElectronRangePower(1.78), // sigma = ElectronRange*(eEnery/ElectronRangeEnergy)**ElectronRangePower
134  NoOfSectors(24),
135  NoOfPads(182),
136  NoOfTimeBins(__MaxNumberOfTimeBins__),
137  mCutEle(1e-5)
138 {
139  memset(beg, 0, end-beg+1);
140  m_Mode = 0;
141  // SETBIT(m_Mode,kHEED);
142  SETBIT(m_Mode,kBICHSEL); // Default is Bichsel
143  SETBIT(m_Mode,kdEdxCorr);
144  SETBIT(m_Mode,kDistortion);
145 }
146 //________________________________________________________________________________
147 StTpcRSMaker::~StTpcRSMaker() {
148  Finish();
149 }
150 //________________________________________________________________________________
151 Int_t StTpcRSMaker::Finish() {
152  // SafeDelete(fTree);
153  if (m_SignalSum) {free(m_SignalSum); m_SignalSum = 0;}
154  for (Int_t io = 0; io < 2; io++) {// Inner/Outer
155  for (Int_t sec = 0; sec < NoOfSectors; sec++) {
156  if (mShaperResponses[io][sec] && !mShaperResponses[io][sec]->TestBit(kNotDeleted)) {SafeDelete(mShaperResponses[io][sec]);}
157  if (mChargeFraction[io][sec] && !mChargeFraction[io][sec]->TestBit(kNotDeleted)) {SafeDelete(mChargeFraction[io][sec]);}
158  if (mPadResponseFunction[io][sec] && !mPadResponseFunction[io][sec]->TestBit(kNotDeleted)) {SafeDelete(mPadResponseFunction[io][sec]);}
159  if (mAltro[io][sec] && !mAltro[io][sec]->TestBit(kNotDeleted)) {SafeDelete(mAltro[io][sec]);}
160  }
161  SafeDelete(mPolya[io]);
162  }
163  if (m_TpcdEdxCorrection && m_TpcdEdxCorrection->TestBit(kCanDelete)) delete m_TpcdEdxCorrection;
164  m_TpcdEdxCorrection = 0;
165  return StMaker::Finish();
166 }
167 //________________________________________________________________________________
168 Int_t StTpcRSMaker::InitRun(Int_t /* runnumber */) {
169  SetAttr("minSector",1);
170  SetAttr("maxSector",24);
171  SetAttr("minRow",1);
172  SetAttr("maxRow",St_tpcPadConfigC::instance()->numberOfRows(20));
173  if (!gStTpcDb) {
174  LOG_ERROR << "Database Missing! Can't initialize TpcRS" << endm;
175  return kStFatal;
176  }
177  mCutEle = GetCutEle();
178  if (mCutEle > 0) {
179  LOG_INFO << "StTpcRSMaker::InitRun: mCutEle set to = " << mCutEle << " from geant \"" << TpcMedium.Data() << "\" parameters" << endm;
180  } else {
181  mCutEle = 1e-4;
182  LOG_ERROR << "StTpcRSMaker::InitRun: mCutEle has not been found in GEANT3 for \"" << TpcMedium.Data() << "\" parameters."
183  << "Probably due to missing Set it to default " << mCutEle << endm;
184  }
185  // Distortions
186  if (TESTBIT(m_Mode,kdEdxCorr)) {
187  LOG_INFO << "StTpcRSMaker:: use Tpc dE/dx correction from calibaration" << endm;
188  Long_t Mask = -1; // 64 bits
189  CLRBIT(Mask,StTpcdEdxCorrection::kAdcCorrection);
190  CLRBIT(Mask,StTpcdEdxCorrection::kAdcCorrectionC);
191  CLRBIT(Mask,StTpcdEdxCorrection::kEdge);
192  CLRBIT(Mask,StTpcdEdxCorrection::kAdcCorrectionMDF);
193  CLRBIT(Mask,StTpcdEdxCorrection::kAdcCorrection3MDF);
194  CLRBIT(Mask,StTpcdEdxCorrection::kAdcCorrection4MDF);
195  CLRBIT(Mask,StTpcdEdxCorrection::kAdcCorrection5MDF);
196  CLRBIT(Mask,StTpcdEdxCorrection::kAdcCorrection6MDF);
197  CLRBIT(Mask,StTpcdEdxCorrection::kEtaCorrection);
198  CLRBIT(Mask,StTpcdEdxCorrection::knTbk);
199  CLRBIT(Mask,StTpcdEdxCorrection::kzCorrectionC);
200  CLRBIT(Mask,StTpcdEdxCorrection::kzCorrection);
201 #if 0
202  CLRBIT(Mask,StTpcdEdxCorrection::kdXCorrection);
203  CLRBIT(Mask,StTpcdEdxCorrection::kTpcPadTBins);
204  CLRBIT(Mask,StTpcdEdxCorrection::kTanL);
205  CLRBIT(Mask,StTpcdEdxCorrection::kAdcI);
206  CLRBIT(Mask,StTpcdEdxCorrection::knPad);
207  CLRBIT(Mask,StTpcdEdxCorrection::kdZdY);
208  CLRBIT(Mask,StTpcdEdxCorrection::kdXdY);
209 #endif
210  m_TpcdEdxCorrection = new StTpcdEdxCorrection(Mask, Debug());
211  m_TpcdEdxCorrection->SetSimulation();
212  }
213  if (TESTBIT(m_Mode,kDistortion)) {
214  LOG_INFO << "StTpcRSMaker:: use Tpc distortion correction" << endm;
215  }
216  Double_t samplingFrequency = 1.e6*gStTpcDb->Electronics()->samplingFrequency(); // Hz
217  Double_t TimeBinWidth = 1./samplingFrequency;
218  /*
219 select firstInnerSectorAnodeWire,lastInnerSectorAnodeWire,numInnerSectorAnodeWires,firstOuterSectorAnodeWire,lastOuterSectorAnodeWire,numOuterSectorAnodeWires from Geometry_tpc.tpcWirePlanes;
220 +---------------------------+--------------------------+--------------------------+---------------------------+--------------------------+--------------------------+
221 | firstInnerSectorAnodeWire | lastInnerSectorAnodeWire | numInnerSectorAnodeWires | firstOuterSectorAnodeWire | lastOuterSectorAnodeWire | numOuterSectorAnodeWires |
222 +---------------------------+--------------------------+--------------------------+---------------------------+--------------------------+--------------------------+
223 | 53.2000000000 | 120.8000000000 | 170 | 122.7950000000 | 191.1950000000 | 172 |
224 +---------------------------+--------------------------+--------------------------+---------------------------+--------------------------+--------------------------+
225  */
226  numberOfInnerSectorAnodeWires = gStTpcDb->WirePlaneGeometry()->numberOfInnerSectorAnodeWires ();
227  firstInnerSectorAnodeWire = gStTpcDb->WirePlaneGeometry()->firstInnerSectorAnodeWire();
228  lastInnerSectorAnodeWire = gStTpcDb->WirePlaneGeometry()->lastInnerSectorAnodeWire ();
229  numberOfOuterSectorAnodeWires = gStTpcDb->WirePlaneGeometry()->numberOfOuterSectorAnodeWires ();
230  firstOuterSectorAnodeWire = gStTpcDb->WirePlaneGeometry()->firstOuterSectorAnodeWire();
231  lastOuterSectorAnodeWire = gStTpcDb->WirePlaneGeometry()->lastOuterSectorAnodeWire ();
232  anodeWirePitch = gStTpcDb->WirePlaneGeometry()->anodeWirePitch ();
233  anodeWireRadius = gStTpcDb->WirePlaneGeometry()->anodeWireRadius();
234  if (St_tpcPadConfigC::instance()->iTPC(1)) { // iTpc for all TPC sectors
235  NoOfPads = St_tpcPadConfigC::instance()->numberOfPadsAtRow(1,72);
236  }
237  Float_t BFieldG[3];
238  Float_t xyz[3] = {0,0,0};
239  StMagF::Agufld(xyz,BFieldG);
240  // Shapers
241  const Char_t *Names[2] = {"I","O"};
242  Double_t CathodeAnodeGap[2] = {0.2, 0.4};
243  for (Int_t sector = 1; sector <= 24; sector++) {
244  innerSectorAnodeVoltage[sector-1] = outerSectorAnodeVoltage[sector-1] = 0;
245  Int_t nAliveInner = 0;
246  Int_t nAliveOuter = 0;
247  for (Int_t row = 1; row <= St_tpcPadConfigC::instance()->numberOfRows(sector); row++) {
248  if (St_tpcPadConfigC::instance()->IsRowInner(sector,row)) {
249  nAliveInner++;
250  innerSectorAnodeVoltage[sector-1] += St_tpcAnodeHVavgC::instance()->voltagePadrow(sector,row);
251  } else {
252  nAliveOuter++;
253  outerSectorAnodeVoltage[sector-1] += St_tpcAnodeHVavgC::instance()->voltagePadrow(sector,row);
254  }
255  }
256  if (! nAliveInner && ! nAliveOuter) {
257  LOG_INFO << "Illegal date/time. Tpc sector " << sector << " Anode Voltage is not set to run condition: AliveInner: " << nAliveInner
258  << "\tAliveOuter: " << nAliveOuter
259  << "\tStop the run" << endm;
260  assert(nAliveInner || nAliveOuter);
261  } else {
262  if (nAliveInner > 1) innerSectorAnodeVoltage[sector-1] /= nAliveInner;
263  if (nAliveOuter > 1) outerSectorAnodeVoltage[sector-1] /= nAliveOuter;
264  }
265  if (GetTFile()) GetTFile()->cd();
266 #ifdef __CHECK_RDOMAP_AND_VOLTAGE__
267  static TH3F *AlivePads = 0;
268  if (! AlivePads) {
269  Int_t nrows = St_tpcPadConfigC::instance()->numberOfRows(20);
270  AlivePads = new TH3F("AlivePads","Active pads from RDO map, tpcGainPadT0, and Tpc Anode Voltage:sector:row:pad",24,0.5,24.5,nrows,0.5,nrows+.5,NoOfPads,0.5,NoOfPads+0.5);
271  }
272  for(Int_t row = 1; row <= St_tpcPadConfigC::instance()->numberOfRows(sector); row++) {
273  Int_t noOfPadsAtRow = St_tpcPadConfigC::instance()->numberOfPadsAtRow(sector,row);
274  if ( ! St_tpcAnodeHVavgC::instance()->livePadrow(sector,row)) continue;
275  for(Int_t pad = 1; pad<=noOfPadsAtRow; pad++) {
276  Int_t iRdo = StDetectorDbTpcRDOMasks::instance()->rdoForPadrow(sector,row,pad);
277  if ( ! StDetectorDbTpcRDOMasks::instance()->isOn(sector,iRdo)) {
278  continue;
279  } else {
280  Double_t gain = St_tpcPadGainT0BC::instance()->Gain(sector,row,pad);
281  if (gain <= 0.0) {
282  continue;
283  } else {
284  AlivePads->Fill(sector, row, pad, gain);
285  }
286  }
287  }
288  }
289 #endif /* __CHECK_RDOMAP_AND_VOLTAGE__ */
290  for (Int_t io = 0; io < 2; io++) {// In/Out
291  if (io == 0) {
292  if (sector > 1 && TMath::Abs(innerSectorAnodeVoltage[sector-1] - innerSectorAnodeVoltage[sector-2]) < 1) {
293  InnerAlphaVariation[sector-1] = InnerAlphaVariation[sector-2];
294  } else {
295  LOG_INFO << "Inner Sector " << sector << " ======================" << endm;
296  InnerAlphaVariation[sector-1] = InducedCharge(anodeWirePitch,
297  CathodeAnodeGap[io],
298  anodeWireRadius,
299  innerSectorAnodeVoltage[sector-1], t0IO[io]);
300  }
301  }
302  else {
303  if (sector > 1 && TMath::Abs(outerSectorAnodeVoltage[sector-1] - outerSectorAnodeVoltage[sector-2]) < 1) {
304  OuterAlphaVariation[sector-1] = OuterAlphaVariation[sector-2];
305  } else {
306  LOG_INFO << "Outer Sector " << sector << " ======================" << endm;
307  OuterAlphaVariation[sector-1] = InducedCharge(anodeWirePitch,
308  CathodeAnodeGap[io],
309  anodeWireRadius,
310  outerSectorAnodeVoltage[sector-1], t0IO[io]);
311  }
312  }
313  } // IO
314  } // sector
315 
316  for (Int_t io = 0; io < 2; io++) {// In/Out
317  // mPolya = new TF1F("Polya;x = G/G_0;signal","sqrt(x)/exp(1.5*x)",0,10); // original Polya
318  // mPolya = new TF1F("Polya;x = G/G_0;signal","pow(x,0.38)*exp(-1.38*x)",0,10); // Valeri Cherniatin
319  // mPoly = new TH1D("Poly","polyaAvalanche",100,0,10);
320  // TF1F *func = new TF1F("funcP","x*sqrt(x)/exp(2.5*x)",0,10);
321  // see http://www4.rcf.bnl.gov/~lebedev/tec/polya.html
322  // Gain fluctuation in proportional counters follows Polya distribution.
323  // x = G/G_0
324  // P(m) = m(m(x)**(m-1)*exp(-m*x)/Gamma(m);
325  // original Polya m = 1.5 (R.Bellazzini and M.A.Spezziga, INFN PI/AE-94/02).
326  // Valeri Cherniatin (cherniat@bnlarm.bnl.gov) recomends m=1.38
327  // Trs uses x**1.5/exp(x)
328  // tss used x**0.5/exp(1.5*x)
329  Double_t gamma;
330  if (!io ) gamma = St_TpcResponseSimulatorC::instance()->PolyaInner();
331  else gamma = St_TpcResponseSimulatorC::instance()->PolyaOuter();
332  if (gamma <= 0) gamma = 1.38;
333  mPolya[io] = new TF1F(io == 0 ? "PolyaInner;x = G/G_0;signal" : "PolyaOuter;x = G/G_0;signal",polya,0,10,3);
334  mPolya[io]->SetParameters(gamma, 0., 1./gamma);
335  Double_t params3[9] = {t0IO[io],
336  St_TpcResponseSimulatorC::instance()->tauF(),
337  St_TpcResponseSimulatorC::instance()->tauP(),
338  St_TpcResponseSimulatorC::instance()->tauIntegration(),
339  TimeBinWidth, 0, (Double_t ) io, 1.,
340  St_TpcResponseSimulatorC::instance()->tMax()[io]
341  };
342  Double_t params0[7] = {t0IO[io], St_TpcResponseSimulatorC::instance()->tauX()[io], TimeBinWidth, St_TpcResponseSimulatorC::instance()->tauC()[io], (Double_t ) io, 1.,
343  St_TpcResponseSimulatorC::instance()->tMax()[io]
344  };
345  if (! fgTimeShape3[io]) {// old electronics, intergation + shaper alltogether
346  fgTimeShape3[io] = new TF1(Form("TimeShape3%s",Names[io]),
347  shapeEI3,timeBinMin*TimeBinWidth,timeBinMax*TimeBinWidth,9);
348  fgTimeShape3[io]->SetParNames("t0","tauF","tauP", "tauI","width","tauC","io","norm","tMax");
349  params3[7] = 1.0;
350  fgTimeShape3[io]->SetParameters(params3);
351  params3[7] = fgTimeShape3[io]->Integral(timeBinMin*TimeBinWidth,timeBinMax*TimeBinWidth);
352  fgTimeShape3[io]->SetTitle(fgTimeShape3[io]->GetName());
353  fgTimeShape3[io]->GetXaxis()->SetTitle("time (secs)");
354  fgTimeShape3[io]->GetYaxis()->SetTitle("signal");
355  fgTimeShape3[io]->SetNpx(200);
356  fgTimeShape3[io]->Save(timeBinMin*TimeBinWidth,timeBinMax*TimeBinWidth, 0,0,0,0);
357  if (GetTFile()) fgTimeShape3[io]->Write();
358  }
359  if (! fgTimeShape0[io]) {// new electronics only integration
360  fgTimeShape0[io] = new TF1(Form("TimeShape0%s",Names[io]),
361  shapeEI,timeBinMin*TimeBinWidth,timeBinMax*TimeBinWidth,7);
362  fgTimeShape0[io]->SetParNames("t0","tauI","width","tauC","io","norm","tMax");
363  params0[5] = 1.0;
364  fgTimeShape0[io]->SetParameters(params0);
365  params0[5] = fgTimeShape0[io]->Integral(0,timeBinMax*TimeBinWidth);
366  fgTimeShape0[io]->SetTitle(fgTimeShape0[io]->GetName());
367  fgTimeShape0[io]->GetXaxis()->SetTitle("time (secs)");
368  fgTimeShape0[io]->GetYaxis()->SetTitle("signal");
369  fgTimeShape3[io]->SetNpx(200);
370  fgTimeShape0[io]->Save(timeBinMin*TimeBinWidth,timeBinMax*TimeBinWidth, 0,0,0,0);
371  if (GetTFile()) fgTimeShape0[io]->Write();
372  }
373  // w h s a l i
374  // Double_t paramsI[6] = {0.2850, 0.2000, 0.4000, 0.0010, 1.1500, 0};
375  // Double_t paramsO[6] = {0.6200, 0.4000, 0.4000, 0.0010, 1.1500, 0};
376  Double_t xmaxP = 4.5;//4.5*St_tpcPadConfigC::instance()->innerSectorPadWidth(sector);// 4.5
377  Double_t xminP = -xmaxP;
378  Double_t paramsPad[6];
379  for (Int_t sector = 1; sector <= NoOfSectors; sector++) {
380  if (! mPadResponseFunction[io][sector-1]) {
381  if (! io) {
382  paramsPad[0] = St_tpcPadConfigC::instance()->innerSectorPadWidth(sector); // w = width of pad
383  paramsPad[1] = gStTpcDb->WirePlaneGeometry()->innerSectorAnodeWirePadPlaneSeparation(); // h = Anode-Cathode gap
384  paramsPad[2] = gStTpcDb->WirePlaneGeometry()->anodeWirePitch(); // s = wire spacing
385  paramsPad[3] = St_TpcResponseSimulatorC::instance()->K3IP();
386  paramsPad[4] = 0;
387  paramsPad[5] = St_tpcPadConfigC::instance()->innerSectorPadPitch(sector);
388  } else {
389  paramsPad[0] = St_tpcPadConfigC::instance()->outerSectorPadWidth(sector); // w = width of pad
390  paramsPad[1] = gStTpcDb->WirePlaneGeometry()->outerSectorAnodeWirePadPlaneSeparation();// h = Anode-Cathode gap
391  paramsPad[2] = gStTpcDb->WirePlaneGeometry()->anodeWirePitch(); // s = wire spacing
392  paramsPad[3] = St_TpcResponseSimulatorC::instance()->K3OP();
393  paramsPad[4] = 0;
394  paramsPad[5] = St_tpcPadConfigC::instance()->outerSectorPadPitch(sector);
395  }
396  xmaxP = 4.5;//4.5*St_tpcPadConfigC::instance()->innerSectorPadWidth(sector);// 4.5
397  xminP = -xmaxP;
398  for (Int_t sec = 1; sec < sector; sec++) {
399  if (mPadResponseFunction[io][sec-1]) {
400  if (! memcmp(paramsPad, mPadResponseFunction[io][sec-1]->GetParameters(),sizeof(paramsPad))) {
401  mPadResponseFunction[io][sector-1] = mPadResponseFunction[io][sec-1];
402  break;
403  }
404  }
405  }
406  }
407  if (! mPadResponseFunction[io][sector-1]) {
408  mPadResponseFunction[io][sector-1] = new TF1F(Form("%s_%02i",io == 0 ? "PadResponseFunctionInner" : "PadResponseFunctionOuter",sector),StTpcRSMaker::PadResponseFunc,xminP,xmaxP,6);
409  mPadResponseFunction[io][sector-1]->SetParameters(paramsPad);
410  mPadResponseFunction[io][sector-1]->SetParNames("PadWidth","Anode-Cathode gap","wire spacing","K3OP","CrossTalk","PadPitch");
411  mPadResponseFunction[io][sector-1]->SetTitle(mPadResponseFunction[io][sector-1]->GetName());
412  mPadResponseFunction[io][sector-1]->GetXaxis()->SetTitle("pads");
413  mPadResponseFunction[io][sector-1]->GetYaxis()->SetTitle("Signal");
414  // Cut tails
415  Double_t x = 2.5;//cm xmaxP;
416 #if 0
417  Double_t ymax = mPadResponseFunction[io][sector-1]->Eval(0);
418  for (; x > 1.5; x -= 0.05) {
419  Double_t r = mPadResponseFunction[io][sector-1]->Eval(x)/ymax;
420  if (r > 1e-2) break;
421  }
422 #endif
423  mPadResponseFunction[io][sector-1]->SetRange(-x,x);
424  mPadResponseFunction[io][sector-1]->Save(xminP,xmaxP,0,0,0,0);
425  if (GetTFile()) mPadResponseFunction[io][sector-1]->Write();
426  }
427  // check if the function has been created
428  Double_t paramsRow[6] = {0};
429  if (! mChargeFraction[io][sector-1]) {
430  xmaxP = 2.5;//5*St_tpcPadConfigC::instance()->innerSectorPadLength(sector); // 1.42
431  xminP = - xmaxP;
432  if (! io) {
433  paramsRow[0] = St_tpcPadConfigC::instance()->innerSectorPadLength(sector);
434  paramsRow[1] = gStTpcDb->WirePlaneGeometry()->innerSectorAnodeWirePadPlaneSeparation(); // h = Anode-Cathode gap
435  paramsRow[2] = gStTpcDb->WirePlaneGeometry()->anodeWirePitch(); // s = wire spacing
436  paramsRow[3] = St_TpcResponseSimulatorC::instance()->K3IR();
437  paramsRow[4] = 0;
438  paramsRow[5] = 1.;
439  } else {
440  paramsRow[0] = St_tpcPadConfigC::instance()->outerSectorPadLength(sector);
441  paramsRow[1] = gStTpcDb->WirePlaneGeometry()->innerSectorAnodeWirePadPlaneSeparation(); // h = Anode-Cathode gap
442  paramsRow[2] = gStTpcDb->WirePlaneGeometry()->anodeWirePitch(); // s = wire spacing
443  paramsRow[3] = St_TpcResponseSimulatorC::instance()->K3OR();
444  paramsRow[4] = 0;
445  paramsRow[5] = 1.;
446  }
447  for (Int_t sec = 1; sec < sector; sec++) {
448  if (mChargeFraction[io][sec-1]) {
449  if (! memcmp(paramsRow, mChargeFraction[io][sec-1]->GetParameters(),sizeof(paramsRow))) {
450  mChargeFraction[io][sector-1] = mChargeFraction[io][sec-1];
451  break;
452  }
453  }
454  }
455  }
456  if (! mChargeFraction[io][sector-1]) {
457  mChargeFraction[io][sector-1] = new TF1F(Form("%s_%02i", io == 0 ? "ChargeFractionInner" : "ChargeFractionOuter", sector),
458  StTpcRSMaker::PadResponseFunc,xminP,xmaxP,6);
459  mChargeFraction[io][sector-1]->SetParameters(paramsRow);
460  mChargeFraction[io][sector-1]->SetParNames("PadLength","Anode-Cathode gap","wire spacing","K3IR","CrossTalk","RowPitch");
461  mChargeFraction[io][sector-1]->SetTitle(mChargeFraction[io][sector-1]->GetName());
462  mChargeFraction[io][sector-1]->GetXaxis()->SetTitle("Distance (cm)");
463  mChargeFraction[io][sector-1]->GetYaxis()->SetTitle("Signal");
464  // Cut tails
465  Double_t x = xmaxP;
466  Double_t ymax = mChargeFraction[io][sector-1]->Eval(0);
467  for (; x > 1.5; x -= 0.05) {
468  Double_t r = mChargeFraction[io][sector-1]->Eval(x)/ymax;
469  if (r > 1e-2) break;
470  }
471  mChargeFraction[io][sector-1]->SetRange(-x,x);
472  mChargeFraction[io][sector-1]->Save(xminP,xmaxP,0,0,0,0);
473  if (GetTFile()) mChargeFraction[io][sector-1]->Write();
474  }
475 #if 0
476  memset(mLocalYDirectionCoupling[io][sector-1], 0, sizeof(mLocalYDirectionCoupling[io][sector-1]));
477  for (Int_t j = 0; j < 7; j++) {
478  mLocalYDirectionCoupling[io][sector-1][j] = mChargeFraction[io][sector-1]->Eval(anodeWirePitch*j);
479  }
480 #endif
481  Int_t l = 0;
482  if (St_tpcAltroParamsC::instance()->Table()->GetNRows() > sector) l = sector - 1;
483  if (io == 0 && St_tpcAltroParamsC::instance()->Table()->GetNRows() > 24 + sector) l += 24;
484 #if 1
485  // Check that Shaper has initialized before
486  tpcAltroParams_st *Ssector = St_tpcAltroParamsC::instance()->Struct(l);
487  for (Int_t sec = 1; sec < sector; sec++) {
488  Int_t lc = 0;
489  if (St_tpcAltroParamsC::instance()->Table()->GetNRows() > sec) lc = sec - 1;
490  if (io == 0 && St_tpcAltroParamsC::instance()->Table()->GetNRows() > 24 + sector) lc += 24;
491  if (mShaperResponses[io][sec-1]) {
492  tpcAltroParams_st *Ssec = St_tpcAltroParamsC::instance()->Struct(lc);
493  if (! memcmp(Ssector,Ssec,sizeof(tpcAltroParams_st))) {
494  mShaperResponses[io][sector-1] = mShaperResponses[io][sec-1];
495  mAltro[io][sector-1] = mAltro[io][sec-1];
496  break;
497  }
498  }
499  }
500  if (mShaperResponses[io][sector-1]) continue;
501 #endif
502  if (St_tpcAltroParamsC::instance()->N(l) < 0) {// old TPC
503  mShaperResponses[io][sector-1] = new TF1F(Form("ShaperFunc_%s_S%02i",Names[io],sector),
504  StTpcRSMaker::shapeEI3_I,timeBinMin,timeBinMax,9);
505  mShaperResponses[io][sector-1]->SetParameters(params3);
506  mShaperResponses[io][sector-1]->SetParNames("t0","tauF","tauP", "tauI", "width","tauC","io","norm");
507  mShaperResponses[io][sector-1]->SetTitle(mShaperResponses[io][sector-1]->GetName());
508  mShaperResponses[io][sector-1]->GetXaxis()->SetTitle("time (buckets)");
509  mShaperResponses[io][sector-1]->GetYaxis()->SetTitle("signal");
510  assert(! mAltro[io][sector-1]);
511  } else if (St_tpcAltroParamsC::instance()->N(l) >= 0) {// Altro/Sampa
512  mShaperResponses[io][sector-1] = new TF1F(Form("ShaperFunc_%s_S%02i",Names[io],sector),
513  StTpcRSMaker::shapeEI_I,timeBinMin,timeBinMax,7);
514  mShaperResponses[io][sector-1]->SetParameters(params0);
515  mShaperResponses[io][sector-1]->SetParNames("t0","tauI", "width","tauC","io","norm");
516  mShaperResponses[io][sector-1]->SetTitle(mShaperResponses[io][sector-1]->GetName());
517  mShaperResponses[io][sector-1]->GetXaxis()->SetTitle("time (buckets)");
518  mShaperResponses[io][sector-1]->GetYaxis()->SetTitle("signal");
519  // Altro/Sampa
520  mAltro[io][sector-1] = new Altro(__MaxNumberOfTimeBins__,mADCs);
521  if (St_tpcAltroParamsC::instance()->N(l) == 0) {// no tail cancellation
522  mAltro[io][sector-1] = new Altro(__MaxNumberOfTimeBins__,mADCs);
523  mAltro[io][sector-1]->ConfigAltro(0,0,0,1,1);
524  } else { // Altro/Sampa with shaping parameters
525  // ConfigAltro(ONBaselineCorrection1, ONTailcancellation, ONBaselineCorrection2, ONClipping, ONZerosuppression)
526  mAltro[io][sector-1]->ConfigAltro( 0, 1, 0, 1, 1);
527  // ConfigBaselineCorrection_1(int mode, int ValuePeDestal, int *PedestalMem, int polarity)
528  //altro->ConfigBaselineCorrection_1(4, 0, PedestalMem, 0); // Tonko 06/25/08
529  mAltro[io][sector-1]->ConfigTailCancellationFilter(St_tpcAltroParamsC::instance()->K1(l),
530  St_tpcAltroParamsC::instance()->K2(l),
531  St_tpcAltroParamsC::instance()->K3(l), // K1-3
532  St_tpcAltroParamsC::instance()->L1(l),
533  St_tpcAltroParamsC::instance()->L2(l),
534  St_tpcAltroParamsC::instance()->L3(l));// L1-3
535  }
536  if (mAltro[io][sector-1]) {
537  mAltro[io][sector-1]->ConfigZerosuppression(St_tpcAltroParamsC::instance()->Threshold(l),
538  St_tpcAltroParamsC::instance()->MinSamplesaboveThreshold(l),
539  0,0);
540  mAltro[io][sector-1]->PrintParameters();
541  }
542  }
543  // Cut tails
544  Double_t t = timeBinMax;
545  Double_t ymax = mShaperResponses[io][sector-1]->Eval(0.5);
546  for (; t > 5; t -= 1) {
547  Double_t r = mShaperResponses[io][sector-1]->Eval(t)/ymax;
548  if (r > 1e-2) break;
549  }
550  mShaperResponses[io][sector-1]->SetRange(timeBinMin,t);
551  mShaperResponses[io][sector-1]->Save(timeBinMin,t,0,0,0,0);
552  if (GetTFile()) mShaperResponses[io][sector-1]->Write();
553  } // sector
554  } // io
555  if (Debug()) Print();
556  if (ClusterProfile && GetTFile()) GetTFile()->cd();
557 #if 0
558  StMagUtilities::SetDoDistortionT(gFile);
559  StMagUtilities::SetUnDoDistortionT(gFile);
560 #endif
561  mHeed = fEc(St_TpcResponseSimulatorC::instance()->W());
562  if ( ClusterProfile) {
563  Int_t color = 1;
564  struct Name_t {
565  const Char_t *Name;
566  const Char_t *Title;
567  };
568  const Name_t InOut[6] = {
569  {"Inner","Inner old electronics or iTPC"},
570  {"Outer","Outer old electronics or ITPC"},
571  {"InnerX","Inner new electronics"},
572  {"OuterX","Outer new electronics"},
573  {"I","Inner"},
574  {"O","Outer"}
575  };
576  const Name_t PadTime[4] = {
577  {"Pad","Pad - <PadMC> ; Z"},
578  {"Time","Time - <TimeMC> ; Z"},
579  {"PixelPad","Pad ; row ; pad"},
580  {"PixelTime","Time ; row ; time "},
581  };
582  for (Int_t io = 2; io < 4; io++) {
583  for (Int_t pt = 0; pt < 2; pt++) {
584  TString Name(InOut[io].Name); Name += PadTime[pt].Name; Name += "Mc";
585  TString Title(InOut[io].Title); Title += PadTime[pt].Title; Title += " Mc";
586  hist[io][pt] = (TProfile2D *) gDirectory->Get(Name);
587  if (! hist[io][pt]) {
588  hist[io][pt] = new TProfile2D(Name,Title,nx[pt],xmin[pt],xmax[pt],nz,zmin,zmax,"");
589  hist[io][pt]->SetMarkerStyle(20);
590  hist[io][pt]->SetMarkerColor(color++);
591  }
592  }
593  }
594  GainHist[0] = new TProfile2D("dEdxCorSecRow","dEdx correction ; sector ; row",
595  NoOfSectors,0.5,NoOfSectors+0.5,
596  St_tpcPadConfigC::instance()->numberOfRows(20),0.5,St_tpcPadConfigC::instance()->numberOfRows(20)+0.5,"");
597  GainHist[1] = new TProfile2D("GainSecRow","Overall gain ; sector ; row",
598  NoOfSectors,0.5,NoOfSectors+0.5,
599  St_tpcPadConfigC::instance()->numberOfRows(20),0.5,St_tpcPadConfigC::instance()->numberOfRows(20)+0.5,"");
600  for (Int_t pt = 0; pt < 2; pt++) {
601  TString Name(PadTime[pt+2].Name);
602  TString Title(PadTime[pt+2].Title); Title += " Mc";
603  PixelHist[pt] = new TProfile2D(Name,Title, 72, 0.5, 72.5, nx[pt+2],xmin[pt+2],xmax[pt+2],"");
604  }
605  const Name_t Checks[22] = {
606  {"dEGeant","dE in Geant"}, // 0
607  {"dSGeant","ds in Geant"}, // 1
608  {"Gain","Gas Gain after Voltage"}, // 2
609  {"GainMc","Gas Gain after MC correction"}, // 3
610  {"dEdxCor","correction of dEdx"}, // 4
611  {"lgam","lgam"}, // 5
612  {"NPGEANT","no. of primary electros from GEANT"}, // 6
613  {"NP","no. of primary electros"}, // 7
614  {"Nt","total no. of electors per cluster"}, // 8
615  {"Qav","Gas gain flactuations"}, // 9
616  {"localYDirectionCoupling","localYDirectionCoupling"}, //10
617  {"n0","No. electrons per primary interaction"}, //11
618  {"padGain","padGain"}, // 12
619  {"localXDirectionCoupling","localXDirectionCoupling"}, // 13
620  {"XYcoupling","XYcoupling"}, //14
621  {"dE","dE"}, // 15
622  {"dS","dS"}, // 16
623  {"adc","adc"},// 17
624  {"NE","Total no. of generated electors"}, // 18
625  {"dECl","Total log(signal/Nt) in a cluster ; Wire Index"}, // 19
626  {"nPdT","log(Total no. of conducting electrons) - log(no. of primary one) ; log(no. primary electrons)"}, // 20
627  {"bgVsbg","log10(bg_from_gkin) ; log10(bg_from_mom"} // 21
628  };
629  const Int_t Npbins = 151;
630  const Int_t NpbinsL = 10;
631  const Double_t Xmax = 1e5;
632  Double_t dX = TMath::Log(Xmax/10)/(Npbins - NpbinsL);
633  Double_t *pbins = new Double_t[Npbins];
634  Double_t *pbinsL = new Double_t[Npbins];
635  pbins[0] = 0.5;
636  pbinsL[0] = TMath::Log(pbins[0]);
637  for (Int_t bin = 1; bin < Npbins; bin++) {
638  if (bin <= NpbinsL) {
639  pbins[bin] = pbins[bin-1] + 1;
640  } else if (bin == Npbins - 1) {
641  pbins[bin] = 1e5;
642  } else {
643  Int_t nM = 0.5*(pbins[NpbinsL-2] + pbins[NpbinsL-1])*TMath::Exp(dX*(bin-NpbinsL));
644  Double_t dbin = TMath::Nint(nM - pbins[bin-1]);
645  if (dbin < 1.0) dbin = 1.0;
646  pbins[bin] = pbins[bin-1] + dbin;
647  }
648  pbinsL[bin] = TMath::Log(pbins[bin]);
649  }
650  for (Int_t io = 0; io < 2; io++) {
651  for (Int_t i = 0; i < nChecks; i++) {
652  TString Name(Checks[i].Name); Name += InOut[4+io].Name;
653  TString Title(Checks[i].Title); Title += InOut[4+io].Title;
654  if (i == 11) checkList[io][i] = new TH2D(Name,Title,nz,zmin,zmax,100,-0.5,99.5);
655  else if (i == 19) checkList[io][i] = new TH2D(Name,Title,173,-.5,172.5,200,-10,10);
656  // else if (i == 20) checkList[io][i] = new TH2D(Name,Title,Npbins-1,pbinsL,Npbins-1,pbinsL);
657  else if (i == 20) checkList[io][i] = new TH2D(Name,Title,Npbins-1,pbinsL,500,-2.0,8.0);
658  else if (i == 21) checkList[io][i] = new TH2D(Name,Title,100,-5,5,100,-5,5);
659  else checkList[io][i] = new TProfile(Name,Title,nz,zmin,zmax,"");
660  }
661  }
662  delete [] pbins;
663  delete [] pbinsL;
664  }
665  return kStOK;
666 }
667 //________________________________________________________________________________
668 Int_t StTpcRSMaker::Make(){ // PrintInfo();
669  static Int_t minSector = IAttr("minSector");
670  static Int_t maxSector = IAttr("maxSector");
671  static Int_t minRow = IAttr("minRow");
672  static Int_t maxRow = IAttr("maxRow");
673  TArrayF rs(1000);
674  // constants
675  static Int_t iBreak = 0;
676 #ifdef __DEBUG__
677  if (Debug()%10) {
678  gBenchmark->Reset();
679  gBenchmark->Start("TpcRS");
680  LOG_INFO << "\n -- Begin TpcRS Processing -- \n";
681  }
682 #endif
683  Double_t vminI = St_tpcGainCorrectionC::instance()->Struct(1)->min;
684  Double_t vminO = St_tpcGainCorrectionC::instance()->Struct(0)->min;
685  St_g2t_tpc_hit *g2t_tpc_hit = (St_g2t_tpc_hit *) GetDataSet("geant/g2t_tpc_hit");
686  if (!g2t_tpc_hit) return kStWarn;
687  Int_t no_tpc_hits = g2t_tpc_hit->GetNRows(); if (no_tpc_hits<1) return kStOK;
688  if (Debug() > 1) g2t_tpc_hit->Print(0,10);
689  St_g2t_track *g2t_track = (St_g2t_track *) GetDataSet("geant/g2t_track"); // if (!g2t_track) return kStWarn;
690  Int_t NoTpcTracks = 0;
691  if (g2t_track) NoTpcTracks = g2t_track->GetNRows();
692  mNoTpcHitsAll = TArrayI(NoTpcTracks+1);
693  mNoTpcHitsReal = TArrayI(NoTpcTracks+1);
694  g2t_track_st *tpc_track = 0;
695  if (g2t_track) tpc_track = g2t_track->GetTable();
696  St_g2t_vertex *g2t_ver = (St_g2t_vertex *) GetDataSet("geant/g2t_vertex");// if (!g2t_ver) return kStWarn;
697  g2t_vertex_st *gver = 0;
698  Int_t NV = 0;
699  fgTriggerT0 = 0;
700  Double_t mTimeBinWidth = 1./StTpcDb::instance()->Electronics()->samplingFrequency();
701  Double_t driftVelocity = StTpcDb::instance()->DriftVelocity(1);
702  if (g2t_ver) {
703  gver = g2t_ver->GetTable();
704  NV = g2t_ver->GetNRows();
705  StEvent* pEvent = dynamic_cast<StEvent*> (GetInputDS("StEvent"));
706  if (pEvent && StTpcBXT0CorrEPDC::instance()->nrows() && pEvent->epdCollection() ) {
707  int TAC = 0;
708  int maxTAC = -1;
709  StEpdCollection * epdCol = pEvent->epdCollection();
710  if (epdCol) {
711  StSPtrVecEpdHit &epdHits = epdCol->epdHits();
712  int nEpdHits = epdHits.size();
713 
714  for(int i = 0; i < nEpdHits; i++) {
715  StEpdHit * epdHit = dynamic_cast<StEpdHit*>(epdHits[i]);
716  TAC = 0;
717  if (epdHit->tile() > 9) continue; // only tiles 1 - 9 have timing info
718  // if (epdHit->id() < 0) ew = -1; // tile is on the east
719  // else ew = 1;
720  if (epdHit->adc() < 100) continue;
721  TAC = epdHit->tac(); // this is the timing
722  if (TAC > maxTAC) maxTAC = TAC;
723  }
724  }
725  fgTriggerT0 = - StTpcBXT0CorrEPDC::instance()->getCorrection(maxTAC, driftVelocity, mTimeBinWidth)*mTimeBinWidth*1e-6;
726  } else if (pEvent && IAttr("EbyET0")) {
727  // StEbyET0 returns microsec, will need it in seconds
728  fgTriggerT0 = - StEbyET0::Instance()->getT0(pEvent)*1e-6;
729 #if 0
730  } else if (g2t_ver->GetNRows() > 0) {
731  Double_t beta = St_beamInfoC::instance()->BetaYellow();
732  fgTriggerT0 = gver->ge_x[2]/(beta*TMath::Ccgs()); // swap sign 01/11/21
733 #endif
734  }
735  }
736  g2t_tpc_hit_st *tpc_hit_begin = g2t_tpc_hit->GetTable();
737  g2t_tpc_hit_st *tpc_hit = tpc_hit_begin;
738  if (m_TpcdEdxCorrection) {
739  St_tpcGainCorrectionC::instance()->Struct(0)->min = -500;
740  St_tpcGainCorrectionC::instance()->Struct(1)->min = -500;
741  if (Debug()) {
742  LOG_INFO << "Reset min for gain Correction to I/O\t"
743  << St_tpcGainCorrectionC::instance()->Struct(1)->min
744  << "\t"
745  << St_tpcGainCorrectionC::instance()->Struct(0)->min
746  << " (V)" << endm;
747  }
748  }
749  mNSplittedHits = 0;
750  // sort
751  TTableSorter sorter(g2t_tpc_hit,&SearchT,&CompareT);//, 0, no_tpc_hits);
752  Int_t sortedIndex = 0;
753  tpc_hit = tpc_hit_begin;
754  for (Int_t sector = minSector; sector <= maxSector; sector++) {
755  Int_t NoHitsInTheSector = 0;
756  free(m_SignalSum); m_SignalSum = 0;
757  ResetSignalSum(sector);
758  // it is assumed that hit are ordered by sector, trackId, pad rows, and track length
759  for (; sortedIndex < no_tpc_hits; sortedIndex++) {
760  Int_t indx = sorter.GetIndex(sortedIndex);
761  if (indx < 0) break;
762  tpc_hit = tpc_hit_begin + indx;
763  Int_t volId = tpc_hit->volume_id%10000;
764  Int_t iSector = volId/100;
765  Int_t row = volId%100;
766  if (row < minRow || row > maxRow) continue;
767  if (iSector != sector) {
768  if (! ( iSector > sector ) ) {
769  LOG_ERROR << "StTpcRSMaker::Make: g2t_tpc_hit table has not been ordered by sector no. " << sector << " and iSector = " << iSector << ". Skip the hit." << endm;
770  g2t_tpc_hit->Print(indx,1);
771  continue;
772  // assert( iSector > sector );
773  }
774  break;
775  }
776  if (tpc_hit->volume_id <= 0 || tpc_hit->volume_id > 1000000) continue;
777  Int_t Id = tpc_hit->track_p;
778  if (Id <= 0 || Id > NoTpcTracks) {
779  LOG_ERROR << "StTpcRSMaker::Make: g2t_tpc_hit table does not matched with g2t_track: track Id = " << Id << " and NoTpcTracks = " << NoTpcTracks << " skip the hit" << endm;
780  g2t_tpc_hit->Print(indx,1);
781  continue;
782  }
783  Int_t id3 = 0, ipart = 8, charge = 1;
784  Double_t mass = 0;
785  if (tpc_track) {
786  id3 = tpc_track[Id-1].start_vertex_p;
787  // assert(id3 > 0 && id3 <= NV);
788  if (id3 <= 0 || id3 > NV) {
789  LOG_ERROR << "StTpcRSMaker::Make: g2t_tpc_hit table does not matched with g2t_vertex: id3 = " << id3 << " and NV = " << NV << " skip the hit" << endm;
790  g2t_tpc_hit->Print(indx,1);
791  if (Id > 0) g2t_track->Print(Id-1,1);
792  continue;
793  }
794  ipart = tpc_track[Id-1].ge_pid;
795  charge = (Int_t) tpc_track[Id-1].charge;
796 
797  }
798  if (ipart == Laserino || ipart == Chasrino) {
799  charge = 0;
800  } else {
801  if (ipart == 1) {// gamma => electron
802  ipart = 3;
803  charge = -1;
804  }
805  if (charge == 0) {
806  continue;
807  }
808  } // special treatment for electron/positron
809  // Track segment to propagate
810  enum {NoMaxTrackSegmentHits = 100};
811  static HitPoint_t TrackSegmentHits[NoMaxTrackSegmentHits];
812  msMin = 9999;
813  msMax = -9999;
814  Int_t nSegHits = 0;
815  Int_t sIndex = sortedIndex;
816  if (Debug() > 13) cout << "sortedIndex = " << sortedIndex << "\tno_tpc_hits = " << no_tpc_hits << endl;
817  Int_t ID = 0;
818  Double_t pOld = 0;
819  Double_t zOld = -999;
820  Int_t TrackDirection = 0; // 0 - increase no of row, 1 - decrease no of. row.
821  for (nSegHits = 0, sIndex = sortedIndex;
822  sIndex < no_tpc_hits && nSegHits < NoMaxTrackSegmentHits; sIndex++) {
823  indx = sorter.GetIndex(sIndex);
824  g2t_tpc_hit_st *tpc_hitC = tpc_hit_begin + indx;
825  if ((tpc_hitC->volume_id%10000)/100 != sector) break;
826  if (ID > 0 && ID != tpc_hitC->track_p) break;
827  // separater delta electons
828  Double_t pNew = TMath::Sqrt(tpc_hitC->p[0]*tpc_hitC->p[0] + tpc_hitC->p[1]*tpc_hitC->p[1] + tpc_hitC->p[2]*tpc_hitC->p[2]);
829  Double_t zNew = tpc_hitC->x[2];
830  if (pOld > 0) {
831  Double_t rNewOld = pNew/pOld;
832  if (rNewOld < 0.1 || rNewOld > 10) break;
833  if (zOld > -999 && TMath::Abs(zNew - zOld) > 10) break;
834  }
835  pOld = pNew;
836  zOld = zNew;
837  ID = tpc_hitC->track_p;
838 
839  if (nSegHits == 1) { // No Loopers !
840  if (TrackSegmentHits[nSegHits-1].tpc_hitC->volume_id%100 <= tpc_hitC->volume_id%100) {
841  TrackDirection = 0;
842  } else {
843  TrackDirection = 1;
844  }
845  } else if (nSegHits > 1) {
846  if ((! TrackDirection && TrackSegmentHits[nSegHits-1].tpc_hitC->volume_id%100 > tpc_hitC->volume_id%100) ||
847  ( TrackDirection && TrackSegmentHits[nSegHits-1].tpc_hitC->volume_id%100 < tpc_hitC->volume_id%100))
848  break;
849  }
850  if (Debug() > 13) cout << "sIndex = " << sIndex << "\tindx = " << indx << "\ttpc_hitC = " << tpc_hitC << endl;
851  TrackSegmentHits[nSegHits].indx = indx;
852  TrackSegmentHits[nSegHits].s = tpc_hitC->length;
853  if (tpc_hitC->length == 0 && nSegHits > 0) {
854  TrackSegmentHits[nSegHits].s = TrackSegmentHits[nSegHits-1].s + tpc_hitC->ds;
855  }
856  TrackSegment2Propagate(tpc_hitC, &gver[id3-1],TrackSegmentHits[nSegHits]);
857  if (TrackSegmentHits[nSegHits].Pad.timeBucket() < 0 || TrackSegmentHits[nSegHits].Pad.timeBucket() > NoOfTimeBins) continue;
858  nSegHits++;
859  }
860  if (! nSegHits) continue;
861  if (Debug() >= 10) {
862  PrPP(Make,nSegHits);
863  for (Int_t s = 0; s < nSegHits; s++) {
864  cout << "Seg[" << Form("%2i",s) << "]\tId " << TrackSegmentHits[s].TrackId << "\ts = " << TrackSegmentHits[s].s
865  << "\tvolumeID :" << Form("%6i",TrackSegmentHits[s].tpc_hitC->volume_id) <<"\t" << TrackSegmentHits[s].Pad
866  << "\ts1/s2 = " << TrackSegmentHits[s].tpc_hitC->length - TrackSegmentHits[s].tpc_hitC->ds/2
867  << "\t" << TrackSegmentHits[s].tpc_hitC->length + TrackSegmentHits[s].tpc_hitC->ds/2 << "\tds = " << TrackSegmentHits[s].tpc_hitC->ds
868  << endl;
869  }
870  }
871  sortedIndex = sIndex-1; // Irakli 05/06/19, reduce extra step in for loop
872  Double_t s = msMin;
873  memset (rowsdE, 0, sizeof(rowsdE));
874  Double_t zIntDr = - TMath::Log(gRandom->Rndm());
875  for (Int_t iSegHits = 0; iSegHits < nSegHits && s < msMax; iSegHits++) {
876  memset (rowsdEH, 0, sizeof(rowsdEH));
877  g2t_tpc_hit_st *tpc_hitC = TrackSegmentHits[iSegHits].tpc_hitC;
878  Double_t dStep = TMath::Abs(tpc_hitC->ds);
879  tpc_hitC->adc = 0;
880  tpc_hitC->dESum = 0;
881  tpc_hitC->dSSum = 0;
882  tpc_hitC->np = 0;
883  tpc_hitC->ne = 0;
884  memset (tpc_hitC->adcs, 0, sizeof(tpc_hitC->adcs));
885  volId = tpc_hitC->volume_id%100000;
886 
887  Int_t row = TrackSegmentHits[iSegHits].coorLS.fromRow();
888  Int_t io = (row <= St_tpcPadConfigC::instance()->numberOfInnerRows(sector)) ? 0 : 1;
889  // switch between Inner / Outer Sector paramters
890  // Extra correction for simulation with respect to data
891  Int_t iowe = 0;
892  if (sector > 12) iowe += 4;
893  if (io) iowe += 2;
894  Float_t *AdditionalMcCorrection = St_TpcResponseSimulatorC::instance()->SecRowCor();
895  Float_t *AddSigmaMcCorrection = St_TpcResponseSimulatorC::instance()->SecRowSig();
896  // Generate signal
897  Double_t sigmaJitterT = St_TpcResponseSimulatorC::instance()->SigmaJitterTI();
898  Double_t sigmaJitterX = St_TpcResponseSimulatorC::instance()->SigmaJitterXI();
899  if(io) { // Outer
900  sigmaJitterT = St_TpcResponseSimulatorC::instance()->SigmaJitterTO();
901  sigmaJitterX = St_TpcResponseSimulatorC::instance()->SigmaJitterXO();
902  }
903  // Generate signal
904  Double_t Gain = St_tss_tssparC::instance()->gain(sector,row);
905  if (ClusterProfile) {
906  checkList[io][2]->Fill(TrackSegmentHits[iSegHits].xyzG.position().z(),Gain);
907  }
908  Double_t GainXCorrectionL = AdditionalMcCorrection[iowe] + row*AdditionalMcCorrection[iowe+1];
909  Gain *= TMath::Exp(-GainXCorrectionL);
910  Double_t GainXSigma = AddSigmaMcCorrection[iowe] + row*AddSigmaMcCorrection[iowe+1];
911  if (GainXSigma > 0) Gain *= TMath::Exp(gRandom->Gaus(0.,GainXSigma));
912  if (ClusterProfile) {
913  checkList[io][3]->Fill(TrackSegmentHits[iSegHits].xyzG.position().z(),Gain);
914  }
915  // dE/dx correction
916  Double_t dEdxCor = 1.0;
917  if (TrackSegmentHits[iSegHits].coorLS.position().z() < -0.6) {// prompt hist
918  dEdxCor = 1.0;
919  } else {
920  dEdxCor = dEdxCorrection(TrackSegmentHits[iSegHits]);
921  }
922 #ifdef __DEBUG__
923  if (TMath::IsNaN(dEdxCor)) {
924  iBreak++;
925  }
926 #endif
927  if (dEdxCor < minSignal) continue;
928  if (ClusterProfile) {
929  checkList[io][4]->Fill(TrackSegmentHits[iSegHits].xyzG.position().z(),dEdxCor);
930  GainHist[0]->Fill(TrackSegmentHits[iSegHits].Pad.sector(),TrackSegmentHits[iSegHits].Pad.row(),dEdxCor);
931  }
932  // Initialize propagation
933  Float_t BField[3] = {(Float_t ) TrackSegmentHits[iSegHits].BLS.position().x(),
934  (Float_t ) TrackSegmentHits[iSegHits].BLS.position().y(),
935  (Float_t ) TrackSegmentHits[iSegHits].BLS.position().z()};
936  StPhysicalHelixD track(TrackSegmentHits[iSegHits].dirLS.position(),
937  TrackSegmentHits[iSegHits].coorLS.position(),
938  BField[2]*kilogauss,charge);
939  StThreeVectorD unit = TrackSegmentHits[iSegHits].dirLS.position().unit();
940  Double_t *cxyz = unit.xyz();
941  TRMatrix L2L(3,3,
942  cxyz[2], - cxyz[0]*cxyz[2] , cxyz[0],
943  cxyz[0], - cxyz[1]*cxyz[2] , cxyz[1],
944  0.0 , cxyz[0]*cxyz[0] + cxyz[1]*cxyz[1], cxyz[2]);
945 #ifdef __DEBUG__
946  if (Debug() > 11) PrPP(Make,track);
947 #endif
948  static StThreeVectorD normal(0,1,0);
949  static StTpcCoordinateTransform transform(gStTpcDb);
950  StThreeVectorD rowPlane(0,transform.yFromRow(TrackSegmentHits[iSegHits].Pad.sector(),TrackSegmentHits[iSegHits].Pad.row()),0);
951  Double_t sR = track.pathLength(rowPlane,normal);
952  if (sR < 1e10) {
953  PrPP(Maker,sR);
954  PrPP(Make,TrackSegmentHits[iSegHits].coorLS);
955  StThreeVectorD xyzP = track.at(sR);
956  TrackSegmentHits[iSegHits].coorLS.setPosition(xyzP); PrPP(Make,TrackSegmentHits[iSegHits].coorLS);
957  // don't useT0, don't useTau
958  PrPP(Make,TrackSegmentHits[iSegHits].Pad);
959  transform(TrackSegmentHits[iSegHits].coorLS,TrackSegmentHits[iSegHits].Pad,kFALSE,kFALSE); // don't use T0, don't use Tau
960  PrPP(Make,TrackSegmentHits[iSegHits].Pad);
961  }
962  Int_t ioH = io;
963  if (St_tpcAltroParamsC::instance()->N(sector-1) >= 0) ioH += 2;
964  TotalSignal = 0;
965  Double_t lgam = tpc_hitC->lgam;
966  if (ClusterProfile) {
967  checkList[io][5]->Fill(TrackSegmentHits[iSegHits].xyzG.position().z(),lgam);
968  }
969  Double_t gamma = TMath::Power(10.,lgam) + 1;
970  Double_t betaGamma = TMath::Sqrt(gamma*gamma - 1.);
971  StThreeVectorD pxyzG(tpc_hitC->p[0],tpc_hitC->p[1],tpc_hitC->p[2]);
972  Double_t bg = 0;
973  static const Double_t m_e = .51099907e-3;
974  if (mass > 0) bg = pxyzG.mag()/mass;
975  if (ClusterProfile && bg > 0) {
976  checkList[io][21]->Fill(TMath::Log10(betaGamma),TMath::Log10(bg));
977  }
978  if (bg > betaGamma) betaGamma = bg;
979  Double_t bg2 = betaGamma*betaGamma;
980  gamma = TMath::Sqrt(bg2 + 1.);
981  Double_t Tmax;
982  if (mass < 2*m_e) {
983  if (charge > 0) Tmax = m_e*(gamma - 1);
984  else Tmax = 0.5*m_e*(gamma - 1);
985  } else {
986  Double_t r = m_e/mass;
987  Tmax = 2*m_e*bg2/(1 + 2*gamma*r + r*r);
988  }
989  if (Tmax > mCutEle) Tmax = mCutEle;
990  Double_t padH = TrackSegmentHits[iSegHits].Pad.pad();
991  Double_t tbkH = TrackSegmentHits[iSegHits].Pad.timeBucket();
992  tpc_hitC->pad = padH;
993  tpc_hitC->timebucket = tbkH;
994  pad0 = TMath::Max(1,TMath::Nint(padH + xmin[0]));
995  tbk0 = TMath::Max(0,TMath::Nint(tbkH + xmin[1]));
996  Double_t OmegaTau = St_TpcResponseSimulatorC::instance()->OmegaTau()*
997  TrackSegmentHits[iSegHits].BLS.position().z()/5.0;// from diffusion 586 um / 106 um at B = 0/ 5kG
998  Double_t NP = TMath::Abs(tpc_hitC->de)/(St_TpcResponseSimulatorC::instance()->W()*eV*
999  St_TpcResponseSimulatorC::instance()->Cluster()); // from GEANT
1000  if (ClusterProfile) {
1001  checkList[io][6]->Fill(TrackSegmentHits[iSegHits].xyzG.position().z(),NP);
1002  }
1003  Double_t driftLength = TMath::Abs(TrackSegmentHits[iSegHits].coorLS.position().z());
1004  Double_t D = 1. + OmegaTau*OmegaTau;
1005  Double_t SigmaT = St_TpcResponseSimulatorC::instance()->transverseDiffusion()* TMath::Sqrt( driftLength/D);
1006  Double_t SigmaL = St_TpcResponseSimulatorC::instance()->longitudinalDiffusion()*TMath::Sqrt( driftLength );
1007  if (St_tpcPadConfigC::instance()->IsRowInner(sector,row)) {
1008  if ( St_TpcResponseSimulatorC::instance()->transverseDiffusionI() > 0.0)
1009  SigmaT = St_TpcResponseSimulatorC::instance()->transverseDiffusionI()* TMath::Sqrt( driftLength/D);
1010  if (St_TpcResponseSimulatorC::instance()->longitudinalDiffusionI() > 0.0)
1011  SigmaL = St_TpcResponseSimulatorC::instance()->longitudinalDiffusionI()*TMath::Sqrt( driftLength );
1012  }
1013  // Double_t SigmaL = St_TpcResponseSimulatorC::instance()->longitudinalDiffusion()*TMath::Sqrt(2*driftLength );
1014  if (sigmaJitterX > 0) {SigmaT = TMath::Sqrt(SigmaT*SigmaT + sigmaJitterX*sigmaJitterX);}
1015  Double_t NoElPerAdc = St_TpcResponseSimulatorC::instance()->NoElPerAdc();
1016  if (NoElPerAdc <= 0) {
1017  if (St_tpcPadConfigC::instance()->iTPC(sector) && St_tpcPadConfigC::instance()->IsRowInner(sector,row)) {
1018  NoElPerAdc = St_TpcResponseSimulatorC::instance()->NoElPerAdcX(); // iTPC
1019  } else if (St_tpcPadConfigC::instance()->IsRowInner(sector,row)) {
1020  NoElPerAdc = St_TpcResponseSimulatorC::instance()->NoElPerAdcI(); // inner TPX
1021  } else {
1022  NoElPerAdc = St_TpcResponseSimulatorC::instance()->NoElPerAdcO(); // outer TPX
1023  }
1024  }
1025 #ifndef __NO_1STROWCORRECTION__
1026  if (row == 1) dEdxCor *= TMath::Exp(St_TpcResponseSimulatorC::instance()->FirstRowC());
1027 #endif /* __NO_1STROWCORRECTION__ */
1028  mGainLocal = Gain/dEdxCor/NoElPerAdc; // Account dE/dx calibration
1029  // end of dE/dx correction
1030  // generate electrons: No. of primary clusters per cm
1031  NP = StdEdxModel::instance()->dNdx(betaGamma,charge); // per cm
1032  if (ClusterProfile) {
1033  checkList[io][7]->Fill(TrackSegmentHits[iSegHits].xyzG.position().z(),NP);
1034  }
1035  memset (padsdE, 0, sizeof(padsdE));
1036  memset (tbksdE, 0, sizeof(tbksdE));
1037  Float_t dEr = 0;
1038  tpc_hitC->dNdx = NP;
1039  do {// Clusters
1040  Float_t dS = 0;
1041  Float_t dE = 0;
1042  NP = tpc_hitC->dNdx;
1043  if (charge) {
1044  dS = zIntDr/NP;
1045  dE = StdEdxModel::instance()->dNdE();
1046  } else { // charge == 0 geantino
1047  // for Laserino assume dE/dx = 25 keV/cm;
1048  dE = 10; // eV
1049  dS = dE*eV/(TMath::Abs(mLaserScale*tpc_hitC->de/tpc_hitC->ds));
1050  }
1051  Double_t E = dE*eV;
1052  Double_t step = dStep - tpc_hitC->dSSum;
1053  if (dS > step) {
1054  zIntDr -= step*NP;
1055  tpc_hitC->dSSum += step;
1056 #ifdef __DEBUG__
1057  if (Debug() > 12) {
1058  LOG_INFO << "dESum = " << tpc_hitC->dESum << " /\tdSSum " << tpc_hitC->dSSum << " /\tds " << tpc_hitC->ds << endm;
1059  }
1060 #endif
1061  break;
1062  } else {
1063  tpc_hitC->dSSum += dS;
1064  zIntDr = - TMath::Log(gRandom->Rndm());
1065  }
1066  if (dE < St_TpcResponseSimulatorC::instance()->W()/2 || E > Tmax) continue;
1067  tpc_hitC->dESum += E;
1068  tpc_hitC->np++;
1069 #ifdef __DEBUG__
1070  if (Debug() > 12) {
1071  LOG_INFO << "dESum = " << tpc_hitC->dESum << " /\tdSSum " << tpc_hitC->dSSum << " /\tds " << tpc_hitC->ds << endm;
1072  }
1073 #endif
1074  Double_t xRange = 0;
1075  if (dE > ElectronRangeEnergy) xRange = ElectronRange*TMath::Power((dE+dEr)/ElectronRangeEnergy,ElectronRangePower);
1076  Int_t Nt = 0; // HeedCsize(dE, dEr,rs);
1077  Float_t dET = dE + dEr;
1078  dEr = dET;
1079  Float_t EC;
1080  Int_t Nr = 0;
1081  if (xRange > 0) {Nr = rs.GetSize();}
1082  while ((EC = mHeed->GetRandom()) < dEr) {
1083  dEr -= EC;
1084  if (Nr) {
1085  if (Nr <= Nt) {Nr = 2*Nt+1; rs.Set(Nr); }
1086  rs[Nt] = 1 - dEr/dET;
1087  }
1088  Nt++;
1089  }
1090  if (! Nt) continue;
1091  if (ClusterProfile) {
1092  checkList[io][8]->Fill(TrackSegmentHits[iSegHits].xyzG.position().z(),Nt);
1093  checkList[io][11]->Fill(TrackSegmentHits[iSegHits].xyzG.position().z(),Nt);
1094  }
1095  StThreeVectorD xyzC = track.at(tpc_hitC->dSSum - dStep/2);
1096  Double_t phiXY = 2*TMath::Pi()*gRandom->Rndm();
1097  Double_t rX = TMath::Cos(phiXY);
1098  Double_t rY = TMath::Sin(phiXY);
1099  Double_t sigmaT = SigmaT;
1100  Double_t sigmaL = SigmaL;
1101  TotalSignalInCluster = 0;
1102  Int_t WireIndex = 0;
1103  for (Int_t ie = 0; ie < Nt; ie++) {
1104  // transport to wire
1105  gRandom->Rannor(rX,rY);
1106  StTpcLocalSectorCoordinate xyzE(xyzC.x()+rX*sigmaT,
1107  xyzC.y()+rY*sigmaT,
1108  xyzC.z()+gRandom->Gaus(0,sigmaL), sector, row);
1109  if (xRange > 0) {
1110  Double_t xR = rs[ie]*xRange;
1111  TRVector xyzRangeL(3, xR*rX, xR*rY, 0.);
1112  TRVector xyzR(L2L,TRArray::kAxB,xyzRangeL);
1113 #ifdef __DEBUG__
1114  if (Debug() > 13) {
1115  LOG_INFO << "xyzRangeL: " << xyzRangeL << endm;
1116  LOG_INFO << "L2L: " << L2L << endm;
1117  LOG_INFO << "xyzR: " << xyzR << endm;
1118  }
1119 #endif
1120  TCL::vadd(xyzE.position().xyz(),xyzR.GetArray(),xyzE.position().xyz(),3);
1121  }
1122  // Check Gatting Grid
1123  Double_t zGG = xyzE.position().z();
1124  Double_t GGTransperency = 1.0;
1125  if (zGG > -0.6) { // non Prompt hits before Cathode wire plane
1126  Double_t driftTime = 1e6*zGG/driftVelocity; // usec
1127  driftTime -= St_starTriggerDelayC::instance()->TrigT0GG(io); // trigger delay + Gating Grid cables
1128  Double_t lGGTransperency = St_GatingGridBC::instance()->CalcCorrection(0,driftTime);
1129  if (lGGTransperency < -9) continue;
1130  GGTransperency = TMath::Exp(lGGTransperency);
1131  if (GGTransperency < minSignal) continue;
1132  }
1133  tpc_hitC->ne++;
1134  QAv = GGTransperency*mPolya[io]->GetRandom();
1135  Double_t y = xyzE.position().y();
1136  Double_t alphaVariation = InnerAlphaVariation[sector-1];
1137  // Transport to wire
1138  if (y <= lastInnerSectorAnodeWire) {
1139  WireIndex = TMath::Nint((y - firstInnerSectorAnodeWire)/anodeWirePitch) + 1;
1140 #ifndef __NO_1STROWCORRECTION__
1141  if (St_tpcPadConfigC::instance()->iTPC(sector)) {// two first and two last wires are removed, and 3rd wire is fat wiere
1142  if (WireIndex <= 3 || WireIndex >= numberOfInnerSectorAnodeWires - 3) continue;
1143  } else { // old TPC the first and last wires are fat ones
1144  if (WireIndex <= 1 || WireIndex >= numberOfInnerSectorAnodeWires) continue;
1145  }
1146 #else /* __NO_1STROWCORRECTION__ */
1147  if (WireIndex <= 1 || WireIndex >= numberOfInnerSectorAnodeWires) continue; // to check the 1-st pad row effect
1148 #endif /* ! __NO_1STROWCORRECTION__ */
1149  yOnWire = firstInnerSectorAnodeWire + (WireIndex-1)*anodeWirePitch;
1150  } else { // the first and last wires are fat ones
1151  WireIndex = TMath::Nint((y - firstOuterSectorAnodeWire)/anodeWirePitch) + 1;
1152  if (WireIndex <= 1 || WireIndex >= numberOfOuterSectorAnodeWires) continue;
1153  yOnWire = firstOuterSectorAnodeWire + (WireIndex-1)*anodeWirePitch;
1154  alphaVariation = OuterAlphaVariation[sector-1];
1155  }
1156  Double_t distanceToWire = y - yOnWire; // Calculated effective distance to wire affected by Lorentz shift
1157  xOnWire = xyzE.position().x();
1158  zOnWire = xyzE.position().z();
1159  // Grid plane (1 mm spacing) focusing effect + Lorentz angle in drift volume
1160  Int_t iGridWire = (Int_t ) TMath::Abs(10.*distanceToWire);
1161  Double_t dist2Grid = TMath::Sign(0.05 + 0.1*iGridWire, distanceToWire); // [cm]
1162  // Ground plane (1 mm spacing) focusing effect
1163  Int_t iGroundWire = (Int_t ) TMath::Abs(10.*dist2Grid);
1164  Double_t distFocused = TMath::Sign(0.05 + 0.1*iGroundWire, dist2Grid);
1165  // OmegaTau near wires taken from comparison with data
1166  Double_t tanLorentz = 0;
1167  if (y < firstOuterSectorAnodeWire) {
1168  if (St_TpcResponseSimulatorC::instance()->OmegaTauScaleI() > 0) tanLorentz = OmegaTau/St_TpcResponseSimulatorC::instance()->OmegaTauScaleI();
1169  } else {
1170  if (St_TpcResponseSimulatorC::instance()->OmegaTauScaleO() > 0) tanLorentz = OmegaTau/St_TpcResponseSimulatorC::instance()->OmegaTauScaleO();
1171  }
1172  xOnWire += distFocused*tanLorentz; // tanLorentz near wires taken from comparison with data
1173  zOnWire += TMath::Abs(distFocused);
1174  if (! iGroundWire ) QAv *= TMath::Exp( alphaVariation);
1175  else QAv *= TMath::Exp(-alphaVariation);
1176  if (ClusterProfile) {
1177  checkList[io][9]->Fill(TrackSegmentHits[iSegHits].xyzG.position().z(),QAv);
1178  }
1179  Double_t dY = mChargeFraction[io][sector-1]->GetXmax();
1180  Double_t yLmin = yOnWire - dY;
1181  Double_t yLmax = yLmin + 2*dY;
1182  Int_t rowMin = transform.rowFromLocalY(yLmin,sector);
1183  Int_t rowMax = transform.rowFromLocalY(yLmax,sector);
1184  Double_t yRmin = transform.yFromRow(sector,rowMin) - St_tpcPadConfigC::instance()->PadLengthAtRow(sector,rowMin)/2;
1185  Double_t yRmax = transform.yFromRow(sector,rowMax) + St_tpcPadConfigC::instance()->PadLengthAtRow(sector,rowMax)/2;
1186  if (yRmin > yLmax || yRmax < yLmin) {
1187  iBreak++; continue;
1188  }
1189  GenerateSignal(TrackSegmentHits[iSegHits],sector,rowMin, rowMax, sigmaJitterT, sigmaJitterX);
1190  } // electrons in Cluster
1191  if (ClusterProfile) {
1192  if (TotalSignalInCluster > 0 && checkList[io][19]) {
1193  checkList[io][19]->Fill(WireIndex,TMath::Log(TotalSignalInCluster));
1194 
1195  }
1196  }
1197  TotalSignal += TotalSignalInCluster;
1198  } while (kTRUE); // Clusters
1199  //03.14.2022 tpc_hitC->adc = -99;
1200  if (tpc_hitC->dESum > 0 && tpc_hitC->dSSum > 0) {
1201 #ifdef __DEBUG__
1202  if (Debug() > 12) {
1203  LOG_INFO << "sIndex = " << sIndex << " volId = " << volId
1204  << " tpc_hitC->de/eV = " << tpc_hitC->de/eV << " /\tdSSum " << tpc_hitC->dSSum << " /\t TotalSignal " << TotalSignal << endm;
1205  }
1206 #endif
1207  if (row > 1) tpc_hitC->adcs[0] += rowsdEH[row-2];
1208  tpc_hitC->adcs[1] += rowsdEH[row-1];
1209  if (row <= kRowMax) tpc_hitC->adcs[2] += rowsdEH[row];
1210  if (ClusterProfile) {
1211  if (TotalSignal > 0) {
1212  if (hist[ioH][0]) {
1213  for (Int_t p = 0; p < kPadMax; p++) {
1214  hist[ioH][0]->Fill((p+pad0)-padH,TrackSegmentHits[iSegHits].xyzG.position().z(),padsdE[p]/TotalSignal);
1215  PixelHist[0]->Fill(row, p+pad0,padsdE[p]/TotalSignal);
1216  }
1217  }
1218  if (hist[ioH][1]) {
1219  for (Int_t t = 0; t < kTimeBacketMax; t++) {
1220  hist[ioH][1]->Fill((t+tbk0+0.5)-tbkH,TrackSegmentHits[iSegHits].xyzG.position().z(),tbksdE[t]/TotalSignal);
1221  PixelHist[1]->Fill(row, (t+tbk0+0.5),tbksdE[t]/TotalSignal);
1222  }
1223  }
1224  }
1225  checkList[io][15]->Fill(TrackSegmentHits[iSegHits].xyzG.position().z(),tpc_hitC->de);
1226  checkList[io][16]->Fill(TrackSegmentHits[iSegHits].xyzG.position().z(),tpc_hitC->ds);
1227  checkList[io][18]->Fill(TrackSegmentHits[iSegHits].xyzG.position().z(),tpc_hitC->ne);
1228  if (tpc_hitC->np > 0 && tpc_hitC->ne > 0)
1229  checkList[io][20]->Fill(TMath::Log(tpc_hitC->np),TMath::Log(tpc_hitC->ne) - TMath::Log(tpc_hitC->np));
1230  }
1231  }
1232  NoHitsInTheSector++;
1233  } // end do loop over segments for a given particle
1234  for (Int_t iSegHits = 0; iSegHits < nSegHits; iSegHits++) {
1235  g2t_tpc_hit_st *tpc_hitC = TrackSegmentHits[iSegHits].tpc_hitC;
1236  if (tpc_hitC->volume_id > 10000) continue;
1237  Int_t row = tpc_hitC->volume_id%100;
1238  //03.14.2022
1239  tpc_hitC->adc += rowsdE[row-1];
1240  Int_t io = (row <= St_tpcPadConfigC::instance()->numberOfInnerRows(sector)) ? 0 : 1;
1241  if (checkList[io][17])
1242  checkList[io][17]->Fill(TrackSegmentHits[iSegHits].xyzG.position().z(),tpc_hitC->adc);
1243  }
1244  } // hits in the sector
1245  if (NoHitsInTheSector) {
1246  StTpcDigitalSector *digitalSector = DigitizeSector(sector);
1247  if (Debug()) LOG_INFO << "StTpcRSMaker: Done with sector\t" << sector << " total no. of hit = " << NoHitsInTheSector << endm;
1248  if (Debug() > 2) digitalSector->Print();
1249  }
1250  } // sector
1251  if (m_SignalSum) {free(m_SignalSum); m_SignalSum = 0;}
1252  if (Debug()%10) gBenchmark->Show("TpcRS");
1253  if (m_TpcdEdxCorrection) {
1254  St_tpcGainCorrectionC::instance()->Struct(1)->min = vminI;
1255  St_tpcGainCorrectionC::instance()->Struct(0)->min = vminO;
1256  if (Debug()) {
1257  LOG_INFO << "Reset min for gain Correction to I/O\t"
1258  << St_tpcGainCorrectionC::instance()->Struct(1)->min
1259  << "\t"
1260  << St_tpcGainCorrectionC::instance()->Struct(0)->min
1261  << " (V)" << endm;
1262  }
1263  }
1264  if (g2t_track) {
1265  // Reset no. Tpc hits in g2t_track
1266  tpc_track = g2t_track->GetTable();
1267  for (Int_t i = 0; i < NoTpcTracks; i++, tpc_track++) {
1268  Int_t Id = tpc_track->id;
1269  tpc_track->n_tpc_hit = (mNoTpcHitsReal[Id-1] << 8) + (0xff & mNoTpcHitsAll[Id-1]);
1270  }
1271  }
1272  return kStOK;
1273 }
1274 //________________________________________________________________________________
1275 Double_t StTpcRSMaker::ShaperFunc(Double_t *x, Double_t *par) {
1276  Double_t tau = par[0];
1277  Double_t width = par[1];
1278  Double_t p = par[2];
1279  Double_t t = x[0]*width/tau;
1280  Double_t Delta = width/tau;
1281  Double_t t1 = t - Delta/2.;
1282  Double_t t2 = t1 + Delta;
1283  Double_t val = TMath::Gamma(p,t2) - TMath::Gamma(p,t1);
1284  return val;
1285 }
1286 //________________________________________________________________________________
1287 Double_t StTpcRSMaker::PadResponseFunc(Double_t *x, Double_t *par) {
1288  Double_t CrossTalk = 0;
1289  Double_t Value = 0;
1290  Double_t X = par[5]*x[0];
1291  if (CrossTalk > 0) {
1292  for (Int_t i = -1; i <= 1; i++) {
1293  Double_t xx = X + par[5]*i;
1294  if (i == 0) Value += (1. - 2.*CrossTalk)*Gatti(&xx,par);
1295  else Value += CrossTalk *Gatti(&xx,par);
1296  }
1297  } else Value = Gatti(&X,par);
1298  return Value;
1299 }
1300 //________________________________________________________________________________
1301 Double_t StTpcRSMaker::Gatti(Double_t *x, Double_t *par) {
1302  /***********************************************x[*************************
1303  * Function : generates the cathode signal using *
1304  * the single-parameter Gatti formula: *
1305  * 1 - tanh(K2 * lambda)**2 *
1306  * GFunc(lambda) = K1 * ------------------------------- *
1307  * 1 + K3 * tanh (K2 *lambda)**2 *
1308  * lambda = x/h, h is anode cathode spacing *
1309  * *
1310  * K2 = pi/2*(1 - 0.5*sqrt(K3)) *
1311  * *
1312  * K2*sqrt(K3) *
1313  * K1 = ------------------- *
1314  * 4 * atan(sqrt(K3)) *
1315  * *
1316  * References : E.Gatti, A.Longoni, NIM 163 (1979) 82-93. *
1317  * Authors : V.Balagura,V.Cherniatin,A.Chikanian *
1318  ************************************************************************/
1319  Double_t y = x[0]; // distance to center of strip [cm]
1320  Double_t w = par[0]; // w = width of pad
1321  Double_t h = par[1]; // h = Anode-Cathode gap
1322  Double_t K3 = par[3];
1323  Double_t lambda = y/h;
1324  Double_t K2 = TMath::PiOver2()*(1. - 0.5*TMath::Sqrt(K3));
1325  // Double_t K1 = K2*TMath::Sqrt(K3)/(2*TMath::ATan(TMath::Sqrt(K3)));
1326  Double_t sqK3 = TMath::Sqrt(K3);
1327  Double_t ATsqK3 = 0.5/TMath::ATan(sqK3);
1328  Double_t Y1 = lambda - w/h/2;
1329  Double_t Y2 = Y1 + w/h;
1330  Double_t X1 = K2*Y1;
1331  Double_t X2 = K2*Y2;
1332  Double_t Z1 = sqK3*TMath::TanH(X1);
1333  Double_t Z2 = sqK3*TMath::TanH(X2);
1334  Double_t val = ATsqK3*(TMath::ATan(Z2) - TMath::ATan(Z1));
1335  return val;
1336 }
1337 //________________________________________________________________________________
1338 void StTpcRSMaker::Print(Option_t */* option */) const {
1339  PrPP(Print, NoOfSectors);
1340  PrPP(Print, St_tpcPadConfigC::instance()->numberOfRows(1));
1341  PrPP(Print, St_tpcPadConfigC::instance()->numberOfRows(20));
1342  PrPP(Print, St_tpcPadConfigC::instance()->numberOfInnerRows(20));
1343  PrPP(Print, NoOfPads);
1344  PrPP(Print, St_TpcResponseSimulatorC::instance()->W());// = 26.2);//*eV
1345  PrPP(Print, St_TpcResponseSimulatorC::instance()->Cluster());
1346  PrPP(Print, St_TpcResponseSimulatorC::instance()->longitudinalDiffusion());
1347  PrPP(Print, St_TpcResponseSimulatorC::instance()->longitudinalDiffusionI());
1348  PrPP(Print, St_TpcResponseSimulatorC::instance()->transverseDiffusion());
1349  PrPP(Print, St_TpcResponseSimulatorC::instance()->transverseDiffusionI());
1350  // PrPP(Print, Gain);
1351  PrPP(Print, NoOfTimeBins);
1352  PrPP(Print, numberOfInnerSectorAnodeWires);
1353  PrPP(Print, firstInnerSectorAnodeWire);
1354  PrPP(Print, lastInnerSectorAnodeWire);
1355  PrPP(Print, numberOfOuterSectorAnodeWires);
1356  PrPP(Print, firstOuterSectorAnodeWire);
1357  PrPP(Print, lastOuterSectorAnodeWire);
1358  PrPP(Print, anodeWirePitch);
1359  PrPP(Print,St_TpcResponseSimulatorC::instance()->OmegaTau()); // tan of Lorentz angle
1360  PrPP(Print, St_TpcResponseSimulatorC::instance()->NoElPerAdcI());
1361  PrPP(Print, St_TpcResponseSimulatorC::instance()->NoElPerAdcO());
1362  PrPP(Print, St_TpcResponseSimulatorC::instance()->NoElPerAdcX());
1363  PrPP(Print, anodeWireRadius);
1364  PrPP(Print, St_TpcResponseSimulatorC::instance()->AveragePedestal());
1365  PrPP(Print, St_TpcResponseSimulatorC::instance()->AveragePedestalRMS());
1366  PrPP(Print, St_TpcResponseSimulatorC::instance()->AveragePedestalRMSX());
1367  PrPP(Print, St_TpcResponseSimulatorC::instance()->FanoFactor());
1368  for (Int_t sector = 1; sector <= 24; sector++) {
1369  PrPP(Print, innerSectorAnodeVoltage[sector-1]);
1370  PrPP(Print, outerSectorAnodeVoltage[sector-1]);
1371  }
1372  PrPP(Print, St_TpcResponseSimulatorC::instance()->K3IP());
1373  PrPP(Print, St_TpcResponseSimulatorC::instance()->K3IR());
1374  PrPP(Print, St_TpcResponseSimulatorC::instance()->K3OP());
1375  PrPP(Print, St_TpcResponseSimulatorC::instance()->K3OR());
1376  PrPP(Print, St_TpcResponseSimulatorC::instance()->SigmaJitterTI());
1377  PrPP(Print, St_TpcResponseSimulatorC::instance()->SigmaJitterTO());
1378 }
1379 //________________________________________________________________________________
1380 StTpcDigitalSector *StTpcRSMaker::DigitizeSector(Int_t sector){
1381  static Int_t iBreak = 0;
1382  static Int_t AdcCut = 500;
1383  // static Int_t PedestalMem[__MaxNumberOfTimeBins__];
1384  TDataSet *event = GetData("Event");
1385  StTpcRawData *data = 0;
1386  if (! event) {
1387  data = new StTpcRawData(NoOfSectors);
1388  event = new TObjectSet("Event", data);
1389  AddData(event);
1390  } else data = (StTpcRawData *) event->GetObject();
1391  assert(data);
1392  SignalSum_t *SignalSum = GetSignalSum(sector);
1393  Double_t ped = 0;
1394  Int_t adc = 0;
1395  Int_t index = 0;
1396  Double_t gain = 1;
1397  Int_t row, pad, bin;
1398  Int_t Sector = TMath::Abs(sector);
1399  StTpcDigitalSector *digitalSector = data->GetSector(Sector);
1400  if (! digitalSector) {
1401  digitalSector = new StTpcDigitalSector(Sector);
1402  data->setSector(Sector,digitalSector);
1403  } else
1404  digitalSector->clear();
1405  for (row = 1; row <= St_tpcPadConfigC::instance()->numberOfRows(sector); row++) {
1406  Int_t noOfPadsAtRow = St_tpcPadConfigC::instance()->St_tpcPadConfigC::instance()->numberOfPadsAtRow(sector,row);
1407  Double_t pedRMS = St_TpcResponseSimulatorC::instance()->AveragePedestalRMS();
1408  Int_t ioA = 1; // Outer
1409  if ( St_tpcPadConfigC::instance()->IsRowInner(sector,row) ) ioA = 0; // Inner
1410  if (St_tpcAltroParamsC::instance()->N(sector-1) > 0) {
1411  if (! (St_tpcPadConfigC::instance()->iTPC(sector) && St_tpcPadConfigC::instance()->IsRowInner(sector,row))) {
1412  pedRMS = St_TpcResponseSimulatorC::instance()->AveragePedestalRMSX();
1413  }
1414  }
1415 #ifdef __DEBUG__
1416  Float_t AdcSumBeforeAltro = 0, AdcSumAfterAltro = 0;
1417 #endif /* __DEBUG__ */
1418  for (pad = 1; pad <= noOfPadsAtRow; pad++) {
1419  gain = St_tpcPadGainT0BC::instance()->Gain(Sector,row,pad);
1420  if (gain <= 0.0) continue;
1421  ped = St_TpcResponseSimulatorC::instance()->AveragePedestal();
1422 #ifdef __TFG__VERSION__
1423  static Int_t IDTs[__MaxNumberOfTimeBins__];
1424 #else /* ! __TFG__VERSION__ */
1425  static UShort_t IDTs[__MaxNumberOfTimeBins__];
1426 #endif /* __TFG__VERSION__ */
1427  Double_t pedRMSpad = pedRMS;
1428  if (pedRMSpad < 0) {
1429  Double_t NPads = St_tpcPadConfigC::instance()->numberOfPadsAtRow(sector,row);
1430  Double_t x = pad/NPads - 0.5;
1431  pedRMSpad = St_TpcPadPedRMSC::instance()->CalcCorrection(1-ioA, x);
1432  ped = gRandom->Gaus(St_TpcPadPedRMSC::instance()->a(1-ioA)[3],St_TpcPadPedRMSC::instance()->a(1-ioA)[4]);
1433  }
1434  memset(mADCs, 0, sizeof(mADCs));
1435  memset(IDTs, 0, sizeof(IDTs));
1436  Int_t NoTB = 0;
1437  index = NoOfTimeBins*((row-1)*NoOfPads+pad-1);
1438  for (bin = 0; bin < NoOfTimeBins; bin++,index++) {
1439  // Int_t index= NoOfTimeBins*((row-1)*NoOfPads+pad-1)+bin;
1440  // Digits : gain + ped
1441  // GG TF1F *ff = new TF1F("ff","TMath::Sqrt(4.76658e+01*TMath::Exp(-2.87987e-01*(x-1.46222e+01)))",21,56)
1442  Double_t pRMS = pedRMSpad;
1443 #if 0
1444  if (bin >= 21 && bin <= 56) {
1445  pRMS = TMath::Sqrt(pedRMSpad*pedRMSpad + 4.76658e+01*TMath::Exp(-2.87987e-01*(bin-1.46222e+01)));
1446  }
1447 #endif
1448  if (pRMS > 0) {
1449  adc = (Int_t) (SignalSum[index].Sum/gain + gRandom->Gaus(ped,pRMS));
1450  adc = adc - (Int_t) ped;
1451  }
1452  else adc = (Int_t) (SignalSum[index].Sum/gain);
1453  if (adc > 1023) adc = 1023;
1454  if (adc < 1) continue;
1455  SignalSum[index].Adc = adc;
1456  NoTB++;
1457  mADCs[bin] = adc;
1458  IDTs[bin] = SignalSum[index].TrackId;
1459 #ifdef __DEBUG__
1460  if (adc > 3*pRMS) AdcSumBeforeAltro += adc;
1461  if (Debug() > 11 && SignalSum[index].Sum > 0) {
1462  LOG_INFO << "digi R/P/T/I = " << row << " /\t" << pad << " /\t" << bin << " /\t" << index
1463  << "\tSum/Adc/TrackId = " << SignalSum[index].Sum << " /\t"
1464  << SignalSum[index].Adc << " /\t" << SignalSum[index].TrackId << endm;
1465  }
1466 #endif
1467  }
1468  if (! NoTB) continue;
1469  if (mAltro[ioA][sector-1]) {
1470  //#define PixelDUMP
1471 #ifdef PixelDUMP
1472  static Short_t ADCsSaved[__MaxNumberOfTimeBins__];
1473  memcpy(ADCsSaved, mADCs,sizeof(ADCsSaved));
1474 #endif
1475  mAltro[ioA][sector-1]->RunEmulation();
1476 #ifdef PixelDUMP
1477  ofstream *out = new ofstream("digi.dump",ios_base::app);
1478  for (Int_t i = 0; i < __MaxNumberOfTimeBins__; i++) {
1479  if (ADCsSaved[i] > 0 || mADCs[i] > 0) {
1480  LOG_INFO << Form("s %2i r %i p %3i t %3i: %10i => %10i keep %10i",sector,row,pad,i,ADCsSaved[i],mADCs[i],mAltro[ioA][sector-1]->ADCkeep[i]) << endm;
1481  *out << Form("s %2i r %i p %3i t %3i: %10i => %10i keep %10i",sector,row,pad,i,ADCsSaved[i],mADCs[i],mAltro[ioA][sector-1]->ADCkeep[i]) << endl;
1482  }
1483  }
1484  delete out;
1485 #endif
1486  NoTB = 0;
1487  Int_t ADCsum = 0;
1488  for (Int_t i = 0; i < __MaxNumberOfTimeBins__; i++) {
1489  if (mADCs[i] && ! mAltro[ioA][sector-1]->ADCkeep[i]) {mADCs[i] = 0;}
1490  if (mADCs[i]) {
1491  NoTB++;
1492  ADCsum += mADCs[i];
1493 #ifdef __DEBUG__
1494  if (mADCs[i] > 3*pedRMSpad) AdcSumAfterAltro += mADCs[i];
1495  if (Debug() > 12) {
1496  LOG_INFO << "Altro R/P/T/I = " << row << " /\t" << pad << " /\t" << i
1497  << "\tAdc/TrackId = " << mADCs[i] << " /\t" << IDTs[i] << endm;
1498  }
1499 #endif
1500  } else {IDTs[i] = 0;}
1501  }
1502 #ifdef __DEBUG__
1503  if (ADCsum > AdcCut) {
1504  iBreak++;
1505  }
1506 #endif
1507  }
1508  else {
1509  if (St_tpcAltroParamsC::instance()->N(sector-1) < 0) NoTB = AsicThresholds();
1510  }
1511  if (NoTB > 0 && digitalSector) {
1512  digitalSector->putTimeAdc(row,pad,mADCs,IDTs);
1513  }
1514  } // pads
1515 #ifdef __DEBUG__
1516  if (Debug() > 10) {
1517  LOG_INFO << "row = " << row << "\tAdcSumBeforeAltro = " << AdcSumBeforeAltro << "\tAdcSumAfterAltro = " << AdcSumAfterAltro << endm;
1518  }
1519 #endif /* __DEBUG__ */
1520  } // row
1521  return digitalSector;
1522 }
1523 //________________________________________________________________________________
1524 Int_t StTpcRSMaker::AsicThresholds() {
1525  Int_t t1 = 0;
1526  Int_t nSeqLo = 0;
1527  Int_t nSeqHi = 0;
1528  Int_t noTbleft = 0;
1529  for (UInt_t tb = 0; tb < __MaxNumberOfTimeBins__; tb++) {
1530  if (mADCs[tb] <= St_asic_thresholdsC::instance()->thresh_lo()) {
1531  if (! t1) mADCs[tb] = 0;
1532  else {
1533  if (nSeqLo <= St_asic_thresholdsC::instance()->n_seq_lo() ||
1534  nSeqHi <= St_asic_thresholdsC::instance()->n_seq_hi())
1535  for (UInt_t t = t1; t <= tb; t++) mADCs[t] = 0;
1536  else noTbleft += nSeqLo;
1537  }
1538  t1 = nSeqLo = nSeqHi = 0;
1539  }
1540  nSeqLo++;
1541  if (! t1) t1 = tb;
1542  if (mADCs[tb] > St_asic_thresholdsC::instance()->thresh_hi()) {nSeqHi++;}
1543  }
1544  return noTbleft;
1545 }
1546 //________________________________________________________________________________
1547 Double_t StTpcRSMaker::InducedCharge(Double_t s, Double_t h, Double_t ra, Double_t Va, Double_t &t0) {
1548  // Calculate variation of induced charge due to different arrived angles
1549  // alpha = -26 and -70 degrees
1550  LOG_INFO << "wire spacing = " << s << " cm"
1551  << "\tcathode anode gap = " << h << " cm"
1552  << "\tanode wire radius = " << ra << " cm"
1553  << "\tpotential on anode wire = " << Va << " V" << endm;
1554  const Double_t B = 30e-3; // 1/V
1555  const Double_t E0 = 20e3; // V/cm
1556  const Double_t mu = 2.26; // cm**2/V/sec CH4+ mobility
1557  // const Double_t mu = 1.87; // cm**2/V/sec Ar+ mobility
1558  Double_t alpha[2] = {-26., -70.};
1559  Double_t pi = TMath::Pi();
1560  // E.Mathieson (3.2b), V.Chernyatin said that it should be used this (Weber ) approximation 07/09/08
1561  Double_t rc = s/(2*pi)*TMath::Exp(pi*h/s); LOG_INFO << "rc(Cylinder approx) = " << rc << " cm" << endm;
1562  // Double_t rc = 4*h/pi; LOG_INFO << "rc = " << rc << " cm" << endm; // E.Mathieson (4.3), no valid for our case
1563  Double_t C = 1./(2*TMath::Log(rc/ra)); LOG_INFO << "C = " << C << endm;
1564  Double_t E = 2*pi*C*Va/s; LOG_INFO << "E = " << E << " V/cm" << endm;
1565  // Gain variation: M = M0*(1 - k*cos(2*alpha))
1566  Double_t k = 2*B/3.*TMath::Power((pi/E0/s),2)*TMath::Power(C*Va,3); LOG_INFO << "k = " << k << endm;
1567  // Induced charge variation
1568  t0 = ra*ra/(4*mu*C*Va);
1569  LOG_INFO << "t0 = " << 1e9*t0 << " ns" << endm; // E.Mathieson (2.10)
1570  Double_t Tav = t0*h/s/(2*pi*C); LOG_INFO << "Tav = " << 1e9*Tav << " ns" << endm;
1571  // Double_t t = 5*55e-9; LOG_INFO << "t = " << 1e9*t << " ns" << endm;
1572  Double_t t = 180e-9; LOG_INFO << "t = " << 1e9*t << " ns" << endm;
1573  Double_t rp = TMath::Sqrt(1. + t/t0); LOG_INFO << "r' = " << rp << endm;
1574  // qc = rp*ra*sin(alpha)/(2*h) + C/2*log(1 + t/t0) = A*sin(alpha) + B
1575  Double_t Aconstant = rp*ra/(2*h); LOG_INFO << "Aconstant = " << Aconstant << endm;
1576  Double_t Bconstant = C/2*TMath::Log(1 + t/t0); LOG_INFO << "Bconstant = " << Bconstant << endm;
1577  Double_t Gains[2];
1578  for (Int_t i = 0; i < 2; i++) {
1579  Gains[i] = Aconstant*TMath::Sin(pi/180*alpha[i]) + Bconstant;
1580  LOG_INFO << "Gain = " << Gains[i] << " at alpha = " << alpha[i] << " degree" << endm;
1581  }
1582  Double_t GainsAv = TMath::Sqrt(Gains[0]*Gains[1]);
1583  Double_t r = 0;
1584  for (Int_t i = 0; i < 2; i++) {
1585  r = TMath::Log(Gains[i]/GainsAv); LOG_INFO << "Relative gain " << r << " at alpha = " << alpha[i] << endm;
1586  }
1587  return r;
1588 }
1589 //________________________________________________________________________________
1590 Int_t StTpcRSMaker::SearchT(const void *elem1, const void **elem2) {
1591  g2t_tpc_hit_st *value1 = (g2t_tpc_hit_st *) elem1;
1592  g2t_tpc_hit_st *value2 = (g2t_tpc_hit_st *) *elem2;
1593  // sectors
1594  if ((value1->volume_id%100000)/100 != (value2->volume_id%100000)/100)
1595  return (value1->volume_id%100000)/100 - (value2->volume_id%100000)/100;
1596  // track id
1597  if (value1->track_p != value2->track_p) return value1->track_p - value2->track_p;
1598  // pad rows
1599  // if (value1->volume_id%100 != value2->volume_id%100) return value1->volume_id%100 - value2->volume_id%100;
1600  // track length
1601  return (Int_t) 100*(value1->length - value2->length);
1602 }
1603 //________________________________________________________________________________
1604 Int_t StTpcRSMaker::CompareT(const void **elem1, const void **elem2) {
1605  return SearchT(*elem1, elem2);
1606 }
1607 #if 0
1608 //________________________________________________________________________________
1609 Double_t StTpcRSMaker::DriftLength(Double_t x, Double_t y) {
1610  static const Double_t Step = 5e-2;
1611  Double_t r = TMath::Sqrt(x*x + y*y);
1612  if (r < 0.25) return r;
1613  x = TMath::Abs(x);
1614  y = TMath::Abs(y);
1615  Int_t Nstep = 0;
1616  while (x > Step || y > Step) {
1617  Double_t Slope = TMath:SinH(TMath::Pi()*y/s)/TMath:Sin(TMath::Pi()*x/s);
1618  Double_t Co2 = 1./(1. + Slope*Slope);
1619  Double_t Si = TMath::Sqrt(1. - Co2);
1620  Double_t Co = TMath::Sqrt(Co2);
1621  x = TMath::Abs(x - Step*Co);
1622  y = TMath::Abs(y - Step*Si);
1623  NStep++;
1624  }
1625  return NStep*Step;
1626 }
1627 #endif
1628 //________________________________________________________________________________
1629 TF1 *StTpcRSMaker::Fei() {
1630  TF1 *f = (TF1 *) gROOT->GetListOfFunctions()->FindObject("Fei");
1631  if (! f) {
1632  f = new TF1("Fei", feiFunc, 0, 5e-6, 2);
1633  f->SetParNames("t0","T");
1634  f->SetParameters(1.11582e-09, 6e-08);
1635  }
1636  return f;
1637 }
1638 //________________________________________________________________________________
1639 Double_t StTpcRSMaker::fei(Double_t t, Double_t t0, Double_t T) {
1640  static const Double_t xmaxD = 100;// XC11 , XC15
1641  if (t < t0) return 0;
1642  Double_t t01 = xmaxD, t11 = xmaxD;
1643  if (T > 0) {
1644  t11 = (t+t0)/T;
1645  t01 = t0/T;
1646  }
1647  if (t11 > xmaxD) t11 = xmaxD;
1648  if (t01 > xmaxD) t01 = xmaxD;
1649  if (t01 >= t11) return 0;
1650  Double_t ex1 = ROOT::Math::expint(t11);
1651  Double_t ex0 = ROOT::Math::expint(t01);
1652  // return TMath::Exp(-t11)*(ex1 - ex0);
1653  return TMath::Exp(-t11)*(ex1 - ex0);
1654 }
1655 //________________________________________________________________________________
1656 Double_t StTpcRSMaker::shapeEI(Double_t *x, Double_t *par) {// does not work. It is needed to 1/s
1657  Double_t t = x[0];
1658  // static Double_t tmax = 6e-6; // Y
1659  // static Double_t tmax = 2e-6; // Z
1660  // static Double_t tmax = 1e-6; // ZZ
1661  // static Double_t tmax = 1.5e-6; // ZZ1.5
1662  Double_t tmax = par[6];
1663  if (tmax <= 0.0) tmax = 1.5e-6;
1664  Double_t value = 0;
1665  if (t <= 0 || t > tmax) return value;
1666  Double_t t0 = par[0];
1667  Double_t T1 = par[1]; // tau_X
1668  Double_t T2 = par[3]; // tau_C
1669  if (TMath::Abs((T1-T2)/(T1+T2)) < 1e-7) {
1670  return TMath::Max(0.,(t + t0)/T1*fei(t,t0,T1) + TMath::Exp(-t/T1) - 1);
1671  }
1672  if (T2 <= 0) return fei(t,t0,T1);
1673  if (T1 <= 0) return 0;
1674  return T1/(T1 - T2)*(fei(t,t0,T1) - fei(t,t0,T2));
1675 }
1676 //________________________________________________________________________________
1677 Double_t StTpcRSMaker::shapeEI3(Double_t *x, Double_t *par) {// does not work. It is needed to 1/s
1678  Double_t t = x[0];
1679  Double_t value = 0;
1680  if (t <= 0) return value;
1681  Double_t t0 = par[0];
1682  Double_t tau_F = par[1];
1683  Double_t tau_P = par[2];
1684  Double_t tau_I = par[3];
1685  Double_t tau_C = par[5];
1686  Double_t d = 1./tau_P;
1687  Double_t a[3] = {- 1./tau_I, - 1./tau_F, 0};
1688  Double_t A[3] = {(a[0]+d)/(a[0]-a[1]), (a[1]+d)/(a[1]-a[0]), 0};
1689  Int_t N = 2;
1690  if (tau_C > 0) {
1691  N = 3;
1692  a[2] = -1./tau_C;
1693  A[0] = (a[0] + d)/a[0]/(a[0] - a[1])/(a[0] - a[2]);
1694  A[1] = (a[1] + d)/a[1]/(a[1] - a[0])/(a[1] - a[2]);
1695  A[2] = (a[2] + d)/a[2]/(a[2] - a[0])/(a[2] - a[1]);
1696  }
1697  for (Int_t i = 0; i < N; i++) {
1698  value += A[i]*TMath::Exp(a[i]*(t+t0))*(ROOT::Math::expint(-a[i]*(t+t0))-ROOT::Math::expint(-a[i]*t0));
1699  }
1700  return value;
1701 }
1702 //________________________________________________________________________________
1703 Double_t StTpcRSMaker::shapeEI_I(Double_t *x, Double_t *par) { //Integral of shape over time bin
1704  static Double_t sqrt2 = TMath::Sqrt(2.);
1705  Double_t TimeBinWidth = par[2];
1706  static Double_t norm = 1;
1707  static Double_t params0[5] = {0};
1708  Int_t io = (Int_t) par[4];
1709  assert(io >= 0 && io <= 1);
1710  Int_t ok = memcmp(par,params0, sizeof(params0));
1711  if (ok) {
1712  fgTimeShape0[io]->SetParameters(par);
1713  norm = fgTimeShape0[io]->Integral(TimeBinWidth*timeBinMin,TimeBinWidth*timeBinMax);
1714  memcpy(params0, par, sizeof(params0));
1715  }
1716  Double_t t1 = TimeBinWidth*(x[0] - 0.5);
1717  Double_t t2 = t1 + TimeBinWidth;
1718  return sqrt2*fgTimeShape0[io]->Integral(t1,t2)/norm;
1719 }
1720 //________________________________________________________________________________
1721 Double_t StTpcRSMaker::shapeEI3_I(Double_t *x, Double_t *par) { //Integral of shape over time bin
1722  static Double_t sqrt2 = TMath::Sqrt(2.);
1723  Double_t TimeBinWidth = par[4];
1724  static Double_t norm = 1;
1725  static Double_t params3[7] = {0};
1726  Int_t io = (Int_t) par[4];
1727  assert(io >= 0 && io <= 1);
1728  Int_t ok = memcmp(par,params3, sizeof(params3));
1729  if (ok) {
1730  fgTimeShape0[io]->SetParameters(par);
1731  norm = fgTimeShape0[io]->Integral(TimeBinWidth*timeBinMin,TimeBinWidth*timeBinMax);
1732  memcpy(par, params3, sizeof(params3));
1733  }
1734  Double_t t1 = TimeBinWidth*(x[0] - 0.5);
1735  Double_t t2 = t1 + TimeBinWidth;
1736  return sqrt2*fgTimeShape3[io]->Integral(t1,t2)/norm;
1737 }
1738 //________________________________________________________________________________
1739 SignalSum_t *StTpcRSMaker::GetSignalSum(Int_t sector) {
1740  if (! m_SignalSum)
1741  m_SignalSum = (SignalSum_t *) malloc(St_tpcPadConfigC::instance()->numberOfRows(sector)*NoOfPads*NoOfTimeBins*sizeof(SignalSum_t));
1742  return m_SignalSum;
1743 }
1744 //________________________________________________________________________________
1745 SignalSum_t *StTpcRSMaker::ResetSignalSum(Int_t sector) {
1746  GetSignalSum(sector);
1747  memset (m_SignalSum, 0, St_tpcPadConfigC::instance()->numberOfRows(sector)*NoOfPads*NoOfTimeBins*sizeof(SignalSum_t));
1748  return m_SignalSum;
1749 }
1750 //________________________________________________________________________________
1751 Double_t StTpcRSMaker::polya(Double_t *x, Double_t *par) {
1752  return TMath::GammaDist(x[0],par[0],par[1],par[2]);
1753 }
1754 //________________________________________________________________________________
1755 Double_t StTpcRSMaker::Ec(Double_t *x, Double_t *p) {
1756  if (x[0] < p[0]/2 || x[0] > 3.064*p[0]) return 0;
1757  if (x[0] < p[0]) return 1;
1758  return TMath::Power(p[0]/x[0],4);
1759 }
1760 //________________________________________________________________________________
1761 TF1 *StTpcRSMaker::StTpcRSMaker::fEc(Double_t w) {
1762  TF1 *f = new TF1("Ec",Ec,0,3.064*w,1);
1763  f->SetParameter(0,w);
1764  return f;
1765 }
1766 //________________________________________________________________________________
1767 #ifndef WIN32
1768 # define gcomad gcomad_
1769 #else
1770 # define gcomad GCOMAD
1771 #endif
1772 //______________________________________________________________________
1773 extern "C"
1774 {
1775  void type_of_call gcomad(DEFCHARD, Int_t*& DEFCHARL);
1776 }
1777 Float_t StTpcRSMaker::GetCutEle() {
1778  //----------GCBANK
1779  // COMMON/GCBANK/NZEBRA,GVERSN,ZVERSN,IXSTOR,IXDIV,IXCONS,FENDQ(16)
1780  // + ,LMAIN,LR1,WS(KWBANK)
1781  struct Gcbank_t {
1782  Int_t nzebra;
1783  Float_t gversn;
1784  Float_t zversn;
1785  Int_t ixstor;
1786  Int_t ixdiv;
1787  Int_t ixcons;
1788  Float_t fendq[16];
1789  Int_t lmain;
1790  Int_t lr1;
1791  };
1792  Gcbank_t *fGcbank;
1793  gcomad(PASSCHARD("GCBANK"),(int*&) fGcbank PASSCHARL("GCBANK"));
1794 //----------GCLINK
1795 // COMMON/GCLINK/JDIGI ,JDRAW ,JHEAD ,JHITS ,JKINE ,JMATE ,JPART
1796 // + ,JROTM ,JRUNG ,JSET ,JSTAK ,JGSTAT,JTMED ,JTRACK,JVERTX
1797 // + ,JVOLUM,JXYZ ,JGPAR ,JGPAR2,JSKLT
1798 typedef struct {
1799  Int_t jdigi;
1800  Int_t jdraw;
1801  Int_t jhead;
1802  Int_t jhits;
1803  Int_t jkine;
1804  Int_t jmate;
1805  Int_t jpart;
1806  Int_t jrotm;
1807  Int_t jrung;
1808  Int_t jset;
1809  Int_t jstak;
1810  Int_t jgstat;
1811  Int_t jtmed;
1812  Int_t jtrack;
1813  Int_t jvertx;
1814  Int_t jvolum;
1815  Int_t jxyz;
1816  Int_t jgpar;
1817  Int_t jgpar2;
1818  Int_t jsklt;
1819 } Gclink_t;
1820  Gclink_t *fGclink;
1821  gcomad(PASSCHARD("GCLINK"),(int*&) fGclink PASSCHARL("GCLINK"));
1822 
1823  //----------GCCUTS
1824  // COMMON/GCCUTS/CUTGAM,CUTELE,CUTNEU,CUTHAD,CUTMUO,BCUTE,BCUTM
1825  // + ,DCUTE ,DCUTM ,PPCUTM,TOFMAX,GCUTS(5)
1826  struct Gccuts_t {
1827  Float_t cutgam;
1828  Float_t cutele;
1829  Float_t cutneu;
1830  Float_t cuthad;
1831  Float_t cutmuo;
1832  Float_t bcute;
1833  Float_t bcutm;
1834  Float_t dcute;
1835  Float_t dcutm;
1836  Float_t ppcutm;
1837  Float_t tofmax;
1838  Float_t gcuts[5];
1839  };
1840  Gccuts_t *fGccuts;
1841  gcomad(PASSCHARD("GCCUTS"),(int*&) fGccuts PASSCHARL("GCCUTS"));
1842  //----------GCNUM
1843  // COMMON/GCNUM/NMATE ,NVOLUM,NROTM,NTMED,NTMULT,NTRACK,NPART
1844  // + ,NSTMAX,NVERTX,NHEAD,NBIT
1845  struct Gcnum_t {
1846  Int_t nmate;
1847  Int_t nvolum;
1848  Int_t nrotm;
1849  Int_t ntmed;
1850  Int_t ntmult;
1851  Int_t ntrack;
1852  Int_t npart;
1853  Int_t nstmax;
1854  Int_t nvertx;
1855  Int_t nhead;
1856  Int_t nbit;
1857  };
1858  Gcnum_t *fGcnum;
1859  gcomad(PASSCHARD("GCNUM"), (int*&) fGcnum PASSCHARL("GCNUM"));
1860  Int_t *addr;
1861  // Variables for ZEBRA store
1862  gcomad(PASSCHARD("IQ"), addr PASSCHARL("IQ"));
1863  Int_t *fZiq = addr;
1864  gcomad(PASSCHARD("LQ"), addr PASSCHARL("LQ"));
1865  Int_t *fZlq = addr;
1866  Float_t *fZq = (float*)fZiq;
1867  Int_t ITPAR = 2; // IF(CHPAR.EQ.'CUTELE')ITPAR=2
1868  Int_t JTMED = fGclink->jtmed;
1869  for (Int_t i = 1; i <= fGcnum->ntmed; i++) {
1870  Int_t JTM = fZlq[JTMED-i];
1871  if (! JTM) continue;
1872  TString Medium((Char_t *)(&fZiq[JTM+1]),20);
1873  if (!Medium.BeginsWith(TpcMedium)) continue;
1874  Int_t JTMN = fZlq[JTM];
1875  if (! JTMN) continue;
1876  Float_t cutele = fZq[JTMN+ITPAR];
1877  return cutele;
1878  }
1879  LOG_INFO << "StTpcRSMaker::GetCutEle: specific CutEle for medium \"" << TpcMedium.Data() << "\" has not been found. Use default." << endm;
1880  return fGccuts->cutele;
1881 }
1882 //________________________________________________________________________________
1883 Bool_t StTpcRSMaker::TrackSegment2Propagate(g2t_tpc_hit_st *tpc_hitC, g2t_vertex_st *gver, HitPoint_t &TrackSegmentHits) {
1884  static Int_t iBreak = 0;
1885  if (! tpc_hitC) {
1886  iBreak++;
1887  }
1888  Int_t Id = tpc_hitC->track_p;
1889  mNoTpcHitsAll[Id-1]++;
1890  if (tpc_hitC->volume_id < 10000) {
1891  // Account hits which can be splitted
1892  mNoTpcHitsReal[Id-1]++;
1893  }
1894  if (tpc_hitC->de > 0) {
1895  mNSplittedHits = 0;
1896  } else if (! mNSplittedHits) {
1897  mNSplittedHits++;
1898  }
1899  Int_t volId = tpc_hitC->volume_id%10000;
1900  Int_t sector = volId/100;
1901  static StGlobalCoordinate coorG; // ideal
1902  TrackSegmentHits.xyzG =
1903  StGlobalCoordinate(tpc_hitC->x[0],tpc_hitC->x[1],tpc_hitC->x[2]); PrPP(Make,TrackSegmentHits.xyzG);
1904  coorG = TrackSegmentHits.xyzG;
1905  static StTpcLocalCoordinate coorLT; // before do distortions
1906  static StTpcLocalDirection dirLT, BLT;
1907  // calculate row
1908  static StTpcLocalSectorCoordinate coorS;
1909  static StTpcCoordinateTransform transform(gStTpcDb);
1910  transform(coorG, coorS,sector,0); PrPP(Make,coorS);
1911  Int_t row = coorS.fromRow();
1912  transform(coorG, coorLT,sector,row); PrPP(Make,coorLT);
1913  Int_t io = (row <= St_tpcPadConfigC::instance()->numberOfInnerRows(sector)) ? 0 : 1;
1914  TrackSegmentHits.TrackId = Id;
1915  TrackSegmentHits.tpc_hitC = tpc_hitC;
1916 
1917  if (ClusterProfile) {
1918  checkList[io][0]->Fill(TrackSegmentHits.tpc_hitC->x[2],TMath::Abs(TrackSegmentHits.tpc_hitC->de));
1919  checkList[io][1]->Fill(TrackSegmentHits.tpc_hitC->x[2], TrackSegmentHits.tpc_hitC->ds );
1920  }
1921  TrackSegmentHits.sMin = TrackSegmentHits.s - TrackSegmentHits.tpc_hitC->ds;
1922  TrackSegmentHits.sMax = TrackSegmentHits.s;
1923  if (TrackSegmentHits.sMin < msMin) msMin = TrackSegmentHits.sMin;
1924  if (TrackSegmentHits.sMax > msMax) msMax = TrackSegmentHits.sMax;
1925  // move up, calculate field at center of TPC
1926  static Float_t BFieldG[3];
1927  StMagF::Agufld(tpc_hitC->x,BFieldG);
1928  // distortion and misalignment
1929  // replace pxy => direction and try linear extrapolation
1930  StThreeVectorD pxyzG(tpc_hitC->p[0],tpc_hitC->p[1],tpc_hitC->p[2]);
1931  StGlobalDirection dirG(pxyzG.unit()); PrPP(Make,dirG);
1932  StGlobalDirection BG(BFieldG[0],BFieldG[1],BFieldG[2]); PrPP(Make,BG);
1933  transform( dirG, dirLT,sector,row); PrPP(Make,dirLT);
1934  transform( BG, BLT,sector,row); PrPP(Make,BLT);
1935  // Distortions
1936  if (TESTBIT(m_Mode, kDistortion) && StMagUtilities::Instance()) {
1937  Float_t pos[3] = {(Float_t ) coorLT.position().x(), (Float_t ) coorLT.position().y(), (Float_t ) coorLT.position().z()};
1938  Float_t posMoved[3];
1939  StMagUtilities::Instance()->DoDistortion(pos,posMoved,sector); // input pos[], returns posMoved[]
1940  StThreeVector<double> position(posMoved[0],posMoved[1],posMoved[2]);
1941  coorLT.setPosition(position); // after do distortions
1942  transform(coorLT,TrackSegmentHits.xyzG); PrPP(Make,coorLT);
1943  }
1944  // end of distortion
1945  transform(coorLT,TrackSegmentHits.coorLS); PrPP(Make,TrackSegmentHits.coorLS);
1946  transform( dirLT, TrackSegmentHits.dirLS); PrPP(Make,TrackSegmentHits.dirLS);
1947  transform( BLT, TrackSegmentHits.BLS); PrPP(Make,TrackSegmentHits.BLS);
1948  Double_t tof = gver->ge_tof + fgTriggerT0;
1949  // if (! TESTBIT(m_Mode, kNoToflight))
1950  tof += tpc_hitC->tof;
1951  Double_t driftLength = TrackSegmentHits.coorLS.position().z() + tof*gStTpcDb->DriftVelocity(sector); // ,row);
1952 #if 0 /* don't use sign swap for hits behind zGG */
1953  if (driftLength > -1.0 && driftLength <= 0) {
1954  if ((row > St_tpcPadConfigC::instance()->numberOfInnerRows(sector) && driftLength > -gStTpcDb->WirePlaneGeometry()->outerSectorAnodeWirePadPlaneSeparation()) ||
1955  (row <= St_tpcPadConfigC::instance()->numberOfInnerRows(sector) && driftLength > -gStTpcDb->WirePlaneGeometry()->innerSectorAnodeWirePadPlaneSeparation()))
1956  driftLength = TMath::Abs(driftLength);
1957  }
1958 #endif
1959  TrackSegmentHits.coorLS.position().setZ(driftLength); PrPP(Make,TrackSegmentHits.coorLS);
1960  // useT0, don't useTau
1961  transform(TrackSegmentHits.coorLS,TrackSegmentHits.Pad,kFALSE,kFALSE); // don't use T0, don't use Tau
1962  PrPP(Make,TrackSegmentHits.Pad);
1963  return kTRUE;
1964 }
1965 //________________________________________________________________________________
1966 void StTpcRSMaker::GenerateSignal(HitPoint_t &TrackSegmentHits, Int_t sector, Int_t rowMin, Int_t rowMax, Double_t sigmaJitterT, Double_t sigmaJitterX) {
1967  static StTpcCoordinateTransform transform(gStTpcDb);
1968  SignalSum_t *SignalSum = GetSignalSum(sector);
1969  for(Int_t row = rowMin; row <= rowMax; row++) {
1970  // if (St_tpcPadConfigC::instance()->numberOfRows(sector) == 45) { // ! iTpx
1971  if ( ! St_tpcAnodeHVavgC::instance()->livePadrow(sector,row)) continue;
1972  // }
1973  Int_t io = (row <= St_tpcPadConfigC::instance()->numberOfInnerRows(sector)) ? 0 : 1;
1974  StTpcLocalSectorCoordinate xyzW(xOnWire, yOnWire, zOnWire, sector, row);
1975  TF1F *ShaperResponse = mShaperResponses[io][sector-1];
1976 
1977  static StTpcPadCoordinate Pad;
1978  transform(xyzW,Pad,kFALSE,kFALSE); // don't use T0, don't use Tau
1979  Float_t bin = Pad.timeBucket();//L - 1; // K
1980  Int_t binT = TMath::Nint(bin); //L bin;//K TMath::Nint(bin);// J bin; // I TMath::Nint(bin);
1981  if (binT < 0 || binT >= NoOfTimeBins) continue;
1982  Double_t dT = bin - binT + St_TpcResponseSimulatorC::instance()->T0offset();
1983  dT += (row <= St_tpcPadConfigC::instance()->numberOfInnerRows(sector)) ?
1984  St_TpcResponseSimulatorC::instance()->T0offsetI() :
1985  St_TpcResponseSimulatorC::instance()->T0offsetO();
1986  if (sigmaJitterT) dT += gRandom->Gaus(0,sigmaJitterT); // #1
1987 #if 1
1988  Double_t dely = {transform.yFromRow(sector,row)-yOnWire};
1989  Double_t localYDirectionCoupling = mChargeFraction[io][sector-1]->GetSaveL(&dely);
1990 #else
1991  Int_t idWire = TMath::Abs(TMath::Nint((transform.yFromRow(sector,row)-yOnWire)/anodeWirePitch));
1992  if (idWire > 6) continue;
1993  Double_t localYDirectionCoupling = mLocalYDirectionCoupling[io][sector-1][idWire];
1994 #endif
1995  if (ClusterProfile) {
1996  checkList[io][10]->Fill(TrackSegmentHits.xyzG.position().z(),localYDirectionCoupling);
1997  }
1998  if(localYDirectionCoupling < minSignal) continue;
1999  Float_t padX = Pad.pad();
2000  Int_t CentralPad = TMath::Nint(padX);
2001  if (CentralPad < 1) continue;
2002  Int_t PadsAtRow = St_tpcPadConfigC::instance()->numberOfPadsAtRow(sector,row);
2003  if(CentralPad > PadsAtRow) continue;
2004  Int_t DeltaPad = TMath::Nint(mPadResponseFunction[io][sector-1]->GetXmax()) + 1;
2005  Int_t padMin = TMath::Max(CentralPad - DeltaPad ,1);
2006  Int_t padMax = TMath::Min(CentralPad + DeltaPad ,PadsAtRow);
2007  Int_t Npads = TMath::Min(padMax-padMin+1, kPadMax);
2008  Double_t xPadMin = padMin - padX;
2009  static Double_t XDirectionCouplings[kPadMax];
2010  static Double_t TimeCouplings[kTimeBacketMax];
2011  mPadResponseFunction[io][sector-1]->GetSaveL(Npads,xPadMin,XDirectionCouplings);
2012  // Double_t xPad = padMin - padX;
2013  for(Int_t pad = padMin; pad <= padMax; pad++) {
2014  if ( ! StDetectorDbTpcRDOMasks::instance()->isRowOn(sector,row,pad)) continue;
2015  Double_t gain = QAv*mGainLocal;
2016  Double_t dt = dT;
2017  // if (St_tpcPadConfigC::instance()->numberOfRows(sector) ==45 && ! TESTBIT(m_Mode, kGAINOAtALL)) {
2018  if (! TESTBIT(m_Mode, kGAINOAtALL)) {
2019  Double_t GC = St_tpcPadGainT0BC::instance()->Gain(sector,row,pad);
2020  if (GC <= 0.0) continue;
2021  gain /= GC;
2022  dt -= St_tpcPadGainT0BC::instance()->T0(sector,row,pad);
2023  }
2024  if (ClusterProfile) {
2025  checkList[io][12]->Fill(TrackSegmentHits.xyzG.position().z(),gain);
2026  GainHist[1]->Fill(sector,row,gain);
2027  }
2028  // Double_t localXDirectionCoupling = localXDirectionCouplings[pad-padMin];
2029  Double_t localXDirectionCoupling = gain*XDirectionCouplings[pad-padMin];
2030  if (localXDirectionCoupling < minSignal) continue;
2031  if (ClusterProfile) {
2032  checkList[io][13]->Fill(TrackSegmentHits.xyzG.position().z(),localXDirectionCoupling);
2033  }
2034  Double_t XYcoupling = localYDirectionCoupling*localXDirectionCoupling;
2035  if (ClusterProfile) {
2036  checkList[io][14]->Fill(TrackSegmentHits.xyzG.position().z(),XYcoupling);
2037  }
2038  if(XYcoupling < minSignal) continue;
2039  Int_t bin_low = TMath::Max(0 ,binT + TMath::Nint(dt+ShaperResponse->GetXmin()-0.5));
2040  Int_t bin_high = TMath::Min(NoOfTimeBins-1,binT + TMath::Nint(dt+ShaperResponse->GetXmax()+0.5));
2041  Int_t index = NoOfTimeBins*((row-1)*NoOfPads+pad-1)+bin_low;
2042  Int_t Ntbks = TMath::Min(bin_high-bin_low+1, kTimeBacketMax);
2043  Double_t tt = -dt + (bin_low - binT);
2044  ShaperResponse->GetSaveL(Ntbks,tt,TimeCouplings);
2045  for(Int_t itbin=bin_low;itbin<=bin_high;itbin++, index++){
2046  Double_t signal = XYcoupling*TimeCouplings[itbin-bin_low];
2047  if (signal < minSignal) continue;
2048 #ifdef __DEBUG__
2049  static Int_t iBreak = 0;
2050  if (TMath::IsNaN(signal) || TMath::IsNaN(SignalSum[index].Sum)) {
2051  iBreak++;
2052  }
2053 #endif
2054  TotalSignalInCluster += signal;
2055  SignalSum[index].Sum += signal;
2056  if (ClusterProfile) {
2057  if (pad >= pad0 && pad < pad0 + kPadMax &&
2058  itbin >= tbk0 && itbin < tbk0 + kTimeBacketMax) {
2059  padsdE[pad-pad0] += signal;
2060  tbksdE[itbin-tbk0] += signal;
2061  }
2062  }
2063  rowsdE[row-1] += signal;
2064  rowsdEH[row-1] += signal;
2065  if ( TrackSegmentHits.TrackId ) {
2066  if (! SignalSum[index].TrackId ) SignalSum[index].TrackId = TrackSegmentHits.TrackId;
2067  else // switch TrackId, works only for 2 tracks, more tracks ?
2068  if ( SignalSum[index].TrackId != TrackSegmentHits.TrackId && SignalSum[index].Sum < 2*signal)
2069  SignalSum[index].TrackId = TrackSegmentHits.TrackId;
2070  }
2071 #ifdef __DEBUG__
2072  if (Debug() > 13 && (SignalSum[index].Sum > 0 || ! TMath::Finite(SignalSum[index].Sum)) ) {
2073  LOG_INFO << "simu row = " << TrackSegmentHits.tpc_hitC->volume_id%100 << "\tR/P/T/I = " << row << " /\t" << pad << " /\t" << itbin << " /\t" << index
2074  << "\tSum/Adc/TrackId = " << SignalSum[index].Sum << " /\t"
2075  << SignalSum[index].Adc << " /\t" << SignalSum[index].TrackId
2076  << "\tsignal = " << signal
2077  << "\trow Min/Max = " << rowMin << "/" << rowMax
2078  << endm;
2079  if (! TMath::Finite(SignalSum[index].Sum)) {
2080  LOG_INFO << "Not Finite" << endm;
2081  }
2082  }
2083 #endif /* __DEBUG__ */
2084  } // time
2085  } // pad limits
2086  } // row limits
2087 }
2088 //________________________________________________________________________________
2089 Double_t StTpcRSMaker::dEdxCorrection(HitPoint_t &TrackSegmentHits) {
2090  Double_t dEdxCor = 1;
2091  if (m_TpcdEdxCorrection) {
2092  // dEdxCor = -1;
2093  Double_t dStep = TMath::Abs(TrackSegmentHits.tpc_hitC->ds);
2094  dEdxY2_t CdEdx;
2095  memset (&CdEdx, 0, sizeof(dEdxY2_t));
2096  CdEdx.DeltaZ = 5.2;
2097  CdEdx.QRatio = -2;
2098  CdEdx.QRatioA = -2.;
2099  CdEdx.QSumA = 0;
2100  CdEdx.sector = TrackSegmentHits.Pad.sector();
2101  CdEdx.row = TrackSegmentHits.Pad.row();
2102  CdEdx.pad = TMath::Nint(TrackSegmentHits.Pad.pad());
2103  CdEdx.edge = CdEdx.pad;
2104  if (CdEdx.edge > 0.5*St_tpcPadConfigC::instance()->numberOfPadsAtRow(CdEdx.sector,CdEdx.row))
2105  CdEdx.edge += 1 - St_tpcPadConfigC::instance()->numberOfPadsAtRow(CdEdx.sector,CdEdx.row);
2106  CdEdx.F.dE = 1;
2107 #if 0
2108  CdEdx.dCharge = 0;
2109  Int_t p1 = tpcHit->minPad();
2110  Int_t p2 = tpcHit->maxPad();
2111  Int_t t1 = tpcHit->minTmbk();
2112  Int_t t2 = tpcHit->maxTmbk();
2113  CdEdx.rCharge= 0.5*m_TpcdEdxCorrection->Adc2GeV()*TMath::Pi()/4.*(p2-p1+1)*(t2-t1+1);
2114  if (TESTBIT(m_Mode, kEmbeddingShortCut) &&
2115  (tpcHit->idTruth() && tpcHit->qaTruth() > 95)) CdEdx.lSimulated = tpcHit->idTruth();
2116 #endif
2117  CdEdx.F.dx = dStep;
2118  CdEdx.xyz[0] = TrackSegmentHits.coorLS.position().x();
2119  CdEdx.xyz[1] = TrackSegmentHits.coorLS.position().y();
2120  CdEdx.xyz[2] = TrackSegmentHits.coorLS.position().z();
2121  Double_t probablePad = St_tpcPadConfigC::instance()->numberOfPadsAtRow(CdEdx.sector,CdEdx.row)/2;
2122  Double_t pitch = (CdEdx.row <= St_tpcPadConfigC::instance()->numberOfInnerRows(CdEdx.sector)) ?
2123  St_tpcPadConfigC::instance()->innerSectorPadPitch(CdEdx.sector) :
2124  St_tpcPadConfigC::instance()->outerSectorPadPitch(CdEdx.sector);
2125  Double_t PhiMax = TMath::ATan2(probablePad*pitch, St_tpcPadConfigC::instance()->radialDistanceAtRow(CdEdx.sector,CdEdx.row));
2126  CdEdx.PhiR = TMath::ATan2(CdEdx.xyz[0],CdEdx.xyz[1])/PhiMax;
2127  CdEdx.xyzD[0] = TrackSegmentHits.dirLS.position().x();
2128  CdEdx.xyzD[1] = TrackSegmentHits.dirLS.position().y();
2129  CdEdx.xyzD[2] = TrackSegmentHits.dirLS.position().z();
2130  CdEdx.ZdriftDistance = CdEdx.xyzD[2];
2131  CdEdx.zG = CdEdx.xyz[2];
2132  if (St_trigDetSumsC::instance()) CdEdx.Zdc = St_trigDetSumsC::instance()->zdcX();
2133  CdEdx.ZdriftDistance = TrackSegmentHits.coorLS.position().z(); // drift length
2134  St_tpcGas *tpcGas = m_TpcdEdxCorrection->tpcGas();
2135  if (tpcGas)
2136  CdEdx.ZdriftDistanceO2 = CdEdx.ZdriftDistance*(*tpcGas)[0].ppmOxygenIn;
2137  Int_t iok = m_TpcdEdxCorrection->dEdxCorrection(CdEdx);
2138  if (! iok) {
2139  dEdxCor = CdEdx.F.dE;
2140  } else {
2141  dEdxCor = 0; // reject hits with wrong correction
2142  }
2143  }
2144  return dEdxCor;
2145 }
2146 //________________________________________________________________________________
2147 #undef PrPP
StMaker::m_Mode
Int_t m_Mode
counters
Definition: StMaker.h:81
StMaker::AddData
virtual void AddData(TDataSet *data, const char *dir=".data")
User methods.
Definition: StMaker.cxx:332
Altro::ConfigAltro
void ConfigAltro(int ONBaselineCorrection1, int ONTailcancellation, int ONBaselineCorrection2, int ONClipping, int ONZerosuppression)
Configures which modules of the Altro should be on.
Definition: Altro.cxx:72
TF1F
Definition: TF1F.h:5
g2t_tpc_hit
Definition: g2t_tpc_hit.hh:34
TTableSorter::GetIndex
Int_t GetIndex(UInt_t sortedIndex) const
returns the original index of the row by its sorted index
Definition: TTableSorter.cxx:643
tof
Definition: tof.h:15
SignalSum_t
Definition: StTpcRSMaker.h:33
StEpdHit::adc
int adc() const
ADC value [0,4095].
Definition: StEpdHit.h:149
dEdxY2_t
Definition: StTpcdEdxCorrection.h:141
HitPoint_t
Definition: StTpcRSMaker.cxx:70
Altro::ConfigTailCancellationFilter
void ConfigTailCancellationFilter(int K1, int K2, int K3, int L1, int L2, int L3)
Configures the Tail Cancellation Filter (TCF) Module.
Definition: Altro.cxx:118
StMagUtilities::DoDistortion
virtual void DoDistortion(const Float_t x[], Float_t Xprime[], Int_t Sector=-1)
Main Entry Point for requests to DO the E and B field distortions (for simulations) ...
Definition: StMagUtilities.cxx:1418
Altro.h
This the header File for the Altro class.
StTpcRSMaker::Make
virtual Int_t Make()
Definition: StTpcRSMaker.cxx:668
StEpdHit
Stores information for tiles in STAR Event Plane Detector.
Definition: StEpdHit.h:43
Altro::PrintParameters
void PrintParameters()
Prints the set Parameters, if module is configured.
Definition: Altro.cxx:182
StMaker::GetName
virtual const char * GetName() const
special overload
Definition: StMaker.cxx:237
kStWarn
Definition: Stypes.h:42
kStOK
Definition: Stypes.h:40
Altro::RunEmulation
void RunEmulation()
Runs the emulation of all configured Modules.
Definition: Altro.cxx:238
StEpdHit::tile
int tile() const
tile on the supersector [1,31]
Definition: StEpdHit.h:148
StTpcRSMaker::Finish
virtual Int_t Finish()
Definition: StTpcRSMaker.cxx:151
StMaker::Finish
virtual Int_t Finish()
Definition: StMaker.cxx:776
Altro
Definition: Altro.h:22
Step
Definition: AgUStep.h:71
Altro::ConfigZerosuppression
void ConfigZerosuppression(int Threshold, int MinSamplesaboveThreshold, int Presamples, int Postsamples)
Configures the Zero Suppression Module (ZSU)
Definition: Altro.cxx:165
StEpdHit::tac
int tac() const
TAC value [0,4095].
Definition: StEpdHit.h:150
Generated by   doxygen 1.8.5