////////////////////////////////////////////////////////////////////////////// // // $Id: $ // //*-- Author : Laura Fabbietti //*-- Revised : Martin Jurkovic 2010 // //_HADES_CLASS_DESCRIPTION ////////////////////////////////////////////////////////////////////////////// // // HRichDigitizer // // For each digitized pad the corresponding track numbers of the // direct hits and photons are stored in a linear category // (catRichTrack<-HRichTrack) that is made sortable by the pad address. // Each pad object of the HRichCalSim class has 2 members, nTrack1, nTrack2, // that give the indexes of the first and the last track numbers in the linear // category (catRichTrack) corresponding to a pad. // This method allows to store as many particle track numbers as the pad // " sees ". // here are listed the values of some useful variables // # means non costant value. // Each value refers to the 1 fired pad. // // realID=gHades->getEmbeddingRealTrackId() // // track Number Flag energy // Cheren. Phot. # 0 # // Feedback Phot. -5 0 0 // Direct Hits # 1 0 // Noise Hits -1 0 0 // REAL Hits (embedding) realID 0 0 //--------------------------------------------------------- // // // ----------------------- ------------- // | HRichUnpacker | | HGeantRich | // | (embedding mode) | \ ------------- // | | \ || // ----------------------- \ || input to digitizer // \ \/ // ------------- read real ------------------ // | HRichCalSim | ------------> | | // ------------- <----------- | | // | HRichDigitizer | // ------------- | | // | HRichTrack | <----------- | | // ------------- ------------------ // write output // // In the case of TRACK EMBEDDING of simulated tracks into // experimental data the real data are written by the HRichUnpacker into // HRichCalSim category. The real hits are taken into // account by the digitizer (adding of charges, adding track numbers to HRichTrack). // The embedding mode is recognized automatically by analyzing the // gHades->getEmbeddingMode() flag. // Mode ==0 means no embedding // ==1 realistic embedding (first real or sim hit makes the game) // ==2 keep GEANT tracks (as 1, but GEANT track numbers will always // win against real data. besides the tracknumber the output will // be the same as in 1) // // Needed parameter containers: // - HRichCalPar (used for noise/threshold calculation) // - HRichGeometryPar // - HRichDigitisationPar // //////////////////////////////////////////////////////////////////////////// #include "TF1.h" #include "TMath.h" #include "TRandom.h" #include "TVector.h" #include "hades.h" #include "hcategory.h" #include "hdebug.h" #include "hevent.h" #include "heventheader.h" #include "hgeantrich.h" #include "hlinearcatiter.h" #include "hmatrixcatiter.h" #include "hparset.h" #include "hrichcalpar.h" #include "hrichcalparcell.h" #include "hrichcalsim.h" #include "hrichdetector.h" #include "hrichdigitisationpar.h" #include "hrichdigitizer.h" #include "hrichgeometrypar.h" #include "hrichpad.h" #include "hrichpadcorner.h" #include "hrichpadfilter.h" #include "hrichtrack.h" #include "hrichwiresignal.h" #include "hruntimedb.h" #include "hspectrometer.h" #include "richdef.h" #include #include #include using namespace std; ClassImp(HRichDigitizer) //---------------------------------------------------------------------------- HRichDigitizer::HRichDigitizer(const Text_t* name, const Text_t* title, Bool_t kNoise, Float_t slope, Bool_t oem) : HReconstructor(name, title) { setDefaults(); isOEM = oem; isActiveNoise = kNoise; fSlopeCorrection = slope; } //============================================================================ //---------------------------------------------------------------------------- HRichDigitizer::HRichDigitizer() { setDefaults(); } //============================================================================ //---------------------------------------------------------------------------- void HRichDigitizer::setDefaults() { catRichPhoton = NULL; catRichDirect = NULL; catTrack = NULL; catCal = NULL; iterRichPhoton = NULL; iterRichDirect = NULL; iterRichTrack = NULL; iterRichCal = NULL; pDigitisationPar = NULL; pGeometryPar = NULL; pCalPar = NULL; ga = NULL; isOEM = kFALSE; isActiveNoise = kFALSE; countFBphot = 0; countNoisePad = 0; fWiresNr = 0; fWiresDist = 0.0; distWirePads = 0.0; fYShift = 0.0; binsQE = 0; fChargePerChannel = 0.0; fChargeScaling = 1.0; fElectronsNr = 0.0; fFactor1 = 0.0; fFactor1Sig = 0.0; fFactor2 = 0.0; fFactor2Sig = 0.0; fIncreaseNoise = 1.0; fParam1 = 0.0; fParam2 = 0.0; fQupper = 0.0; fSigmaValue = 0.0; fExpSlope .Set(6); fExpSlope .Reset(0.0); fSlopeCorrection = 1.0; photlength .Set(18); photlength .Reset(0.0); photeffic .Set(18); photeffic .Reset(0.0); for (Int_t sec = 0; sec < 6; ++sec) { correction[sec].Set(18); correction[sec].Reset(0.0); } noiseProb = 0.0; for (Int_t i = 0; i < 9; ++i) fDigitPadMatrix[i] = 0.0; for (Int_t i = 0; i < 9; ++i) fDigitPadMatrix[i] = 0.0; loc1 .set(1, 0); locat.set(3, 0, 0, 0); loc .set(3, 0, 0, 0); } //============================================================================ //---------------------------------------------------------------------------- HRichDigitizer::~HRichDigitizer() { if (iterRichPhoton) { delete iterRichPhoton; iterRichPhoton = NULL; } if (iterRichDirect) { delete iterRichDirect; iterRichDirect = NULL; } if (iterRichTrack) { delete iterRichTrack; iterRichTrack = NULL; } if (iterRichCal) { delete iterRichCal; iterRichCal = NULL; } fChargeOnWireList.Delete(); fTrackNumberList.Delete(); if (ga) { delete ga; ga = NULL; } } //============================================================================ //---------------------------------------------------------------------------- Bool_t HRichDigitizer::init() { HRichDetector* pRichDet = static_cast(gHades->getSetup() ->getDetector("Rich")); HRuntimeDb* rtdb = gHades->getRuntimeDb(); if (kTRUE == isActiveNoise) { Info("init", "------------------------------------------------------"); Info("init", "RICH digitizer NOISE MODE is ON"); Info("init", "------------------------------------------------------"); } else { Info("init", "------------------------------------------------------"); Info("init", "RICH digitizer NOISE is OFF"); Info("init", "------------------------------------------------------"); } // Set pointer to RICH calibration parameters pCalPar = static_cast(rtdb->getContainer("RichCalPar")); if (NULL == pCalPar) { Error("init", "Initialization of calibration parameters failed, returning..."); return kFALSE; } // Set pointer to RICH geometry parameters pGeometryPar = static_cast(rtdb->getContainer("RichGeometryParameters")); if (NULL == pGeometryPar) { pGeometryPar = new HRichGeometryPar; rtdb->addContainer(pGeometryPar); } if (NULL == pGeometryPar) { Error("init", "Initialization of geometry parameters failed, returning..."); return kFALSE; } // Set pointer to RICH digitization parameters pDigitisationPar = static_cast(rtdb->getContainer("RichDigitisationParameters")); if (NULL == pDigitisationPar) { pDigitisationPar = new HRichDigitisationPar; rtdb->addContainer(pDigitisationPar); } if (NULL == pDigitisationPar) { Error("init", "Initialization of digitization parameters failed, returning..."); return kFALSE; } // Initialize geant rich cherenkov photon category and set appropriate iterator catRichPhoton = gHades->getCurrentEvent()->getCategory(catRichGeantRaw); if (NULL == catRichPhoton) { Error("init", "Initializatin of Cherenkov photon category failed, returning..."); return kFALSE; } iterRichPhoton = static_cast(catRichPhoton->MakeIterator()); if (NULL == iterRichPhoton) { Error("init", "Can not create catRichPhoton iterator, returning..."); return kFALSE; } // Initialize geant category for direct hits and set appropriate iterator catRichDirect = gHades->getCurrentEvent()->getCategory(catRichGeantRaw + 1); if (NULL == catRichDirect) { Error("init", "Initialization of geant category for direct hits failed, returning..."); return kFALSE; } iterRichDirect = static_cast(catRichDirect->MakeIterator()); if (NULL == iterRichDirect) { Error("init", "Can not create catRichDirect iterator, returning..."); return kFALSE; } // Initialize output category for RICH HGeant tracks and its appropriate iterator catTrack = gHades->getCurrentEvent()->getCategory(catRichTrack); if (NULL == catTrack) { catTrack = pRichDet->buildCategory(catRichTrack); if (NULL == catTrack) { Error("init", "Can not build output category catRichTrack, returning..."); return kFALSE; } else gHades->getCurrentEvent()->addCategory(catRichTrack, catTrack, "Rich"); } iterRichTrack = static_cast(catTrack->MakeIterator()); if (NULL == iterRichTrack) { Error("init", "Can not create catRichTrack iterator, returning..."); return kFALSE; } // Initialize output category catRichCalSim and its appropriate iterator catCal = gHades->getCurrentEvent()->getCategory(catRichCal); if (NULL == catCal) { catCal = pRichDet->buildMatrixCat("HRichCalSim", 1); if (NULL == catCal) { Error("init", "Can not build output category catRichCalSim, returning..."); return kFALSE; } else gHades->getCurrentEvent()->addCategory(catRichCal, catCal, "Rich"); } iterRichCal = static_cast(catCal->MakeIterator()); if (NULL == iterRichCal) { Error("init", "Can not create catRichCalSim iterator, returning..."); return kFALSE; } ga = new TF1("ga", "[0]*exp([1]*(x-[2])*(x-[2]))", -10, 10); return kTRUE; } //============================================================================ //---------------------------------------------------------------------------- Bool_t HRichDigitizer::reinit() { if (NULL == pDigitisationPar || NULL == pGeometryPar) { Error("reinit", "One of the needed pointers is NULL, returning..."); return kFALSE; } // Get geometry parameters fWiresNr = pGeometryPar->getWiresPar()->getNrWires(); fWiresDist = pGeometryPar->getWiresPar()->getDistWire(); distWirePads = pGeometryPar->getDistanceWiresPads(); fYShift = pGeometryPar->getSectorShift(); fYShift = fYShift / cos(20. * TMath::DegToRad()); // Get digitization parameters fChargePerChannel = pDigitisationPar->getChargePerChannel(); fChargeScaling = pDigitisationPar->getChargeScaling(); fElectronsNr = pDigitisationPar->getElectronsNr(); fFactor1 = pDigitisationPar->getFactor1(); fFactor1Sig = pDigitisationPar->getFactor1Sig(); fFactor2 = pDigitisationPar->getFactor2(); fFactor2Sig = pDigitisationPar->getFactor2Sig(); fIncreaseNoise = pDigitisationPar->getIncreaseNoise(); fParam1 = pDigitisationPar->getParameter1(); fParam2 = pDigitisationPar->getParameter2(); fQupper = pDigitisationPar->getQupper(); fSigmaValue = pDigitisationPar->getSigmaValue(); binsQE = pDigitisationPar->getQEBinsNr(); photlength.Set(binsQE, pDigitisationPar->getPhotonLenArray()); photeffic .Set(binsQE, pDigitisationPar->getPhotonEffArray()); fExpSlope .Set(6, pDigitisationPar->getExpSlope()); for (Int_t sec = 0; sec < 6; ++sec) correction[sec].Set(binsQE, pDigitisationPar->getCorrectionParams(sec)); // this variable is the probability that a given pad will have // a signal induced by the electronic noise above threshold ga->SetRange(-10.0, 10.0); ga->SetParameters(1.0 / ((TMath::Sqrt(2. * TMath::Pi())) * fIncreaseNoise), -1.0 / (2. * fIncreaseNoise * fIncreaseNoise)); noiseProb = 0.5 - ga->Integral(0.0, fSigmaValue); return kTRUE; } //============================================================================ //---------------------------------------------------------------------------- Int_t HRichDigitizer::execute() { // At the end of the execute() the TClonesArray that contains // the track numbers is sorted by the pad addreses. Then // every CalSim object is retrieved from the catRichCal // container to assign to each pad the corresponding nTrack1, nTrack2. // These 2 numbers are the indexes of the track nbs. in the TClonesArray // HRichTrack. Finally the noise is calculated. vector index; mapTracks.clear(); HGeantRichPhoton* pCherenkovHit = NULL; HGeantRichDirect* pDirectHit = NULL; HRichCalSim* calsimobj = NULL; HRichPadFilter padFilter; Bool_t wasReal = kFALSE; Int_t fCerNr = 0; Int_t fDirNr = 0; countNoisePad = 0; fChargeOnWireList.Delete(); fTrackNumberList.Delete(); if (NULL == pCalPar) { Error("execute", "Pointer to HRichCalPar is NULL, needed for THR, noise,... simulation!!!"); return kSkipEvent; } //--------------------------------------------------------------------- // In embedding mode the trackNumber of the // real data has to be set // remember the address of pads fired by real data // for later vector < Int_t > addr_Real(200, 0); // alocate for 200 addr_Real.clear(); // actual size = 0 if (gHades->getEmbeddingMode() > 0) { HRichCalSim* cal = NULL; iterRichCal->Reset(); while ((cal = static_cast(iterRichCal->Next()))) { addr_Real.push_back(cal->getAddress()); // remember which objects came from real data cal->setEventNr(gHades->getCurrentEvent()->getHeader()->getEventSeqNumber()); cal->setEnergy(0); // realID Float_t dummy[2]; dummy[0] = gHades->getEmbeddingRealTrackId(); dummy[1] = 0; // flag: realID==0 TVector rTrack(1, 2, dummy); loc[0] = cal->getSector(); loc[1] = cal->getRow(); loc[2] = cal->getCol(); updateTrack(cal, loc, &rTrack); } } //--------------------------------------------------------------------- //------------------------------------------------------------------- // digitization of GEANT data starts here // if one wants to have pure simulation in embedding mode, HRichCalSim // category has to be cleared (filled with real data in unpacker task) // catCal->Clear(); iterRichPhoton->Reset(); while ((pCherenkovHit = static_cast(iterRichPhoton->Next()))) { ++fCerNr; // for the oem measurement I have to increase the track number // of each produced photon by myself (fCerNr). This number // will be the track number od the photon. if (kTRUE == isOEM) digitiseCherenkovHits(pCherenkovHit, fCerNr); else digitiseCherenkovHits(pCherenkovHit, 0); } iterRichDirect->Reset(); while ((pDirectHit = static_cast(iterRichDirect->Next()))) { digitiseDirectHits(pDirectHit); fDirNr++; } digitisePads(); //------------------------------------------------------------------- //------------------------------------------------------------------- // generate noise and add the charge to the already fired pads iterRichCal->Reset(); while ((calsimobj = static_cast(iterRichCal->Next()))) { //loop on cal sim objects. First the calsim container // contains only those pads fired by the photons. // if the electronic noise is switched a charge fluctuation // dues to the noise is added to the pad charge // and then all those pads that are below threshold // get a 0 charge wasReal = kFALSE; if (find(addr_Real.begin(), addr_Real.end(), calsimobj->getAddress()) != addr_Real.end()) { wasReal = kTRUE; } //do not add correlated noise to real data if (kTRUE == isActiveNoise && kFALSE == wasReal && NULL != pCalPar) { addNoiseToCharge(calsimobj); } if (NULL != pCalPar) { if (0 == checkPad(calsimobj) && kFALSE == wasReal) { // if charge is below threshold and // pad does not contain real data calsimobj->setCharge(0.);//LF } } // in the updateCharge member for each fired pad the energy of // the "hitting" photon is added. At the same time the number // of the photons hits is stored, so that after having digizited // all pads the average photon energy for each pad can be calculated // as follows. if (0 != calsimobj->getNTrack1()) calsimobj->setEnergy(calsimobj->getEnergy() / calsimobj->getNTrack1()); //------------------------------------------------------------------- // set index range } // end of loop on calsim object, pads fired by photons //------------------------------------------------------------------- //------------------------------------------------------------------- // after the charge propagation due tot he photons is completed // the electronic noise is produced. This part is skipped in embedding. if (kTRUE == isActiveNoise && gHades->getEmbeddingMode() == 0 && NULL != pCalPar) { makeNoiseOnPads(); } //------------------------------------------------------------------- //------------------------------------------------------------------- // last action: sort track category and fill // the index range to the cal object catTrack->sort(); // sort track objects iterRichCal->Reset(); while ((calsimobj = static_cast(iterRichCal->Next()))) { Int_t addPad = calsimobj->getAddress(); //calsimobj->setNTrack1(0); map >::iterator iter = mapTracks.find(addPad); // pad addr -> list of HRichTrack if (iter != mapTracks.end()) { // found address vector & list = iter->second; UInt_t n = list.size(); if (n > 1) { // more than 1 track hit the pad index.clear(); for (UInt_t i = 0; i < n; i ++) { // store the indices in HRichTrack cat for each HRichTrack index.push_back(catTrack->getIndex(list[i])); } std::sort(index.begin(), index.end()); // sort indices in ascending order calsimobj->setNTrack1(index[0]); calsimobj->setNTrack2(index[n - 1]); } else { // only 1 track hit the pad calsimobj->setNTrack1(catTrack->getIndex(list[0])); // both values are the same calsimobj->setNTrack2(calsimobj->getNTrack1()); } } else { Error("execute()", "address of pad =%i not found in map!",addPad); } } //------------------------------------------------------------------- //------------------------------------------------------------------- // all the pads with charge 0 are removed from the calsim container catCal->filter(padFilter); //------------------------------------------------------------------- return 0; } //============================================================================ //---------------------------------------------------------------------------- Float_t HRichDigitizer::GaussFun(const Float_t mean, const Float_t sigma) { Double_t x = 0.0; Double_t y = 0.0; Double_t z = 0.0; do { y = gRandom->Rndm(); } while (!y); z = gRandom->Rndm(); x = z * 2 * TMath::Pi(); return mean + sigma * static_cast(sin(x) * sqrt(-2.0 * log(y))); } //============================================================================ //---------------------------------------------------------------------------- Bool_t HRichDigitizer::calcQE(const Float_t photlen, const Int_t sec) { Int_t i = 0; Float_t fDraw = 0.0; if (100 == photlen) return kTRUE; fDraw = gRandom->Rndm(); for (i = 0; i < binsQE; ++i) { if (i == 0) { if (photlen < ((photlength[i] + photlength[i+1]) / 2)) { if (fSlopeCorrection * fDraw > photeffic[i] * correction[sec][i]) { return kFALSE; } else { return kTRUE; } } } else if (i == binsQE - 1) { if (photlen >= ((photlength[i-1] + photlength[i]) / 2)) { if (fSlopeCorrection * fDraw > photeffic[i] * correction[sec][i]) { return kFALSE; } else { return kTRUE; } } } else if (((photlength[i-1] + photlength[i]) / 2) <= photlen && ((photlength[i] + photlength[i+1]) / 2) > photlen) { if (fSlopeCorrection * fDraw > photeffic[i] * correction[sec][i]) { return kFALSE; } else { return kTRUE; } } } return kTRUE; } //============================================================================ //---------------------------------------------------------------------------- Int_t HRichDigitizer::getWireNr(const Float_t xhit) { Int_t i = 0; Int_t nWireNr = -1; Float_t Xwir = 0.0; for (i = 0; i < fWiresNr; ++i) { Xwir = pGeometryPar->getWiresPar()->getWire(i)->getWireX(); if ((xhit >= Xwir - fWiresDist) && (xhit < Xwir + fWiresDist)) { nWireNr = i; break; } } return nWireNr; } //============================================================================ //---------------------------------------------------------------------------- Float_t HRichDigitizer::calcChargeOnWire(const Int_t sector, const Float_t xhit, const Float_t yhit, const Float_t nTrack, const Float_t nFlag, Float_t ene) { // this function calculates the charge on each wire and saves it in a TList // moreover, it saves the track number and the corresponding flag // (cf HRichTrack) in a second TList. // The charge is calculated according to one exponential function // obtained fitting the total charge distribution of the photon candidates // belonging to class one ( OEM data). // That means that the so calculated charge isn't in reality the // original charge deposited on the wire but the total amount of the charge // coupled to the pads. Therefore it isn't any longer necessary to // use an additional coupling factor. ( 75%) Bool_t endLoop = kFALSE; Int_t fHitedWireNr = 0; Float_t fX = 0.0; Float_t fY = 0.0; Float_t fQ = 0.0; Float_t charge = 0.0; HRichFrame* pFrame = NULL; TVector* gTrack1 = NULL; Float_t dummy[2]; if (NULL == (pFrame = pGeometryPar->getFramePar())) { Error("calcChargeOnWire", "Can not get RICH frame parameters"); return -1; } if (-1 == (fHitedWireNr = getWireNr(xhit))) return -1; dummy[0] = nTrack; dummy[1] = nFlag; gTrack1 = new TVector(1, 2, dummy); fX = pGeometryPar->getWiresPar()->getWire(fHitedWireNr)->getWireX(); fY = yhit; while (kFALSE == endLoop) { charge = 1.0 / fExpSlope[sector] * log(gRandom->Rndm() * (exp(fExpSlope[sector] * fQupper) - 1.0) + 1.0); if (charge <= fQupper) endLoop = kTRUE; } fQ = charge * fChargePerChannel; if (ene == 100) { ene = 0.; } if (fQ > 0. && kFALSE == pFrame->isOut(fX, fY)) { fChargeOnWireList.Add(new HRichWireSignal(sector, fHitedWireNr, fQ, fX, fY, ene)); fTrackNumberList.Add(gTrack1); } else { delete gTrack1; } return fQ; } //============================================================================ //---------------------------------------------------------------------------- void HRichDigitizer::digitiseCherenkovHits(HGeantRichPhoton *pCerHit, const Int_t count) { // for every photon hit on a pad the resulting charge on a wire // is calculated and the track number of the photon parent is stored. // (cf. calcChargeOnWire). Float_t xHitcm = (pCerHit->getX() / 10.0); //HGeant output is in mm!!! Float_t yHitcm = (pCerHit->getY() / 10.0 + fYShift); //HGeant output is in mm!!! countFBphot = 0; if (0 == count) processPhoton(pCerHit->getEnergy(), xHitcm, yHitcm, pCerHit->getTrack(), pCerHit->getSector()); else processPhoton(pCerHit->getEnergy(), xHitcm, yHitcm, count, pCerHit->getSector()); } //============================================================================ //---------------------------------------------------------------------------- void HRichDigitizer::processPhoton(const Float_t ene, const Float_t xPos, const Float_t yPos, const Int_t track, const Int_t sector) { Float_t eneFeedBack = 0.0; Float_t xHitcmFB = 0.0; Float_t yHitcmFB = 0.0; Float_t chargeOnWire = 0.0; HRichFrame* pFrame = NULL; if (NULL == (pFrame = pGeometryPar->getFramePar())) { Error("processPhoton", "Can not get RICH frame parameters"); return; } if (kTRUE == calcQE(ene, sector) && kFALSE == pFrame->isOut(xPos, yPos)) { chargeOnWire = calcChargeOnWire(sector, xPos, yPos, track, 0, ene); if (calcFeedBack(sector, xPos, yPos, eneFeedBack, xHitcmFB, yHitcmFB, chargeOnWire) && countFBphot < 7) { // -5 is the tracknumber assigned to all the feedback photons. processPhoton(eneFeedBack, xHitcmFB, yHitcmFB, -5, sector); } } } //============================================================================ //---------------------------------------------------------------------------- Bool_t HRichDigitizer::calcFeedBack(const Int_t sec, const Float_t xhit, const Float_t yhit, Float_t &ene, Float_t &xhittFB, Float_t &yhittFB, const Float_t charge) { // We assume that the feed back photon is produced on the anodic wire, // the number of the feed back photon is proportional to the value A0. countFBphot++; Float_t thetaDir = 0.0; Float_t phiDir = 0.0; Float_t r = 0.0; Float_t fDraw1 = 0.0; Float_t fDraw2 = 0.0; if (gRandom->Rndm() < (charge * 7.8e-6)) { fDraw1 = gRandom->Rndm(); fDraw2 = gRandom->Rndm(); thetaDir = fDraw1 * 90; phiDir = fDraw1 * 360; r = distWirePads * TMath::Tan(thetaDir * TMath::DegToRad()); xhittFB = xhit + r * TMath::Cos(phiDir * TMath::DegToRad()); yhittFB = yhit + r * TMath::Sin(phiDir * TMath::DegToRad()); // the Alice TRD gives as VUV photon energies the following // values : 156nm,166nm,193nm with weights 30%,57%,13% // they correspond to :7.94 eV,7.46 eV,6.42eV. if (fDraw2 < 0.13) ene = 6.42; else if (fDraw2 < 0.43) ene = 7.94; else ene = 7.46; return kTRUE; } else { ene = 0; xhittFB = 0; yhittFB = 0; return kFALSE; } } //============================================================================ //---------------------------------------------------------------------------- void HRichDigitizer::digitiseDirectHits(HGeantRichDirect *pDirHit) { // for every direct hit (charge particle hitting the RICH or ionizing // the gas near the surface of the photon detector) the // resulting charge on the wires is calculated and the track number // of the charged particle hitting the pad is stored. // (cf. calcChargeOnWire) Int_t i = 0; Int_t fStepsNr = 0; Int_t dTrack = -1; Int_t dId = -1; Float_t fStepLength = 0.0; Float_t fNewX = 0.0; Float_t fNewY = 0.0; HRichFrame* pFrame = NULL; Float_t xHitcm = pDirHit->getX() / 10.0F; //HGeant output is in mm!!! Float_t yHitcm = pDirHit->getY() / 10.0F + fYShift ; //HGeant output is in mm!!! if (NULL == (pFrame = pGeometryPar->getFramePar())) { Error("digitiseDirectHits", "Can not get RICH frame parameters"); return; } if (pDirHit->getEnergyLoss() > 0.0F && // one day it will be additional condition if it is mirror or detector kFALSE == pFrame->isOut(xHitcm, yHitcm)) { // Number of electrons is given in keV while the energy loss is given // in MeV, therefore a factor 1000 is needed to calculate the correct // number of steps fStepsNr = static_cast(1 + fElectronsNr * 1000 * pDirHit->getEnergyLoss()); if (fStepsNr > 5000) { #ifdef HRICH_DEBUGMODE Warning("digitiseDirectHits", "Evt %d: number of maximum steps (5000) exceeded by %d!", gHades->getCurrentEvent()->getHeader()->getEventSeqNumber(), fStepsNr); #endif fStepsNr = 5000; } // units: trackLength is in mm and coordinates in cm! fStepLength = 0.1 * pDirHit->getTrackLengthInPhotDet() / static_cast(fStepsNr); for (i = 0; i < fStepsNr; ++i) { fNewX = xHitcm + (i + 1) * fStepLength * sin(TMath::DegToRad() * pDirHit->getTheta()) * cos(TMath::DegToRad() * pDirHit->getPhi()); fNewY = yHitcm + (i + 1) * fStepLength * sin(TMath::DegToRad() * pDirHit->getTheta()) * sin(TMath::DegToRad() * pDirHit->getPhi()); pDirHit->getTrack(dTrack, dId); // the energy of the direct hits is set to 0 // not to modify the average energy deposited by one or more // photons on one pad. (this is very important for the OEM // analysis). calcChargeOnWire(pDirHit->getSector(), fNewX, fNewY, dTrack, 1, 0); } } } //============================================================================ //---------------------------------------------------------------------------- HRichPad* HRichDigitizer::translateCorners(HRichPad *pPad, const Float_t dx, const Float_t dy) { // check what you want to translate - probably real corners Int_t i; HRichPad* fpTranslatedPad = new HRichPad(*pPad); if (NULL == fpTranslatedPad) { Error("translateCorners", "Can not get pointer to HRichPad"); return NULL; } for (i = 0; i < fpTranslatedPad->getCornersNr(); ++i) { fpTranslatedPad->getCorner(i)->setX(fpTranslatedPad->getCorner(i)->getX() - dx); fpTranslatedPad->getCorner(i)->setY(fpTranslatedPad->getCorner(i)->getY() - dy); } fpTranslatedPad->setXmin(fpTranslatedPad->getXmin() - dx); fpTranslatedPad->setYmin(fpTranslatedPad->getYmin() - dy); fpTranslatedPad->setXmax(fpTranslatedPad->getXmax() - dx); fpTranslatedPad->setYmax(fpTranslatedPad->getYmax() - dy); // the ownership of this pointer is passed to outer variable-pointer return fpTranslatedPad; } //=========================================================================== //---------------------------------------------------------------------------- void HRichDigitizer::makeNoiseOnPads() { // nbNoisePads is the number of pads on which the electronic noise // will produce a non zero signal. This number is extimated taken // into account the threshold taht has been used for the pedestals. HRichCalSim* pCalSim = NULL; HRichCalParCell* pCell = NULL; const Int_t maxCols = pGeometryPar->getColumns(); const Int_t maxPads = pGeometryPar->getPadsNr(); Int_t nbNoisePads = 0; Int_t padNr = 0; Int_t iCol = 0; Int_t iRow = 0; Float_t fCharge = 0.0; Float_t fSigmaPad = 0.0; Float_t fMeanPad = 0.0; Float_t fFloatMean = 0.0; // here the number of pads that will have a noise induced signal above // threshold is calculated. The number countNoisePad/6 represents // the number of pads that had already a noise signal above threshold // in each sector. nbNoisePads = static_cast(TMath::Floor(4800 * noiseProb)); nbNoisePads = gRandom->Poisson(nbNoisePads) - static_cast(countNoisePad / 6.); if (nbNoisePads < 0) { nbNoisePads = 0; } // a Gauss distributed signal is induced on nbNoisePads pads // for each sector. for (Int_t sec = 0; sec < 6; ++sec) { if (pGeometryPar->getSectorActive(sec) > 0) { for (Int_t j = 0; j < nbNoisePads; ++j) { do { padNr = static_cast(TMath::Floor(maxPads * gRandom->Rndm())); iRow = pGeometryPar->calcRow(padNr); iCol = pGeometryPar->calcCol(padNr); // the following condition is necessary to test that the pad // exists and is connected to the electronic if (padNr >= maxPads) { Error("makeNoiseOnPads", "Index %d out of bounds...", maxPads); continue; } } while (!pGeometryPar->getPadsPar()->getPad(padNr)->getPadActive() && iCol < maxCols); loc[0] = sec; loc[1] = iRow; loc[2] = iCol; pCell = static_cast(pCalPar->getObject(loc)); if (NULL != pCell) { fSigmaPad = pCell->getSigma(); } else { Error("makeNoiseOnPads", "Did not get any valid HRichCalParCell..."); fSigmaPad = 0.0; } if (fSigmaPad > 0.4) { fMeanPad = pCell->getOffset(); fFloatMean = fMeanPad - (TMath::Floor(fMeanPad)); fCharge = calcNoiseOnPad(fSigmaPad, fFloatMean); // the noise contribution is calculated with a gauss dist. // centred in 0, hence we have to add the float part of the // mean. Hence we cast the charge to an integer, then we // subtract the float part of the mean value and finally // we add a random number between 0 and 1 to have a float. pCalSim = static_cast(catCal->getObject(loc)); // if the pad has already been fired no additional // charge is added and a new pad is drawn. if (NULL == pCalSim && fCharge > 0) { pCalSim = static_cast(catCal->getSlot(loc)); if (NULL != pCalSim) { pCalSim = new(pCalSim) HRichCalSim; pCalSim->setEventNr(gHades->getCurrentEvent()->getHeader()->getEventSeqNumber()); pCalSim->setSector(loc[0]); pCalSim->setCol(loc[2]); pCalSim->setRow(loc[1]); pCalSim->setCharge(fCharge); // Now update HRichTrack Float_t dummy[2]; dummy[0] = -1; // noise pads have tracknumber == -1 dummy[1] = 0; // flag == 0 TVector rTrack(1, 2, dummy); updateTrack(pCalSim, loc, &rTrack); } } else j--; } else j--; } } } } //============================================================================ //---------------------------------------------------------------------------- Float_t HRichDigitizer::calcNoiseOnPad(const Float_t fSigmaPad, const Float_t fFloatMean) { // the charge induced on the pad by the electronic noise // is calculated according to a gauss function with // sigma = fSigmaPad*fIncreaseNoise and mean 0. // Anyway only values above threshold // (fSigmaValue * fSigmaPad) // are drawn const Int_t cutLoop = 10; const Int_t iThreshold = static_cast(fFloatMean + fSigmaValue * fSigmaPad); Int_t count = 0; Int_t chan = 0; Float_t fCharge = 0.0; if (fSigmaValue < 2.0) return static_cast(TMath::Sqrt(2.) * fSigmaPad * fIncreaseNoise * noiseCharge[0] + fFloatMean) - fFloatMean + gRandom->Rndm(); do { // the InverseErf function is defined in the array noiseCharge // with 1000 bin in the range (0.95,1). This range should be enough // to calculate the noise contribution up to a 2 sigma cut if the // variable fIncreaseNoise is 1. If the noise is increased ( e.g. // fSigmaNoise=1.9) the values are consistent with a cut of 3.9 sigma. chan = static_cast(((1. - 2. * noiseProb * gRandom->Rndm()) - 0.95) / (0.05 / 1000.)); if (chan >= 0 && chan < 1000) { fCharge = TMath::Sqrt(2.) * fSigmaPad * fIncreaseNoise * noiseCharge[chan] + fFloatMean; } else fCharge = 0.; count++; } while ((static_cast(fCharge) <= iThreshold || fCharge > 1000.) && count < cutLoop); // the value of sigma can be either too big or too small // we try to avoid extreme cases. // If the sigma of the pad is too big (>10) or too small(<0.4) // then the array used to calculate the value of the charge // has not the required granularity and the previous loop // is neverending. Therefore we allow only 10 tries, that // should be sufficient for a reasonable sigma. if ((static_cast(fCharge) <= iThreshold || fCharge > 1000.) && count == cutLoop) return 0.; return static_cast(fCharge) - fFloatMean + gRandom->Rndm(); } //=========================================================================== //---------------------------------------------------------------------------- void HRichDigitizer::addNoiseToCharge(HRichCalSim* calSim) { // An additional charge is added on pads that have already // been "fired" by a photon. // This charge is drawn according to a Gauss distribution // with sigma = fSigmaPad*fIncreaseNoise and mean = 0. HRichCalParCell* pCell = NULL; Float_t fCharge = 0.0; Float_t fSigmaPad = 0.0; loc[0] = calSim->getSector(); loc[1] = calSim->getRow(); loc[2] = calSim->getCol(); pCell = static_cast(pCalPar->getObject(loc)); if (NULL != pCell) { fSigmaPad = pCell->getSigma(); } else { Error("addNoiseToCharge", "Did not get any valid HRichCalParCell..."); fSigmaPad = 1.0; } fCharge = GaussFun(0, fSigmaPad * fIncreaseNoise); if (fCharge > fSigmaPad * fSigmaValue) { countNoisePad++; } calSim->addCharge(fCharge); } //============================================================================ //---------------------------------------------------------------------------- Int_t HRichDigitizer::checkPad(HRichCalSim *calSim) { // this function checks if the pads fired by the photons // are above threshold. // The procedure used to calculate the threshols is the // same used for real data. The charge on the pad is composed // of the signal induced by the photon plus the noise fluctuation. // The float part of the mean value of the noise distribution // is added to the charge then // the charge value is casted to an integer and is // compared to the threshold. // If the charge is above threshold, the float part of the mean is // subtracted from the total charge and a random number between 0 and 1 // is added to the integer produce a float. HRichCalParCell* pCell = NULL; Float_t fCharge = calSim->getCharge(); Int_t iRow = calSim->getRow(); Int_t iCol = calSim->getCol(); Int_t iSec = calSim->getSector(); Int_t padNr = iCol + iRow * pGeometryPar->getColumns(); Int_t iCharge = 0; Int_t iThreshold = 0; Float_t fSigmaPad = 0.0; Float_t fMeanPad = 0.0; Float_t fFloatMean = 0.0; loc[0] = iSec; loc[1] = iRow; loc[2] = iCol; pCell = static_cast(pCalPar->getObject(loc)); if (NULL != pCell) { fSigmaPad = pCell->getSigma(); fMeanPad = pCell->getOffset(); } else { Error("checkPad", "Did not get any valid HRichCalParCell..."); Error("checkPad", "Is pad active??? %d", pGeometryPar->getPadsPar()->getPad(padNr)->getPadActive()); fSigmaPad = 1.0; } // check, if current pad has valid calibration parameters if (0 == fSigmaPad && 0 == fMeanPad) { return 0; } fFloatMean = fMeanPad - static_cast(fMeanPad); fCharge += fFloatMean; iCharge = static_cast(fCharge); iThreshold = static_cast(fFloatMean + fSigmaValue * fSigmaPad); if (iCharge > iThreshold) { fCharge = iCharge - fFloatMean + gRandom->Rndm(); calSim->setCharge(fCharge); return iCharge; } else return 0; } //============================================================================ //---------------------------------------------------------------------------- void HRichDigitizer::digitisePads() { // for each wire on which some charge has been deposited // the corresponding coupled pads are calculated. // The function updateCharge creates a HRichCalSim obj if // the charge on the pad is greater than zero. // A cut off threshold is applied to each pad in the execute func. // after the digitilasation of all the hits. // If the pad is hit twice in the same event, the charges // corresponding to the 2 different hits are added. // The particle track number is passed to the function // updateCharge, too. UInt_t xcenter = 0; UInt_t ycenter = 0; Int_t sector = -1; Int_t iPadInd = 0; Float_t energyPhot = 0.0; Float_t yRel = 0.0; Float_t charge4 = 0.0; Float_t param11 = 0.0; Float_t param21 = 0.0; HRichPadTab* pPadsPar = NULL; HRichPad* pPad = NULL; HRichPad* pTmpPad = NULL; TVector* nTrackTemp = NULL; HRichWireSignal* pCharge = NULL; if (NULL == (pPadsPar = pGeometryPar->getPadsPar())) { Error("digitisePads", "Can not get RICH pad parameters"); return; } for (Int_t i = 0; i < fChargeOnWireList.GetSize(); ++i) { for (Int_t j = 0; j < 9; ++j) fDigitPadMatrix[j] = 0.; pCharge = static_cast(fChargeOnWireList.At(i)); nTrackTemp = static_cast(fTrackNumberList.At(i)); energyPhot = pCharge->getEnergy(); sector = pCharge->getSector(); pPad = pPadsPar->getPad(pCharge->getX(), pCharge->getY()); xcenter = pPad->getPadX(); ycenter = pPad->getPadY(); // the position of the photon hit on the central pad is calculated. yRel = TMath::Abs((pCharge->getY() - pPad->getYmin()) / (pPad->getYmax() - pPad->getYmin())); pTmpPad = translateCorners(pPad, pCharge->getX(), pCharge->getY()); // do not forget to delete it at the end of loop // ------------------------ central pad (0,0) charge4 = 0.; if (2 == pTmpPad->getPadFlag()) { // First we calculate the charge on the central pad taking as // contributions on the side pads, 0 for the weak coupling and t // he mean for the strong one. // Then we recalculate the coupling parameter for the // strong coupling drawing a number on the landau // distribution centered in the mean value. param11 = gRandom->Landau(fFactor1, fFactor1Sig); param21 = GaussFun(fFactor2, fFactor2Sig); if (param21 < 0.) param21 = 0.; charge4 = q4Calc(pCharge->getCharge(), yRel, param11, param21); if (charge4 < 0.) fDigitPadMatrix[4] = 0.; else { updateCharge(sector, pPad, charge4, nTrackTemp, energyPhot); fDigitPadMatrix[4] = charge4; } } //------------------------ loop on other pads iPadInd = 0; for (Int_t n = -1; n < 2; ++n) { for (Int_t m = -1; m < 2; ++m) { if ((m != 0 || n != 0) && (static_cast(xcenter) + m) >= 0 && (static_cast(ycenter) + n) >= 0) { if (pPadsPar->getPad(xcenter + m, ycenter + n) && pPadsPar->getPad(xcenter + m, ycenter + n)->getPadActive()) { delete pTmpPad; pTmpPad = translateCorners(pPadsPar->getPad(xcenter + m, ycenter + n), pCharge->getX(), pCharge->getY()); if (2 == pTmpPad->getPadFlag()) { fDigitPadMatrix[iPadInd] = calcIndCharge(yRel, fDigitPadMatrix[4], iPadInd, pCharge->getWireNr(), param11, param21); if (fDigitPadMatrix[iPadInd] <= 0.) fDigitPadMatrix[iPadInd] = 0.; else updateCharge(sector, pTmpPad, fDigitPadMatrix[iPadInd], nTrackTemp, energyPhot); } } } iPadInd ++; } } delete pTmpPad; } } //============================================================================ //---------------------------------------------------------------------------- Float_t HRichDigitizer::calcIndCharge(const Float_t yCharge, const Float_t q4, const Int_t iPdaIndex, const Int_t iWireNr, const Float_t param11, const Float_t param21) { // iWireNr is used to determine whether the wire is // the left or the right one. // Even wire numbers correspond to right. // for the side pads a coupling costant is multiplied // by the charge on the central pad, the upper an lower // pad get an amount of charge that is position dependent. switch (iPdaIndex) { case 1: return (qX(1 - yCharge) * q4); case 7: return (qX(yCharge) * q4); case 3: { if (iWireNr % 2 == 0) return (param21 * q4); else return (param11 * q4); } case 5: { if (iWireNr % 2 == 0) return (param11 * q4); else return (param21 * q4); } case 2: { if (iWireNr % 2 == 0) return (1.37 * param11 * q4 * qX(1 - yCharge)); else return (1.37 * param21 * q4 * qX(1 - yCharge)); } case 8: { if (iWireNr % 2 == 0) return (1.37 * param11 * q4 * qX(yCharge)); else return (1.37 * param21 * q4 * qX(yCharge)); } case 6: { if (iWireNr % 2 == 0) return (1.37 * param21 * q4 * qX(yCharge)); else return (1.37 * param11 * q4 * qX(yCharge)); } case 0: { if (iWireNr % 2 == 0) return (1.37 * param21 * q4 * qX(1 - yCharge)); else return (1.37 * param11 * q4 * qX(1 - yCharge)); } } return 0; } //============================================================================ //---------------------------------------------------------------------------- Float_t HRichDigitizer::qX(const Float_t pos) { // This function calculate the charge ratio between the upper/lower // pad and the central one Float_t charge = 0.0; Float_t q0 = -0.03; Float_t integral = 0.0; integral = 1.0 / (fParam1 * fParam1 * q0 + fParam1 * fParam2) - 1.0 / (fParam1 * fParam1 + fParam1 * fParam2); charge = ((q0 + fParam2 / fParam1) / (1.0 - (q0 + fParam2 / fParam1) * fParam1 * fParam1 * integral * pos)) - fParam2 / fParam1; if (charge < 0.) charge = 0.; return charge; } //============================================================================ //---------------------------------------------------------------------------- Float_t HRichDigitizer::q4Calc(const Float_t charge, const Float_t pos, const Float_t par1, const Float_t par2) { // This function calculate the charge induced on the central pad. // the coupling factor (75%) from the wire to the pad has // been discarded since the charge distribution on the wire // has been calculated fitting pads responses. Float_t gX = 1 + qX(pos) + qX(1 - pos); Float_t gX1 = 1 + 1.37 * qX(pos) + 1.37 * qX(1 - pos); if (0 != gX) return charge / ((par1 + par2) * gX1 + gX); else return 0.0; } //============================================================================ //---------------------------------------------------------------------------- void HRichDigitizer::updateCharge(const Int_t sector, HRichPad* pPad, const Float_t charge, TVector* rTrack, const Float_t ene) { // This function creates an HRichCalSim obj that corresponds // to a fired pad, it calls the function updateTrack, // that stores the corresponding track numbers. Int_t fCol = 0; Int_t fRow = 0; HRichCalSim* pCalSim = NULL; pPad->getPadXY(&fCol, &fRow); locat[0] = sector; locat[1] = fRow; locat[2] = fCol; pCalSim = static_cast(catCal->getObject(locat)); if (NULL == pCalSim) { pCalSim = static_cast(catCal->getSlot(locat)); if (NULL != pCalSim) { pCalSim = new(pCalSim) HRichCalSim; pCalSim->setEventNr(gHades->getCurrentEvent()->getHeader()->getEventSeqNumber()); pCalSim->setSector(locat[0]); pCalSim->setCol(locat[2]); pCalSim->setRow(locat[1]); updateTrack(pCalSim, locat, rTrack); } } if (NULL != pCalSim) { pCalSim->addCharge(charge * fChargeScaling / fChargePerChannel); // the feed back photon and the direct hits are assigned energy // equal to 0. Infact the energy of the other photons emitted in the // oem experiment depends on the polar emission angle theta // (dispersion in the solid radiators) and this is not true // for the feed back photons and for the direct hits. // There fore their energy isnt added to the pad they fire // and isnt taken into account in the calculation theta->Energy // because it would falsify the distribution. // cut-off in ADC if (pCalSim->getCharge() > 800.0) pCalSim->setCharge(800.0); if (ene > 0) { pCalSim->addHit(); pCalSim->addEnergy(ene); } // the track object is created for all the photon and direct hits // the feedback photons have track number -5 and flag 0. // Only the pad fired just because of electronic noise // dont correspond to any track objects. updateTrack(pCalSim, locat, rTrack); } } //============================================================================ //--------------------------------------------------------------------------- void HRichDigitizer::updateTrack(HRichCalSim* pCal, HLocation& loc, TVector* rTrack) { // this function stores the track numbers of parent // particles of photons and of direct hits in a linear // category (HRichTrack). This category is set sortable. HRichTrack *pRichTrack = NULL; Int_t Ad = pCal->getAddress(); map >::iterator iter = mapTracks.find(Ad); if (iter != mapTracks.end()) { // found address vector & list = iter->second; UInt_t n = list.size(); for (UInt_t i = 0; i < n; i++) { pRichTrack = list[i]; if (pRichTrack->getTrack() == (*rTrack)(1) && pRichTrack->getFlag() == (*rTrack)(2) ) return; // already registered } } pRichTrack = static_cast(catTrack->getNewSlot(loc1)); if (NULL != pRichTrack) { pRichTrack = new(pRichTrack) HRichTrack; pRichTrack->setEventNr(pCal->getEventNr()); //usage of TVector because of used in TList, explicit cast to integer is needed since TVector converts all arguments into double. pRichTrack->setTrack((Int_t)((*rTrack)(1))); pRichTrack->setFlag((Int_t)((*rTrack)(2))); pRichTrack->setAddress(Ad); if (iter != mapTracks.end()) { // old address vector & list = iter->second; list.push_back(pRichTrack); } else { // new address vector list; list.push_back(pRichTrack); mapTracks[Ad] = list; } } } //============================================================================ //--------------------------------------------------------------------------- Bool_t HRichDigitizer::finalize() { return kTRUE; } //============================================================================ //--------------------------------------------------------------------------- // This table contains values of the inverse error function in the range // [0.95;1.] needed to create a set of random numbers following this // funcion, used in makeNoiseOnPads const Float_t HRichDigitizer::noiseCharge[] = {1.3859, 1.38621, 1.38651, 1.38681, 1.38712, 1.38742, 1.38772, 1.38803, 1.38833, 1.38864, 1.38894, 1.38925, 1.38955, 1.38986, 1.39016, 1.39047, 1.39078, 1.39108, 1.39139, 1.3917, 1.392, 1.39231, 1.39262, 1.39293, 1.39324, 1.39355, 1.39386, 1.39416, 1.39447, 1.39478, 1.39509, 1.3954, 1.39572, 1.39603, 1.39634, 1.39665, 1.39696, 1.39727, 1.39759, 1.3979, 1.39821, 1.39852, 1.39884, 1.39915, 1.39946, 1.39978, 1.40009, 1.40041, 1.40072, 1.40104, 1.40135, 1.40167, 1.40199, 1.4023, 1.40262, 1.40294, 1.40325, 1.40357, 1.40389, 1.40421, 1.40453, 1.40485, 1.40516, 1.40548, 1.4058, 1.40612, 1.40644, 1.40676, 1.40708, 1.40741, 1.40773, 1.40805, 1.40837, 1.40869, 1.40901, 1.40934, 1.40966, 1.40998, 1.41031, 1.41063, 1.41096, 1.41128, 1.41161, 1.41193, 1.41226, 1.41258, 1.41291, 1.41323, 1.41356, 1.41389, 1.41422, 1.41454, 1.41487, 1.4152, 1.41553, 1.41586, 1.41619, 1.41651, 1.41684, 1.41717, 1.4175, 1.41784, 1.41817, 1.4185, 1.41883, 1.41916, 1.41949, 1.41983, 1.42016, 1.42049, 1.42083, 1.42116, 1.42149, 1.42183, 1.42216, 1.4225, 1.42283, 1.42317, 1.4235, 1.42384, 1.42418, 1.42451, 1.42485, 1.42519, 1.42553, 1.42587, 1.4262, 1.42654, 1.42688, 1.42722, 1.42756, 1.4279, 1.42824, 1.42858, 1.42892, 1.42927, 1.42961, 1.42995, 1.43029, 1.43064, 1.43098, 1.43132, 1.43167, 1.43201, 1.43236, 1.4327, 1.43305, 1.43339, 1.43374, 1.43408, 1.43443, 1.43478, 1.43512, 1.43547, 1.43582, 1.43617, 1.43652, 1.43687, 1.43722, 1.43757, 1.43792, 1.43827, 1.43862, 1.43897, 1.43932, 1.43967, 1.44002, 1.44038, 1.44073, 1.44108, 1.44144, 1.44179, 1.44215, 1.4425, 1.44286, 1.44321, 1.44357, 1.44392, 1.44428, 1.44464, 1.44499, 1.44535, 1.44571, 1.44607, 1.44643, 1.44679, 1.44715, 1.44751, 1.44787, 1.44823, 1.44859, 1.44895, 1.44931, 1.44967, 1.45004, 1.4504, 1.45076, 1.45113, 1.45149, 1.45185, 1.45222, 1.45259, 1.45295, 1.45332, 1.45368, 1.45405, 1.45442, 1.45479, 1.45515, 1.45552, 1.45589, 1.45626, 1.45663, 1.457, 1.45737, 1.45774, 1.45811, 1.45848, 1.45886, 1.45923, 1.4596, 1.45997, 1.46035, 1.46072, 1.4611, 1.46147, 1.46185, 1.46222, 1.4626, 1.46297, 1.46335, 1.46373, 1.46411, 1.46448, 1.46486, 1.46524, 1.46562, 1.466, 1.46638, 1.46676, 1.46714, 1.46753, 1.46791, 1.46829, 1.46867, 1.46906, 1.46944, 1.46982, 1.47021, 1.47059, 1.47098, 1.47136, 1.47175, 1.47214, 1.47253, 1.47291, 1.4733, 1.47369, 1.47408, 1.47447, 1.47486, 1.47525, 1.47564, 1.47603, 1.47642, 1.47681, 1.47721, 1.4776, 1.47799, 1.47839, 1.47878, 1.47918, 1.47957, 1.47997, 1.48036, 1.48076, 1.48116, 1.48156, 1.48195, 1.48235, 1.48275, 1.48315, 1.48355, 1.48395, 1.48435, 1.48475, 1.48516, 1.48556, 1.48596, 1.48636, 1.48677, 1.48717, 1.48758, 1.48798, 1.48839, 1.4888, 1.4892, 1.48961, 1.49002, 1.49043, 1.49083, 1.49124, 1.49165, 1.49206, 1.49247, 1.49289, 1.4933, 1.49371, 1.49412, 1.49454, 1.49495, 1.49536, 1.49578, 1.49619, 1.49661, 1.49703, 1.49744, 1.49786, 1.49828, 1.4987, 1.49912, 1.49954, 1.49996, 1.50038, 1.5008, 1.50122, 1.50164, 1.50207, 1.50249, 1.50291, 1.50334, 1.50376, 1.50419, 1.50461, 1.50504, 1.50547, 1.50589, 1.50632, 1.50675, 1.50718, 1.50761, 1.50804, 1.50847, 1.5089, 1.50934, 1.50977, 1.5102, 1.51064, 1.51107, 1.5115, 1.51194, 1.51238, 1.51281, 1.51325, 1.51369, 1.51413, 1.51457, 1.51501, 1.51545, 1.51589, 1.51633, 1.51677, 1.51721, 1.51765, 1.5181, 1.51854, 1.51899, 1.51943, 1.51988, 1.52033, 1.52077, 1.52122, 1.52167, 1.52212, 1.52257, 1.52302, 1.52347, 1.52392, 1.52437, 1.52483, 1.52528, 1.52573, 1.52619, 1.52665, 1.5271, 1.52756, 1.52802, 1.52847, 1.52893, 1.52939, 1.52985, 1.53031, 1.53077, 1.53123, 1.5317, 1.53216, 1.53262, 1.53309, 1.53355, 1.53402, 1.53449, 1.53495, 1.53542, 1.53589, 1.53636, 1.53683, 1.5373, 1.53777, 1.53824, 1.53871, 1.53919, 1.53966, 1.54014, 1.54061, 1.54109, 1.54156, 1.54204, 1.54252, 1.543, 1.54348, 1.54396, 1.54444, 1.54492, 1.5454, 1.54589, 1.54637, 1.54685, 1.54734, 1.54783, 1.54831, 1.5488, 1.54929, 1.54978, 1.55027, 1.55076, 1.55125, 1.55174, 1.55223, 1.55273, 1.55322, 1.55372, 1.55421, 1.55471, 1.55521, 1.5557, 1.5562, 1.5567, 1.5572, 1.5577, 1.55821, 1.55871, 1.55921, 1.55972, 1.56022, 1.56073, 1.56123, 1.56174, 1.56225, 1.56276, 1.56327, 1.56378, 1.56429, 1.5648, 1.56532, 1.56583, 1.56635, 1.56686, 1.56738, 1.56789, 1.56841, 1.56893, 1.56945, 1.56997, 1.57049, 1.57102, 1.57154, 1.57206, 1.57259, 1.57312, 1.57364, 1.57417, 1.5747, 1.57523, 1.57576, 1.57629, 1.57682, 1.57735, 1.57789, 1.57842, 1.57896, 1.57949, 1.58003, 1.58057, 1.58111, 1.58165, 1.58219, 1.58273, 1.58328, 1.58382, 1.58437, 1.58491, 1.58546, 1.58601, 1.58655, 1.5871, 1.58765, 1.58821, 1.58876, 1.58931, 1.58987, 1.59042, 1.59098, 1.59154, 1.59209, 1.59265, 1.59321, 1.59378, 1.59434, 1.5949, 1.59547, 1.59603, 1.5966, 1.59717, 1.59773, 1.5983, 1.59887, 1.59945, 1.60002, 1.60059, 1.60117, 1.60174, 1.60232, 1.6029, 1.60348, 1.60406, 1.60464, 1.60522, 1.6058, 1.60639, 1.60697, 1.60756, 1.60815, 1.60874, 1.60933, 1.60992, 1.61051, 1.6111, 1.6117, 1.6123, 1.61289, 1.61349, 1.61409, 1.61469, 1.61529, 1.61589, 1.6165, 1.6171, 1.61771, 1.61832, 1.61892, 1.61953, 1.62015, 1.62076, 1.62137, 1.62199, 1.6226, 1.62322, 1.62384, 1.62446, 1.62508, 1.6257, 1.62632, 1.62695, 1.62757, 1.6282, 1.62883, 1.62946, 1.63009, 1.63072, 1.63136, 1.63199, 1.63263, 1.63327, 1.6339, 1.63454, 1.63519, 1.63583, 1.63647, 1.63712, 1.63777, 1.63841, 1.63906, 1.63972, 1.64037, 1.64102, 1.64168, 1.64233, 1.64299, 1.64365, 1.64431, 1.64498, 1.64564, 1.64631, 1.64697, 1.64764, 1.64831, 1.64898, 1.64966, 1.65033, 1.65101, 1.65168, 1.65236, 1.65304, 1.65372, 1.65441, 1.65509, 1.65578, 1.65647, 1.65716, 1.65785, 1.65854, 1.65924, 1.65993, 1.66063, 1.66133, 1.66203, 1.66273, 1.66344, 1.66414, 1.66485, 1.66556, 1.66627, 1.66698, 1.6677, 1.66841, 1.66913, 1.66985, 1.67057, 1.67129, 1.67202, 1.67274, 1.67347, 1.6742, 1.67493, 1.67567, 1.6764, 1.67714, 1.67788, 1.67862, 1.67936, 1.68011, 1.68085, 1.6816, 1.68235, 1.6831, 1.68386, 1.68461, 1.68537, 1.68613, 1.68689, 1.68766, 1.68842, 1.68919, 1.68996, 1.69073, 1.69151, 1.69228, 1.69306, 1.69384, 1.69462, 1.6954, 1.69619, 1.69698, 1.69777, 1.69856, 1.69936, 1.70015, 1.70095, 1.70175, 1.70256, 1.70336, 1.70417, 1.70498, 1.70579, 1.7066, 1.70742, 1.70824, 1.70906, 1.70988, 1.71071, 1.71154, 1.71237, 1.7132, 1.71404, 1.71487, 1.71571, 1.71656, 1.7174, 1.71825, 1.7191, 1.71995, 1.72081, 1.72166, 1.72252, 1.72339, 1.72425, 1.72512, 1.72599, 1.72686, 1.72774, 1.72862, 1.7295, 1.73038, 1.73127, 1.73216, 1.73305, 1.73394, 1.73484, 1.73574, 1.73664, 1.73755, 1.73846, 1.73937, 1.74028, 1.7412, 1.74212, 1.74304, 1.74397, 1.7449, 1.74583, 1.74677, 1.7477, 1.74865, 1.74959, 1.75054, 1.75149, 1.75244, 1.7534, 1.75436, 1.75532, 1.75629, 1.75726, 1.75823, 1.75921, 1.76019, 1.76117, 1.76216, 1.76315, 1.76415, 1.76514, 1.76614, 1.76715, 1.76816, 1.76917, 1.77018, 1.7712, 1.77223, 1.77325, 1.77428, 1.77532, 1.77635, 1.7774, 1.77844, 1.77949, 1.78054, 1.7816, 1.78266, 1.78373, 1.7848, 1.78587, 1.78695, 1.78803, 1.78912, 1.79021, 1.7913, 1.7924, 1.7935, 1.79461, 1.79572, 1.79684, 1.79796, 1.79909, 1.80022, 1.80135, 1.80249, 1.80363, 1.80478, 1.80594, 1.80709, 1.80826, 1.80942, 1.8106, 1.81178, 1.81296, 1.81415, 1.81534, 1.81654, 1.81774, 1.81895, 1.82017, 1.82139, 1.82261, 1.82384, 1.82508, 1.82632, 1.82757, 1.82882, 1.83008, 1.83135, 1.83262, 1.83389, 1.83518, 1.83646, 1.83776, 1.83906, 1.84037, 1.84168, 1.843, 1.84433, 1.84566, 1.847, 1.84835, 1.8497, 1.85106, 1.85243, 1.8538, 1.85518, 1.85657, 1.85796, 1.85937, 1.86078, 1.86219, 1.86362, 1.86505, 1.86649, 1.86794, 1.86939, 1.87086, 1.87233, 1.87381, 1.8753, 1.87679, 1.8783, 1.87981, 1.88133, 1.88286, 1.8844, 1.88595, 1.88751, 1.88908, 1.89065, 1.89224, 1.89383, 1.89544, 1.89705, 1.89868, 1.90031, 1.90196, 1.90361, 1.90528, 1.90696, 1.90864, 1.91034, 1.91205, 1.91377, 1.9155, 1.91725, 1.919, 1.92077, 1.92255, 1.92434, 1.92615, 1.92796, 1.92979, 1.93163, 1.93349, 1.93536, 1.93724, 1.93914, 1.94105, 1.94297, 1.94491, 1.94687, 1.94884, 1.95082, 1.95282, 1.95484, 1.95687, 1.95892, 1.96098, 1.96306, 1.96516, 1.96728, 1.96941, 1.97156, 1.97373, 1.97592, 1.97813, 1.98036, 1.98261, 1.98487, 1.98716, 1.98947, 1.9918, 1.99415, 1.99653, 1.99892, 2.00135, 2.00379, 2.00626, 2.00875, 2.01127, 2.01381, 2.01638, 2.01898, 2.02161, 2.02426, 2.02694, 2.02965, 2.0324, 2.03517, 2.03797, 2.04081, 2.04368, 2.04658, 2.04952, 2.0525, 2.05551, 2.05856, 2.06164, 2.06477, 2.06794, 2.07115, 2.0744, 2.0777, 2.08105, 2.08444, 2.08788, 2.09137, 2.09491, 2.09851, 2.10216, 2.10587, 2.10963, 2.11346, 2.11735, 2.1213, 2.12533, 2.12942, 2.13358, 2.13783, 2.14214, 2.14654, 2.15103, 2.1556, 2.16027, 2.16503, 2.16989, 2.17486, 2.17993, 2.18512, 2.19044, 2.19587, 2.20144, 2.20716, 2.21301, 2.21903, 2.22521, 2.23156, 2.2381, 2.24484, 2.25179, 2.25896, 2.26637, 2.27404, 2.28199, 2.29023, 2.2988, 2.30773, 2.31703, 2.32675, 2.33694, 2.34763, 2.35889, 2.37078, 2.38339, 2.39679, 2.41112, 2.42652, 2.44316, 2.46127, 2.48115, 2.50322, 2.52803, 2.5564, 2.58961, 2.62974, 2.68069, 2.75106, 2.86776 };