00001
00002
00003
00004 #define HISTOGRAM 1
00005 #define UNPACK_ALL
00006 #define ALL_OF_IT
00007
00008
00010 #include <Stiostream.h>
00011 #include <unistd.h>
00012
00013 #include <string>
00014 #include <vector>
00015 #include <list>
00016 #include <utility>
00017 #include <algorithm>
00018
00019
00020 #include "StGlobals.hh"
00021 #include "StMatrix.hh"
00022 #include "Randomize.h"
00023 #ifdef HISTOGRAM
00024 #include "StHbook.hh"
00025 #endif
00026
00027 #include "StCoordinates.hh"
00028 #include "StTpcCoordinateTransform.hh"
00029
00030
00031
00032 #include "StTpcSimpleGeometry.hh"
00033 #include "StTpcSimpleSlowControl.hh"
00034 #include "StTpcSimpleElectronics.hh"
00035 #include "StSimpleMagneticField.hh"
00036 #include "StTrsDeDx.hh"
00037
00038
00039 #include "StTrsFastChargeTransporter.hh"
00040 #include "StTrsSlowAnalogSignalGenerator.hh"
00041 #include "StTrsFastDigitalSignalGenerator.hh"
00042
00043
00044 #include "StTrsChargeSegment.hh"
00045 #include "StTrsMiniChargeSegment.hh"
00046 #include "StTrsAnalogSignal.hh"
00047 #include "StTrsWireBinEntry.hh"
00048 #include "StTrsWireHistogram.hh"
00049
00050 #include "StTrsSector.hh"
00051 #include "StTrsDigitalSector.hh"
00052
00053
00054 #include "StTrsRawDataEvent.hh"
00055 #include "StTrsUnpacker.hh"
00056 #include "StSequence.hh"
00057
00058
00059 #include "g2t_tpc_hit.hh"
00060
00061
00062
00063 void whichSector(int volId, int* isDet, int* sector, int* padrow){
00064
00065
00066
00067 *isDet = (volId/100000);
00068
00069 volId -= (*isDet)*100000;
00070 *sector = volId/100;
00071
00072 volId -= (*sector)*100;
00073 *padrow = volId;
00074
00075 }
00076
00077 int main()
00078 {
00079 #ifdef HISTOGRAM
00080 cout << "Histogram" << endl;
00081
00082
00083
00084
00085 string fname = "hbook";
00086 StHbookFile *hbookFile =
00087 new StHbookFile(fname.c_str());
00088
00089
00090
00091 const int tupleSize1 = 5;
00092 float tuple[tupleSize1];
00093 StHbookTuple *theTuple =
00094 new StHbookTuple("data", tupleSize1);
00095 *theTuple << "sec" << "row" << "pad" << "amp" << "t" << book;
00096
00097
00098 const int tupleSize2 = 5;
00099 float tuple2[tupleSize2];
00100 StHbookTuple *secTuple =
00101 new StHbookTuple("dedx", tupleSize2);
00102
00103 *secTuple << "sec" << "row" << "x" << "y" << "z" << book;
00104 #endif
00105
00106
00107
00108
00109
00110 cout << "*********************************************" << endl;
00111 cout << "Try Initialize" << endl;
00112 string geoFile("../run/TPCgeo.conf");
00113 if (access(geoFile.c_str(),R_OK)) {
00114 cerr << "ERROR:\n" << geoFile.c_str() << " cannot be opened" << endl;
00115
00116 cerr << "Exitting..." << endl;
00117 exit(1);
00118 }
00119
00120 string scFile = "../run/sc.conf";
00121 if (access(scFile.c_str(),R_OK)) {
00122 cerr << "ERROR:\n" << scFile.c_str() << " cannot be opened" << endl;
00123 cerr << "Exitting..." << endl;
00124 exit(1);
00125 }
00126
00127 string electronicsFile = "../run/electronics.conf";
00128 if (access(electronicsFile.c_str(),R_OK)) {
00129 cerr << "ERROR:\n" << electronicsFile.c_str() << " cannot be opened" << endl;
00130 cerr << "Exitting..." << endl;
00131 exit(1);
00132 }
00133
00134 string magFile = "../run/example.conf";
00135 if (access(magFile.c_str(),R_OK)) {
00136 cerr << "ERROR:\n" << magFile.c_str() << " cannot be opened" << endl;
00137 cerr << "Exitting..." << endl;
00138 exit(1);
00139 }
00140
00141
00142
00143
00144 StTpcGeometry *mGeometryDb =
00145 StTpcSimpleGeometry::instance(geoFile.c_str());
00146
00147
00148 StTpcSlowControl *mSlowControlDb =
00149 StTpcSimpleSlowControl::instance(scFile.c_str());
00150
00151 StMagneticField *mMagneticFieldDb =
00152 StSimpleMagneticField::instance(magFile.c_str());
00153
00154 StTpcElectronics *mElectronicsDb =
00155 StTpcSimpleElectronics::instance(electronicsFile.c_str());
00156
00157
00158 string gas =("Ar");
00159 StTrsDeDx *mGasDb = new StTrsDeDx(gas);
00160
00161
00162
00163
00164
00165
00166
00167 StTrsWireHistogram *mWireHistogram =
00168 StTrsWireHistogram::instance(mGeometryDb, mSlowControlDb);
00169
00170
00171 mWireHistogram->setGasGainInnerSector(2000);
00172 mWireHistogram->setGasGainOuterSector(1000);
00173 mWireHistogram->setDoTimeDelay(false);
00174
00175
00176
00177 StTrsSector *mSector =
00178 new StTrsSector(mGeometryDb);
00179
00180
00181
00182 StTrsUnpacker *mUnPacker = new StTrsUnpacker;
00183
00184 StTrsRawDataEvent *mAllTheData =
00185 new StTrsRawDataEvent();
00186
00187
00188
00189 StTrsChargeTransporter *mChargeTransporter =
00190 StTrsFastChargeTransporter::instance(mGeometryDb, mSlowControlDb, mGasDb, mMagneticFieldDb);
00191
00192 mChargeTransporter->setChargeAttachment(true);
00193 mChargeTransporter->setGatingGridTransparency(false);
00194
00195
00196 mChargeTransporter->setExB(false);
00197
00198
00199 StTrsAnalogSignalGenerator *mAnalogSignalGenerator =
00200 StTrsSlowAnalogSignalGenerator::instance(mGeometryDb, mSlowControlDb, mElectronicsDb, mSector);
00201
00202
00203
00204 static_cast<StTrsSlowAnalogSignalGenerator*>(mAnalogSignalGenerator)->
00205 setChargeDistribution(StTrsSlowAnalogSignalGenerator::endo);
00206
00207
00208
00209
00210
00211 static_cast<StTrsSlowAnalogSignalGenerator*>(mAnalogSignalGenerator)->
00212 setElectronicSampler(StTrsSlowAnalogSignalGenerator::symmetricGaussianApproximation);
00213
00214
00215
00216
00217 mAnalogSignalGenerator->setDeltaRow(0);
00218 mAnalogSignalGenerator->setDeltaPad(2);
00219 mAnalogSignalGenerator->setSignalThreshold(1.*(.001*volt));
00220 mAnalogSignalGenerator->setSuppressEmptyTimeBins(true);
00221 mAnalogSignalGenerator->addNoise(false);
00222 mAnalogSignalGenerator->generateNoiseUnderSignalOnly(false);
00223 mAnalogSignalGenerator->setNoiseRMS(900.);
00224
00225 StTrsDigitalSignalGenerator *mDigitalSignalGenerator =
00226 StTrsFastDigitalSignalGenerator::instance(mElectronicsDb, mSector);
00227
00229
00230
00231
00232
00233 cout << "Processing" << endl;
00234
00235 int bisdet, bsectorOfHit, bpadrow;
00236
00237
00238
00239
00240
00241
00242
00243
00244
00245 ifstream ifs("/data/stix/event_14.txt",ios::in);
00246
00247 g2t_tpc_hit tpc_hit;
00248
00249
00250
00251 int firstSectorToProcess = 1;
00252 int lastSectorToProcess = 1;
00253 PR(firstSectorToProcess);
00254 PR(lastSectorToProcess);
00255 bool low = false;
00256 bool high = false;
00257
00258
00259
00260 int currentSectorProcessed = firstSectorToProcess;
00261 int lastSectorProcessed = currentSectorProcessed;
00262
00263
00264
00265 int numberOfPointsInWirePlane = 0;
00266
00267 int no_tpc_hits = 400000;
00268
00269
00270
00271 int i=0;
00272 do {
00273 i++;
00274 cout << "--> tpc_hit: " << i << endl;
00275 ifs >> tpc_hit.volume_id
00276 >> tpc_hit.de
00277 >> tpc_hit.ds
00278 >> tpc_hit.x
00279 >> tpc_hit.p;
00280 if(ifs.eof() || ifs.bad()) break;
00281
00282 whichSector(tpc_hit.volume_id, &bisdet, &bsectorOfHit, &bpadrow);
00283
00284
00285
00286
00287 if(bsectorOfHit >= firstSectorToProcess)
00288 low = true;
00289 else
00290 low = false;
00291
00292 if(bsectorOfHit <= lastSectorToProcess)
00293 high = true;
00294 else
00295 high = false;
00296
00297
00298 if(bisdet) {
00299
00300 if(i != no_tpc_hits) continue;
00301 }
00302
00303
00304 if( low &&
00305 high &&
00306 i != no_tpc_hits &&
00307 (bsectorOfHit == currentSectorProcessed)) {
00308
00309 int pID = 1;
00310
00311
00312
00313
00314
00315
00316
00317 double GEANTDriftLength = 208.55119*centimeter;
00318 double GEANTOffSet = mGeometryDb->frischGrid() - GEANTDriftLength;
00319
00320
00321
00322
00323
00324
00325 double zPosition = tpc_hit.x[2]*centimeter;
00326
00327
00328
00329 StThreeVector<double> hitPosition(tpc_hit.x[0]*centimeter,
00330 tpc_hit.x[1]*centimeter,
00331 tpc_hit.x[2]*centimeter);
00332
00333
00334
00335
00336
00337
00338
00339
00340
00341
00342
00343
00344
00345
00346
00347
00348 double beta = (bsectorOfHit>12) ?
00349 -bsectorOfHit*M_PI/6 :
00350 bsectorOfHit*M_PI/6 ;
00351
00352 double xp = hitPosition.x()*cos(beta) - hitPosition.y()*sin(beta);
00353 double yp = hitPosition.x()*sin(beta) + hitPosition.y()*cos(beta);
00354
00355 StThreeVector<double>
00356 sector12Coordinate(xp*centimeter,
00357 yp*centimeter,
00358 zPosition);
00359
00360
00361 StThreeVector<double> hitMomentum(tpc_hit.p[0]*GeV,
00362 tpc_hit.p[1]*GeV,
00363 tpc_hit.p[2]*GeV);
00364
00365
00366 StTrsChargeSegment aSegment(sector12Coordinate,
00367 hitMomentum,
00368 (fabs(tpc_hit.de*GeV)),
00369 tpc_hit.ds*centimeter,
00370 pID);
00371
00372 #ifdef HISTOGRAM
00373 tuple2[0] = static_cast<float>(bsectorOfHit);
00374 tuple2[1] = static_cast<float>(bpadrow);
00375 tuple2[2] = static_cast<float>(sector12Coordinate.x());
00376 tuple2[3] = static_cast<float>(sector12Coordinate.y());
00377 tuple2[4] = static_cast<float>(sector12Coordinate.z());
00378 secTuple->fill(tuple2);
00379 #endif
00380
00381
00382
00383
00384 #ifndef ST_NO_TEMPLATE_DEF_ARGS
00385 vector<int> all[3];
00386 #else
00387 vector<int,allocator<int> > all[3];
00388 #endif
00389
00390 #ifndef ST_NO_TEMPLATE_DEF_ARGS
00391 list<StTrsMiniChargeSegment> comp;
00392 list<StTrsMiniChargeSegment>::iterator iter;
00393 #else
00394 list<StTrsMiniChargeSegment,allocator<StTrsMiniChargeSegment> > comp;
00395 list<StTrsMiniChargeSegment,allocator<StTrsMiniChargeSegment> >::iterator iter;
00396 #endif
00397 int breakNumber = 1;
00398
00399
00400
00401 aSegment.split(mGasDb, mMagneticFieldDb, breakNumber, &comp);
00402
00403
00404 #ifndef ST_NO_TEMPLATE_DEF_ARGS
00405
00406
00407 #endif
00408
00409
00410
00411
00412 for(iter = comp.begin();
00413 iter != comp.end();
00414 iter++) {
00415
00416
00417
00418
00419 mChargeTransporter->transportToWire(*iter);
00420
00421
00422
00423
00424
00425
00426 #ifndef __sun // Bug in the sun iterators. Must Explicitly dereference!
00427 StTrsWireBinEntry anEntry(iter->position(), iter->charge());
00428
00429 #else
00430 StTrsWireBinEntry anEntry((*iter).position(), (*iter).charge());
00431 #endif
00432 mWireHistogram->addEntry(anEntry);
00433
00434 }
00435
00436
00437
00438
00439
00440
00441
00442
00443
00444
00445 numberOfPointsInWirePlane++;
00446 continue;
00447 }
00448
00449
00450 if(numberOfPointsInWirePlane) {
00451 cout << "Processing Sector: " << currentSectorProcessed << endl;
00452
00453
00454
00455
00456 cout << "--->inducedChargeOnPad()..." << endl;
00457 mAnalogSignalGenerator->inducedChargeOnPad(mWireHistogram);
00458
00459 cout << "--->sampleAnalogSignal()..." << endl;
00460 mAnalogSignalGenerator->sampleAnalogSignal();
00461
00462
00463
00464
00465 StTrsDigitalSector* aDigitalSector =
00466 new StTrsDigitalSector(mGeometryDb);
00467
00468
00469
00470 mDigitalSignalGenerator->fillSector(aDigitalSector);
00471
00472
00473
00474 cout << "--->digitizeSignal()..." << endl;
00475 mDigitalSignalGenerator->digitizeSignal();
00476
00477
00478
00479
00480
00481
00482
00483
00484 mAllTheData->mSectors[(currentSectorProcessed-1)] = aDigitalSector;
00485
00486
00487
00488
00489 mWireHistogram->clear();
00490 mSector->clear();
00491 numberOfPointsInWirePlane = 0;
00492
00493
00494 PR(currentSectorProcessed);
00495 currentSectorProcessed = bsectorOfHit;
00496
00497 PR(currentSectorProcessed);
00498 PR(bsectorOfHit);
00499
00500
00501
00502
00503
00504 }
00505
00506
00507
00508 } while(true);
00509
00510
00511
00512 #ifdef UNPACK_ALL
00513
00514
00515
00516
00517 cout << "Try Access the data!" << endl;
00518
00519
00520 PR(mAllTheData->mSectors.size());
00521 for(int isector=1; isector<=24; isector++) {
00522 int getSectorStatus =
00523 mUnPacker->getSector(isector,
00524 static_cast<StTpcRawDataEvent*>(mAllTheData));
00525
00526
00527
00528 if(getSectorStatus) continue;
00529
00530
00531 unsigned char* padList;
00532 for(int irow=1; irow<=45; irow++) {
00533 #ifdef ALL_OF_IT
00534 PR(irow);
00535 #endif
00536 int numberOfPads = mUnPacker->getPadList(irow, &padList);
00537 PR(numberOfPads);
00538
00539 if(!numberOfPads)
00540 continue;
00541
00542 for(int ipad = 0; ipad<numberOfPads; ipad++) {
00543 #ifdef ALL_OF_IT
00544
00545 PR(static_cast<int>(padList[ipad]));
00546 #endif
00547 int nseq;
00548
00549
00550 StSequence* listOfSequences;
00551 int getSequencesStatus =
00552 mUnPacker->getSequences(irow, padList[ipad], &nseq, &listOfSequences);
00553
00554 PR(getSequencesStatus);
00555 #ifdef ALL_OF_IT
00556 PR(nseq);
00557 #endif
00558 for(int kk=0; kk<nseq; kk++) {
00559 #ifdef ALL_OF_IT
00560 PR(listOfSequences[kk].length);
00561 PR(listOfSequences[kk].startTimeBin);
00562
00563 for(int zz=0; zz<listOfSequences[kk].length; zz++) {
00564 cout << " " << kk << " " << zz << '\t'
00565 << '\t' << static_cast<int>(*(listOfSequences[kk].firstAdc)) << endl;
00566 tuple[0] = static_cast<float>(isector);
00567 tuple[1] = static_cast<float>(irow);
00568 tuple[2] = static_cast<float>(padList[ipad]);
00569 tuple[3] = static_cast<float>(static_cast<int>(*(listOfSequences[kk].firstAdc)));
00570 tuple[4] = static_cast<float>(listOfSequences[kk].startTimeBin+zz);
00571 theTuple->fill(tuple);
00572
00573 listOfSequences[kk].firstAdc++;
00574 }
00575 #endif
00576
00577 }
00578
00579 }
00580
00581
00582 mUnPacker->clear();
00583
00584 }
00585
00586 }
00587 #endif
00588
00589 cout << "Got to the end of the maker" << endl;
00590
00591 #ifdef HISTOGRAM
00592 cout << "Save and close " << endl;
00593 hbookFile->saveAndClose();
00594 #endif
00595 return 0;
00596 }
00597
00598
00599
00600
00601
00602
00603
00604
00605
00606
00607
00608
00609
00610
00611
00612
00613
00614
00615
00616
00617
00618
00619
00620