//*-- Author : 08/06/2006 D. Gonzalez-Diaz //*-- Modified : 12/06/2006 D. Gonzalez-Diaz //*-- Modified : 13/12/2009 D. Gonzalez-Diaz //*-- Modified : 24/11/2010 P. Cabanelas DigiPar container moved to //*-- Modified : condition style. //*-- Modified : 22/06/2012 A.Mangiarotti General code cleanup. //*-- Modified : Geometry parameters used in several places //*-- Modified : are now stored as local variables (see D, //*-- Modified : geaXLoc). Reading of the digitiser //*-- Modified : parameters, which are constant, moved out //*-- Modified : of the hit processing loop (see vprop, //*-- Modified : sigma_el, sigma_T, t_offset, Qmean, gap)! //*-- Modified : The inefficiency per gap for normally //*-- Modified : incident particles ineff_gap_n is //*-- Modified : calculated once for all outside the hit //*-- Modified : loop. Namespace TMath used consistently //*-- Modified : everywhere. The procedure has been //*-- Modified : simplified. First the ideal left and //*-- Modified : right times left->gtime and right->gtime //*-- Modified : are calculated and after the resolutions //*-- Modified : are added to obtain left->time and //*-- Modified : right->time. //*-- Modified : 20/09/2012 A.Mangiarotti Special treatment of //*-- Modified : particles produced in the gap. //*-- Modified : 20/11/2012 A.Mangiarotti modified to include the new //*-- Modified : model II with parametrised dependences on //*-- Modified : the beta of the particles. Efficiency and //*-- Modified : time of flight resolution employ a //*-- Modified : sigmoid function. The charge distribution //*-- Modified : is now a Landau with parameters depending //*-- Modified : on beta with second degree polynomials. //*-- Modified : An attempt to correct on average for the //*-- Modified : number of active gaps is included. The //*-- Modified : old digitizer is now available as model //*-- Modified : I. If the old parameter container is //*-- Modified : found, model I is used automatically. //*-- Modified : 27/11/2012 A.Mangiarotti added new cell wise mode of //*-- Modified : operation controllable with a flag //*-- Modified : (fMode=1). If the new mode is selected, //*-- Modified : the gaps are first reorganised into //*-- Modified : cells before applying the response with //*-- Modified : model I or model II. Two new internal //*-- Modified : structs celltrack and celldat added to //*-- Modified : handle the reorganisation of the gaps //*-- Modified : into cells. If model II is used in cell //*-- Modified : wise mode, no further corrections are //*-- Modified : needed for the average number of active //*-- Modified : gaps. //*-- Modified : 27/02/2012 A.Mangiarotti A small bug has been //*-- Modified : corrected: the information about the //*-- Modified : local track length normalised to the gap //*-- Modified : was not propagated from gaps to cells. //*-- Modified : The access to the HGeantRpc is done with //*-- Modified : the new optimised getHitDigi to avoid //*-- Modified : retrieving information that is not used. //*-- Modified : 03/03/2013 A.Mangiarotti Two new (private) methods //*-- Modified : have been introduced to streamline the //*-- Modified : digitiser operation and prepare for the //*-- Modified : different options available with the new //*-- Modified : HGeantRpc structure. //*-- Modified : 10/03/2013 A.Mangiarotti Digitizer adapted to the //*-- Modified : new HGeantRpc structure: all the old //*-- Modified : options remain available for the old //*-- Modified : HGeantRpc to have full backward //*-- Modified : compatibility. With the new HGeantRpc //*-- Modified : only the cell wise mode is available. //*-- Modified : 25/03/2013 J.Markert Code clean up. The particle //*-- Modified : mass is now consistently retrieved from //*-- Modified : HPhysicsConstants. //*-- Modified : 09/07/2013 A.Mangiarotti counters nHGeantRpc and //*-- Modified : nDigitizedRpc introduced to flag the //*-- Modified : events with no generated hits and no //*-- Modified : digitised hits, respectively. This //*-- Modified : allows to avoid checking all //*-- Modified : cellobjects when nHGeantRpc=0 and all //*-- Modified : rpcobjects when nDigitizedRpc=0. //*-- Modified : Warning: when nHGeantRpc=0, the //*-- Modified : information on the HGeantRpc version //*-- Modified : in HGeantRpc_version is not available. //*-- Modified : 04/12/2013 G.Kornakov Added the embedding mode ///////////////////////////////////////////////////////////////////// // HRpcDigitizer digitizes HGeantRpc data, puts output values into // HRpcCalSim data category. ///////////////////////////////////////////////////////////////////// #include "hrpcdigitizer.h" #include "hades.h" #include "hphysicsconstants.h" #include "hruntimedb.h" #include "hdebug.h" #include "hspectrometer.h" #include "hrpcdetector.h" #include "hevent.h" #include "hcategory.h" #include "hiterator.h" #include "hlocation.h" #include "rpcdef.h" #include "hgeantrpc.h" #include "hgeantkine.h" #include "hrpccalsim.h" #include "hrpccal.h" #include "hrpcdigipar.h" #include "hrpcgeomcellpar.h" #include "TRandom.h" #include #include using namespace std; HRpcDigitizer::HRpcDigitizer(void) { initVars(); } HRpcDigitizer::HRpcDigitizer(const Text_t *name, const Text_t *title) : HReconstructor(name,title) { initVars(); } HRpcDigitizer::~HRpcDigitizer(void) { if(iterGeantRpc) delete iterGeantRpc; if(iterRpcCal) delete iterRpcCal; clearObjects(); cellobjects.clear(); rpcobjects.clear(); } void HRpcDigitizer::initVars(){ fGeantRpcCat = 0; fCalCat = 0; fKineCat = 0; fGeomCellPar = 0; fDigiPar = 0; fLoc.set(3,0,0,0); iterGeantRpc = 0; iterRpcCal = 0; } Bool_t HRpcDigitizer::initParContainer() { fGeomCellPar=(HRpcGeomCellPar *)gHades->getRuntimeDb()->getContainer("RpcGeomCellPar"); if(!fGeomCellPar){ Error("HRpcDigitizer::init()","No RpcGeomCellPar Parameters"); return kFALSE; } fDigiPar =(HRpcDigiPar *) gHades->getRuntimeDb()->getContainer("RpcDigiPar"); if(!fDigiPar){ Error("HRpcDigitizer::init()","No RpcDigiPar Parameters"); return kFALSE; } return kTRUE; } Bool_t HRpcDigitizer::init(void) { HRpcDetector* rpc = (HRpcDetector*)(gHades->getSetup()->getDetector("Rpc")); if(!rpc){ Error("HRpcDigitizer::init()","No Rpc Detector found"); return kFALSE; } //--------------------------------------------------------------- // Working arrays maxCol = rpc->getMaxColumns(); maxCell = rpc->getMaxCells(); rpcobjects.resize(6*maxCol*maxCell, 0 ); // Size is constant over run time cellobjects.resize(6*maxCol*maxCell, 0 ); // Size is constant over run time //--------------------------------------------------------------- //--------------------------------------------------------------- // Parameters if(!initParContainer()) return kFALSE; //--------------------------------------------------------------- //--------------------------------------------------------------- // Data fKineCat = gHades->getCurrentEvent()->getCategory(catGeantKine); if(!fKineCat) { Error("HRpcDigitizer::init()","HGeant kine input missing"); return kFALSE; } fGeantRpcCat = gHades->getCurrentEvent()->getCategory(catRpcGeantRaw); if(!fGeantRpcCat) { Error("HRpcDigitizer::init()","HGeant RPC input missing"); return kFALSE; } // Build the Calibration category //Adding the embedding mode! if(gHades->getEmbeddingMode()>0) { //Getting the temporary category for calibrated rpc hits. fCalCat = gHades->getCurrentEvent()->getCategory(catRpcCalTmp); if(!fCalCat) { fCalCat=rpc->buildMatrixCategory("HRpcCalSim",0.5); gHades->getCurrentEvent()->addCategory(catRpcCalTmp,fCalCat,"Rpc"); } } else { fCalCat = gHades->getCurrentEvent()->getCategory(catRpcCal); if(!fCalCat) { fCalCat=rpc->buildMatrixCategory("HRpcCalSim",0.5); gHades->getCurrentEvent()->addCategory(catRpcCal,fCalCat,"Rpc"); } } iterGeantRpc = (HIterator *)fGeantRpcCat->MakeIterator("native"); iterRpcCal = (HIterator *)fCalCat ->MakeIterator("native"); //--------------------------------------------------------------- return kTRUE; } Int_t HRpcDigitizer::execute(void) { HGeantRpc* geantrpc = 0; HRpcCalSim* cal = 0; Float_t mass; Int_t RefL=-1,trackNumber,gpid,mode,ngap,HGeantRpc_version=-1; Int_t nHGeantRpc=0,nDigitizedRpc=0; // Load the relevant parameters vprop = fDigiPar->getVprop(); // [mm/ns] sigma_el = (fDigiPar->getSigmaX()) / vprop; // [ns] sigma0_T = (fDigiPar->getSigmaT()) / 1000.; // [ns] sigma1_T = (fDigiPar->getSigmaT1()) / 1000.; // [ns] sigma2_T = (fDigiPar->getSigmaT2()); // pure number sigma3_T = (fDigiPar->getSigmaT3()); // pure number t_offset = (fDigiPar->getToff()) / 1000.; // [ns] Qmean0 = (fDigiPar->getQmean()); // [pC] Qmean1 = (fDigiPar->getQmean1()); // [pC] Qmean2 = (fDigiPar->getQmean2()); // [pC] Qwid0 = (fDigiPar->getQwid()); // [pC] Qwid1 = (fDigiPar->getQwid1()); // [pC] Qwid2 = (fDigiPar->getQwid2()); // [pC] Eff0 = (fDigiPar->getEff()); // pure number Eff1 = (fDigiPar->getEff1()); // pure number Eff2 = (fDigiPar->getEff2()); // pure number Eff3 = (fDigiPar->getEff3()); // pure number gap = fDigiPar->getGap(); // [mm] mode = fDigiPar->getMode(); // mode flag // if((sigma1_T<=0.)||(sigma2_T<=0.)||(sigma3_T<=0.)) { // Digitizing model I. sigma_el *= TMath::Sqrt2(); } // //if((mode==0)&&(HGeantRpc_version==5)) { // Error("HRpcDigitizer::init()","HGeantRpc version 5 does not allow gap mode"); // return kFALSE; //} //--------------------------------------------------------------- // Loop over the HGeantrpc objects, i.e, over the gaps/cells // and fill temporary working objects clearObjects(); iterGeantRpc->Reset(); nHGeantRpc=0; nDigitizedRpc=0; while ((geantrpc=(HGeantRpc *)iterGeantRpc->Next()) != 0) { geantrpc->setVersion(5); // Skip hit in virtual volume EBOX. if(geantrpc->getDetectorID() < 0) continue; // Count number of HGeantRpc objects excluding hits // in virtual volume EBOX. nHGeantRpc++; // Initialise HGeantRpc version and efficiency caluclation. if(HGeantRpc_version==-1) { HGeantRpc_version = geantrpc->getVersion(); if(HGeantRpc_version<=4) { // Initialise the inefficency of a gap. calc_eff_hit(-1); } else { // Initialise the inefficency of a cell. calc_eff_hit(-2); } } // Retrive hit. // Units are [mm], [ns], [MeV/c], [mm] in module ref system. if(HGeantRpc_version<=4) { // Get simulation output of a gap. geantrpc->getHitDigi(geaX, geaTof, geaMom, geaLocLen); } else { // Get simulation output of a cell. geantrpc->getCellAverageDigi(gap, geaX, geaTof, geaMom, geaLocLen); } // If the gap enetering was missed in Geant, assume the gap as local track length. if(geaLocLen<0.) { geaLocLen=gap; geaLocLenNorm=1.; } else { geaLocLenNorm=geaLocLen/gap; } // Determine the beta of the particle. trackNumber = geantrpc->getTrack(); HGeantKine* kine = (HGeantKine*)fKineCat->getObject(trackNumber-1); if(!kine) { Warning("HRpcDigitizer::execute()", "missing kine entry for track %i, assuming beta=0",trackNumber); beta=0.; } else { gpid = kine->getID(); // From the Geant particle type find the mass (only charged // particles need to be considered). mass = HPhysicsConstants::mass(gpid); if(mass <= 0 ){ if(mass < 0 ) Warning("HRpcDigitizer::execute()", "unknown particle id %i, assuming proton mass",gpid); mass = HPhysicsConstants::mass(14); // unknown + gamma } beta=mass/geaMom; beta=1./TMath::Sqrt(1.+beta*beta); } fLoc[0] = geantrpc->getSector(); fLoc[1] = geantrpc->getColumn(); fLoc[2] = geantrpc->getCell(); // Apply efficiency. if(HGeantRpc_version<=4) { // Calculate the efficiency of a gap. calc_eff_hit(1); } else { // Calculate the efficiency of a cell. calc_eff_hit(2); } // This may depend a bit on the digitization algorithm. if(gRandom->Uniform(1) > eff_hit) continue; if((mode==0)||(HGeantRpc_version==5)) { ///////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////// //////////// BEGIN GAP WISE OPEREATION (MODE=0) ///////////// //////////// WITH OLD HGeantRpc before Jan 2013 ///////////// //////////// BEGIN CELL WISE OPEREATION ///////////// //////////// WITH NEW HGeantRpc after Jan 2013 ///////////// ///////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////// //--------------------------------------------------------- // Collect data // objects are stored linear Int_t ind = fLoc[0] * maxCol * maxCell + fLoc[1] * maxCell + fLoc[2]; if(rpcobjects[ind] == 0) { // First time we have to create the object rpcobjects[ind] = new rpcdat; //rpcobjects[ind]->linIndex= ind; rpcobjects[ind]->sec = fLoc[0]; rpcobjects[ind]->col = fLoc[1]; rpcobjects[ind]->cell = fLoc[2]; } rpcdat* dat = rpcobjects[ind]; //--------------------------------------------------------- // Create new objects for left/right of the gap // and add them to the list of gaps for this cell gaptrack* left = new gaptrack; gaptrack* right = new gaptrack; left ->reset(); right->reset(); dat->left .push_back(left); dat->right.push_back(right); //--------------------------------------------------------- left ->refDgtr = fGeantRpcCat->getIndex(geantrpc); left ->trackDgtr = ((HGeantRpc*)fGeantRpcCat->getObject(left->refDgtr))->getTrack(); // Right is the same as left right->refDgtr = left->refDgtr; right->trackDgtr = left->trackDgtr; RefL = findMother(left->refDgtr); // Right is the same as left , do it only once if(RefL < 0){ // Mother not found. Very seldom. Keep info of daughter. left ->ref = left ->refDgtr; left ->track = left ->trackDgtr; left ->isAtBox = kFALSE; right->ref = left ->refDgtr; right->track = left ->trackDgtr; right->isAtBox = kFALSE; Warning("HRpcDigitizer::execute()", "mother of track not found in RPC box! will use RefLDgtr=%i, TrackLDgtr=%i", left->refDgtr,left->trackDgtr); } else { // Mother found left ->track = ((HGeantRpc*)fGeantRpcCat->getObject(RefL))->getTrack(); left ->ref = RefL; left ->isAtBox = kTRUE; right->track = left->track; right->ref = left->ref; right->isAtBox = kTRUE; } // Cell length in mm D = fGeomCellPar->getLength(fLoc[1],fLoc[2]); if(HGeantRpc_version<=4) { // Apply the response of a gap. digitize_one_hit(left,right,0); } else { // Apply the response of a cell. digitize_one_hit(left,right,1); } nDigitizedRpc++; //--------------------------------------------------------- ///////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////// //////////// END GAP WISE OPEREATION (MODE=0) ///////////// //////////// WITH OLD HGeantRpc before Jan 2013 ///////////// //////////// END CELL WISE OPEREATION ///////////// //////////// WITH NEW HGeantRpc after Jan 2013 ///////////// ///////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////// } else { ///////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////// //////////// BEGIN CELL WISE OPEREATION (MODE=1) //////////// //////////// WITH OLD HGeantRpc before Jan 2013 //////////// ///////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////// //--------------------------------------------------------- // Collect data // objects are stored linear Int_t ind = fLoc[0] * maxCol * maxCell + fLoc[1] * maxCell + fLoc[2]; if(cellobjects[ind] == 0) { // First time we have to create the object cellobjects[ind] = new celldat; //rpcobjects[ind]->linIndex= ind; cellobjects[ind]->sec = fLoc[0]; cellobjects[ind]->col = fLoc[1]; cellobjects[ind]->cell = fLoc[2]; } celldat* cdat = cellobjects[ind]; //--------------------------------------------------------- // Create new objects for the track crossing the cell // and add them to the list of tracks for this cell celltrack* celltr; // CASE I: no tacks in this cell yet if(cdat->celltr.size()==0) { celltr = new celltrack; celltr->reset(); cdat->celltr.push_back(celltr); // CASE II: this cell has already tracks associated // with it } else { // Look if current track is already present in the list for(UInt_t i=0;icelltr.size();i++){ celltr = cdat->celltr[i]; // CASE II a: the current track is already associated, // so select it if(celltr->track == trackNumber) break; } // CASE II b: the current track is not associated, // so add it if(celltr->track != trackNumber) { celltr = new celltrack; celltr->reset(); cdat->celltr.push_back(celltr); } } //--------------------------------------------------------- ngap = geantrpc->getGap(); // The condition on the time takes care of possible multiple // re-entries of the same track. if(celltr->gaptime[ngap]==0.) { celltr->gaptime[ngap] = geaTof; celltr->gappos[ngap] = geaX; celltr->gapltln[ngap] = geaLocLenNorm; celltr->gapbeta[ngap] = beta; celltr->track = trackNumber; celltr->geantrpcIndex = fGeantRpcCat->getIndex(geantrpc); } else if(celltr->gaptime[ngap]gaptime[ngap] = geaTof; celltr->gappos[ngap] = geaX; celltr->gapltln[ngap] = geaLocLenNorm; celltr->gapbeta[ngap] = beta; celltr->track = trackNumber; celltr->geantrpcIndex = fGeantRpcCat->getIndex(geantrpc); Warning("HRpcDigitizer::execute()", "reentry of track=%i in Sector=%i Module=%i Cell=%i", trackNumber,fLoc[0],fLoc[1],fLoc[2]); } ///////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////// //////////// END CELL WISE OPEREATION (MODE=1) ///////////// //////////// WITH OLD HGeantRpc before Jan 2013 ///////////// ///////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////// } // End of algorithm } // End of HGeantRpc loop //--------------------------------------------------------------- // If in cell mode and version of HGeantRpc is after Jan 2013, // process cells if((mode>0)&&(HGeantRpc_version<=4)&&(nHGeantRpc>0)) { for(UInt_t i = 0; i < cellobjects.size(); i ++) { celldat* cdat = cellobjects[i]; if(cdat){ fLoc[0] = cdat->sec; fLoc[1] = cdat->col; fLoc[2] = cdat->cell; // Change to cell coordinate system (with origin at the cell center) // Note that fGeomCellPar->getX(column,cell) corresponds to the left-down corner // Cell length in mm D = fGeomCellPar->getLength(fLoc[1],fLoc[2]); // Objects are stored linear Int_t ind = fLoc[0] * maxCol * maxCell + fLoc[1] * maxCell + fLoc[2]; if(rpcobjects[ind] == 0) { // First time we have to create the object rpcobjects[ind] = new rpcdat; //rpcobjects[ind]->linIndex= ind; rpcobjects[ind]->sec = fLoc[0]; rpcobjects[ind]->col = fLoc[1]; rpcobjects[ind]->cell = fLoc[2]; } rpcdat* dat = rpcobjects[ind]; for(UInt_t j = 0; j < cdat->celltr.size(); j ++) { celltrack* celltr = cdat->celltr[j]; celltr->calcMeans(); //------------------------------------------------------- // Create new objects for left/right of the gap // and add them to the list of gaps for this cell gaptrack* left = new gaptrack; gaptrack* right = new gaptrack; left ->reset(); right->reset(); dat->left .push_back(left); dat->right.push_back(right); //------------------------------------------------------- left ->refDgtr = celltr->geantrpcIndex; left ->trackDgtr = ((HGeantRpc*)fGeantRpcCat->getObject(left->refDgtr))->getTrack(); // Right is the same as left right->refDgtr = left->refDgtr; right->trackDgtr = left->trackDgtr; RefL = findMother(left->refDgtr); // Right is the same as left , do it only once if(RefL < 0){ // Mother not found. Very seldom. Keep info of daughter. left ->ref = left ->refDgtr; left ->track = left ->trackDgtr; left ->isAtBox = kFALSE; right->ref = left ->refDgtr; right->track = left ->trackDgtr; right->isAtBox = kFALSE; Warning("HRpcDigitizer::execute()", "mother of track not found in RPC box! will use RefLDgtr=%i, TrackLDgtr=%i", left->refDgtr,left->trackDgtr); } else { // Mother found left ->track = ((HGeantRpc*)fGeantRpcCat->getObject(RefL))->getTrack(); left ->ref = RefL; left ->isAtBox = kTRUE; right->track = left->track; right->ref = left->ref; right->isAtBox = kTRUE; } // Geant Position geaX = celltr->pos; // Geant Time geaTof = celltr->time; // Geant local track length normalised to the gap geaLocLenNorm = celltr->ltln; // Geant Beta beta = celltr->beta; // Apply response of a cell digitize_one_hit(left,right,1); nDigitizedRpc++; //------------------------------------------------------- } // End of celltrack loop } } // End of cellobjects loop } //--------------------------------------------------------------- // Fill output if(nDigitizedRpc>0) { Int_t Track [10]; Int_t TrackDgtr[10]; Bool_t isAtBox [10]; vector tracklist; tracklist.reserve(20); // Minimum size for(UInt_t i = 0; i < rpcobjects.size(); i ++) { rpcdat* dat = rpcobjects[i]; if(dat){ fLoc[0] = dat->sec; fLoc[1] = dat->col; fLoc[2] = dat->cell; //------------------------------------------------------- // Test if cell in use. Needed in embedding mode: //THIS OPTION IS NO MORE EXISTING: its used the tmp cat in the embedding mode // Neither in embedding mode nor in the standard this part should not be used! cal = (HRpcCalSim*) fCalCat->getObject(fLoc); if(cal) { // Add track from real data for sorting procedure gaptrack* left = new gaptrack; gaptrack* right = new gaptrack; left ->reset(); right->reset(); dat->left .push_back(left); dat->right.push_back(right); left ->time = cal->getLeftTime(); left ->gtime = cal->getLeftTime(); left ->charge = cal->getLeftCharge(); left ->track = gHades->getEmbeddingRealTrackId(); left ->trackDgtr = gHades->getEmbeddingRealTrackId(); left ->isAtBox = kTRUE; right->time = cal->getRightTime(); right->gtime = cal->getRightTime(); right->charge = cal->getRightCharge(); right->track = gHades->getEmbeddingRealTrackId(); right->trackDgtr = gHades->getEmbeddingRealTrackId(); right->isAtBox = kTRUE; } else { cal = (HRpcCalSim*) fCalCat->getSlot(fLoc); if(!cal) Error("execute()","Error: could not allocate a new slot in HRpcCalSim!"); cal = new(cal) HRpcCalSim; cal->setAddress(fLoc[0], fLoc[1], fLoc[2]); } //------------------------------------------------------- // Sort by GEANT time and select the fastest hit dat->sortTime(kTRUE); // Right dat->sortTime(kFALSE); // Left if(((mode==0)&&((sigma1_T<=0.)||(sigma2_T<=0.)||(sigma3_T<=0.))) ||(mode>0)||(HGeantRpc_version==5)) { // Digitizing model I and II (except gap-wise) // Take the fastest in the left and in the right. cal->setRightTime(dat->getSmallestTof(kTRUE )); // Right cal->setLeftTime (dat->getSmallestTof(kFALSE)); // Left } else { // Digitizing model II gap-wise. // Take the average in the left and in the right. cal->setRightTime(dat->getMeanTof(kTRUE )); // Right cal->setLeftTime (dat->getMeanTof(kFALSE)); // Left } cal->setRightCharge(dat->getSumCharge(kTRUE )); // Right cal->setLeftCharge (dat->getSumCharge(kFALSE)); // Left cal->setRefR (dat->right[0]->ref); cal->setRefRDgtr(dat->right[0]->refDgtr); cal->setRefL (dat->left [0]->ref); cal->setRefLDgtr(dat->left [0]->refDgtr); //------------------------------------------------------- //------------------------------------------------------- // Fill list of tracks // and remove doubles //------------------------------------------------------- // Right side tracklist.clear(); Int_t ct = 0; for(UInt_t j = 0; j < 10; j ++) { Track[j] = TrackDgtr[j] = -999; isAtBox[j] = kFALSE; } for(UInt_t j = 0; j < dat->right.size() && ct < 10; j ++ ){ Int_t tr = dat->right[j]->track; if(find(tracklist.begin(),tracklist.end(),tr ) == tracklist.end()){ // Track was not stored before Track [ct] = tr; TrackDgtr[ct] = dat->right[j]->trackDgtr; isAtBox [ct] = dat->right[j]->isAtBox; // Remember used track tracklist.push_back(tr); ct ++; } } cal->setTrackRArray (Track); cal->setTrackRDgtrArray(TrackDgtr); cal->setRisAtBoxArray (isAtBox); cal->setNTracksR(ct); //------------------------------------------------------- //------------------------------------------------------- // Left side tracklist.clear(); ct = 0; for(UInt_t j = 0; j < 10; j ++) { Track[j] = TrackDgtr[j] = -999; isAtBox[j] = kFALSE; } for(UInt_t j = 0; j < dat->left.size() && ct < 10 ; j ++ ){ Int_t tr = dat->left[j]->track; if(find(tracklist.begin(),tracklist.end(),tr) == tracklist.end()){ // Track was not stored before Track [ct] = tr; TrackDgtr[ct] = dat->left[j]->trackDgtr; isAtBox [ct] = dat->left[j]->isAtBox; // Remember used track tracklist.push_back(tr); ct ++; } } cal->setTrackLArray (Track); cal->setTrackLDgtrArray(TrackDgtr); cal->setLisAtBoxArray (isAtBox); cal->setNTracksL(ct); //------------------------------------------------------- //------------------------------------------------------- } } } //--------------------------------------------------------------- return 0; } void HRpcDigitizer::clearObjects(){ // Delete objects in working array and // set pointer to 0. vector is still // not cleared. for(UInt_t i = 0; i < rpcobjects.size(); i ++){ if(rpcobjects[i]){ rpcobjects[i]->reset(); delete rpcobjects[i]; rpcobjects[i] = 0; } } for(UInt_t i = 0; i < cellobjects.size(); i ++){ if(cellobjects[i]){ cellobjects[i]->reset(); delete cellobjects[i]; cellobjects[i] = 0; } } } Int_t HRpcDigitizer::findMother(Int_t Ref_initial) { // Algorithm for finding mother. When not found, -1 is returned // in order to tag the track. It very rarely happens and is // always associated to neutral particles at the box. Int_t detectorID = 0, track_mother = -1, Ref_final = -1; HGeantRpc* geantrpc = 0; HGeantKine* kine = 0; track_mother = ((HGeantRpc*)fGeantRpcCat->getObject(Ref_initial))->getTrack(); kine = (HGeantKine*)fKineCat ->getObject(track_mother - 1); while(kine && detectorID >= 0){ if(track_mother == 0) { //End of loop. Mother not found Ref_final = -999; break; } kine = (HGeantKine*)fKineCat->getObject(track_mother - 1); if(kine && kine->getNRpcHits() == 0) { //Mother has no hit in RPC cells track_mother = kine->getParentTrack(); continue; } Ref_final = kine->getFirstRpcHit(); //Look for the first hit (time-wise) of the track geantrpc = (HGeantRpc*)fGeantRpcCat->getObject(Ref_final); detectorID = geantrpc->getDetectorID(); track_mother = kine->getParentTrack(); } return Ref_final; } void HRpcDigitizer::calc_eff_hit(Int_t mode) { // Algorithm to calculate the inefficiency of a gap. The efficiency // in the digitizer parameters is ALWAYS per CELL. // If mode=-1 the efficiency/inefficiency per hit are initialised as per gap, // If mode=-2 the efficiency/inefficiency per hit are initialised as per cell, // If mode=1 the efficiency/inefficiency per hit are calculated as per gap, // If mode=2 the efficiency/inefficiency per hit are calculated as per cell. if(mode<0) { // // Initialisation call // if((Eff1==0.)||(Eff2<=0.)||(Eff3<=0.)) { // Digitizing model I. if(mode==-1) { // HIT == GAP ineff_hit_n = TMath::Sqrt(TMath::Sqrt(1. - fDigiPar->getEff())); } else { // HIT == CELL ineff_hit_n = 1. - fDigiPar->getEff(); } } else { ineff_hit_n = 0.; } } else { // // Calculation call // if((Eff1==0.)||(Eff2<=0.)||(Eff3<=0.)) { // Digitizing model I. ineff_hit = TMath::Power(ineff_hit_n,geaLocLenNorm); // correction for angle of incidence } else { // Digitizing model II. eff_hit = Eff0+Eff1/(1.+TMath::Exp(-Eff2*(beta-Eff3))); // efficeny of the cell if(mode==1) { // HIT == GAP ineff_hit = TMath::Sqrt(TMath::Sqrt(1. - eff_hit)); // convert to inefficency of a gap } else { // HIT == CELL ineff_hit = 1. - eff_hit; } ineff_hit = TMath::Power(ineff_hit,(geaLocLenNorm)); // correction for angle of incidence } // Efficiency per gap eff_hit = 1. - ineff_hit; } } void HRpcDigitizer::digitize_one_hit(gaptrack* left,gaptrack* right,Int_t mode) { // Algorithm to calculate the response of a gap. // If mode=0 in MODEL II a corection for the average number of active // gaps is applyed. Float_t geaXLoc = 0.; Float_t beta2,corr_ngap,sigma_T,sigma_el_here,X_smearing,T_smearing; Float_t Qmean,Qwid,Qgap; // Change to cell coordinate system (with origin at the cell center) // Note that fGeomCellPar->getX(column,cell) corresponds to the left-down corner geaXLoc = fGeomCellPar->getX(fLoc[1],fLoc[2]) - geaX; // Position in cell frame // Calculate GEANT time left/right without smearing for sorting left->gtime = geaTof + geaXLoc / vprop; // Geant Time at left-end right->gtime = geaTof + (D-geaXLoc) / vprop; // Geant Time at right-end // Common to left and right!. Intrinsic RPC fluctuations. if((sigma1_T<=0.)||(sigma2_T<=0.)||(sigma3_T<=0.)) { // Digitizing model I. // Common to left and right!. Intrinsic RPC fluctuations. T_smearing = gRandom->Gaus(t_offset, sigma0_T); left ->time = left ->gtime + T_smearing + gRandom->Gaus(0.,sigma_el); right->time = right->gtime + T_smearing + gRandom->Gaus(0.,sigma_el); } else { // Digitizing model II. sigma_T = sigma0_T+sigma1_T/(1.+TMath::Exp(-sigma2_T*(beta-sigma3_T))); //cout << sigma0_T << " " << sigma1_T << " " << sigma2_T << " " << sigma3_T << " " // << beta << " " << sigma_T << endl; sigma_el_here=sigma_el; // Attempt to approximately correct for the number of active gaps. if(mode==0) { Float_t norm=0.; corr_ngap=0.; for(Int_t i = 1; i < 5; i ++) { Float_t w4i=TMath::Binomial(4,i)*TMath::Power(eff_hit,i)*TMath::Power(ineff_hit,4-i); norm += w4i; corr_ngap += w4i/TMath::Sqrt((Float_t)i); } // The correction factor for the sigma of the Gaussian is the // inverse of corr_ngap/norm, to have a multiplicative correction // it is calculated here once for all. corr_ngap = norm/corr_ngap; //cout << beta << " " << Eff0+Eff1/(1.+TMath::Exp(-Eff2*(beta-Eff3))) // << " " << eff_hit << " " << corr_ngap << " " << sigma_T << endl; sigma_T *= corr_ngap; sigma_el_here *= corr_ngap; } T_smearing = gRandom->Gaus(t_offset,sigma_T); // The position and time resolutions are decoupled. X_smearing = gRandom->Gaus(0.,sigma_el_here); //X_smearing = 0.; left ->time = left ->gtime + T_smearing + X_smearing; right->time = right->gtime + T_smearing - X_smearing; //cout << fLoc[0] << " " << fLoc[1] << " " << fLoc[2] << " " // << geantrpc->getGap() << " " // << left ->gtime << " " << left ->time << " " // << right->gtime << " " << right->gtime << " " << endl; } // Generated charge. if(Qwid0<=0.) { // Digitizing model I. // Charge distribution assumed to be flat per gap. At present slewing // correction is not included. Also charges in left and in right are // *exactly* the same by definition. A smearing in charge to account // for the resolution of the Q2W method can be easily included. Qgap = (gRandom->Uniform(1)) * (2. * Qmean0) / 4. * geaLocLenNorm; // Correction for angle of incidence } else { // Digitizing model II. // Use Landau distribution with a mean and width parametrised with // second order polynomials in beta. beta2 = beta*beta; Qmean = Qmean0+Qmean1*beta+Qmean2*beta2; Qwid = Qwid0 +Qwid1*beta +Qwid2*beta2; // Attempt to approximately correct for the number of active gaps. if(mode==0) { Float_t norm=0.; corr_ngap=0.; for(Int_t i = 1; i < 5; i ++) { Float_t w4i=TMath::Binomial(4,i)*TMath::Power(eff_hit,i)*TMath::Power(ineff_hit,4-i); norm += w4i; corr_ngap += w4i*TMath::Sqrt((Float_t)i); } // The correction factor for the sigma of the Gaussian is the // inverse of corr_ngap/norm, to have a multiplicative correction // it is calculated here once for all. Note that this is very // approximate for a Landau distribution. corr_ngap = norm/corr_ngap; // Per gap. Qmean /= (4.*eff_hit); // Per gap. Qwid *= corr_ngap; } // Correction for angle of incidence. Qmean *= geaLocLenNorm; // Correction for angle of incidence. Qwid *= geaLocLenNorm; Qgap = gRandom->Landau(Qmean,Qwid); } left ->charge = Qgap; right->charge = Qgap; } ClassImp(HRpcDigitizer)