using namespace std; #include "TRandom.h" #include #include "hwalldigitizer.h" #include "walldef.h" #include "hdebug.h" #include "hades.h" #include "hiterator.h" #include "hruntimedb.h" #include "hspectrometer.h" #include "hwalldetector.h" #include "hwalldigipar.h" #include "hwallgeompar.h" #include "hgeantwall.h" #include "hwallrawsim.h" #include "hevent.h" #include "hcategory.h" #include "hlocation.h" #include "hwallrawsimfilter.h" #include #include //*-- Author : D.Vasiliev //*-- Modified: 07/04/2005 F.Krizek //*-- Modified: 27/10/2002 D.Zovinec //*-- Modified: 13/02/2002 by D.Zovinec //*-- Modified: 30/11/2000 by R.Holzmann //*-- Modified: 16/12/99 by J.Bielcik //*-- Modified: 9/12/99 by V.Charriere //*-- Modified: 8/11/99 by R.Holzmann //*-- Modified: 24/10/99 by D.Vasiliev // Modified by M.Golubeva 01.11.2006 //_HADES_CLASS_DESCRIPTION ///////////////////////////////////////////////////////////////////// // // HTofDigitizer digitizes data, puts output values into // raw data category // // Need to implement quenching? // ///////////////////////////////////////////////////////////////////// void HWallDigitizer::initParContainer() { pWallDigiPar=(HWallDigiPar *)gHades->getRuntimeDb()->getContainer("WallDigiPar"); } HWallDigitizer::HWallDigitizer(void) { fGeantCat = 0; fRawCat = 0; pWallDigiPar = 0; pWallGeomPar = 0; fLoc.set(1,0); iterGeant = 0; iterWallRaw = 0; } HWallDigitizer::HWallDigitizer(const Text_t *name,const Text_t *title) : HReconstructor(name,title) { fGeantCat = 0; fRawCat = 0; pWallDigiPar = 0; pWallGeomPar = 0; fLoc.set(1,0); iterGeant = 0; iterWallRaw = 0; } HWallDigitizer::~HWallDigitizer(void) { } Bool_t HWallDigitizer::init(void) { time_t curtime; initParContainer(); fGeantCat = gHades->getCurrentEvent()->getCategory(catWallGeantRaw); if (!fGeantCat) { fGeantCat = gHades->getSetup()->getDetector("Wall")->buildCategory(catWallGeantRaw); if (!fGeantCat) return kFALSE; else gHades->getCurrentEvent()->addCategory(catWallGeantRaw,fGeantCat,"Wall"); } fRawCat = gHades->getCurrentEvent()->getCategory(catWallRaw); if (!fRawCat) { HWallDetector* wall=(HWallDetector*)(gHades->getSetup()->getDetector("Wall")); if(wall){ fRawCat=wall->buildMatrixCategory("HWallRawSim",0.5F); } if (!fRawCat) return kFALSE; else gHades->getCurrentEvent()->addCategory(catWallRaw,fRawCat,"Wall"); } iterGeant = (HIterator *)fGeantCat->MakeIterator("native"); iterWallRaw = (HIterator *)fRawCat ->MakeIterator("native"); time(&curtime); return kTRUE; } Int_t HWallDigitizer::execute(void) { const Float_t quantEff = 0.24; // PMT quantum efficiency. const Float_t photYield = 11100.0; // photon yield in scintillator (phot/MeV). const Float_t timeResWall[3] = {0.42, 0.52, 0.67}; //values based on apr07 analysis //{0.35, 0.42, 0.60}; // time resolution in the centres of the // small(40x40mm), middle (80x80mm) and // large (160x160mm) cells [ns]. const Float_t relAmplResol = 0.08; // sigma of Gaus distribution. const Float_t minEnerRelease = 1.8; // minimum energy release (MeV/cm) HGeantWall* geant = 0; HWallRawSim* raw = 0; HWallRawSimFilter fRawFilter; TString cellName,str; Float_t slopeTDC; Int_t thrCFD, thrADC; Int_t cell = -1; fLoc.set(1,-1); Int_t numTrack, numTrack1 = -1, numTrack2 = -1; //track numbers Float_t trackLen; Float_t time, charge; Int_t timeCh, chargeCh; Float_t prevTime, prevCharge; Float_t timeResol, amplResol, chargeRef; Float_t geaTof = 0.; Float_t geaTof1[383]; Float_t geaTof2[383]; Float_t geaEner = 0.; Float_t geaX = 0.; Float_t geaY = 0.; // used by Tofino only Float_t geaMom = 0.; iterGeant->Reset(); // this works only in split mode=2 // (for 0 and 1 the iterator must be recreated) while ((geant = (HGeantWall *)iterGeant->Next()) != 0) { fLoc[0] = geant->getCell(); raw = (HWallRawSim*) fRawCat->getObject(fLoc); // test if cell in use if(raw) { raw->incNHit(); // increment counter for additional hits numTrack1 = raw->getNTrack1(); numTrack2 = raw->getNTrack2(); prevTime = raw->getTime(); prevCharge = raw->getCharge(); } else { prevTime = 100000.; prevCharge = 0.; raw = (HWallRawSim*) fRawCat->getNewSlot(fLoc); // get new slot if(!raw) continue; raw = new(raw) HWallRawSim; raw->setNHit(1); } cell = fLoc[0]; slopeTDC = pWallDigiPar->getTDC_Slope(cell); geant->getHit(geaEner,geaX,geaY,geaTof,geaMom,trackLen); numTrack = geant->getTrack(); cellName = pWallGeomPar->getCellName(fLoc[0]); str=cellName.Replace(0,3,"",0); str.Replace(1,3,"",0); time = geaTof; timeResol = timeResWall[atoi(str)-1]; time = gRandom->Gaus(time,timeResol); timeCh = (Int_t) (time/slopeTDC); if (timeCh < 0) timeCh = 0; if (timeCh > 4095) timeCh = 4095; time = ( Float_t ) (timeCh); if(raw->getNHit() > 1){ if(geaTof < geaTof1[fLoc[0]]){ numTrack2 = numTrack1; numTrack1 = numTrack; geaTof2[fLoc[0]] = geaTof1[fLoc[0]]; geaTof1[fLoc[0]] = geaTof; } else { if(geaTof < geaTof2[fLoc[0]]){ numTrack2 = numTrack; geaTof2[fLoc[0]] = geaTof; } } } if(time > prevTime) time = prevTime; //charge = geaEner*photYield*quantEff*0.5*(1 - cos(ar*deg2rad))*exp(-(hl-geaX)/al); charge = geaEner*photYield*quantEff; amplResol = charge*relAmplResol; charge = gRandom->Gaus(charge,amplResol); chargeRef = photYield*quantEff*minEnerRelease; chargeCh = (Int_t) ((charge/chargeRef)*256. + prevCharge); if (chargeCh < 0) chargeCh = 0; if (chargeCh > 4095) chargeCh = 4095; charge = (Float_t)chargeCh; raw->setTime(time); raw->setCharge(charge); raw->setCell((Int_t) fLoc[0]); if(raw->getNHit() > 1){ raw->setNTrack1(numTrack1); raw->setNTrack2(numTrack2); } else { raw->setNTrack1(numTrack); raw->setNTrack2(-1); geaTof1[fLoc[0]] = geaTof; geaTof2[fLoc[0]] = 100000.; } } // Exclusion of hits with charge less than ADC threshold. // Time of hits with charge less than CFD threshold is set // to zero. These hits are excluded later in the hitfinder. iterWallRaw->Reset(); while ( (raw = (HWallRawSim *)iterWallRaw->Next()) != NULL) { fLoc[0] = raw->getCell(); time = raw->getTime(); charge = raw->getCharge(); // overflow suppression if(((Int_t)time) >= 4095) raw->setTime(0.0); cell=fLoc[0]; thrCFD = (Int_t)pWallDigiPar->getCFD_Threshold(cell); thrADC = (Int_t)pWallDigiPar->getADC_Threshold(cell); // CFD and ADC thresholds if(((Int_t)charge) < thrCFD){ raw->setTime(0.0); if(((Int_t)charge) < thrADC){ raw->setCharge(0.0); } } } fRawCat->filter(fRawFilter); return 0; } ClassImp(HWallDigitizer)