|
//____________________________________________________________________ // // BrSiDigModule // // Digitisation Module class for Silicon wafers. You have the // option to use calibrations rather than the ASCII file parameters. // Set UseCalibrations to kTRUE. // //____________________________________________________________________ // // $Id: BrSiDigModule.cxx,v 1.15 2002/04/22 16:58:35 cholm Exp $ // $Author: cholm $ // $Date: 2002/04/22 16:58:35 $ // $Copyright: (c) 2001 BRAHMS Collaboration <brahmlib@rcf.rhic.bnl.gov> // #include <iomanip.h> #include <fstream.h> #ifndef ROOT_TRandom #include "TRandom.h" #endif #ifndef ROOT_TSystem #include "TSystem.h" #endif #ifndef ROOT_TMath #include "TMath.h" #endif #ifndef BRAT_BrGeometryDbManager #include "BrGeometryDbManager.h" #endif #ifndef BRAT_BrParameterDbManager #include "BrParameterDbManager.h" #endif #ifndef BRAT_BrCalibrationManager #include "BrCalibrationManager.h" #endif #ifndef BRAT_BrDetectorList #include "BrDetectorList.h" #endif #ifndef BRAT_BrDetectorVolume #include "BrDetectorVolume.h" #endif #ifndef BRAT_BrSiParameters #include "BrSiParameters.h" #endif #ifndef BRAT_BrDataTable #include "BrDataTable.h" #endif #ifndef BRAT_BrTableNames #include "BrTableNames.h" #endif #ifndef BRAT_BrGeantHit #include "BrGeantHit.h" #endif #ifndef BRAT_BrSiDigModule #include "BrSiDigModule.h" #endif #ifndef BRAT_BrSiDig #include "BrSiDig.h" #endif #ifndef BRAT_BrPathManager #include "BrPathManager.h" #endif #ifndef BRAT_BrSiCalibration #include "BrSiCalibration.h" #endif #ifndef ROOT_TDirectory #include "TDirectory.h" #endif ClassImp(BrSiDigModule); //__________________________________________________________________ BrSiDigModule::BrSiDigModule(void) { // default Constructor. Don't use this one. Only here for IO SetState(kSetup); fVolumeParams = 0; fDetectorParams = 0; fCalibration = 0; UseCalibrations(); } //____________________________________________________________________ BrSiDigModule::BrSiDigModule(const Char_t *name, const Char_t *title) :BrModule(name,title) { // Normal Named constructor SetState(kSetup); fVolumeParams = 0; fDetectorParams = 0; fCalibration = 0; UseCalibrations(); } //___________________________________________________________________________ BrSiDigModule::~BrSiDigModule(void) { // Destructor } //___________________________________________________________________________ void BrSiDigModule::Init(void) { // Initialize entry point. This is not called at creation time but // this would be a good place to setup histograms, statistics // variables etc. This could also be done at the time of the // constructor. This might likely be more usefull since one has the // opportunity to setup variables before this call (e.g. detector // parameters or alike). // // open a file for detector location to adc location conversion SetState(kInit); #if 0 BrGeometryDbManager *gGeomDb = BrGeometryDbManager::Instance(); fVolumeParams = (BrDetectorVolume*) gGeomDb->GetDetectorVolume(const_cast<Char_t*>("BrDetectorVolume"), const_cast<Char_t*>(GetName())); #endif BrParameterDbManager *gParamDb = BrParameterDbManager::Instance(); fDetectorParams = (BrSiParameters*) gParamDb->GetDetectorParameters(const_cast<Char_t*>("BrSiParameters"), GetName()); Int_t ActiveDetTable[6][6]; Int_t side, loc; Int_t i,j; // initialize ActiveDetTable to unrialistic number for (i = 0; i<6; i++) { for (j = 0; j<6; j++) { ActiveDetTable[i][j]=499; } } const Char_t *bratsys = BrPathManager::Instance()->GetDataDir(); // const Char_t *bratsys = gSystem->Getenv("BRATSYS"); char multdetector_file[64]; sprintf(multdetector_file,"%s/params/ActiveMultDet.dat",bratsys); ifstream activeDet(multdetector_file); if (!activeDet) { cout << "ActiveMultDet.dat file could not be found." << endl; cout << "Use the default setting for coverting Geant channel" << endl << " to Brat data table channel number." << endl; cout << "User need to check the output file carefully." << endl; Int_t k=0; for (i = 0; i<6; i++) { for (j = 0; j<6; j++) { ActiveDetTable[i][j]=k; k=k+8; } } } else { Int_t k; Char_t HistName[32], mname[5]; while (activeDet >> HistName >> k) { if (HistName[0]=='s') { sscanf(HistName,"%2s%d_%d",mname,&side,&loc); ActiveDetTable[side][(loc-1)]=k; } } } // Now with Active Detector Table, convet the geant sub volume number to // actual ADC Channel number for (i = 0; i<36; i++) { side = (Int_t)floor((i/6)); loc = (Int_t)fmod(i,6); GeantFERATable[i]=ActiveDetTable[side][loc]; } //cout << "end of Init" << endl; if (DebugLevel() > 9) Print("T"); if (fUseCalibrations) { // Calibrations #ifdef BR_MULT_CAL_TMP fCalibration = BrSiTmpCalibration::Instance(); #else BrCalibrationManager *manager = BrCalibrationManager::Instance(); fCalibration = static_cast<BrSiCalibration*> (manager->Register("BrSiCalibration", BrDetectorList::GetDetectorName(kBrSi))); #endif if (!fCalibration) Failure("Init", "didn't get calibration from manager"); fCalibration->Use("pedestal"); fCalibration->Use("pedestalWidth"); // I think we need to flag this parameter, but I'm not really sure. fCalibration->Use("adcGap"); fCalibration->Use("adcGain"); fCalibration->Use("conversionFunc"); // These are geometry things proper, but put here for now, until a // concensus of how to do the geometry database. fCalibration->Use("ringMap"); fCalibration->Use("rowMap"); fCalibration->Use("ringPostion"); fCalibration->Use("tilt"); } } //___________________________________________________________________________ void BrSiDigModule::DefineHistograms() { if (GetState() != kInit) { Stop("DefineHistograms", "Must be called after Init"); return; } if (!fDetectorParams) { Failure("DefineHistograms", "no detector parameters"); return; } // Get some parameters Int_t nSingles = fDetectorParams->GetNoModules(); TDirectory* saveDir = gDirectory; TDirectory* histDir = gDirectory->mkdir(GetName()); histDir->cd(); // fSingleVsHits = new TH2D("singleVsHits", "Single Hits", nSingles, 1, nSingles+1, 75, 0,75); fSingleVsEnergy = new TH2D("singleVsEnergy", "Single energies", nSingles, 1, nSingles+1, 100, 0, 0.03); fSingleVsEnergyPerHit = new TH2D("singleVsEnergyPerHit", "Single energies per hit", nSingles, 1, nSingles+1, 100, 0, 0.03); fSingleVsAdc = new TH2D("singleVsAdc", "Single ADCs", nSingles, 1, nSingles+1, 100, 0, 2500); fSingleVsAdcPerHit = new TH2D("singleVsAdcPerHit", "Single ADCs per hit", nSingles, 1, nSingles+1, 100, 0, 250); gDirectory = saveDir; } //___________________________________________________________________________ void BrSiDigModule::Event(BrEventNode* inNode, BrEventNode* outNode) { // Event method to be called once per event. // This performes the actual slow simulation // of the hits on the MultSi. // // Description of method: // Each tile can in fact be hit by multiple hits // 1) First a loop over all hits is made to see how many // hits there are per waifer. If the is only a single (charged) // hit The adc value is calculated from the de/dx. // 2) For multiple hits each end of the tube is dealt with seperately. // The fastest signals is found. If the following hits are within // the resolving time / gate width (100n sec) the adc values are // added to those from the first hit. // SetState(kEvent); Debug(5, "Event", " Inside Digitize of %s", GetName()); // inNode->ListObjects(); // // Get simulation parameters if(!fDetectorParams) { Stop("Event","parameters are not set"); return ; } if(DebugLevel() > 14) fDetectorParams->ListParameters(); // Numbers for pulser correction. Double_t l1 = 25; // Double_t l2 = 400; // Year 2000 data Double_t l2 = 250; // Year 2001 data // Get some of the parameters from the file Int_t nWafers = fDetectorParams->GetNoWafers(); Int_t nRings = fDetectorParams->GetNoRings(); //cout << "noWafer = " << noWafers << endl; // // Prepare the output table // BrSiDig* digit = new BrSiDig(BRTABLENAMES kDigSi,"sim"); outNode->AddObject(digit); // // Decode the inputtable. // Char_t tableName[7][32]; sprintf(tableName[0], "%s STRA", BRTABLENAMES kGeantHits); sprintf(tableName[1], "%s STRB", BRTABLENAMES kGeantHits); sprintf(tableName[2], "%s STRC", BRTABLENAMES kGeantHits); sprintf(tableName[3], "%s STRD", BRTABLENAMES kGeantHits); sprintf(tableName[4], "%s STRE", BRTABLENAMES kGeantHits); sprintf(tableName[5], "%s STRF", BRTABLENAMES kGeantHits); sprintf(tableName[6], "%s STRG", BRTABLENAMES kGeantHits); // Variables for use in loop over hits Float_t badTdc = 1000000; Int_t nHits[nWafers]; Float_t tdc[nWafers]; Float_t adc[nWafers]; // Loop over Table by each strip stra ---> strg. Int_t j; for (Int_t i = 0; i < 7; i++) { Info(10, "Event", "Inspecting table # %d: %s", i, tableName[i]); // Get the table from the BrEvent object. BrDataTable* hits = (BrDataTable*)inNode->GetObject(tableName[i]); // Test if there were any hits at all if(!hits) { Abort("Event", "No GeantHit Table for %s - Should not happen", tableName[i]); return; } // Get total number of hits in this table. (each strip) Int_t noHits = hits->Entries(); Info(10 ,"Event", " Found %d hits in %s", noHits, tableName[i]); // Initialize the arrays. Notice that the TObjArrays are // initaliazed with noHits/noWafers to minimize the number of new // calls for (j = 0; j < nWafers; j++) { nHits[j] = 0; tdc[j] = badTdc; adc[j] = 0; } // First find the earliest TDC signal, so we can emulate a time // gate. Debug(10, "Event", "* Loop over the %d hits (TDC)", noHits); for (j = 0; j < noHits; j++) { BrGeantHit* hit = (BrGeantHit*)(hits->At(j)); // parictular strip number Int_t id = hit->Isub() - 1; Int_t rid = GeantFERATable[id] + i; if (GeantFERATable[id] == 499) { Debug(20, "Event", "** Active map says wafer %d (FERA # %d) isn't on", id, GeantFERATable[id]); continue; } if (fUseCalibrations && (fCalibration->GetRingMap(rid) < 0 || fCalibration->GetRingMap(rid) > nRings || fCalibration->GetRingMap(rid) > 100)) { Debug(20, "Event", "Ring map says si %d isn't on", rid); continue; } Debug(20, "Event", "** Found a hit at %d (channel # %d)", id, rid); // Find the earliest time tdc[id] = TMath::Min(tdc[id - 1], hit->Tof()); } // Now loop over the hitTables to get the actual hit ADC value to // store. The time gate is assumed to be 100 nsec. wide. Debug(10, "Event", "* Loop over the %d hits (ADC)", noHits); for (j = 0; j < noHits; j++) { BrGeantHit* hit = (BrGeantHit*)hits->At(j); // parictular strip number Int_t id = hit->Isub() - 1; Int_t rid = GeantFERATable[id] + i; if (GeantFERATable[id] == 499) { Debug(20, "Event", "** Active map says wafer %d (FERA # %d) isn't on", id, GeantFERATable[id]); continue; } if (fUseCalibrations && (fCalibration->GetRingMap(rid) < 0 || fCalibration->GetRingMap(rid) > nRings || fCalibration->GetRingMap(rid) > 100)) { Debug(20, "Event", "Ring map says si %d isn't on", rid); continue; } // If the TOF for this hit is outside the time gate, we go // onto the next hit, silently ignoring this hit. if (hit->Tof() >= tdc[id] + 100) continue; Debug(20, "Event", "** Found a hit at %d (channel # %d)", id, rid); // Otherwise, we sum up the ADC value, and increment counter adc[id] += hit->Dedx(); nHits[id]++; } // Having accumelated all the (relevant) ADCs we can now make the // digitised outputs for (j = 0; j < nWafers; j++) { // id is the number of channel in RAW File ADC Channel // i is to give right channel number for each strip stra --> strg Int_t id = GeantFERATable[j] + i; if (GeantFERATable[j] >= 499) { Debug(20, "Event", "** Active map says wafer %d (FERA # %d) isn't on", id, GeantFERATable[id]); continue; } Debug(10, "Event", "* Ch: %3d W: %3d H: %4d E: %8f T: %8f", id, j + 1, nHits[j], adc[j], tdc[j]); if (HistOn()) { fSingleVsHits->Fill(id + 1,nHits[j]); fSingleVsEnergy->Fill(id + 1,adc[j]); if (nHits[j] != 0) fSingleVsEnergyPerHit->Fill(id + 1, adc[j] / nHits[j]); } if (!fUseCalibrations) { adc[j] *= fDetectorParams->GetAdcGain(); adc[j] += fDetectorParams->GetAdcOffset(); } else { if (fCalibration->GetRingMap(id) < 0 || fCalibration->GetRingMap(id) > nRings || fCalibration->GetRingMap(id) > 100) { Debug(20, "Even", "** Channel %d isn't on", id); continue; } if (nHits[j] != 0) { Debug(15, "Event", "*** Ch: %3d W: %3d Old E: %8f Old A: %8f", id, j+1, adc[j] * fDetectorParams->GetAdcGain(), adc[j] * fDetectorParams->GetAdcGain() + fDetectorParams->GetAdcOffset()); adc[j] *= 1000. / fCalibration->GetAdcGain(id); Debug(10, "Event", "* Ch: %3d W: %3d Gain A: %8f", id, j+1, adc[j]); Double_t p0 = fCalibration->GetPulserFuncPar(id,0); Double_t p1 = fCalibration->GetPulserFuncPar(id,1); Double_t p2 = fCalibration->GetPulserFuncPar(id,2); Double_t p3 = fCalibration->GetPulserFuncPar(id,3); Double_t a = p0 + p1 * l1 + p2 * l1 * l1; Double_t b = p1 + 2 * p2 * l1; Double_t c = b + 2 * p3 * (l2 - l1); Double_t d = a + b * (l2 - l1) + p3 * (l2 - l1) * (l2 - l1); Debug(15, "Event", "* Ch: %3d W: %3d Puls a=%f b=%f c=%f d=%f", id, j + 1, a, b, c, d); if (adc[j] < a) { // f(x) = p0 + p1 * x + p2 * x^2 // // -p1 + sqrt(p1*p1 - 4 * p2 * (p0 - y)) // g(y) = ------------------------------------- // 2 * p2 Double_t det = p1 * p1 - 4 * p2 * (p0 - adc[j]); Debug(15, "Event", "D (%4f^2 - 4 * %4f * (%4f - %4f) = %4f", float(p2), float(p1), float(p0), adc[j], float(det)); adc[j] = (-p1 + TMath::Sqrt(det)) / 2 / p2; } else if (adc[j] < d) { // f(x) = a + b * (x -l1) + p2 * (x - l1)^2 // // -b + sqrt(b*b - 4 * p3 * (a - y)) // g(y) = l1 + ------------------------------------- // 2 * p3 Double_t det = b * b - 4 * p3 * (a - adc[j]); Debug(15, "Event", "D (%4f^2 - 4 * %4f * (%4f - %4f) = %4f", float(b), float(p3), float(a), adc[j], float(det)); adc[j] = l1 + (-b + TMath::Sqrt(det)) / 2 / p3; } else // f(x) = d + c * (x - l2) // // x - d // g(y) = l2 + ----- // c adc[j] = l2 + (adc[j] - d) / c; Debug(10, "Event", "* Ch: %3d W: %3d Puls A: %8f", id, j+1, adc[j]); } adc[j] += gRandom->Gaus(fCalibration->GetPedestal(id), fCalibration->GetPedestalWidth(id)); } if (nHits[j] != 0) Debug(10, "Event", "* C: %3d W: %3d A: %8f", id, j+1, adc[j]); if (HistOn()) { fSingleVsAdc->Fill(id + 1, adc[j]); if (nHits[j] != 0) fSingleVsAdcPerHit->Fill(id + 1, adc[j] / nHits[j]); } digit->SetAdc(id + 1,(Int_t)adc[j]); // digit->SetId(id + 1, id + 1); } // End of loop over hits in table star,...,strg } // end of loop over each strip type star--> strg // // List the contents of the added datatable, if the debug level // is high enough. // if(DebugLevel() > 19) { cout << "Output from BrSiDigModule:" << endl; digit->Print("R"); } } //______________________________________________________________________ void BrSiDigModule::ListDetectorParameters() const { // // List the current value of the digitization parameters on // standard out. // if( fDetectorParams != NULL){ cout << "The name is " << GetName() << endl; fDetectorParams->ListParameters(); } else cout << "No parameters set for this Digitisation module" << endl; } //____________________________________________________________________________ void BrSiDigModule::Print(Option_t* option) const { // Standard information printout. // // Options: See BrModule::Print // TString opt(option); opt.ToLower(); BrModule::Print(option); if (opt.Contains("d")) cout << endl << " Original author: Hironori Ito, Version 0.2, Jul 15 2000" << endl << " $Author: cholm $" << endl << " $Date: 2002/04/22 16:58:35 $" << endl << " $Revision: 1.15 $ " << endl << endl << "-------------------------------------------------" << endl << " " << (!fUseCalibrations ? "Not using" : "Using") << " real calibrations" << endl; if (opt.Contains("t")) { cout << " BRAG subid to FERA ADC channel conversion: " << endl; for (Int_t i = 0; i < 36 ; i++) { cout << " " << setw (3) << i + 1<< ": " << setw(3) << GeantFERATable[i] << flush; if ((i + 1) % 6 == 0) cout << endl; } cout << endl; } } // $Log: BrSiDigModule.cxx,v $ // Revision 1.15 2002/04/22 16:58:35 cholm // Some fixes for the modes using the calibrations rahter than // parameters. Also changed a lot of // // if (DebuLevel() > ...) // cout << .... // // to // // Debug(..., ..., ...) // // and similar for Info. // // Revision 1.14 2002/01/03 19:51:34 cholm // Prepared to use BrTableNames class (or perhaps BrDetectorList) for table names // // Revision 1.13 2001/12/14 15:40:52 cholm // Updated for new BrModule::Info method // // Revision 1.12 2001/10/08 11:28:50 cholm // Changed to use new DB access classes // // Revision 1.11 2001/09/20 13:40:43 cholm // Added some documentation and such. // // Revision 1.10 2001/09/01 10:13:52 cholm // The digitisation using calibrations rather than the pseudo numbers in // the ASCII parameter file is now fully implmented. It turns out, that // using the calibrations gives far more consistent (with real data) // numbers than using the parameter file. Use the method // UseCalibrations(kTRUE) to turn on this feature. Perhaps this kind of // thing should be implemented for all detectors. // // Revision 1.9 2001/08/30 12:42:50 cholm // Added the option to use calibrations rather than the ASCII file // parameters. This also simulates the non-linear behaviour of the // readout electronics on the silicon detectors. This feature is not // fully operational yet (should be later today though). // // Revision 1.8 2001/08/24 18:50:45 cholm // Better histogram limits // // Revision 1.7 2001/08/21 12:42:44 cholm // Increased some debugging levels // // Revision 1.6 2001/08/12 13:13:51 cholm // Fixed a few things. // // Revision 1.5 2001/08/12 09:55:04 cholm // Moved description in front of CVS keywords (for HTML documentation) // // Revision 1.4 2001/08/12 09:54:17 cholm // Added some debug output, fixed boundaries on histograms, and some // general clean up. // // Revision 1.3 2001/08/10 15:48:09 cholm // Added histograms and state information // // Revision 1.2 2001/08/10 14:22:39 cholm // Cleaned up the Event method. This method was very inefficient (many // nested loops). Also removed some erroneous failures (due to missing // non-existing geometry data). // // Revision 1.1.1.1 2001/06/21 14:55:07 hagel // Initial revision of brat2 // // Revision 1.5 2001/06/04 13:37:11 cholm // Changes to use BrPathManager, BrVersion, and perpare to use BrFileTag. // // Revision 1.4 2001/05/03 18:45:48 cholm // Fixed a few superflous warnings, and uninitialised array in BrTileParameters // // Revision 1.3 2001/03/07 12:20:18 cholm // * Made the method BrModule::Info() const obsolete in favour of // BrModule::Print(Option_t* option="B") const. // // Revision 1.2 2001/02/02 19:50:48 videbaek // Made Warning contingent on Verbose // // Revision 1.1 2001/01/29 20:25:54 cholm // Renamed class BrDigitizeMultSi to BrSiDigModule to fit naming scheme. // // Revision 1.13 2001/01/17 02:01:48 hagel // Fixed float to int warning that made into error // // Revision 1.12 2000/11/22 21:16:44 videbaek // Moved geometry and parameter calls from constructor to Initm method // // Revision 1.11 2000/09/21 20:49:43 cholm // Several changes, since this was really messy code. Most important, moved // the arrays of constant size to pointers. // // Revision 1.10 2000/08/19 21:01:43 hito // Minor update. // // Revision 1.9 2000/08/17 14:50:43 hito // *** empty log message *** // // Revision 1.8 2000/08/17 13:14:25 sanders // Several enhancements to RDO mult modules // // Revision 1.7 2000/07/26 22:08:37 hito // Mapping of digitization for si and tile were changed to reflect the raw data // and the change in gbrahms. // // Revision 1.6 2000/05/24 13:58:05 cholm // Added some new method. Most noticably "BrSiDig::GetNumberOfChannels()", // and "BrDigMultSi::GetName()". Hito has maked "BrDigMultSi" as obsolete, // should probably be put in the Attic. // Also. "BrDitigitizeMultSi" now uses "BrSiDig", similar to what // "BrDigitizeMultTile" does. // Put a body into "BrSiDig::List()". // // Revision 1.5 2000/04/07 20:33:07 cholm // Added targets "realclean" and "distclean", to _really_ clean up // the working directories. Hmm, see the log message in BrGeantDigitize // from the same date as this commit - stupid me! // // Revision 1.4 2000/04/06 20:21:13 cholm // Corrected bugs in BrSiDigModule // // Revision 1.3 2000/03/23 15:44:12 cholm // Bug corrections to Mult classes // // Revision 1.2 2000/03/21 21:22:24 cholm // Several changes: A few hacks where needed to compile on Digital Unix, noticably in my - Christian Holm - classes for the DB and Exceptions. Proberly still need to do something to some of Konstantins stuff. Also, I removed a lot of warnings fromt the compiler, by teddying up the code. Several places, old C-style indecies in for loops were assumed. Corrected. Unused variables have been commented out. // // Revision 1.1 2000/03/17 16:17:57 cholm // Added classes BrSiDigModule and BrDetectorPramMultSi for GBRAHMS simulations // // |
||||||
|
This page automatically generated by script docBrat by Christian Holm |
Copyright ; 2002 BRAHMS Collaboration
<brahmlib@rcf.rhic.bnl.gov>
|