StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StTrsMaker.cxx
1 // $Id: StTrsMaker.cxx,v 1.94 2018/12/07 17:59:39 genevb Exp $
2 //
3 
4 // $Log: StTrsMaker.cxx,v $
5 // Revision 1.94 2018/12/07 17:59:39 genevb
6 // Avoid using 0 as random seed when intent is to fix seed
7 //
8 // Revision 1.93 2018/06/21 22:23:08 perev
9 // TpcGroup fixes
10 //
11 // Revision 1.92 2018/02/20 22:45:53 smirnovd
12 // Revert "Changes from Irakli's directory to make the code compile"
13 //
14 // Revision 1.90 2014/12/16 04:09:08 perev
15 // In Make() SetSeed() for StTrsRand called. Input seed tkane from g2t_event::ge_rndm[2]
16 // It is made for reproducion event after skip previous ones.
17 // Simulation of event is independent now.
18 //
19 // Revision 1.89 2013/03/26 15:56:00 genevb
20 // Replace agufld(x,b) with direct call to StarMagField::Instance()->BField(x,b)
21 //
22 // Revision 1.88 2013/02/18 16:31:13 fisyak
23 // gufld => agufld
24 //
25 // Revision 1.87 2011/01/18 14:40:15 fisyak
26 // Clean up TpcDb interfaces and Tpc coordinate transformation
27 //
28 // Revision 1.86 2010/01/27 21:33:08 perev
29 // Account Prompt hits
30 //
31 // Revision 1.85 2009/11/03 14:34:19 fisyak
32 // Remove default in zFromTB
33 //
34 // Revision 1.84 2009/07/28 14:40:46 fisyak
35 // Comment out cut on TPC fiducial volume
36 //
37 // Revision 1.83 2008/10/13 19:56:09 fisyak
38 // Account that Z-offset is sector dependent
39 //
40 // Revision 1.82 2008/06/20 15:00:57 fisyak
41 // move from StTrsData to StTpcRawData
42 //
43 // Revision 1.81 2007/07/12 20:25:04 fisyak
44 // Use StarLogger, use time of flight, fix cluster shape
45 //
46 // Revision 1.80 2007/04/28 17:57:27 perev
47 // Redundant StChain.h removed
48 //
49 // Revision 1.79 2005/12/12 21:00:11 perev
50 // 3 random generators ==> 1
51 //
52 // Revision 1.78 2005/09/09 22:12:48 perev
53 // Bug fix + IdTruth added
54 //
55 // Revision 1.76 2005/07/19 22:20:52 perev
56 // Bug fix
57 //
58 // Revision 1.75 2004/03/31 19:45:57 jeromel
59 // Initialize data member (valgrind report)
60 //
61 // Revision 1.74 2003/12/24 13:44:42 fisyak
62 // Add (GEANT) track Id information in Trs; propagate it via St_tpcdaq_Maker; account interface change in StTrsZeroSuppressedReaded in StMixerMaker
63 //
64 // Revision 1.73 2003/09/02 17:59:14 perev
65 // gcc 3.2 updates + WarnOff
66 //
67 // Revision 1.72 2003/05/02 23:54:18 hardtke
68 // Allow user to adjust mNormalFactor (i.e. Fudge Factor)
69 //
70 // Revision 1.71 2003/04/30 20:38:56 perev
71 // Warnings cleanup. Modified lines marked VP
72 //
73 // Revision 1.70 2003/04/10 21:30:34 hardtke
74 // only call Init Run once per job
75 //
76 // Revision 1.69 2002/02/05 22:21:27 hardtke
77 // Move Init code to InitRun
78 //
79 // Revision 1.68 2001/11/21 01:49:25 long
80 // add log message for 3/2001 long;
81 // adding in Make():
82 // < //update the gain for very events
83 // < mWireHistogram->setGasGainInnerSector(mSlowControlDb->innerSectorGasGain());
84 // < mWireHistogram->setGasGainOuterSector(mSlowControlDb->outerSectorGasGain());
85 //
86 // Revision 1.67 2001/11/14 22:24:32 jeromel
87 // mAllTheData deleted in maker. Was in maker AND in Event (destroyed
88 // twice) causing a crash.
89 //
90 // Revision 1.66 2001/03/20 01:44:24 perev
91 // drift velocity print added
92 //
93 // Revision 1.65 2001/03/13 22:09:19 long
94 // *** empty log message ***
95 //
96 // Revision 1.64 2001/02/15 21:34:45 perev
97 // clear improved
98 //
99 // Revision 1.63 2000/08/11 02:24:27 long
100 // comment out sampleAnalogSignal();
101 //
102 // Revision 1.62 2000/08/05 18:32:58 long
103 // add check for the No.of input hits. if(no_tpc_hits<1)return kStOK;
104 //
105 // Revision 1.61 2000/08/04 21:03:57 perev
106 // Leaks + Clear() cleanup
107 //
108 // Revision 1.60 2000/08/04 03:33:45 long
109 // mMiniSegmentLength(3.*millimeter)--->mMiniSegmentLength(4.*millimeter)
110 //
111 // Revision 1.59 2000/07/30 02:36:35 long
112 // add dx(d[0]),dy(d[1]),dz calculation instead of just ds
113 //
114 // Revision 1.58 2000/07/21 22:31:11 calderon
115 // Added checks at the beginning of Make() to avoid
116 // seg faults when dereferencing an invalid pointer. This
117 // happens when no g2t tables are present. For the case of
118 // embedding this might happen for the very low multiplicity
119 // events where none of the embedded monte carlo tracks reach
120 // a sensitive volume.
121 //
122 // Revision 1.57 2000/06/25 23:58:28 fisyak
123 // Remove params
124 //
125 // Revision 1.56 2000/03/24 02:43:46 long
126 // comment out "hitMomentum.setZ(-(tpc_hit->p[2]*GeV)); when bsectorOfHit>12"
127 //
128 // Revision 1.55 2000/03/15 23:33:30 calderon
129 // Modified Finish() to properly take care of pointers when
130 // using the mixer chain.
131 //
132 // Revision 1.54 2000/02/24 16:30:41 long
133 // 1) modified for field on case
134 // //2) modified for pileup pp events ---Balewski
135 //
136 // Revision 1.54 2000/02/23 01:34:06 long
137 //1) modified for field on case
138 //2) modified for pileup pp events ---Balewski
139 // Revision 1.53 2000/02/10 01:21:27 calderon
140 // Switch to use StTpcDb.
141 // Coordinates checked for consistency.
142 // Fixed problems with StTrsIstream & StTrsOstream.
143 //
144 // Revision 1.52 2000/01/27 22:56:58 calderon
145 // add Magnetic field Db instantiation when using StTpcDb. The
146 // Magnetic Field is still obtained through gufld, so the call
147 // is the same as for the StSimpleMagneticField case. StTpcDb
148 // is still not used until further tests.
149 //
150 // Revision 1.51 2000/01/10 23:11:07 lasiuk
151 // Include MACROS for compatibility with SUN CC5.0
152 //
153 // Revision 1.50 1999/11/12 01:34:06 long
154 // 1)fix a bug for track with z<0
155 // 2) turn on the switch on pseudo pad row and set the signal threshold to
156 // 0 for the time being
157 //
158 // Revision 1.49 1999/11/11 19:42:24 calderon
159 // Add #ifdef HISTOGRAM for Ntuple Diagnostics.
160 // Use ROOT_DATABASE_PARAMETERS. As soon as Jeff and Dave give Ok,
161 // we will switch to TPC_DATABASE_PARAMETERS.
162 //
163 // Revision 1.48 1999/11/10 15:45:39 calderon
164 // Made changes to reduce timing, including:
165 // Made coordinate transfrom a data member of StTrsAnalogSignalGenerator
166 // Added upper-lower bound instead of symmetric cut.
167 // Revived checking if signal is above threshold.
168 //
169 // Revision 1.47 1999/11/05 22:10:13 calderon
170 // Added Clear() method in StTrsMaker.
171 // Removed StTrsUnpacker from maker.
172 // Added StTrsZeroSuppressedReader and StTrsZeroSuppressedReader
173 // for DAQ type interface to event data.
174 // Made private copy constructor and operator= in classes that needed it.
175 // Renamed the DigitalSignalGenerators: Fast -> Old, Parameterized -> Fast,
176 // and the default to use is the new "Fast" which has the 10-8 bit conversion.
177 // Removed vestigial loop in AnalogSignalGenerator.
178 // Added Time diagnostics.
179 // Added trsinterface.pdf in doc/ directory.
180 // Write version of data format in .trs data file.
181 //
182 // Revision 1.46 1999/10/21 23:50:32 calderon
183 // Changes are
184 // -string for gasDb
185 // -modified mDeltaRow to 0 if using SlowAnalogSignalGenerator
186 //
187 // Revision 1.45 1999/10/13 17:36:53 calderon
188 // Fixed path to find SimpleDb files.
189 //
190 // Revision 1.44 1999/10/12 22:51:17 long
191 // fix a bug in switching from sector to sector
192 //
193 // Revision 1.43 1999/10/11 23:54:31 calderon
194 // Version with Database Access and persistent file.
195 // Not fully tested due to problems with cons, it
196 // doesn't find the local files at compile time.
197 // Yuri suggests forcing commit to work directly with
198 // files in repository.
199 //
200 // Revision 1.42 1999/10/04 16:13:14 long
201 // minor change on how to loop over geant hits
202 //
203 //
204 // Revision 1.40 1999/09/24 01:23:29 fisyak
205 // Reduced Include Path
206 //
207 // Revision 1.39 1999/07/20 02:16:58 lasiuk
208 // bring in line with new options (TSS algorithms)
209 //
210 // Revision 1.38 1999/07/15 13:57:30 perev
211 // cleanup
212 //
213 // Revision 1.37 1999/07/09 03:45:50 lasiuk
214 // set switch for the wireHistogram to determine the use
215 // of a Gaussian or exponential gas gain based on a collection
216 // or single electrons
217 //
218 // Revision 1.36 1999/06/17 19:04:40 lasiuk
219 // Rotate the momentum the same way that the position is rotated
220 //
221 // Revision 1.35 1999/06/16 14:26:49 fisyak
222 // Add flags for egcs on Solaris
223 //
224 // Revision 1.34 1999/05/28 02:55:44 lasiuk
225 // change default settings for testing. Only in the Maker!
226 // Remove histograms
227 //
228 // Revision 1.33 1999/04/29 00:15:10 lasiuk
229 // make sure to clean up pointers that are store
230 // in the StTrsRawDataEvent() because they are allocated
231 // for each event. This is done with StTrsRawDataEvent::clear()
232 // which has been added.
233 //
234 // Revision 1.32 1999/04/27 19:38:37 lasiuk
235 // Bfield units in kG from STAR database
236 //
237 // Revision 1.31 1999/04/23 19:18:08 lasiuk
238 // change magnetic field initialization to use gufld from GEANT
239 //
240 // Revision 1.30 1999/04/07 01:02:41 lasiuk
241 // use gas gain from db
242 //
243 // Revision 1.29 1999/03/28 02:59:11 perev
244 // Put interfaces to .const area
245 //
246 // Revision 1.28 1999/03/24 22:16:24 lasiuk
247 // ROOT deletes the dataSet so you have to make
248 // a new one each time thru Make(). Don't check the
249 // pointers!
250 //
251 // Revision 1.27 1999/03/23 03:37:26 lasiuk
252 // incorporate ROOT dataSets for DB initialization
253 // move construction and destruction of "mAllthedata"
254 //
255 // Revision 1.26 1999/03/20 20:07:56 fisyak
256 // Add access to DataSet with parameters
257 //
258 // Revision 1.25 1999/03/20 03:23:57 lasiuk
259 // setMiniSegmentLength()
260 // setFirstSectorToProcess()
261 // setLastSectorToProcess()
262 //
263 // Revision 1.24 1999/03/19 13:27:11 lasiuk
264 // change sectors to process 1 ONLY!
265 //
266 // Revision 1.23 1999/03/17 17:11:12 lasiuk
267 // comment out data set deletion for SL99d
268 //
269 // Revision 1.22 1999/03/16 02:03:46 lasiuk
270 // Add Finish(); correct breakNumber calculation; Use C++ casts;
271 // Change defaults (again) including P10;
272 // New mechanism for selecting which sectors to process;
273 // Add flag for processing pseudo-pad rows
274 //
275 // Revision 1.21 1999/03/15 02:52:26 perev
276 // new Maker schema
277 //
278 // Revision 1.20 1999/03/02 17:50:39 lasiuk
279 // for testing/defaults/geantPID
280 //
281 // Revision 1.19 1999/02/24 16:59:24 lasiuk
282 // turn off histos
283 //
284 // Revision 1.18 1999/02/24 16:57:02 lasiuk
285 // x --> -x for sectors>12
286 //
287 // Revision 1.17 1999/02/23 19:15:05 lasiuk
288 // Change in the rotation angle for GEANT global coordinates
289 //
290 // Revision 1.16 1999/02/23 01:10:42 lasiuk
291 // 1st production version
292 //
293 // Revision 1.15 1999/02/17 17:02:16 lasiuk
294 // streamline for production
295 // remove debug.
296 // switch for #sectors to be processed
297 //
298 // Revision 1.14 1999/02/16 18:15:40 fisyak
299 // Check in the latest updates to fix them
300 //
301 // Revision 1.13 1999/02/15 03:32:09 lasiuk
302 // coordinate system for input data is global
303 // deltapad(1)
304 //
305 // Revision 1.12 1999/02/12 01:26:51 lasiuk
306 // Limit debug output
307 //
308 // Revision 1.11 1999/02/10 18:01:31 lasiuk
309 // remove debug/sleep
310 // ROOT passing
311 // set defaults
312 //
313 // Revision 1.10 1999/02/10 04:30:02 lasiuk
314 // add unpacker and rawevent as data members/ passed by dataset
315 //
316 // Revision 1.9 1999/02/05 23:08:34 fisyak
317 // Add Valery's update of DataSet
318 //
319 // Revision 1.7 1999/01/28 02:46:09 lasiuk
320 // SUN compile with new GEANT interface
321 //
323 // //
324 // StTrsMaker class for Makers //
325 // //
326 //
327 // You must select a data base initializer method
328 // When using TPC_DATABASE, change also definition in
329 // StTrsChargeSegment.cc
330 // StTrsAnalogSignalGenerator.hh
332 
333 #include "StTrsMaker.h"
334 #include "TDataSetIter.h"
335 #include "TObjectSet.h"
336 
337 #define ST_TRS_RANDOM_SRC
338 #include "StTrsRandom.hh"
339 
340 #include <Stiostream.h>
341 #include <unistd.h> // needed for access()/sleep()
342 #include <math.h>
343 #include <string>
344 #include <algorithm>
345 #include <vector>
346 #include <list>
347 #include <utility> // pair
348 #include <algorithm> // min() max()
349 #include <ctime> // time functions
350 
351 #if defined (__SUNPRO_CC) && __SUNPRO_CC >= 0x500
352 using std::list;
353 using std::vector;
354 using std::string;
355 using std::min;
356 using std::max;
357 #endif
358 
359 #include "StMem.h"
360 // SCL
361 #include "StGlobals.hh"
362 #include "Randomize.h"
363 
364 #include "TNtuple.h"
365 #include "TFile.h"
366 
367 
368 // TRS
369 // DataBase Initialization
370 // Dave's Header file
371 #include "StTpcDb/StTpcDb.h"
372 
373 #include "StTpcDbGeometry.hh"
374 #include "StTpcDbSlowControl.hh"
375 #include "StTpcDbElectronics.hh"
376 //#include "StDbMagneticField.hh" // To be done
377 #include "StarMagField.h"
378 
379 #include "StSimpleMagneticField.hh"
380 #include "StTrsDeDx.hh"
381 
382 #include "electronicsDataSet.h"
383 #include "geometryDataSet.h"
384 #include "slowcontrolDataSet.h"
385 
386 // processes
387 #include "StTrsFastChargeTransporter.hh"
388 #include "StTrsSlowAnalogSignalGenerator.hh"
389 //***************************** SEE the default options for Analog Signal generator
390 #include "StTrsParameterizedAnalogSignalGenerator.hh"
391 //*****************************
392 #include "StTrsFastDigitalSignalGenerator.hh"
393 // #include "StTrsOldDigitalSignalGenerator.hh"
394 
395 // containers
396 #include "StTrsChargeSegment.hh"
397 #include "StTrsMiniChargeSegment.hh"
398 #include "StTrsAnalogSignal.hh"
399 #include "StTrsWireBinEntry.hh"
400 #include "StTrsWireHistogram.hh"
401 
402 #include "StTrsSector.hh"
403 #include "StTrsDigitalSector.hh"
404 
405 // outPut Data--decoder
406 #include "StTrsDetectorReader.hh"
407 #include "StTrsZeroSuppressedReader.hh"
408 #include "StDaqLib/GENERIC/EventReader.hh"
409 #include "StSequence.hh"
410 
411 // g2t tables
412 #include "tables/St_g2t_event_Table.h"
413 #include "tables/St_g2t_tpc_hit_Table.h"
414 #include "tables/St_g2t_track_Table.h"
415 #include "tables/St_g2t_vertex_Table.h"
416 #include "StMessMgr.h"
417 #include "TString.h"
418 #define PILEUP_ON (m_Mode )
419 
420 //#define VERBOSE 1
421 //#define ivb if(VERBOSE)
422 
423 static const char rcsid[] = "$Id: StTrsMaker.cxx,v 1.94 2018/12/07 17:59:39 genevb Exp $";
424 
425 ClassImp(electronicsDataSet)
426 ClassImp(geometryDataSet)
427 ClassImp(slowcontrolDataSet)
428 ClassImp(StTrsMaker)
429 
430 
431 
432 
433 StTrsMaker::StTrsMaker(const char *name):StMaker(name)
434 {
435  memset(mBeg,0,mEnd-mBeg+1); //Zero all
436 
437  mMiniSegmentLength =(4.*millimeter); // test trial,Hui Long
438  mFirstSectorToProcess =(1);
439  mLastSectorToProcess =(24);
440  mUseParameterizedSignalGenerator =(1); // test trial,Hui Long
441  mNormalFactor = 1.25;
442 }
443 
444 StTrsMaker::~StTrsMaker() { /* nopt */ }
445 
446 
447 Int_t StTrsMaker::Init()
448 {
449  mAllTheData = 0;
450  int seed = IAttr("trsInitSeed");
451  if (!seed) seed = 19460510;
452  StTrsRandom::inst().SetSeed(seed);
453  return StMaker::Init();
454 }
455 
456 Int_t StTrsMaker::InitRun(int runnumber)
457 {
458  if (mAllTheData) {gMessMgr->QAInfo() << "StTrsMaker::InitRun Already called" << endm; return kStOK;}
459  // The global pointer to the Db is gStTpcDb and it should be created in the macro.
460 
461  if (!gStTpcDb) {
462  gMessMgr->QAInfo() << "DATABASE MISSING!" << endm;
463  PR(gStTpcDb);
464  gMessMgr->QAInfo() << "Can't initialize TRS" << endm;
465  return kStFatal;
466  }
467  mGeometryDb =
468  StTpcDbGeometry::instance(gStTpcDb);
469  if (Debug()) mGeometryDb->print();
470 
471  // The print statements are done in Make() because the SlowControl DB is only available then.
472 
473  mElectronicsDb =
474  StTpcDbElectronics::instance(gStTpcDb);
475  if (Debug()) mElectronicsDb->print();
476 
477  mSlowControlDb =
478  StTpcDbSlowControl::instance(gStTpcDb);
479  if (Debug()) mSlowControlDb->print();
480 
481 // mMagneticFieldDb =
482 // StTpcDbMagneticField::instance(); // default is .5T field in z direction
483 
484 
485  float x[3] = {0,0,0};
486  float B[3];
487  StarMagField::Instance()->BField(x,B);
488  StThreeVector<double> Bfield(B[0],B[1],B[2]);
489  Bfield*=kilogauss;
490  PR(Bfield/tesla);
491  mMagneticFieldDb =
492  StSimpleMagneticField::instance(Bfield); // default is .5T field in z direction
493 
494 
495  //
496  // Select the gas: Ar, NeCO2, P10 available
497  string gas = "P10";
498  mGasDb = new StTrsDeDx(gas);
499  if (Debug()) mGasDb->print();
500 
501 
502  //
503  // Containers
504  //
505 
506  // A Wire Plane
507  mWireHistogram =
508  StTrsWireHistogram::instance(mGeometryDb, mSlowControlDb ,mGasDb,mMagneticFieldDb);
509  mWireHistogram->setDoGasGain(true);
510  mWireHistogram->setDoGasGainFluctuations(true); // used to be true
511  // mWireHistogram->setDoGasGainFluctuations(false);
512  mWireHistogram->setDoSingleElectronMultiplication(false);
513  mWireHistogram->setGasGainInnerSector(mSlowControlDb->innerSectorGasGain());
514  mWireHistogram->setGasGainOuterSector(mSlowControlDb->outerSectorGasGain());
515  mWireHistogram->setDoTimeDelay(false);
516  // mWireHistogram->setRangeOfWiresForChargeDistribution(0);
517  //
518  // An Analog Sector(for calculation)
519  mSector =
520  new StTrsSector(mGeometryDb);
521 
522  //
523  // Processes
524  //
525  mChargeTransporter =
526  StTrsFastChargeTransporter::instance(mGeometryDb, mSlowControlDb, mGasDb, mMagneticFieldDb);
527  // set status:
528  mChargeTransporter->setChargeAttachment(false); // used to be true
529  mChargeTransporter->setGatingGridTransparency(false);
530  mChargeTransporter->setTransverseDiffusion(true); // used to be true
531  mChargeTransporter->setLongitudinalDiffusion(true); // used to be true
532  mChargeTransporter->setExB(false);
533 
534  mAnalogSignalGenerator = 0;
535  if(mUseParameterizedSignalGenerator) {
536  //**************for the ParameterizedAnalogSignalGenerator
537  mAnalogSignalGenerator = StTrsParameterizedAnalogSignalGenerator::instance(mGeometryDb, mSlowControlDb, mElectronicsDb, mSector);
538  StTrsParameterizedAnalogSignalGenerator *myGen = dynamic_cast<StTrsParameterizedAnalogSignalGenerator*>(mAnalogSignalGenerator);
539 
540  // Any options that are needed can go here:
541  //*******************************************************
542  myGen->setDeltaPad(2);
543  myGen->setSignalThreshold(.1*millivolt);
544  myGen->setSuppressEmptyTimeBins(true);
545  myGen->addNoise(true);
546  myGen->generateNoiseUnderSignalOnly(true);
547  myGen->setNormalFactor(mNormalFactor);
548  }
549  else {
551  (StTrsSlowAnalogSignalGenerator::instance(mGeometryDb, mSlowControlDb, mElectronicsDb, mSector));
552  mAnalogSignalGenerator = myGen;
553  //
554  // Set the function for the induced charge on Pad
555  // -->StTrsSlowAnalogSignalGenerator::endo
556  //-->StTrsSlowAnalogSignalGenerator::gatti
557  //-->StTrsSlowAnalogSignalGeneratommetricGaussianApproximation:
558 // retr::dipole
559  myGen->setChargeDistribution(StTrsSlowAnalogSignalGenerator::endo);
560  //
561  // Set the function for the Analog Electronics signal shape
562  //-->StTrsSlowAnalogSignalGenerator::delta
563  //-->StTrsSlowAnalogSignalGenerator::symmetricGaussianApproximation
564  //-->StTrsSlowAnalogSignalGenerator::symmetricGaussianExact
565  //-->asymmetricGaussianApproximation
566  //-->StTrsSlowAnalogSignalGenerator::realShaper
567  myGen->setElectronicSampler(StTrsSlowAnalogSignalGenerator::symmetricGaussianApproximation);
568  //myGen->setDeltaRow(0);
569  myGen->setDeltaPad(2);
570  myGen->setSignalThreshold(.0*millivolt);
571  myGen->setSuppressEmptyTimeBins(true);
572  myGen->addNoise(true);
573  myGen->generateNoiseUnderSignalOnly(true);
574  myGen->setNoiseRMS(900); // set in #e
575  }
576 
577  mDigitalSignalGenerator =
578  StTrsFastDigitalSignalGenerator::instance(mElectronicsDb, mSector);
579 
580  //
581  // Output is into an StTpcRawDataEvent* vector<StTrsDigitalSector*>
582  // which is accessible via the StTrsUnpacker. Initialize the pointers!
583  //DH now in Init mAllTheData=0;
584 
585  //
586  // Maker Initialization ---now given by default arguments in the constructor
587  // mFirstSectorToProcess = 1;
588  // mLastSectorToProcess = 1;
589 
590 
591  // This should really be a boolean...when ROOT can handle it, change it!
592  mProcessPseudoPadRows = 1; // 0 is no, 1 is yes!
593 
594  //
595  // Construct constant data sets. This is what is passed downstream
596  if (!mAllTheData) {
597  mAllTheData = new StTrsRawDataEvent(mGeometryDb->numberOfSectors());
598  AddConst(new TObjectSet("Event" , mAllTheData, 0));
599  }
600  if ((Debug()/10)%10 > 0) {
601  mTrsNtupleFile = new TFile("TrsOutput.root","RECREATE","Trs Ntuples");
602  mWireNtuple = new TNtuple("WireNtuple", "Wire Histogram Info.", "electrons:wire:sector:id");
603  mContinuousAnalogNtuple = new TNtuple("CAnalogNtuple", "Cont. Analog Sector", "charge:time:pad:row:id");
604  mDiscreteAnalogNtuple = new TNtuple("DAnalogNtuple", "Disc. Analog Sector", "charge:timebin:pad:row:id");
605  mDigitalNtuple = new TNtuple("DigitalSignalNtuple", "Digital Sector", "adc:timebin:pad:row:id");
606  }
607  return kStOK;
608 }
609 
610 //
611 // My Member Functions
612 //
613 void StTrsMaker::whichSector(int volId, int* isDet, int* sector, int* padrow){
614 
615  //gMessMgr->QAInfo() << "StTrsMaker::whichSector()" << endm;
616  *isDet = (volId/100000)%10;
617 
618  volId -= (*isDet)*100000;
619  *sector = volId/100;
620 
621  volId -= (*sector)*100;
622  *padrow = volId;
623 
624 }
626 
627  // PrintInfo();
628  double d[3],absP[3];
629 
630  //Do not use this unless you really know what you are
631  // doing...This is a histogram diagnostic to compare
632  // the GEANT hits to those produced by TRS!
633  time_t trsMakeBegin = time(0);
634  gMessMgr->QAInfo() << "\n -- Begin TRS Processing -- " << endm;
635  gMessMgr->QAInfo() << "Started at: " << ctime(&trsMakeBegin);
636  gMessMgr->QAInfo() << "========= TRS driftVelocity used = "
637  << mSlowControlDb->driftVelocity(13) << "(East) "
638  << mSlowControlDb->driftVelocity(1) << "(West)" << endm;
639  int seed = IAttr("trsMakeSeed");
640  if (seed) StTrsRandom::inst().SetSeed(seed);
641  gMessMgr->QAInfo() << "========= TRS Seed used = " << StTrsRandom::inst().GetSeed()<< endm;
642 
643 
644  int currentSectorProcessed = mFirstSectorToProcess;
645  if (Debug()) {
646  gMessMgr->QAInfo() << "Processing sectors "
647  << mFirstSectorToProcess
648  << "--"
649  << mLastSectorToProcess << endm;
650 
651  //
652  gMessMgr->QAInfo() << "make sure pointer are clean" << endm;
653  }
654  mAllTheData->clear();
655  //
656  //
657  //update the gain for very events
658  mWireHistogram->setGasGainInnerSector(mSlowControlDb->innerSectorGasGain());
659  mWireHistogram->setGasGainOuterSector(mSlowControlDb->outerSectorGasGain());
660 
661  // Normal processing of TRS through GEANT
662  //gMessMgr->QAInfo() << "Make ofstream" << endm;
663  //ofstream ofs("/star/u2b/lasiuk/geantdebug.txt", ios::out);
664  //ofstream raw("/star/u2b/lasiuk/event.txt",ios::out);
665  //
666  // Read the GEANT info
667  // these structures/classes are defined in:
668  // $STAR/pams/sim/idl/g2t_tpc_hit.idl
669  // $STAR/StRoot/base/TDataSet.h & St_Table.h
670  //
671  TDataSetIter geant(GetDataSet("geant"));
672 
673  St_g2t_event *g2tevent = (St_g2t_event *)(geant("g2t_event"));
674  if ( !g2tevent ) return kStWarn;
675  g2t_event_st *event = g2tevent->GetTable();
676  seed = event->ge_rndm[0]^event->ge_rndm[1];
677  if (seed) StTrsRandom::inst().SetSeed(seed);
678  gMessMgr->QAInfo() << "========= TRS Seed used = " << StTrsRandom::inst().GetSeed()<< endm;
679 
680 
681  St_g2t_tpc_hit *g2t_tpc_hit =
682  static_cast<St_g2t_tpc_hit *>(geant("g2t_tpc_hit"));
683  //PR(g2t_tpc_hit);
684 
685  if (!g2t_tpc_hit) return kStWarn;
686 
687  int no_tpc_hits = g2t_tpc_hit->GetNRows();
688  g2t_tpc_hit_st *tpcHit = g2t_tpc_hit->GetTable();
689 
690  St_g2t_track *g2t_track =
691  static_cast<St_g2t_track *>(geant("g2t_track"));
692  if (!g2t_track) return kStWarn;
693 
694  g2t_track_st *tpc_track = g2t_track->GetTable();
695 
696  St_g2t_vertex *g2t_ver=static_cast<St_g2t_vertex *>(geant("g2t_vertex"));
697  if (!g2t_ver) return kStWarn;
698 
699  g2t_vertex_st *gver=g2t_ver->GetTable();
700  assert(gver);
701  if(PILEUP_ON) {
702  gMessMgr->QAInfo() << Form("\n TRS(): Pileup is ON (m_Mode=%d)\n",m_Mode) << endm;
703  }
704  else {
705  gMessMgr->QAInfo() << Form("\n TRS(): Pileup is OFF (m_Mode=%d)\n",m_Mode) << endm;
706  }
707  //int geantPID = tpc_track->ge_pid;
708  //PR(geantPID);
709 
710  // inRange should be a boolean
711  bool inRange,start=true;
712  int bisdet, bsectorOfHit, bpadrow;
713  int numberOfProcessedPointsInCurrentSector = 0;
714  // Limit the processing to a fixed number of segments
715  // no_tpc_hits = 4;
716  if(no_tpc_hits<1)return kStOK;
717  g2t_tpc_hit_st *tpc_hit = tpcHit;
718  for (int i=1; i<=no_tpc_hits; i++,tpc_hit++){
719  int id2=tpc_hit->track_p;
720  int id3=tpc_track[id2-1].start_vertex_p; // "-1" is (Fortran-->C++)
721  whichSector(tpc_hit->volume_id, &bisdet, &bsectorOfHit, &bpadrow);
722  float BunchZoffset=(gver[id3-1].ge_tof+tpc_hit->tof)* mSlowControlDb->driftVelocity(bsectorOfHit);
723  mChargeTransporter->setDriftVelocity(mSlowControlDb->driftVelocity(bsectorOfHit));
724 #if 0
725  float absHitZ=fabs(tpc_hit->x[2]);
726  if(PILEUP_ON)
727  {
728  if(absHitZ - tpc_hit->ds + BunchZoffset<0) continue;//crossed central membrane
729  if(absHitZ + tpc_hit->ds + BunchZoffset> mGeometryDb->frischGrid())
730  continue;//out of TPC
731  } // for piled up events
732 #endif
733 // gMessMgr->QAInfo() << "--> tpc_hit: " << i << endm;
734 // raw << tpc_hit->volume_id << ' '
735 // << tpc_hit->de << ' '
736 // << tpc_hit->ds << ' '
737 // << tpc_hit->x[0] << ' '
738 // << tpc_hit->x[1] << ' '
739 // << tpc_hit->x[2] << ' '
740 // << tpc_hit->p[0] << ' '
741 // << tpc_hit->p[1] << ' '
742 // << tpc_hit->p[2] << ' ' << endm;
743 
744 // PR(bsectorOfHit);
745 // gMessMgr->QAInfo() << mFirstSectorToProcess<<" "<< mLastSectorToProcess<<endm;
746  if(bsectorOfHit >= mFirstSectorToProcess &&
747  bsectorOfHit <= mLastSectorToProcess)
748  inRange = true;
749  else
750  inRange = false;
751 
752  // Save time initially - by not processing pseudo padrows
753  if(bisdet && !mProcessPseudoPadRows) {
754  // gMessMgr->QAInfo() << "Segment in a pseudo-padRow. Skipping..." << endm;
755  // tpc_hit++;
756  if(i != no_tpc_hits) continue;
757  inRange=false;
758  }
759  //sleep(2);
760 
761  //
762  // If not in range AND there are no points processed, skip to next point
763  if(!inRange && !numberOfProcessedPointsInCurrentSector) {
764  //gMessMgr->QAInfo() << "out of range and no points" << endm;
765  // tpc_hit++;
766  continue;
767  }
768  if((inRange) && //HL
769  (bsectorOfHit != currentSectorProcessed&&start==false) &&//HL
770  (i <= no_tpc_hits )) {//HL
771 
772 
773  tpc_hit--;
774  i--;
775  } //HL ,continue to another sector....
776  else if((inRange) &&
777  (bsectorOfHit == currentSectorProcessed||start==true) &&
778  (i <= no_tpc_hits )) {
779  currentSectorProcessed=bsectorOfHit;//HL
780  start=false;//HL,9/9/99
781  // CAREFUL:
782  // GEANT uses: (which is not correct!)
783  //double GEANTDriftLength = 208.55119*centimeter;
784  //double GEANTOffSet = mGeometryDb->frischGrid() - GEANTDriftLength;
785 
786  // Now relative to this, we get the zPosition in coordinates where :
787  // 0 is the membrane, 208+/-dz is the wire grid
788  //double zPosition =
789  //GEANTDriftLength/2. + tpc_hit->x[2]*centimeter + GEANTOffSet;
790  //--->double zPosition = tpc_hit->x[2]*centimeter;
791  //PR(tpc_hit->x[2]*centimeter);
792  //PR(zPosition);
793 
794  // StThreeVector<double> hitPosition(tpc_hit->x[0]*centimeter,
795  // tpc_hit->x[1]*centimeter,
796  // tpc_hit->x[2]*centimeter);
797 // PR(hitPosition);
798 
799 // // Drift Length is calculated with respect to the FG!
800 // double fgOffSet = (bpadrow <= mGeometryDb->numberOfInnerRows()) ?
801 // mGeometryDb->innerSectorFrischGridPadPlaneSeparation() :
802 // mGeometryDb->innerSectorFrischGridPadPlaneSeparation();
803 
804  // In GEANT Global Coordinates we have to rotate
805  // to the sector 12 position
806  // Not using StTpcCoordinateTransform
807  // because it is slower, so this requires that we make sure
808  // the code below follows exactly what is in StTpcCoordinateTransform
809  // It is also in StTrsChargeSegment::rotate()
810  double beta = (bsectorOfHit>12) ?
811  -bsectorOfHit*M_PI/6. :
812  bsectorOfHit*M_PI/6. ; //(30 degrees)
813  double cb = cos(beta);
814  double sb = sin(beta);
815  // double xp = hitPosition.x()*cb - hitPosition.y()*sb;
816  // double yp = hitPosition.x()*sb + hitPosition.y()*cb;
817  double xp = tpc_hit->x[0]*cb -tpc_hit->x[1]*sb;
818  double yp = tpc_hit->x[0]*sb +tpc_hit->x[1]*cb;
819 
821  sector12Coordinate(xp,yp,(tpc_hit->x[2]));
822 
823 
824  // Must rotate the momentum as well, BUT you should
825  // only incur this calculational penalty if you split
826  // the segment into more than 1 mini segement
827  double pxPrime = tpc_hit->p[0]*cb - tpc_hit->p[1]*sb;
828  double pyPrime = tpc_hit->p[0]*sb + tpc_hit->p[1]*cb;
829  absP[0]=fabs(pxPrime);
830  absP[1]=fabs(pyPrime);
831  absP[2]=fabs(tpc_hit->p[2]);
832  StThreeVector<double> hitMomentum(pxPrime*GeV,
833  pyPrime*GeV,
834  tpc_hit->p[2]*GeV);
835 
836  //added by Hui Long ,8/24/99
837  //Modified by M. Calderon 2/5/00
838 
839  if(bsectorOfHit>12)
840  {
841  sector12Coordinate.setZ((tpc_hit->x[2])+mGeometryDb->driftDistance());
842 
843 
844  }
845  else
846  {
847  sector12Coordinate.setX(-xp);
848  sector12Coordinate.setZ(-(tpc_hit->x[2])+mGeometryDb->driftDistance());
849  hitMomentum.setX(-pxPrime*GeV);
850  hitMomentum.setZ(-(tpc_hit->p[2]*GeV));
851  }
852  //
853 
854 // PR(tpc_hit->p[0]*GeV);
855 // PR(pxPrime*GeV);
856 // PR(hitMomentum);
857 
858 
859  // I need PID info here, for the ionization splitting (beta gamma)!
860  // int geantPID = tpc_track[tpc_hit->track_p].ge_pid;
861  int geantPID = tpc_track[tpc_hit->track_p-1].ge_pid;
862 
863  // WARNING: cannot use "abs" (not overloaded (double) for LINUX!
864  StTrsChargeSegment aSegment(sector12Coordinate,
865  hitMomentum,id2,
866  (fabs(tpc_hit->de*GeV)),
867  tpc_hit->ds*centimeter,
868  geantPID);
869 // PR(hitPosition);
870 // PR(sector12Coordinate);
871 // PR(hitMomentum.mag());
872 
873 // ofs << " " << aSegment << endl;
874 // PR(aSegment);
875 
876 
877 #ifndef ST_NO_TEMPLATE_DEF_ARGS
878  list<StTrsMiniChargeSegment> comp;
879  list<StTrsMiniChargeSegment>::iterator iter;
880 #else
881  list<StTrsMiniChargeSegment,allocator<StTrsMiniChargeSegment> > comp;
882  list<StTrsMiniChargeSegment,allocator<StTrsMiniChargeSegment> >::iterator iter;
883 #endif
884  //
885  // Break the segment for diffusion reproduction.
886  // Fast Simulation can use breakNumber = 1
887  double ptot=::sqrt(absP[0]*absP[0]+absP[1]*absP[1]+absP[2]*absP[2]);
888 
889  d[0] =tpc_hit->ds*absP[0]/ptot;
890  d[1] =tpc_hit->ds*absP[1]/ptot;
891  d[2] =tpc_hit->ds*absP[2]/ptot;// approximation
892  // int breakNumber = (int)max(d[0]/mMiniSegmentLength,d[2]/mMiniSegmentLength);
893  int breakNumber = (int)max(aSegment.ds()/mMiniSegmentLength,1.);
894  if(breakNumber<1) breakNumber = 1;
895  breakNumber = min(breakNumber,16); // set limit
896  int numberOfLevels =
897  static_cast<int>(::log(static_cast<double>(breakNumber))/M_LN2 + .999);
898  d[0]/=::pow(2.,numberOfLevels);
899  d[1]/=::pow(2.,numberOfLevels);
900  d[2]/=::pow(2.,numberOfLevels);// approximation
901  numberOfLevels++; // take care of the zero! in aSegment.tssSplit
902 // PR(aSegment.ds()/millimeter);
903 // PR(breakNumber);
904 
905  //
906  // int breakNumber = (int)max(aSegment.ds()/mMiniSegmentLength,1.);
907  // breakNumber = min(breakNumber,16); // set limit
908 // PR(aSegment.ds()/millimeter);
909 // PR(breakNumber);
910 
911 // aSegment.split(mGasDb, mMagneticFieldDb, breakNumber, &comp);
912 
913  // aSegment.tssSplit(mGasDb, mMagneticFieldDb, breakNumber, &comp);
914  aSegment.tssSplit(mGasDb, mMagneticFieldDb, numberOfLevels , &comp);
915 
916 #ifndef ST_NO_TEMPLATE_DEF_ARGS
917 // copy(comp.begin(), comp.end(), ostream_iterator<StTrsMiniChargeSegment>(cout,"\n"));
918 #endif
919 
920  double SigmaL,SigmaT;//HL,added for field on case,02/08/00
921 
922  // Loop over the miniSegments
923  for(iter = comp.begin();
924  iter != comp.end();
925  iter++) {
926 
927  //
928  // TRANSPORT HERE
929  //
930  // mChargeTransporter->transportToWire(*iter);
931  mChargeTransporter->transportToWire(*iter,SigmaL,SigmaT);//HL,08/02/00 for field on
932 
933  if(PILEUP_ON)
934  {
935  StTrsMiniChargeSegment *seg1=&(*iter);
936  float Znew=seg1->position().z()-BunchZoffset;
937  seg1->position().setZ(Znew);
938  }
939 
940 
941 // PR(*iter);
942 
943  //
944  // CHARGE COLLECTION AND AMPLIFICATION
945  //
946 
947  //#if defined(__sun) && !defined(__GNUG__)
948 // Bug in the sun iterators. Must Explicitly dereference!
949  // StTrsWireBinEntry anEntry(iter->position(), iter->charge());
950 // PR(anEntry);
951  //#else
952 
953  // StTrsWireBinEntry anEntry((*iter).position(), (*iter).charge());
954 
955  //#endif
956  // mWireHistogram->addEntry(anEntry);
957 
958  StTrsWireBinEntry anEntry((*iter).position(), (*iter).charge(),SigmaL,SigmaT,d,(*iter).id());//HL,for field on ,08/02/00
959  mWireHistogram->addEntry(anEntry,bsectorOfHit);
960  } // Loop over the list of iterators
961 
962  // tpc_hit++; // increase the pointer to the next hit
963  numberOfProcessedPointsInCurrentSector++;
964 
965  if(i<no_tpc_hits)continue; // don't digitize, you still have data in the same sector to process
966  } // if (currentSector == bsectorOfHit)
967  // Otherwise, do the digitization...
968 
969  if (mWireNtuple) {
970  if (mWireHistogram->minWire() >= 0) {
971  float wireValues[4];
972  // Loop over Wire Histogram
973  for(int jj=mWireHistogram->minWire(); jj<=mWireHistogram->maxWire(); jj++) {
974  aTpcWire currentWire = mWireHistogram->getWire(jj);
975  aTpcWire::iterator iter;
976  for(iter = currentWire.begin();
977  iter != currentWire.end();
978  iter++) {
979  wireValues[0] = iter->numberOfElectrons();
980  wireValues[1] = jj;
981  wireValues[2] = currentSectorProcessed;
982  wireValues[3] = iter->id();
983  mWireNtuple->Fill(wireValues);
984  }
985  }
986  mWireNtuple->Write();
987  }
988  }
989  if (Debug()) PR(currentSectorProcessed);
990 
991  //
992  // Generate the ANALOG Signals on pads
993  //
994  time_t inducedChargeBegin = time(0);
995  if (Debug()) {gMessMgr->QAInfo() << "--->inducedChargeOnPad()..." << endm;}
996  mAnalogSignalGenerator->inducedChargeOnPad(mWireHistogram,currentSectorProcessed);
997  if (Debug()) {
998  time_t inducedChargeEnd= time(0);
999  double inducedChargeTime = difftime(inducedChargeEnd,inducedChargeBegin);
1000  gMessMgr->QAInfo() << "Time to process induced Charge: " << inducedChargeTime << " sec\n" << endm;
1001  }
1002  if (mContinuousAnalogNtuple) {
1003  tpcTimeBins continuousAnalogTimeSequence;
1004  timeBinIterator timeSeqIter;
1005  float cAnalogValues[5];
1006  // Loop over Continuous Analog Sector
1007  for(int jrow=1; jrow<=mGeometryDb->numberOfRows(); jrow++) {
1008  for (int jpad=1; jpad<=mGeometryDb->numberOfPadsAtRow(jrow); jpad++){
1009  continuousAnalogTimeSequence = mSector->timeBinsOfRowAndPad(jrow,jpad);
1010  if(!continuousAnalogTimeSequence.size()) continue;
1011  for(timeSeqIter = continuousAnalogTimeSequence.begin();
1012  timeSeqIter != continuousAnalogTimeSequence.end();
1013  timeSeqIter++) {
1014  cAnalogValues[0] = timeSeqIter->amplitude();
1015  cAnalogValues[1] = timeSeqIter->time();
1016  cAnalogValues[2] = jpad;
1017  cAnalogValues[3] = jrow;
1018  cAnalogValues[4] = timeSeqIter->id();
1019  mContinuousAnalogNtuple->Fill(cAnalogValues);
1020  }
1021  }
1022  }
1023  mContinuousAnalogNtuple->Write();
1024  }
1025 
1026  time_t sampleAnalogSignalBegin = time(0);
1027  if (Debug()) {
1028  gMessMgr->QAInfo() << "--->sampleAnalogSignal()..." << endm;
1029  // mAnalogSignalGenerator->sampleAnalogSignal();
1030  time_t sampleAnalogSignalEnd= time(0);
1031  double sampleAnalogSignalTime = difftime(sampleAnalogSignalEnd,sampleAnalogSignalBegin);
1032  gMessMgr->QAInfo() << "Time to sample Analog Signal: " << sampleAnalogSignalTime << " sec\n" << endm;
1033  }
1034  if (mDiscreteAnalogNtuple) {
1035  tpcTimeBins discreteAnalogTimeSequence;
1036  timeBinIterator timeBinIter;
1037  float dAnalogValues[5];
1038  // Loop over Discrete Analog Sector
1039  for(int drow=1; drow<=mGeometryDb->numberOfRows(); drow++) {
1040  for (int dpad=1; dpad<=mGeometryDb->numberOfPadsAtRow(drow); dpad++){
1041  discreteAnalogTimeSequence = mSector->timeBinsOfRowAndPad(drow,dpad);
1042  if(!discreteAnalogTimeSequence.size()) continue;
1043  for(timeBinIter = discreteAnalogTimeSequence.begin();
1044  timeBinIter != discreteAnalogTimeSequence.end();
1045  timeBinIter++) {
1046  dAnalogValues[0] = timeBinIter->amplitude();
1047  dAnalogValues[1] = timeBinIter->time();
1048  dAnalogValues[2] = dpad;
1049  dAnalogValues[3] = drow;
1050  dAnalogValues[4] = timeBinIter->id();
1051  mDiscreteAnalogNtuple->Fill(dAnalogValues);
1052  }
1053  }
1054  }
1055  mDiscreteAnalogNtuple->Write();
1056  }
1057 
1058  //
1059  // Digitize the Signals
1060  //
1061  // First make a sector where the data can go...
1062  StTrsDigitalSector* aDigitalSector =
1063  new StTrsDigitalSector(20);
1064  aDigitalSector->setSector(currentSectorProcessed);
1065  //
1066  // Point to the object you want to fill
1067  //
1068  mDigitalSignalGenerator->fillSector(aDigitalSector);
1069  mDigitalSignalGenerator->SetSectorNo(currentSectorProcessed);
1070  //
1071  // ...and digitize it
1072  time_t digitizeSignalBegin = time(0);
1073  if (Debug()) {gMessMgr->QAInfo() << "--->digitizeSignal()..." << endm;}
1074  mDigitalSignalGenerator->digitizeSignal();
1075  if (Debug()) {
1076  gMessMgr->QAInfo() <<"--->digitizeSignal() Finished..." << endm;
1077  time_t digitizeSignalEnd= time(0);
1078  double digitizeSignalTime = difftime(digitizeSignalEnd,digitizeSignalBegin);
1079  gMessMgr->QAInfo() << "Time to digitize Signal: " << digitizeSignalTime << " sec\n" << endm;
1080  }
1081  //
1082  // Fill it into the event structure...
1083  // and you better check the sector number!
1084 
1085  mAllTheData->setSector(currentSectorProcessed,aDigitalSector);
1086  // Clear and reset for next sector:
1087  mWireHistogram->clear();
1088  mSector->clear();
1089  if (Debug()) gMessMgr->QAInfo() << endm;
1090 
1091  //
1092  // Go to the next sector --> should be identical to a simple increment
1093  currentSectorProcessed = bsectorOfHit;
1094  numberOfProcessedPointsInCurrentSector = 0;
1095 // you can skip out here if you only want to process a single sector...
1096 // if(currentSectorProcessed>3)
1097 // break; // Finish here
1098 
1099  //
1100  } // loop over all segments: for(int i...
1101 
1102  // The access stuff:
1103  if (Debug() > 2) {
1104  //
1105  // Access the data with
1106  // *mDetectorReader
1107 
1108  string version = "TrsDRv1.0";
1109  StTrsDetectorReader mDetectorReader(mAllTheData, version);
1110  //
1111  // Loop around the sectors: (should be from db, or size of the structure!)
1112  //
1113  // for(int isector=1; isector<=mGeometryDb->numberOfSectors(); isector++) {
1114  for(int isector=mFirstSectorToProcess; isector<=mLastSectorToProcess; isector++) {
1115  StTrsZeroSuppressedReader* zsr = mDetectorReader.getZeroSuppressedReader(isector);
1116  if(!zsr) continue;
1117 // PR(isector);
1118  // otherwise, let's decode it
1119  unsigned char* padList;
1120  if (mDigitalNtuple) {
1121  float digitalValues[5];
1122  for (int irow=1; irow<=mGeometryDb->numberOfRows();irow++) {
1123  int numberOfPads = zsr->getPadList(irow, &padList);
1124  // If there are no pads, go to the next row...
1125  if(!numberOfPads) continue;
1126  for(int ipad = 0; ipad<numberOfPads; ipad++) {
1127  //PR(static_cast<int>(padList[ipad]));
1128  int nseq;
1129 
1130  StSequence* listOfSequences;
1131  //Sequence* listOfSequences;
1132  UShort_t ** listOfIds = 0;
1133  zsr->getSequences(irow,
1134  static_cast<int>(padList[ipad]),
1135  &nseq,
1136  &listOfSequences, &listOfIds);
1137  //PR(getSequencesStatus);
1138  for(int kk=0; kk<nseq; kk++) {
1139  for(int zz=0; zz<listOfSequences[kk].length; zz++) {
1140  if (Debug()%10 > 2 && kk < Debug()) {
1141  gMessMgr->QAInfo() << "sector\t" << isector << "\trow\t" << irow
1142  << "\tpad\t" << ipad << "\tseq\t" << kk
1143  << "\tz\t" << zz + (int)(listOfSequences[kk].startTimeBin)
1144  << "\tADC\t" << (int)(*(listOfSequences[kk].firstAdc))
1145  << "\tid\t" << (int)(*(listOfIds[kk]))
1146  << endm;
1147  }
1148  digitalValues[0] = (int)(*(listOfSequences[kk].firstAdc)) ;
1149  digitalValues[1] = zz + (int)(listOfSequences[kk].startTimeBin);
1150  digitalValues[2] = ipad;
1151  digitalValues[3] = irow;
1152  digitalValues[4] = (int)(*(listOfIds[kk]));
1153  mDigitalNtuple->Fill(digitalValues);
1154 
1155  listOfSequences[kk].firstAdc++;
1156  listOfIds[kk]++;
1157  } // zz
1158  } // Loop kk
1159  } // loop over pads
1160  //
1161  // One would do the data manipulation here!
1162  // Then deallocate the memory
1163  // dynamic_cast<StTrsZeroSuppressedReader*>(zsr);
1164  // zsr->clear();
1165  } // Loop over rows!
1166  mDigitalNtuple->Write();
1167  }
1168  } // Loop over sectors
1169  }
1170 
1171  //gMessMgr->QAInfo() << "Got to the end of the maker" << endm;
1172  // CAUTION: ROOT is resposible for the memory at this point
1173  // ROOT deletes m_DataSet in the chain after every event.
1174  if (Debug()) {
1175  time_t trsMakeEnd = time(0);
1176  gMessMgr->QAInfo() << "\nFinished at: " << ctime(&trsMakeEnd);
1177  double trsMakeTotal = difftime(trsMakeEnd,trsMakeBegin);
1178  gMessMgr->QAInfo() << "Total StTrsMaker::Make() Processing Time: " << trsMakeTotal << " sec" << endm;
1179  }
1180 #if 0
1181  CheckTruth(no_tpc_hits, tpcHit);
1182 #endif
1183  return kStOK;
1184 }
1185 
1186 // *****************************************************************
1187 // Make sure the memory is deallocated!
1188 //
1189 void StTrsMaker::Clear(const char *)
1190 {
1191  if (mAllTheData) mAllTheData ->clear(); //This deletes all the StTpcDigitalSectors in the StTrsRawDataEvent
1192  if (mWireHistogram)mWireHistogram->clear();
1193 
1194  StMaker::Clear();
1195 }
1197 {
1198  //Clean up all the pointers that were initialized in StTrsMaker::Init()
1199  // Don't delete pointers to databases, they're singletons.
1200 // if (mGeometryDb) delete mGeometryDb;
1201 // if (mSlowControlDb) delete mSlowControlDb;
1202 // if (mMagneticFieldDb) delete mMagneticFieldDb;
1203 // if (mElectronicsDb) delete mElectronicsDb;
1204 // if (mGasDb) delete mGasDb;
1205 
1206  StTrsWireHistogram::dropit();
1207  mWireHistogram = 0;
1208  if (mSector) delete mSector;
1209  mSector = 0;
1210  if (mAllTheData) delete mAllTheData;
1211  mAllTheData = 0;
1212  if (mChargeTransporter) delete mChargeTransporter;
1213  mChargeTransporter = 0;
1214  if (mAnalogSignalGenerator) delete mAnalogSignalGenerator;
1215  mAnalogSignalGenerator = 0;
1216  if (mDigitalSignalGenerator)delete mDigitalSignalGenerator;
1217  mDigitalSignalGenerator = 0;
1218  return kStOK;
1219 }
1220 
1221 void StTrsMaker::setNormalFactor(double FudgeFactor) {
1222  mNormalFactor = FudgeFactor;
1223 //VP if (mUseParameterizedSignalGenerator&&mAnalogSignalGenerator) {
1224 //VP mAnalogSignalGenerator->setNormalFactor(mNormalFactor);
1225 //VP }
1226 }
1227 #include "StDbUtilities/StTpcCoordinateTransform.hh"
1228 #include "StarClassLibrary/StSequence.hh"
1229 
1230 void StTrsMaker::CheckTruth(int no_tpc_hits, g2t_tpc_hit_st *tpc_hit)
1231 {
1232  StSequence *seq; int nSeq;
1233  UShort_t **ids;
1234  int sector, row,pad,tim;
1235  TObjectSet *trsEvent=(TObjectSet*)GetDataSet("Event"); if(!trsEvent) return;
1236  StTrsZeroSuppressedReader *mZsr=0;
1237  StTpcRawDataEvent *mEvent=(StTpcRawDataEvent*)(trsEvent->GetObject()); if(!mEvent) return;
1238  StTrsDetectorReader mTdr(mEvent);
1239  StTpcCoordinateTransform transform(gStTpcDb);
1240 
1241  int noData=0,lowTruth=0;
1242  for (int ih=0;ih<no_tpc_hits;ih++) {
1243  UShort_t idtru = (UShort_t) tpc_hit[ih].track_p;
1244  StGlobalCoordinate global(tpc_hit[ih].x[0],tpc_hit[ih].x[1],tpc_hit[ih].x[2]);
1245  sector = (tpc_hit[ih].volume_id/100)%100;
1246  row = tpc_hit[ih].volume_id % 100;
1247  StTpcPadCoordinate Pad; transform(global,Pad,sector,row);
1248  sector = Pad.sector();
1249  pad = (Int_t) Pad.pad();
1250  row = Pad.row();
1251  tim = (Int_t) Pad.timeBucket();
1252  mZsr=mTdr.getZeroSuppressedReader(sector);
1253  if (!mZsr) {
1254  noData++;
1255 // Warning("CheckTruth","%d - No Data in %d.%d.%d",noData,sector,row,pad);
1256  continue;
1257  }
1258  nSeq=0;
1259  mZsr->getSequences(row, pad, &nSeq, &seq, &ids);
1260  if (!nSeq) {
1261  noData++;
1262 // Warning("CheckTruth","%d - No Data in %d.%d.%d",noData,sector,row,pad);
1263  continue;
1264  }
1265 
1266  int nAdc=0,nIds=0;
1267 
1268  for (int jSeq=0;jSeq<nSeq;jSeq++) {
1269  int nnAdc = seq[jSeq].length;
1270  int startTimeBin = seq[jSeq].startTimeBin;
1271  if (tim < startTimeBin ) continue;
1272  if (tim > startTimeBin+nnAdc) continue;
1273  for (int j=0;j<nnAdc;j++) {
1274  if(abs(startTimeBin+j-tim)>5) continue;
1275  nAdc++;
1276  if(ids[jSeq][j]==idtru) nIds++;
1277  }
1278  }
1279  if (!nAdc) nAdc=1;
1280  double pct = double(nIds)/nAdc*100; if (pct){};
1281  if (!nIds) {
1282  lowTruth++;
1283 // Warning("CheckTruth","%d - Low contrib %2.1f%% in %d.%d.%d",lowTruth,pct,sector,row,pad);
1284  }
1285  }
1286  if (lowTruth)
1287  Warning("CheckTruth","nHits=%d lowTruth=%d pct=%2.1f%%",
1288  no_tpc_hits,lowTruth,(lowTruth*100.)/no_tpc_hits);
1289 
1290 }
Int_t m_Mode
counters
Definition: StMaker.h:81
virtual void Clear(Option_t *option="")
User defined functions.
Definition: StMaker.cxx:634
Int_t Make()
Definition: StTrsMaker.cxx:625
Int_t Finish()
Definition: Stypes.h:42
Definition: Stypes.h:40
virtual TObject * GetObject() const
The depricated method (left here for the sake of the backward compatibility)
Definition: TObjectSet.h:56