arich modules

Classes
class	arichBtestModule
	The UserTutorial module. More...
class	ARICHChannelMaskModule
	ARICH channel-mask calibration: data collector. More...
class	ARICHDigitizerModule
	ARICH digitizer module. More...
class	ARICHDQMModule
	Make summary of data quality from reconstruction. More...
class	ARICHFillHitsModule
	Fill ARICHHit collection from ARICHDigits. More...
class	ARICHMCParticlesModule
	Module to match ARICH hits to MCParticles. More...
class	ARICHNtupleModule
	ARICH reconstruction efficiency test module. More...
class	ARICHPackerModule
	Raw data packer for ARICH. More...
class	ARICHRateCalModule
	Fill ARICHHit collection from ARICHDigits. More...
class	ARICHRawUnpackerModule
	Fill ARICHHit collection from ARICHDigits. More...
class	ARICHReconstruction
	Internal ARICH track reconstruction. More...
class	ARICHReconstructorModule
	ARICH subdetector main module. More...
class	ARICHRelateModule
	Creates relations between ARICHAeroHits and ExtHits. More...
class	arichToNtupleModule
	This module extends existing variablesToNtuple to append detailed arich information to candidates in the analysis output ntuple. More...
struct	ARICHRawHeader
	ARICH raw-data header. More...
class	ARICHUnpackerModule
	Raw data unpacker for ARICH. More...

Macros
#define	ARICH_BUFFER_NWORDS   252
	Arich number of words (ints) in buffer: 3 + 33 + 6 * 36 (3 merger header words + 5.5 FEB header words / FEB + 36 data words per / FEB).
#define	ARICHFEB_HEADER_SIZE   10
	FEB header size in bytes.
#define	ARICHRAW_HEADER_SIZE   12
	Raw header size in bytes.

Functions
	REG_MODULE (arichBtest)
	Register the Module.
	REG_MODULE (ARICHDigitizer)
	Register the Module.
TH2 *	moduleHitMap (TH1 *hitMap, int moduleID)
	Make hit map in HAPD view (12*12 channels).
TH2 *	moduleDeadMap (TH1 *hitMap, int moduleID)
	Make chip dead/alive map in HAPD view (2*2 chips).
TH1 *	mergerClusterHitMap1D (TH1 *hitMap, int mergerID)
	Make 1D hit map of specified Merger Board.
TCanvas *	mergerClusterHitMap2D (TH1 *hitMap, int mergerID)
	Make display of 6 HAPDs' 2D hit map of the Merger Board.
TCanvas *	sectorHitMap (TH1 *hitMap, int sector)
	Make display of 70 HAPDs' 2D hit map of the sector.
TCanvas *	sectorDeadMap (TH1 *hitMap, int sector)
	Make display of 70 HAPDs' 2D dead/alive map of the sector.
	REG_MODULE (ARICHDQM)
void	deadPalette ()
	Set palette for sector dead-chip map.
	REG_MODULE (ARICHFillHits)
	Register module.
	REG_MODULE (ARICHMCParticles)
	Register module.
	REG_MODULE (ARICHNtuple)
	Register module.
	REG_MODULE (ARICHPacker)
	REG_MODULE (ARICHRateCal)
	REG_MODULE (ARICHRawUnpacker)
	REG_MODULE (ARICHReconstructor)
	Register the Module.
	REG_MODULE (ARICHRelate)
	Register module.
	REG_MODULE (ARICHUnpacker)
	Register module.
	arichBtestModule ()
	Constructor.
virtual void	initialize () override
	Initialize the Module.
virtual void	beginRun () override
	Called when entering a new run.
int	skipdata (gzFile fp)
	Skip the data part of the record.
void	readmwpc (unsigned int *dbuf, unsigned int len, int print=0)
	Read the MWPC information from the data buffer.
int	readhapd (unsigned int len, unsigned int *data)
	Read the HAPD hits from the data buffer.
int	getTrack (int mask, ROOT::Math::XYZVector &r, ROOT::Math::XYZVector &dir)
	Beamtest Track reconstruction.
int	readdata (gzFile fp, int rec_id, int print)
	Read the data from the file (can be compressed)
virtual void	event () override
	Running over all events.
virtual void	endRun () override
	Is called after processing the last event of a run.
virtual void	terminate () override
	Is called at the end of your Module.
	ARICHDigitizerModule ()
	Constructor.
virtual	~ARICHDigitizerModule ()
	Destructor.
virtual void	initialize () override
	Initialize the Module.
virtual void	beginRun () override
	Called when entering a new run.
virtual void	event () override
	Event processor.
void	magFieldDistorsion (ROOT::Math::XYVector &hit, int copyno)
	Apply correction to hit position due to non-perpendicular component of magnetic field.
	ARICHDQMModule ()
	Constructor.
virtual	~ARICHDQMModule ()
	Destructor.
virtual void	defineHisto () override
	Definition of the histograms.
virtual void	initialize () override
	Initialize the Module.
virtual void	beginRun () override
	Called when entering a new run.
virtual void	event () override
	Event processor.
virtual void	endRun () override
	End-of-run action.
	ARICHFillHitsModule ()
	Constructor.
virtual	~ARICHFillHitsModule ()
	Destructor.
virtual void	initialize () override
	Initialize the Module.
virtual void	event () override
	Event processor.
void	magFieldCorrection (ROOT::Math::XYZVector &hitpos)
	Corrects hit position for distorsion due to non-perpendicular magnetic field component.
	ARICHMCParticlesModule ()
	Constructor.
virtual	~ARICHMCParticlesModule ()
	Destructor.
virtual void	initialize () override
	Initialize the Module.
virtual void	event () override
	Event processor.
	ARICHNtupleModule ()
	Constructor.
virtual	~ARICHNtupleModule ()
	Destructor.
virtual void	initialize () override
	Initialize the Module.
virtual void	event () override
	Event processor.
virtual void	terminate () override
	Termination action.
	ARICHPackerModule ()
	Constructor.
virtual	~ARICHPackerModule ()
	Destructor.
virtual void	initialize () override
	Initialize the Module.
virtual void	event () override
	Event processor.
void	writeHeader (int *buffer, unsigned &ibyte, const ARICHRawHeader &head)
	TODO!
unsigned int	calbyte (const int *buf)
	Get calculated byte.
unsigned int	calword (const int *buf)
	Get calculated word.
	ARICHRateCalModule ()
	Constructor.
virtual	~ARICHRateCalModule ()
	Destructor.
virtual void	defineHisto () override
	Definition of the histograms.
virtual void	initialize () override
	Initialize the Module.
virtual void	beginRun () override
	Called when entering a new run.
virtual void	event () override
	Event processor.
unsigned int	calbyte (const int *buf)
	Read byte with number m_ibyte from the buffer and increase the number by 1.
unsigned int	calword (const int *buf)
	Read word (4 bytes) from the buffer and increase the byte number m_ibyte by 4.
	ARICHRawUnpackerModule ()
	Constructor.
virtual	~ARICHRawUnpackerModule ()
	Destructor.
virtual void	defineHisto () override
	Definition of the histograms.
virtual void	initialize () override
	Initialize the Module.
virtual void	event () override
	Event processor.
	ARICHReconstruction (int storePhotons=0)
	Constructor.
void	initialize ()
	Read geometry parameters from xml and initialize class members.
int	InsideDetector (ROOT::Math::XYZVector a, int copyno)
	Returns 1 if vector "a" lies on "copyno"-th detector active surface of detector and 0 else.
int	smearTrack (ARICHTrack &arichTrack)
	Smears track parameters ("simulate" the uncertainties of tracking).
ROOT::Math::XYZVector	FastTracking (ROOT::Math::XYZVector dirf, ROOT::Math::XYZVector r, double *refind, double *z, int n, int opt)
	Calculates the intersection of the Cherenkov photon emitted from point "r" in "dirf" direction with the detector plane.
ROOT::Math::XYZVector	HitVirtualPosition (const ROOT::Math::XYZVector &hitpos, int mirrorID)
	Returns the hit virtual position, assuming that it was reflected from mirror.
bool	HitsMirror (const ROOT::Math::XYZVector &pos, const ROOT::Math::XYZVector &dir, int mirrorID)
	Returns true if photon at position pos with direction dir hits mirror plate with ID mirrorID.
int	CherenkovPhoton (ROOT::Math::XYZVector r, ROOT::Math::XYZVector rh, ROOT::Math::XYZVector &rf, ROOT::Math::XYZVector &dirf, double *refind, double *z, int n, int mirrorID=0)
	Calculates the direction of photon emission.
int	likelihood2 (ARICHTrack &arichTrack, const StoreArray< ARICHHit > &arichHits, ARICHLikelihood &arichLikelihood)
	Computes the value of identity likelihood function for different particle hypotheses.
void	setTrackPositionResolution (double pRes)
	Sets track position resolution (from tracking).
void	setTrackAngleResolution (double aRes)
	Sets track direction resolution (from tracking).
ROOT::Math::XYZVector	getTrackMeanEmissionPosition (const ARICHTrack &track, int iAero)
	Returns mean emission position of Cherenkov photons from i-th aerogel layer.
ROOT::Math::XYZVector	getTrackPositionAtZ (const ARICHTrack &track, double zout)
	Returns track direction at point with z coordinate "zout" (assumes straight track).
void	transformTrackToLocal (ARICHTrack &arichTrack, bool align)
	Transforms track parameters from global Belle2 to ARICH local frame.
ROOT::Math::XYZVector	getMirrorPoint (int mirrorID)
	Returns point on the mirror plate with id mirrorID.
ROOT::Math::XYZVector	getMirrorNorm (int mirrorID)
	Returns normal vector of the mirror plate with id mirrorID.
void	correctEmissionPoint (int tileID, double r)
	Correct mean emission point z position.
	ARICHReconstructorModule ()
	Constructor.
virtual	~ARICHReconstructorModule ()
	Destructor.
virtual void	initialize () override
	Initialize the Module.
virtual void	beginRun () override
	Called when entering a new run.
virtual void	event () override
	Event processor.
void	printModuleParams ()
	Print module parameters.
	ARICHRelateModule ()
	Constructor.
virtual void	initialize () override
	Initialize the Module.
virtual void	event () override
	Event processor.
unsigned int	calbyte (const int *buf)
	calculate number of bytes in raw Unpacker
unsigned int	cal2byte (const int *buf)
	calculate number of lines (2 bytes) in raw Unpacker
unsigned int	calword (const int *buf)
	calculate number of words in raw Unpacker
	ARICHUnpackerModule ()
	Constructor.
virtual	~ARICHUnpackerModule ()
	Destructor.
virtual void	initialize () override
	Initialize the Module.
virtual void	event () override
	Event processor.
void	readHeader (const int *buffer, unsigned &ibyte, ARICHRawHeader &head)
	Read Merger header.
void	readFEHeader (const int *buffer, unsigned &ibyte, ARICHRawHeader &head)
	Read FE header.
void	printBits (const int *buffer, int bufferSize)
	Unpack raw data given in production format.

Variables
TNtuple *	m_tuple
	ntuple containing hapd hits
TH1F *	hapd [6]
	histogram of hits for each hapd
TH1F *	mwpc_tdc [4][5]
	tdc information from mwpcs
TH1F *	mwpc_diff [4][2]
	tdc difference from mwpcs
TH1F *	mwpc_sum [4][2]
	tdc sum from mwpcs
TH1F *	mwpc_sum_cut [4][2]
	tdc sum from mwpcs, with sum cut applied
TH1F *	mwpc_residuals [4][2]
	residuals from mwpcs
TH2F *	mwpc_xy [4]
	calculated x-y track positions
TH2F *	mwpc_residualsz [4][2]
	z-residuals from mwpcs

Detailed Description

Macro Definition Documentation

ARICH_BUFFER_NWORDS

#define ARICH_BUFFER_NWORDS   252

Arich number of words (ints) in buffer: 3 + 33 + 6 * 36 (3 merger header words + 5.5 FEB header words / FEB + 36 data words per / FEB).

Definition at line 35 of file ARICHPackerModule.cc.

ARICHFEB_HEADER_SIZE

#define ARICHFEB_HEADER_SIZE   10

FEB header size in bytes.

Definition at line 21 of file ARICHRawDataHeader.h.

ARICHRAW_HEADER_SIZE

#define ARICHRAW_HEADER_SIZE   12

Raw header size in bytes.

Definition at line 24 of file ARICHRawDataHeader.h.

Function Documentation

arichBtestModule()

arichBtestModule

(

)

Constructor.

Sets the module parameters.

Definition at line 63 of file arichBtestModule.cc.

                                     : Module(), m_end(0), m_events(0), m_file(NULL), m_timestart(0), m_mwpc(NULL)
  {
    //Set module properties
    setDescription("Module for the ARICH Beamtest data analysis. It creates track form the MWPC hits and reads the HAPD hits");
 
    //Parameter definition
    addParam("outputFileName", m_outFile, "Output Root Filename", std::string("output.root"));
    std::vector<std::string> defaultList;
    addParam("runList", m_runList, "Data Filenames.", defaultList);
    std::vector<int> defaultMask;
    addParam("mwpcTrackMask", m_MwpcTrackMask, "Create track from MWPC layers", defaultMask);
    m_fp = NULL;
    addParam("beamMomentum", m_beamMomentum, "Momentum of the beam [GeV]", 0.0);
 
    for (int i = 0; i < 32; i++) m_tdc[i] = 0;
 
  }

ARICHDigitizerModule()

ARICHDigitizerModule

(

)

Constructor.

Definition at line 51 of file ARICHDigitizerModule.cc.

                                             :  Module(),
    m_maxQE(0)
  {
 
    // Set description()
    setDescription("This module creates ARICHDigits from ARICHSimHits. Here spatial digitization is done, channel-by-channel QE is applied, and readout time window cut is applied.");
 
    // Set property flags
    setPropertyFlags(c_ParallelProcessingCertified);
 
    // Add parameters
    addParam("TimeWindow", m_timeWindow, "Readout time window width in ns", 250.0);
    addParam("TimeWindowStart", m_timeWindowStart, "Shift of the readout time window with respect to the global zero in ns", -50.0);
    addParam("BackgroundHits", m_bkgLevel, "Number of background hits per hapd per readout (electronics noise)", 0.066);
    addParam("MagneticFieldDistorsion", m_bdistort, "Apply image distortion due to non-perp. magnetic field", 0);
  }

ARICHDQMModule()

ARICHDQMModule

(

)

Constructor.

Definition at line 44 of file ARICHDQMModule.cc.

                                 : HistoModule()
  {
    // set module description (e.g. insert text)
    setDescription("Make summary of data quality.");
    setPropertyFlags(c_ParallelProcessingCertified);
    addParam("debug", m_debug, "debug mode", false);
    addParam("UpperMomentumLimit", m_momUpLim, "Upper momentum limit of tracks included in monitoring", 10.0);
    addParam("LowerMomentumLimit", m_momDnLim, "Lower momentum limit of tracks included in monitoring", 2.5);
    addParam("ArichEvents", m_arichEvents, "Include only hits from events where an extrapolated track to arich exists", false);
    addParam("MaxHits", m_maxHits, "Include only events with less than MaxHits hits in ARICH (remove loud events)", 70000);
    addParam("MinHits", m_minHits, "Include only events with more than MinHits hits in ARICH", 0);
  }

ARICHFillHitsModule()

ARICHFillHitsModule

(

)

Constructor.

Definition at line 41 of file ARICHFillHitsModule.cc.

                                           : Module()
  {
    // set module description (e.g. insert text)
    setDescription("Fills ARICHHits collection from ARICHDigits");
    setPropertyFlags(c_ParallelProcessingCertified);
    addParam("bitMask", m_bitMask, "hit bit mask (8 bits/channel)", (uint8_t)0xFF);
    addParam("maxApdHits", m_maxApdHits, "Remove hits with more than MaxApdHits per APD chip", (uint8_t)18);
    addParam("maxHapdHits", m_maxHapdHits, "Remove hits with more than MaxHapdHits per HAPD", (uint8_t)100);
    addParam("MagFieldCorrection", m_bcorrect, "Apply hit position correction due to non-perp. mag. field", 0);
    addParam("fillAll", m_fillall, "Make hits for all active channels (useful for likelihood PDF studies)", 0);
  }

ARICHMCParticlesModule()

ARICHMCParticlesModule

(

)

Constructor.

Definition at line 39 of file ARICHMCParticlesModule.cc.

  {
    // set module description (e.g. insert text)
    setDescription("Creates collection of MCParticles related to tracks that hit ARICH.");
    // Set property flags
    setPropertyFlags(c_ParallelProcessingCertified);
 
  }

ARICHNtupleModule()

ARICHNtupleModule

(

)

Constructor.

Definition at line 54 of file ARICHNtupleModule.cc.

                                       : Module(),
    m_file(0),
    m_tree(0)
  {
    // set module description (e.g. insert text)
    setDescription("The module saves variables needed for performance analysis, such as postion and momentum of the hit, likelihoods for hypotheses and number of photons.");
 
    // Add parameters
    addParam("outputFile", m_outputFile, "ROOT output file name", std::string("ARICHNtuple.root"));
  }

ARICHPackerModule()

ARICHPackerModule

(

)

Constructor.

Definition at line 46 of file ARICHPackerModule.cc.

                                       : Module(),
    m_nonSuppressed(0), m_bitMask(0), m_debug(0)
 
  {
    // set module description (e.g. insert text)
    setDescription("Raw data packer for ARICH");
    setPropertyFlags(c_ParallelProcessingCertified);
    addParam("nonSuppressed", m_nonSuppressed, "Pack in non-suppressed format (store all channels)", unsigned(0));
    addParam("version", m_version, "dataformat version", unsigned(6));
    addParam("bitMask", m_bitMask, "hit bit mask (4 bits/channel)", (unsigned)0xF);
    addParam("debug", m_debug, "print packed bitmap", 0);
    addParam("inputDigitsName", m_inputDigitsName, "name of ARICHDigit store array", std::string(""));
    addParam("outputRawDataName", m_outputRawDataName, "name of RawARICH store array", std::string(""));
 
  }

ARICHRateCalModule()

ARICHRateCalModule

(

)

Constructor.

Definition at line 45 of file ARICHRateCalModule.cc.

                                         : HistoModule()
  {
    // set module description (e.g. insert text)
    setDescription("Fills ARICHHits collection from ARICHDigits");
    setPropertyFlags(c_ParallelProcessingCertified);
    addParam("nrun", m_nrun, "# of scan runs", 100);
    addParam("nevents", m_nevents, "# of events per run", 1000);
    double v = 0;
    addParam("dth", m_dth, "dth", v);
    addParam("th0", m_th0, "th0", v);
    addParam("debug", m_debugmode, "debug mode", false);
    addParam("internal", m_internalmode, "Internal thscan mode", false);
    addParam("daqdb", m_daqdb, "daqdb config name", std::string(""));
  }

ARICHRawUnpackerModule()

ARICHRawUnpackerModule

(

)

Constructor.

Definition at line 36 of file ARICHRawUnpackerModule.cc.

                                                 : HistoModule()
  {
    m_debug = false;
  }

ARICHReconstruction()

ARICHReconstruction ( int  storePhotons = 0 )

explicit

Constructor.

Definition at line 38 of file ARICHReconstruction.cc.

                                                       :
    m_arichgp(),
    m_recPars(),
    m_trackPosRes(0),
    m_trackAngRes(0),
    m_alignMirrors(true),
    m_nAerogelLayers(0),
    m_storePhot(storePhot)
  {
    for (unsigned i = 0; i < c_noOfHypotheses; i++) {p_mass[i] = 0;}
    for (unsigned i = 0; i < c_noOfAerogels; i++) {
      m_refractiveInd[i] = 0;
      m_zaero[i] = 0;
      m_thickness[i] = 0;
      m_transmissionLen[i] = 0;
      m_n0[i] = 0;
      m_anorm[i] = ROOT::Math::XYZVector();
    }
 
  }

ARICHReconstructorModule()

ARICHReconstructorModule

(

)

Constructor.

Definition at line 49 of file ARICHReconstructorModule.cc.

                                                     :
    m_ana(NULL)
  {
    // Set description()
    setDescription("This module calculates the ARICHLikelihood values for all particle id. hypotheses, for all tracks that enter ARICH in the event.");
 
    // Set property flags
    setPropertyFlags(c_ParallelProcessingCertified);
 
    // Add parameters
    addParam("trackPositionResolution", m_trackPositionResolution,
             "Resolution of track position on aerogel plane (for additional smearing of MC tracks)", 1.0 * Unit::mm);
    addParam("trackAngleResolution", m_trackAngleResolution,
             "Resolution of track direction angle on aerogel plane (for additional smearing of MC tracks)", 2.0 * Unit::mrad);
    addParam("inputTrackType", m_inputTrackType, "Input tracks switch: tracking (0), from AeroHits - MC info (1)", 0);
    addParam("storePhotons", m_storePhot, "Set to 1 to store reconstructed photon information (Ch. angle,...)", 0);
    addParam("useAlignment", m_align, "Use ARICH global position alignment constants", true);
    addParam("useMirrorAlignment", m_alignMirrors, "Use ARICH mirror alignment constants", true);
  }

ARICHRelateModule()

ARICHRelateModule

(

)

Constructor.

Definition at line 40 of file ARICHRelateModule.cc.

  {
    // set module description (e.g. insert text)
    setDescription("Creates relations between ARICHAeroHits and ExtHits. Allows to store simulation output without MCParticles");
    // Set property flags
    setPropertyFlags(c_ParallelProcessingCertified);
 
  }

ARICHUnpackerModule()

ARICHUnpackerModule

(

)

Constructor.

Definition at line 51 of file ARICHUnpackerModule.cc.

                                           : Module(), m_bitMask(0), m_debug(0)
  {
    // set module description
    setDescription("Raw data unpacker for ARICH");
    setPropertyFlags(c_ParallelProcessingCertified);
 
    addParam("bitMask", m_bitMask, "hit bit mask (8 bits/channel, only used for unsuppresed format!)", (uint8_t)0xFF);
    addParam("debug", m_debug, "prints debug information", 0);
 
    addParam("inputRawDataName", m_inputRawDataName, "name of RawARICH store array", std::string(""));
    addParam("outputDigitsName", m_outputDigitsName, "name of ARICHDigit store array", std::string(""));
    addParam("outputRawDigitsName", m_outputRawDigitsName, "name of ARICHRawDigit store array", std::string(""));
    addParam("outputarichinfoName", m_outputarichinfoName, "name of ARICHInfo store array", std::string(""));
    addParam("RawUnpackerMode", m_rawmode, "Activate RawUnpacker mode", 0);
    addParam("DisableUnpackerMode", m_disable_unpacker, "Disable Regular Unpacker mode", 0);
 
  }

beginRun() [1/5]

void beginRun ( void  )

overridevirtual

Called when entering a new run.

At the beginning of each run, the function gives you the chance to change run dependent constants like alignment parameters, etc.

Reimplemented from Module.

Definition at line 149 of file arichBtestModule.cc.

  {
 
    B2INFO("arichBtestModule::beginRun()");
 
    StoreObjPtr<EventMetaData> eventMetaDataPtr;
    B2INFO("arichBtestModule::eventMetaDataPtr run:" << eventMetaDataPtr->getRun());
    B2INFO("arichBtestModule::eventMetaDataPtr exp:" << eventMetaDataPtr->getExperiment());
 
    ARICHBtestGeometryPar* _arichbtgp = ARICHBtestGeometryPar::Instance();
    m_mwpc = _arichbtgp->getMwpc();
 
    static int first = 1;
    if (first) {
      m_runCurrent = m_runList.begin();
      first = 0;
      m_mwpc[0].Print();
      m_mwpc[1].Print();
      m_mwpc[2].Print();
      m_mwpc[3].Print();
 
    }
    m_end = 1;
 
    if (m_runCurrent < m_runList.end()) {
      if (m_fp) {
        m_eveList.push_back(m_events);
        gzclose(m_fp);
      }
      m_fp = gzopen((*m_runCurrent).c_str(), "rb");
      if (m_fp == NULL) {
        B2INFO("Cannot open " << *m_runCurrent);
        m_end = 1;
      } else {
        B2INFO("File opened " << *m_runCurrent);
        m_end = 0;
      }
    }
    m_events = 0;
  }

beginRun() [2/5]

void beginRun ( void  )

overridevirtual

Called when entering a new run.

Set run dependent things like run header parameters, alignment, etc.

Reimplemented from Module.

Definition at line 89 of file ARICHDigitizerModule.cc.

  {
    if (m_simPar->getNBkgHits() > 0)  m_bkgLevel = m_simPar->getNBkgHits();
  }

beginRun() [3/5]

void beginRun ( void  )

overridevirtual

Called when entering a new run.

Set run dependent things like run header parameters, alignment, etc.

Reimplemented from HistoModule.

Definition at line 175 of file ARICHDQMModule.cc.

  {
 
    h_chStat->Reset();
    h_aeroStat->Reset();
    h_chDigit->Reset();
    h_chipDigit->Reset();
    h_hapdDigit->Reset();
 
    h_chHit->Reset();
    h_chipHit->Reset();
    h_hapdHit->Reset();
    h_mergerHit->Reset();
    h_bitsPerMergerNorm->Reset();
    h_bitsPerHapdMerger->Reset();
 
    h_aerogelHit->Reset();
    h_bits->Reset();
    h_bitsPerChannel->Reset();
    h_hitsPerTrack2D->Reset();
    h_tracks2D->Reset();
    h_aerogelHits3D->Reset();
 
    h_hitsPerEvent->Reset();
    h_theta->Reset();
    h_hitsPerTrack->Reset();
    h_trackPerEvent->Reset();
    h_flashPerAPD->Reset();
    h_mirrorThetaPhi->Reset();
    h_thetaPhi->Reset();
 
    for (int i = 0; i < 6; i++) {
      h_secTheta[i]->Reset();
      h_secHitsPerTrack[i]->Reset();
    }
 
    h_ARICHOccAfterInjLer->Reset();
    h_ARICHEOccAfterInjLer->Reset();
    h_ARICHOccAfterInjHer->Reset();
    h_ARICHEOccAfterInjHer->Reset();
  }

beginRun() [4/5]

void beginRun ( void  )

overridevirtual

Called when entering a new run.

Set run dependent things like run header parameters, alignment, etc.

Reimplemented from HistoModule.

Definition at line 103 of file ARICHRateCalModule.cc.

  {
    StoreObjPtr<EventMetaData> evtmetadata;
    //int expno = evtmetadata->getExperiment();
    int runno = evtmetadata->getRun();
    if (runno == 100 && !m_internalmode) {
      terminate();
    }
  }

beginRun() [5/5]

void beginRun ( void  )

overridevirtual

Called when entering a new run.

Reimplemented from Module.

Definition at line 114 of file ARICHReconstructorModule.cc.

  {
    m_ana->initialize();
  }

cal2byte()

unsigned int cal2byte ( const int *  buf )

inlineprotected

calculate number of lines (2 bytes) in raw Unpacker

Definition at line 107 of file ARICHUnpackerModule.h.

  {
    return (calbyte(buf) << 8) | calbyte(buf);
  }

calbyte() [1/3]

unsigned int calbyte ( const int *  buf )

inlineprotected

Get calculated byte.

Definition at line 86 of file ARICHRateCalModule.h.

  {
    int shift = (3 - m_ibyte % 4) * 8;
    unsigned int val = 0xff & (buf[m_ibyte / 4] >> shift);
    m_ibyte++;
    return val;
  }

calbyte() [2/3]

unsigned int calbyte ( const int *  buf )

inlineprotected

Read byte with number m_ibyte from the buffer and increase the number by 1.

Parameters

[in]

buf

Buffer.

Definition at line 90 of file ARICHRawUnpackerModule.h.

  {
    int shift = (3 - m_ibyte % 4) * 8;
    unsigned int val = 0xff & (buf[m_ibyte / 4] >> shift);
    m_ibyte++;
    return val;
  }

calbyte() [3/3]

unsigned int calbyte ( const int *  buf )

inlineprotected

calculate number of bytes in raw Unpacker

Definition at line 96 of file ARICHUnpackerModule.h.

  {
    int shift = (3 - m_ibyte % 4) * 8;
    unsigned int val = 0xff & (buf[m_ibyte / 4] >> shift);
    m_ibyte++;
    return val;
  }

calword() [1/3]

unsigned int calword ( const int *  buf )

inlineprotected

Get calculated word.

Definition at line 94 of file ARICHRateCalModule.h.

  {
    return (calbyte(buf) << 24) | (calbyte(buf) << 16)
           | (calbyte(buf) << 8) | calbyte(buf);
  }

calword() [2/3]

unsigned int calword ( const int *  buf )

inlineprotected

Read word (4 bytes) from the buffer and increase the byte number m_ibyte by 4.

Parameters

[in]

buf

Buffer.

Definition at line 98 of file ARICHRawUnpackerModule.h.

  {
    return (calbyte(buf) << 24) | (calbyte(buf) << 16)
           | (calbyte(buf) << 8) | calbyte(buf);
  }

calword() [3/3]

unsigned int calword ( const int *  buf )

inlineprotected

calculate number of words in raw Unpacker

Definition at line 115 of file ARICHUnpackerModule.h.

  {
    return (calbyte(buf) << 24) | (calbyte(buf) << 16)
           | (calbyte(buf) << 8) | calbyte(buf);
  }

CherenkovPhoton()

int CherenkovPhoton ( ROOT::Math::XYZVector  r, ROOT::Math::XYZVector  rh, ROOT::Math::XYZVector &  rf, ROOT::Math::XYZVector &  dirf, double *  refind, double *  z, int  n, int  mirrorID = 0  )

private

Calculates the direction of photon emission.

Given the Cherenkov photon emission point "r" and its position on detector plane "rh" (hit position) this methods calculates the direction "dirf" in which photon was emitted, under the assumption that it was reflected from "mirrorID"-th mirror plate (mirrorID=-1 for no reflection).

Parameters

[in]	r	Vector of photon emission point.
[in]	rh	Photon hit position.
[in]	dirf	Vector of photon emission direction (this is output of method).
[in]	rf	Vector of photon position on aerogel exit.
[in]	refind	Array of layers refractive indices.
[in]	z	Array of z coordinates of borders between layers.
[in]	n	Number of aerogel layers through which photon passes.
[in]	mirrorID	Id of mirror from which the photon was reflected (assumption).

Definition at line 253 of file ARICHReconstruction.cc.

  {
    //
    // Description:
    // The method calculates the direction of the cherenkov photon
    // and intersection with the aerogel exit point
    //
    // Called by: CerenkovAngle
 
 
    // Arguments:
    // Input:
    // dir, r track position
    // rh photon hit
 
 
    // Output:
    // rf aerogel exit from which the  photon was emitted
    // dirf photon direction in aerogel
    static ROOT::Math::XYZVector norm(0, 0, 1); // detector plane normal vector
 
    double dwin    = m_arichgp->getHAPDGeometry().getWinThickness();
    double refractiveInd0 = m_arichgp->getHAPDGeometry().getWinRefIndex();
 
    // iteration is stoped when the difference of photon positions on first aerogel exit
    // between two iterations is smaller than this value.
    const double dmin  = 0.0000001;
    const int    niter = 100; // maximal number of iterations
    ROOT::Math::XYZVector dirw;
    ROOT::Math::XYZVector rh1 = rh - dwin * norm;
 
    std::vector <ROOT::Math::XYZVector > rf0(n + 2);
    //  rf0[0] .. track point
    //  rf0[1] 1. 1st aerogel exit
    //  rf0[n] n.  aerogel exit ...
    std::vector <ROOT::Math::XYZVector > dirf0(n + 2);
    //  dirf0[0] .. 1st aerogel direction
    //  dirf0[1] .. 2nd aerogel direction
    //  dirf0[n] .. direction after aerogels
 
    //  z[0] .. 1st aerogel exit
    //  z[n-1] .. 2nd aerogel exit
    //  z[n]    .. detector position
    //  z[n+1]  .. detector + window
 
    rf0[0] = r;
    rf0[1] = rf;
 
    for (int iter = 0; iter < niter; iter++) {
 
      // direction in the space between aerogels and detector
      // *************************************
      if (iter == 0) dirf0[n] = (rh1 - rf0[1]).Unit();
      else  dirf0[n] = (rh1 - rf0[n]).Unit();
 
      // *************************************
      // n-layers of aerogel // refractiveInd relative refractive index
      for (int a = n - 1; a >= 0 ; a--) {
        double rind = refractiveInd[a] / refractiveInd[a + 1];
        dirf0[a] = Refraction(dirf0[a + 1], rind);
      }
 
      double path = (z[0] - r.Z()) / dirf0[0].Z();
      double x1 = rf0[1].X();
      double y1 = rf0[1].Y();
      for (int a = 0; a < n ; a++) {
        rf0[a + 1] = rf0[a] + dirf0[a] * path;
        path = (z[a + 1] - rf0[a + 1].Z()) / dirf0[a + 1].Z();
      }
 
      Refraction(dirf0[n], norm, refractiveInd0, dirw);
 
      // *************************************
 
      path = dwin / (dirw.Dot(norm));
      rh1 = rh - dirw * path;
 
      double d2 = (rf0[1].X() - x1) * (rf0[1].X() - x1) + (rf0[1].Y() - y1) * (rf0[1].Y() - y1);
 
      if ((d2 < dmin) && (iter)) {
 
        // if mirror hypothesis check if reflection point lies on mirror plate
        if (mirrorID) {
          if (!HitsMirror(rf0[n], dirf0[n], mirrorID)) return -1;
        }
 
        rf = rf0[1];
        dirf = dirf0[0];
        return iter;
      }
    }
    return -1;
  }

correctEmissionPoint()

void correctEmissionPoint	(	int	tileID,
		double	r
	)

Correct mean emission point z position.

Definition at line 768 of file ARICHReconstruction.cc.

  {
 
    double ang = m_tilePars[tileID - 1][0] + m_tilePars[tileID - 1][1] * r;
    m_zaero[0] = m_arichgp->getAerogelPlane().getAerogelZPosition() + m_thickness[0] - ang * 50.;
    m_zaero[1] = m_zaero[0] +  m_thickness[1];
 
  }

deadPalette()

void deadPalette

(

)

Set palette for sector dead-chip map.

Definition at line 214 of file hitMapMaker.cc.

  {
    static Int_t colors[50];
    static Bool_t initialized = kFALSE;
    Double_t Red[3]    = { 1.00, 0.00, 0.00};
    Double_t Green[3]  = { 0.00, 0.00, 0.00};
    Double_t Blue[3]   = { 0.00, 1.00, 1.00};
    Double_t Length[3] = { 0.00, 0.20, 1.00 };
    if (!initialized) {
      Int_t FI = TColor::CreateGradientColorTable(3, Length, Red, Green, Blue, 50);
      for (int i = 0; i < 50; i++) colors[i] = FI + i;
      initialized = kTRUE;
      return;
    }
    gStyle->SetPalette(50, colors);
  }

defineHisto() [1/3]

void defineHisto ( )

overridevirtual

Definition of the histograms.

Reimplemented from HistoModule.

Definition at line 61 of file ARICHDQMModule.cc.

  {
 
    TDirectory* oldDir = gDirectory;
    TDirectory* dirARICHDQM = oldDir->mkdir("ARICH");
    dirARICHDQM->cd();
 
    //Histograms for analysis and statistics
 
    h_chStat = new TH1D("chStat", "Status of channels;Channel serial;Status", 420 * 144, -0.5, 420 * 144 - 0.5);
    h_aeroStat = new TH1D("aeroStat", "Status of aerogels;Aerogel tile serial;Status", 160, -0.5, 160 - 0.5);
    h_chHit = new TH1D("chHit", "Number of hits in each channel;Channel serial;Hits", 420 * 144, -0.5, 420 * 144 - 0.5);
    h_chipHit = new TH1D("chipHit", "Number of hits in each chip;Chip serial;Hits", 420 * 4, -0.5, 420 * 4 - 0.5);
    h_hapdHit = new TH1D("hapdHit", "Number of hits in each HAPD;HAPD serial;Hits", 420, 0.5, 421 - 0.5);
    h_hapdHitPerEvent = new TH2D("hapdHitPerEvent", "Number of hits in each HAPD per Event;HAPD serial;Hits/event", 420, 0.5, 420 + 0.5,
                                 144, -0.5, 143.5);
    h_trackPerEvent = new TH1D("trackPerEvent", "Number of tracks in ARICH per event; # of tracks;Events", 6, -0.5, 5.5);
 
    h_mergerHit = new TH1D("mergerHit", "Number of hits in each merger board;MB ID;Hits", 72, 0.5, 72 + 0.5);
    h_bitsPerMergerNorm = new TH2D("bitsPerMergerNorm", "Normalised number of hits in each bit in each Merger;Bit;MB ID;Hits", 5,
                                   -1 - 0.5,
                                   4 - 0.5, 72, 0.5, 72 + 0.5); // copy of h_bits, normalised to number of connected Hapds and to sum(bit1, bit2)
    h_bitsPerHapdMerger = new TH2D("bitsPerHapdMerger",
                                   "Number of hits in each bit in each Hapd sorted by mergers;Bit;HAPD unsorted;Hits", 5, -1 - 0.5,
                                   4 - 0.5, 432, 1, 432);
 
    h_aerogelHit = new TH1D("aerogelHit", "Number track associated hits in each aerogel tile;Aerogel slot ID;Hits", 125, -0.5,
                            125 - 0.5);
    h_bits = new TH1D("bits", "Number of hits in each bit;Bit;Hits", 5, -1 - 0.5,
                      4 - 0.5); //Bin at -1 is added to set minimum 0 without SetMinimum(0) due to bugs in jsroot on DQM server.
    h_bitsPerChannel = new TH2D("bitsPerChannel", "Number of hits in each bit in each chanenel;Bit;channel #;Hits", 4, - 0.5,
                                4 - 0.5, 420 * 144, -0.5, 420 * 144 - 0.5); // copy of h_bits
    h_hitsPerTrack2D = new TH2D("hitsPerTrack2D", "2D distribution of track associated hits;X[cm];Y[cm];Hits", 230, -115, 115, 230,
                                -115, 115);
    h_tracks2D = new TH2D("tracks2D", "Distribution track positions;X[cm];Y[cm];Tracks", 230, -115, 115, 230, -115, 115);
 
    h_hitsPerEvent = new TH1D("hitsPerEvent", "Number of hit per event;Number of hits;Events", 150, -0.5, 150 - 0.5);
    h_theta = new TH1D("theta", "Cherenkov angle distribution;Angle [rad];Events", 60, 0, M_PI / 6);
    h_hitsPerTrack = new TH1D("hitsPerTrack", "Number of hit per track;Number of hits;Tracks", 150, -0.5, 150 - 0.5);
 
    for (int i = 0; i < 6; i++) {
      h_secTheta[i] = new TH1D(Form("thetaSec%d", i + 1), Form("Cherenkov angle distribution in sector %d;Angle [rad];Events", i + 1),
                               60, 0, M_PI / 6);
      h_secHitsPerTrack[i] = new TH1D(Form("hitsPerTrackSec%d", i + 1),
                                      Form("Number of hit per track in sector %d;Number of hits;Tracks", i + 1), 150, -0.5, 150 - 0.5);
      h_secHapdHit[i] = new TH1D(Form("hapdHit%d", i + 1), Form("Number of hits in each HAPD in sector %d;HAPD serial;Hits", i + 1), 70,
                                 0.5, 71 - 0.5);
    }
 
    h_chDigit = new TH1D("chDigit", "Number of raw digits in each channel;Channel serial;Hits", 420 * 144, -0.5, 420 * 144 - 0.5);
    h_chipDigit = new TH1D("chipDigit", "Number of raw digits in each chip;Chip serial;Hits", 420 * 4, -0.5, 420 * 4 - 0.5);
    h_hapdDigit = new TH1D("hapdDigit", "Number of raw digits in each HAPD;HAPD serial;Hits", 420, 0.5, 421 - 0.5);
 
    h_aerogelHits3D = new TH3D("aerogelHits3D", "Number of track associated hits for each aerogel tile; #phi section; r section", 125,
                               -0.5, 124.5, 20, 0, 20, 20, 0, 20);
    h_mirrorThetaPhi = new TH3D("mirrorThetaPhi",
                                "Cherenkov theta vs Cherenkov phi for mirror reflected photons; mirroID; #phi_{c} [rad]; #theta_{c} [rad]", 18, 0.5, 18.5, 100,
                                -M_PI, M_PI, 100, 0, 0.5);
    h_thetaPhi = new TH2D("thetaPhi", "Cherenkov theta vs phi;#phi [rad];#theta_{c} [rad]", 100, -M_PI, M_PI, 100, 0., 0.5);
 
    h_flashPerAPD = new TH1D("flashPerAPD", "Number of flashes per APD; APD serial; number of flash", 420 * 4, -0.5, 420 * 4 - 0.5);
 
    h_ARICHOccAfterInjLer = new TH1F("ARICHOccInjLER", " ARICHOccInjLER /Time;Time in #mus;Nhits/Time (#mus bins)", 4000, 0, 20000);
    h_ARICHOccAfterInjHer = new TH1F("ARICHOccInjHER", " ARICHOccInjHER /Time;Time in #mus;Nhits/Time (#mus bins)", 4000, 0, 20000);
    h_ARICHEOccAfterInjLer = new TH1F("ARICHEOccInjLER", "ARICHEOccInjLER/Time;Time in #mus;Triggers/Time (#mus bins)", 4000, 0, 20000);
    h_ARICHEOccAfterInjHer = new TH1F("ARICHEOccInjHER", "ARICHEOccInjHER/Time;Time in #mus;Triggers/Time (#mus bins)", 4000, 0, 20000);
 
    //Set the minimum to 0
    h_chDigit->SetMinimum(0);
    h_chipDigit->SetMinimum(0);
    h_hapdDigit->SetMinimum(0);
    h_chHit->SetMinimum(0);
    h_chipHit->SetMinimum(0);
    h_hapdHit->SetMinimum(0);
    h_mergerHit->SetMinimum(0);
    h_bitsPerMergerNorm->SetMinimum(0);
    h_bitsPerHapdMerger->SetMinimum(0);
    h_aerogelHit->SetMinimum(0);
    h_bits->SetMinimum(0);
    h_bitsPerChannel->SetMinimum(0);
    h_hitsPerTrack2D->SetMinimum(0);
    h_tracks2D->SetMinimum(0);
    h_flashPerAPD->SetMinimum(0);
    h_hitsPerEvent->SetMinimum(0);
    h_theta->SetMinimum(0);
    h_hitsPerTrack->SetMinimum(0);
    h_ARICHOccAfterInjLer->SetMinimum(0);
    h_ARICHEOccAfterInjLer->SetMinimum(0);
    h_ARICHOccAfterInjHer->SetMinimum(0);
    h_ARICHEOccAfterInjHer->SetMinimum(0);
 
    for (int i = 0; i < 6; i++) {
      h_secTheta[i]->SetMinimum(0);
      h_secHitsPerTrack[i]->SetMinimum(0);
    }
 
    oldDir->cd();
  }

defineHisto() [2/3]

void defineHisto ( )

overridevirtual

Definition of the histograms.

Reimplemented from HistoModule.

Definition at line 64 of file ARICHRateCalModule.cc.

  {
    /*
    if (m_daqdb.size() > 0) {
      ConfigFile config("slowcontrol");
      PostgreSQLInterface db(config.get("database.host"),
                             config.get("database.dbname"),
                             config.get("database.user"),
                             config.get("database.password"),
                             config.getInt("database.port"));
      DBObject obj = DBObjectLoader::load(db, "arich_th",
                                          StringUtil::replace(m_daqdb, ".", ":"));
      obj.print();
      m_dth = obj.getFloat("dth");
      m_th0 = obj.getFloat("th0");
      m_nrun = obj.getInt("nth");
      db.close();
    }
    */
 
    if (m_dth == 0) {
      B2WARNING("dth is set to 0");
    }
    for (int i = 0; i < 100; i++) {
      h_rate2D[i] = new TH2F(Form("h_rate2D_%d", i), Form("MRG#%d;Channel ID; Vth [mV]", i), 144 * 6, -0.5, -0.5 + 144 * 6,
                             m_nrun, (m_th0 - 0.5 * m_dth) * 1000, (m_th0 + (m_nrun - 0.5)*m_dth) * 1000);
    }
  }

defineHisto() [3/3]

void defineHisto ( )

overridevirtual

Definition of the histograms.

Reimplemented from HistoModule.

Definition at line 45 of file ARICHRawUnpackerModule.cc.

  {
    h_rate_a_all = new TH1F("h_rate_a_all", ";Channel ID", 144 * 6, 0, 144 * 6); //yone
    h_rate_b_all = new TH1F("h_rate_b_all", ";Channel ID", 144 * 6, 0, 144 * 6); //yone
    h_rate_c_all = new TH1F("h_rate_c_all", ";Channel ID", 144 * 6, 0, 144 * 6); //yone
    h_rate_d_all = new TH1F("h_rate_d_all", ";Channel ID", 144 * 6, 0, 144 * 6); //yone
  }

endRun() [1/2]

void endRun ( void  )

overridevirtual

Is called after processing the last event of a run.

Good e.g. for storing stuff, that you want to aggregate over one run.

Reimplemented from Module.

Definition at line 555 of file arichBtestModule.cc.

  {
    B2INFO(" arichBtestModule: End Run !!!");
 
    m_file->Write();
    m_file->Close();
 
    if (m_fp) {
      gzclose(m_fp);
      m_eveList.push_back(m_events);
    }
  }

endRun() [2/2]

void endRun ( void  )

overridevirtual

End-of-run action.

Save run-related stuff, such as statistics.

Reimplemented from HistoModule.

Definition at line 397 of file ARICHDQMModule.cc.

  {
    if (h_theta->GetEntries() < 200) return;
    TF1* f1 = new TF1("arichFitFunc", "gaus(0)+pol1(3)", 0.25, 0.4);
    f1->SetParameters(0.8 * h_theta->GetMaximum(), 0.323, 0.016, 0, 0);
    f1->SetParName(0, "C");
    f1->SetParName(1, "mean");
    f1->SetParName(2, "sigma");
    f1->SetParName(3, "p0");
    f1->SetParName(4, "p1");
    h_theta->Fit(f1, "R");
 
    //Normalise bins in histogram bitsPerMergerNorm
    for (int mergerID = 1; mergerID < 73; ++mergerID) {
      double NHapd = 0;
      for (int febSlot = 1; febSlot < 7; ++febSlot) {
        if (m_arichMergerMap->getModuleID(mergerID, febSlot) > 0) NHapd++;
      }
 
      double bin_value[5];
      for (int i = 1; i <= 5; ++i) {
        // loop over bits and save values
        bin_value[i - 1] = h_bitsPerMergerNorm->GetBinContent(i, mergerID);
      }
 
      for (int i = 1; i <= 5; ++i) {
        // loop over bits again and set bin content
        h_bitsPerMergerNorm->SetBinContent(i, mergerID,
                                           bin_value[i - 1] / (bin_value[3] + bin_value[2]) / NHapd); // normalise with sum of bit1 and bit2, and number of connected HAPDs
      }
    }
 
 
  }

event() [1/12]

void event ( void  )

overridevirtual

Running over all events.

Function is called for each evRunning over all events This means, this function is called very often, and good performance of the code is of strong interest.

Reimplemented from Module.

Definition at line 466 of file arichBtestModule.cc.

  {
 
 
 
 
 
    if (!m_fp) return;
    unsigned int hdr[2];
    EventRec rec;
    BeginRec beginrec;
    BeginRec endrec;
    static char msg[1024];
    int type;
 
    const int sint = sizeof(unsigned int);
    do {
      if (gzread(m_fp, hdr, sint * 2) != 2 * sint  && m_end) {
        B2INFO("[" << m_events << "] End of file: " << *m_runCurrent);
        ++m_runCurrent;
        if (m_runCurrent == m_runList.end()) {
          StoreObjPtr<EventMetaData> eventMetaDataPtr;
          eventMetaDataPtr->setEndOfData();
          return;
        }
        beginRun();
      }
 
    } while (m_end &&  m_runCurrent != m_runList.end());
    if (m_events % 1000 == 0)  {
      time_t m_time;
      time(&m_time);
      B2INFO("neve= [" << m_events << "] in " << (double)(m_time - m_timestart) / 60. << " min (" << int(
               m_time - m_timestart) << "s) from " << *m_runCurrent);
 
    }
    m_events++;
 
    type = hdr[0];
    //   len  = hdr[1];
 
    gzseek(m_fp, 2 * sizeof(unsigned int), SEEK_CUR);
    switch (type) {
      case BEGIN_RECORD_TYPE: {
        gzread(m_fp, &beginrec, sizeof(beginrec));
        time_t t = beginrec.time;
        sprintf(msg, "BeginRec run %u time %s", beginrec.runno, ctime(&t));
        B2INFO(msg);
        break;
      }
      case END_RECORD_TYPE: {
        gzread(m_fp, &endrec, sizeof(endrec));
        time_t t = endrec.time;
        sprintf(msg, "EndRec run %u time %s", endrec.runno, ctime(&t));
        B2INFO(msg);
        break;
      }
      case EVENT_RECORD_TYPE: {
        gzread(m_fp, &rec, sizeof(rec));
        int print = !(rec.evtno % 10000);
        time_t t = rec.time;
        if (print) {
          sprintf(msg, "EventRec run %u evt %u mstime %u, time %s", rec.runno, rec.evtno, rec.mstime, ctime(&t));
          B2INFO(msg);
        }
        /* if you just want to jump to the end */
        int pos = gztell(m_fp);
        /* try to read inside data */
        int buflen = rec.len - sizeof(rec) / sizeof(int);
 
        //if (print) printf ("%d: buflen\n",buflen);
        int n[5] = { 0 };
        int i(0), j(0);
        while (i < buflen && j < 5) {
          n[j] = readdata(m_fp, j, print);
          //      n[j] = skipdata( m_fp );
          i += n[j];
          j++;
        }
        if (gzseek(m_fp, pos + buflen * sizeof(unsigned int), SEEK_SET) < 0) B2ERROR("gzseek returns -1 ");
        break;
      }
      default:
        B2ERROR("IO error unknown record type " << type);
        break;
    }
    if (gzeof(m_fp)) m_end = 1;
  }

event() [2/12]

void event ( void  )

overridevirtual

Event processor.

Convert ARICHSimHits of the event to arichDigits.

Reimplemented from Module.

Definition at line 94 of file ARICHDigitizerModule.cc.

  {
 
    //---------------------------------------------------------------------
    // Convert SimHits one by one to digitizer hits.
    //---------------------------------------------------------------------
 
    // We try to include the effect of opposite-polarity crosstalk among channels
    // on each chip, which depend on number of p.e. on chip
 
    std::map<std::pair<int, int>, int> photoElectrons; // this contains number of photoelectrons falling on each channel
    std::map<std::pair<int, int>, int> chipHits; // this contains number of photoelectrons on each chip
 
    // Loop over all photon hits
    for (const ARICHSimHit& aSimHit : m_ARICHSimHits) {
 
      // check for time window
      double time = aSimHit.getGlobalTime() - m_timeWindowStart;
 
      if (time < 0 || time > m_timeWindow) continue;
 
      ROOT::Math::XYVector locpos = aSimHit.getLocalPosition();
 
      // Get id of module
      int moduleID = aSimHit.getModuleID();
 
      if (m_bdistort) magFieldDistorsion(locpos, moduleID);
 
      // skip if not active
      if (!m_modInfo->isActive(moduleID)) continue;
 
      // Get id number of hit channel
      int chX, chY;
      m_geoPar->getHAPDGeometry().getXYChannel(locpos.X(), locpos.Y(), chX, chY);
 
      if (chX < 0 && chY < 0) continue;
 
      int asicChannel = m_chnMap->getAsicFromXY(chX, chY);
 
 
      // eliminate un-active channels
      if (asicChannel < 0 || !m_chnMask->isActive(moduleID, asicChannel)) continue;
 
      // apply channel dependent QE scale factor
      //double qe_scale =  0.27 / m_maxQE;
      //double qe_scale = m_modInfo->getChannelQE(moduleID, asicChannel) * m_simPar->getColEff() / m_maxQE; // eventually move collection efficiency to here!
      double qe_scale = m_modInfo->getChannelQE(moduleID, asicChannel) / m_maxQE;
 
      if (qe_scale > 1.) B2ERROR("Channel QE is higher than QE applied in SensitiveDetector");
      if (gRandom->Uniform(1.) > qe_scale) continue;
 
      // photon was converted to photoelectron
      chipHits[std::make_pair(moduleID, asicChannel / 36)] += 1;
      photoElectrons[std::make_pair(moduleID, asicChannel)] += 1;
 
    }
 
    // loop over produced photoelectrons. Apply suppression due to the reverse polarization crosstalk
    // among channels on the same chip, and produce hit bitmap (4 bits).
 
    for (std::map<std::pair<int, int>, int>::iterator it = photoElectrons.begin(); it != photoElectrons.end(); ++it) {
 
      std::pair<int, int> modch = it->first;
      double npe = double(it->second);
 
      // reduce efficiency
      npe /= (1.0 + m_simPar->getChipNegativeCrosstalk() * (double(chipHits[std::make_pair(modch.first, modch.second / 36)]) - 1.0));
      if (npe < 1.0 && gRandom->Uniform(1) > npe) continue;
 
      // Make hit bitmap (depends on number of p.e. on channel). For now bitmap is 0001 for single p.e., 0011 for 2 p.e., ...
      // More proper implementation is to be done ...
      uint8_t bitmap = 0;
      for (int i = 0; i < npe; i++) {
        bitmap |= 1 << i;
        if (i == 3) break;
      }
 
      // make new digit!
      m_ARICHDigits.appendNew(modch.first, modch.second, bitmap);
 
    }
 
    //--- if background not overlayed add electronic noise hits
    if (m_bgOverlay) return;
    uint8_t bitmap = 1;
    unsigned nSlots = m_geoPar->getDetectorPlane().getNSlots();
    for (unsigned id = 1; id < nSlots + 1; id++) {
      if (!m_modInfo->isActive(id)) continue;
      int nbkg = gRandom->Poisson(m_bkgLevel);
      for (int i = 0; i < nbkg; i++) {
        m_ARICHDigits.appendNew(id, gRandom->Integer(144), bitmap);
      }
    }
 
  }

event() [3/12]

void event ( void  )

overridevirtual

Event processor.

Reimplemented from HistoModule.

Definition at line 217 of file ARICHDQMModule.cc.

  {
    const ARICHGeoDetectorPlane& arichGeoDec = m_arichGeoConfig->getDetectorPlane();
    const ARICHGeoAerogelPlane& arichGeoAero = m_arichGeoConfig->getAerogelPlane();
 
    setReturnValue(1);
 
    if (m_arichHits.getEntries() < m_minHits || m_arichHits.getEntries() > m_maxHits) { setReturnValue(0); return;}
 
    if (!m_arichLikelihoods.getEntries() && m_arichEvents) { setReturnValue(0); return;}
 
    std::vector<int> apds(420 * 4, 0);
    for (const auto& digit : m_arichDigits) {
      uint8_t bits = digit.getBitmap();
      int moduleID  = digit.getModuleID();
      int channelID = digit.getChannelID();
      int mergerID = m_arichMergerMap->getMergerID(moduleID);
      int febSlot = m_arichMergerMap->getFEBSlot(moduleID);
      unsigned binID = (mergerID - 1) * N_FEB2MERGER + (febSlot);
 
      for (int i = 0; i < 8; i++) {
        if ((bits & (1 << i)) && !(bits & ~(1 << i))) {
          h_bits->Fill(i);
          h_bitsPerChannel->Fill(i, (moduleID - 1) * 144 + channelID);
          h_bitsPerMergerNorm->Fill(i, mergerID);
          h_bitsPerHapdMerger->Fill(i, binID);
        } else if (!bits) {
          h_bits->Fill(8);
          h_bitsPerChannel->Fill(8, (moduleID - 1) * 144 + channelID);
          h_bitsPerMergerNorm->Fill(8, mergerID);
          h_bitsPerHapdMerger->Fill(8, binID);
        }
      }
 
      // fill occupancy histograms for raw data
      h_chDigit  ->Fill((moduleID - 1) * 144 + channelID);
      int chip = (moduleID - 1) * 4   + channelID / 36;
      h_chipDigit->Fill(chip);
      apds[chip] += 1;
      h_hapdDigit->Fill(moduleID);
    }
 
    int iapd = 0;
    for (auto apd : apds) { if (apd > 20) h_flashPerAPD->Fill(iapd); iapd++;}
 
    std::vector<int> hpd(420, 0);
    int nHit = 0;
 
    for (int i = 0; i < m_arichHits.getEntries(); i++) {
 
      int moduleID = m_arichHits[i]->getModule();
      int channelID = m_arichHits[i]->getChannel();
      h_chHit->Fill((moduleID - 1) * 144 + channelID);
      hpd[moduleID - 1]++;
      h_chipHit->Fill((moduleID - 1) * 4 + channelID / 36);
      h_hapdHit->Fill(moduleID);
      if (moduleID > 420) B2INFO("Invalid hapd number " << LogVar("hapd ID", moduleID));
 
      for (int j = 1; j <= 7; j++) {
        int ringStart = (j - 1) * (84 + (j - 2) * 6) / 2 + 1; // The smallest module ID in each ring
        int ringEnd = j * (84 + (j - 1) * 6) / 2; // The biggest module ID in each ring
        if (ringStart <= moduleID && moduleID <= ringEnd) {
          h_secHapdHit[(moduleID - ringStart) / (6 + j)]->Fill((moduleID - ringStart) % (6 + j) + 1 + (ringStart - 1) / 6);
        }
      }
 
      int mergerID = m_arichMergerMap->getMergerID(moduleID);
      h_mergerHit->Fill(mergerID);
      nHit++;
    }
 
    h_hitsPerEvent->Fill(nHit);
    int mmid = 1;
    for (auto hh : hpd) { h_hapdHitPerEvent->Fill(mmid, hh); mmid++;}
 
    int ntrk = 0;
    for (const auto& arichTrack : m_arichTracks) {
 
 
      //Momentum limits are applied
      //if (arichTrack.getPhotons().size() == 0) continue;
      if (arichTrack.getMomentum() > 0.5) ntrk++; // count tracks with momentum larger than 0.5 GeV
      if (arichTrack.getMomentum() < m_momDnLim || arichTrack.getMomentum() > m_momUpLim) continue;
 
      ROOT::Math::XYZVector recPos = arichTrack.getPosition();
      int trSector = 0;
      double dPhi = 0;
      if (recPos.Phi() >= 0) {
        dPhi = recPos.Phi();
      } else {
        dPhi = 2 * M_PI + recPos.Phi();
      }
      while (dPhi > M_PI / 3 && trSector < 5) {
        dPhi -= M_PI / 3;
        trSector++;
      }
      double trR = recPos.Rho();
      double trPhi = 0;
      if (recPos.Phi() >= 0) {
        trPhi = recPos.Phi();
      } else {
        trPhi = 2 * M_PI + recPos.Phi();
      }
 
      int iRing = 0;
      int iAzimuth = 0;
      while (arichGeoAero.getRingRadius(iRing + 1) < trR && iRing < 4) {
        iRing++;
      }
      if (iRing == 0) continue;
      while (arichGeoAero.getRingDPhi(iRing) * (iAzimuth + 1) < trPhi) {
        iAzimuth++;
      }
 
      h_tracks2D->Fill(recPos.X(), recPos.Y());
 
      std::vector<ARICHPhoton> photons = arichTrack.getPhotons();
      for (auto& photon : photons) {
        if (photon.getMirror() == 0) {
          if (trR < 93.) {
            h_thetaPhi->Fill(photon.getPhiCer(), photon.getThetaCer());
            h_theta->Fill(photon.getThetaCer());
          }
          int hitSector = 0;
          double hitPhi = arichGeoDec.getSlotPhi(m_arichHits[photon.getHitID()]->getModule());
          if (hitPhi < 0) hitPhi += 2 * M_PI;
          while (hitPhi > M_PI / 3 && hitSector < 5) {
            hitPhi -= M_PI / 3;
            hitSector++;
          }
          h_secTheta[hitSector]->Fill(photon.getThetaCer());
        } else {
          if (trR > 95.) h_mirrorThetaPhi->Fill(photon.getMirror(), photon.getPhiCer(), photon.getThetaCer());
        }
      }
 
      //Get ARICHLikelihood related to the ARICHTrack
      const ExtHit* extHit = arichTrack.getRelated<ExtHit>();
      const Track* mdstTrack = NULL;
      if (extHit) mdstTrack = extHit->getRelated<Track>();
      const ARICHLikelihood* lkh = NULL;
      if (mdstTrack) lkh = mdstTrack->getRelated<ARICHLikelihood>();
      else lkh = arichTrack.getRelated<ARICHLikelihood>();
 
      if (lkh) {
        if (!lkh->getFlag()) continue; //Fill only when the number of expected photons is more than 0.
        double nphoton = lkh->getDetPhot();
        h_hitsPerTrack->Fill(nphoton);
        h_secHitsPerTrack[trSector]->Fill(nphoton);
        h_hitsPerTrack2D->Fill(recPos.X(), recPos.Y(), nphoton);
        int aeroID = arichGeoAero.getAerogelTileID(recPos.X(), recPos.Y());
        h_aerogelHits3D->Fill(aeroID, (trPhi - arichGeoAero.getRingDPhi(iRing)*iAzimuth) / (arichGeoAero.getRingDPhi(iRing) / 20),
                              (trR - arichGeoAero.getRingRadius(iRing)) / ((arichGeoAero.getRingRadius(iRing + 1) - arichGeoAero.getRingRadius(iRing)) / 20),
                              nphoton);
        h_aerogelHit->Fill(aeroID, nphoton);
      }
    }
 
    h_trackPerEvent->Fill(ntrk);
 
    for (auto& it : m_rawFTSW) {
      B2DEBUG(29, "TTD FTSW : " << std::hex << it.GetTTUtime(0) << " " << it.GetTTCtime(0) << " EvtNr " << it.GetEveNo(0)  << " Type " <<
              (it.GetTTCtimeTRGType(0) & 0xF) << " TimeSincePrev " << it.GetTimeSincePrevTrigger(0) << " TimeSinceInj " <<
              it.GetTimeSinceLastInjection(0) << " IsHER " << it.GetIsHER(0) << " Bunch " << it.GetBunchNumber(0));
      auto difference = it.GetTimeSinceLastInjection(0);
      if (difference != 0x7FFFFFFF) {
        unsigned int nentries = m_arichDigits.getEntries();
        float diff2 = difference / 127.; //  127MHz clock ticks to us, inexact rounding
        if (it.GetIsHER(0)) {
          h_ARICHOccAfterInjHer->Fill(diff2, nentries);
          h_ARICHEOccAfterInjHer->Fill(diff2);
        } else {
          h_ARICHOccAfterInjLer->Fill(diff2, nentries);
          h_ARICHEOccAfterInjLer->Fill(diff2);
        }
      }
    }
 
  }

event() [4/12]

void event ( void  )

overridevirtual

Event processor.

Reimplemented from Module.

Definition at line 69 of file ARICHFillHitsModule.cc.

  {
    StoreObjPtr<EventMetaData> evtMetaData;
 
    StoreArray<ARICHDigit> digits;
    StoreArray<ARICHHit> arichHits;
 
    // calculate number of hits on each apd and on each hapd
    std::vector<uint8_t> apdHits(420 * 4, 0);
    std::vector<uint8_t> hapdHits(420, 0);
    for (const auto& digit : digits) {
      uint8_t bits = digit.getBitmap();
      if (!(bits & m_bitMask)) continue;
 
      int moduleID  = digit.getModuleID();
      if (moduleID > 420 || moduleID < 1) continue;
      moduleID--;
      int channelID = digit.getChannelID();
      if (channelID > 143 || channelID < 0) continue;
      int chipID    = moduleID * 4   + channelID / 36;
      apdHits[chipID]++;
      hapdHits[moduleID]++;
    }
 
    if (!m_fillall) {
      for (const auto& digit : digits) {
        int asicCh = digit.getChannelID();
        int modID = digit.getModuleID();
        if (modID > 420 || modID < 1) continue;
        if (asicCh > 143 || asicCh < 0) continue;
        uint8_t hitBitmap = digit.getBitmap();
        if (!(hitBitmap & m_bitMask)) continue;
 
        // remove hot and dead channels
        if (!m_chnMask->isActive(modID, asicCh)) continue;
 
        int chipID    = (modID - 1) * 4   + asicCh / 36;
 
        if (apdHits[chipID]   > m_maxApdHits) continue;
        if (hapdHits[modID - 1] > m_maxHapdHits) continue;
 
 
        int xCh, yCh;
        if (not m_chnMap->getXYFromAsic(asicCh, xCh, yCh)) {
          B2ERROR("Invalid ARICH hit! This hit will be ignored.");
          continue;
        }
 
        ROOT::Math::XYVector hitpos2D = m_geoPar->getChannelPosition(modID, xCh, yCh);
        ROOT::Math::XYZVector hitpos3D(hitpos2D.X(), hitpos2D.Y(),
                                       m_geoPar->getDetectorZPosition() + m_geoPar->getHAPDGeometry().getWinThickness());
        hitpos3D = m_geoPar->getMasterVolume().pointToGlobal(hitpos3D);
 
        if (m_bcorrect) magFieldCorrection(hitpos3D);
        ARICHHit* hh = arichHits.appendNew(hitpos3D, modID, asicCh);
        hh->addRelationTo(&digit);
      }
    } else {
      for (int imod = 1; imod < 421; imod++) {
        for (int ichn = 0; ichn < 144; ichn++) {
 
          // remove hot and dead channels
          if (!m_chnMask->isActive(imod, ichn)) continue;
          int chipID    = (imod - 1) * 4   + ichn / 36;
 
          if (apdHits[chipID]   > m_maxApdHits) continue;
          if (hapdHits[imod - 1] > m_maxHapdHits) continue;
 
          int xCh, yCh;
          if (not m_chnMap->getXYFromAsic(ichn, xCh, yCh)) {
            B2ERROR("Invalid ARICH hit! This hit will be ignored.");
            continue;
          }
 
          ROOT::Math::XYVector hitpos2D = m_geoPar->getChannelPosition(imod, xCh, yCh);
          ROOT::Math::XYZVector hitpos3D(hitpos2D.X(), hitpos2D.Y(),
                                         m_geoPar->getDetectorZPosition() + m_geoPar->getHAPDGeometry().getWinThickness());
          hitpos3D = m_geoPar->getMasterVolume().pointToGlobal(hitpos3D);
          if (m_bcorrect) magFieldCorrection(hitpos3D);
          arichHits.appendNew(hitpos3D, imod, ichn);
        }
      }
    }
 
  }

event() [5/12]

void event ( void  )

overridevirtual

Event processor.

Reimplemented from Module.

Definition at line 64 of file ARICHMCParticlesModule.cc.

  {
    const Const::EDetector arich = Const::EDetector::ARICH; // arich
    const Const::ChargedStable particleHypothesis = Const::pion;
    const int pdgCode = abs(particleHypothesis.getPDGCode());
 
    for (int itrk = 0; itrk < m_tracks.getEntries(); ++itrk) {
 
      const Track* track = m_tracks[itrk];
      const TrackFitResult* fitResult = track->getTrackFitResultWithClosestMass(particleHypothesis);
      if (!fitResult) {
        B2ERROR("No TrackFitResult for " << particleHypothesis.getPDGCode());
        continue;
      }
 
      const MCParticle* particle = track->getRelated<MCParticle>();
      if (!particle) continue;
 
      RelationVector<ExtHit> extHits = DataStore::getRelationsWithObj<ExtHit>(track);
 
      for (unsigned i = 0; i < extHits.size(); i++) {
        const ExtHit* extHit = extHits[i];
        if (abs(extHit->getPdgCode()) != pdgCode or
            extHit->getDetectorID() != arich or
            extHit->getStatus() != EXT_EXIT or // particles registered at the EXIT of the Al plate
            extHit->getMomentum().Z() < 0.0) continue; // track passes in backward
        if (extHit->getCopyID() == 6789) {
          //MCParticle arichMCP = *particle;
          MCParticle* arichP = m_arichMCPs.appendNew(*particle);
          track->addRelationTo(arichP);
 
        }
 
      }
    }
 
  }

event() [6/12]

void event ( void  )

overridevirtual

Event processor.

Reimplemented from Module.

Definition at line 138 of file ARICHNtupleModule.cc.

  {
 
    StoreObjPtr<EventMetaData> evtMetaData;
    StoreArray<ARICHHit> arichHits;
 
    if (!m_arichTracks.isValid()) return;
 
#ifdef ALIGNMENT_USING_BHABHA
    double E_tot = 0.;
    StoreArray<ECLCluster> eclClusters;
    for (const auto& trk : m_tracks) {
      const ECLCluster* eclTrack = trk.getRelated<ECLCluster>();
      if (eclTrack) {
        if (eclTrack->getEnergy() > 0.1) E_tot += eclTrack->getEnergy();
      }
    }
    for (const auto& cluster : eclClusters) {
      if (!cluster.isNeutral()) continue;
      if (cluster.getHypothesisId() != 5 && cluster.getHypothesisId() != 1)
        continue;
      if (cluster.getEnergy() > 0.1) E_tot += cluster.getEnergy();
    }
#endif
 
    for (const auto& arichTrack : m_arichTracks) {
 
      const ExtHit* extHit = arichTrack.getRelated<ExtHit>();
 
      const Track* mdstTrack = NULL;
      if (extHit) mdstTrack = extHit->getRelated<Track>();
 
      const ARICHAeroHit* aeroHit = arichTrack.getRelated<ARICHAeroHit>();
 
      const ARICHLikelihood* lkh = NULL;
      if (mdstTrack) lkh = mdstTrack->getRelated<ARICHLikelihood>();
      else lkh = arichTrack.getRelated<ARICHLikelihood>();
 
      if (!lkh) continue;
 
      m_arich.clear();
 
      m_arich.evt = evtMetaData->getEvent();
      m_arich.run = evtMetaData->getRun();
      m_arich.exp = evtMetaData->getExperiment();
 
      // set hapd window hit if available
      if (arichTrack.hitsWindow()) {
        m_arich.winHit[0] = arichTrack.windowHitPosition().X();
        m_arich.winHit[1] = arichTrack.windowHitPosition().Y();
      }
 
      if (lkh->getFlag() == 1) m_arich.inAcc = 1;
      m_arich.logL.e = lkh->getLogL(Const::electron);
      m_arich.logL.mu = lkh->getLogL(Const::muon);
      m_arich.logL.pi = lkh->getLogL(Const::pion);
      m_arich.logL.K = lkh->getLogL(Const::kaon);
      m_arich.logL.p = lkh->getLogL(Const::proton);
      m_arich.logL.d = lkh->getLogL(Const::deuteron);
 
      m_arich.expPhot.e = lkh->getExpPhot(Const::electron);
      m_arich.expPhot.mu = lkh->getExpPhot(Const::muon);
      m_arich.expPhot.pi = lkh->getExpPhot(Const::pion);
      m_arich.expPhot.K = lkh->getExpPhot(Const::kaon);
      m_arich.expPhot.p = lkh->getExpPhot(Const::proton);
      m_arich.expPhot.d = lkh->getExpPhot(Const::deuteron);
 
      m_arich.detPhot = lkh->getDetPhot();
 
#ifdef ALIGNMENT_USING_BHABHA
      m_arich.etot = E_tot;
#endif
      m_arich.status = 1;
 
      m_arich.trgtype = 0;
      if (m_arichInfo.getEntries() > 0) {
        auto arichInfo = *m_arichInfo[0];
        m_arich.trgtype = arichInfo.gettrgtype();
      }
 
      const MCParticle* particle = 0;
 
      if (mdstTrack) {
        const TrackFitResult* fitResult = mdstTrack->getTrackFitResultWithClosestMass(Const::pion);
        if (fitResult) {
          m_arich.pValue = fitResult->getPValue();
          ROOT::Math::XYZVector trkPos = fitResult->getPosition();
          m_arich.charge = fitResult->getChargeSign();
          m_arich.z0 = trkPos.Z();
          m_arich.d0 = trkPos.Rho();
          m_arich.nCDC = fitResult->getHitPatternCDC().getNHits();
#ifdef ALIGNMENT_USING_BHABHA
          ROOT::Math::XYZVector trkMom = fitResult->getMomentum();
          const ECLCluster* eclTrack = mdstTrack->getRelated<ECLCluster>();
          if (eclTrack) {
            m_arich.eop = eclTrack->getEnergy() / trkMom.R();
            m_arich.e9e21 = eclTrack->getE9oE21();
          }
#endif
        }
        m_arich.status += 10;
        int fromARICHMCP = 0;
        particle = mdstTrack->getRelated<MCParticle>();
        if (!particle) { particle = mdstTrack->getRelated<MCParticle>("arichMCParticles"); fromARICHMCP = 1;}
        if (particle) {
          m_arich.PDG = particle->getPDG();
          m_arich.primary = particle->getStatus(MCParticle::c_PrimaryParticle);
          m_arich.seen = particle->hasSeenInDetector(Const::ARICH);
          ROOT::Math::XYZVector prodVertex = particle->getProductionVertex();
          m_arich.rhoProd = prodVertex.Rho();
          m_arich.zProd = prodVertex.Z();
          m_arich.phiProd = prodVertex.Phi();
          ROOT::Math::XYZVector decVertex = particle->getDecayVertex();
          m_arich.rhoDec = decVertex.Rho();
          m_arich.zDec = decVertex.Z();
          m_arich.phiDec = decVertex.Phi();
 
          if (!fromARICHMCP) {
            MCParticle* mother = particle->getMother();
            if (mother) m_arich.motherPDG = mother->getPDG();
            std::vector<Belle2::MCParticle*> daughs =  particle->getDaughters();
            for (const auto daugh : daughs) {
              if (daugh == NULL) continue;
              if (daugh->getPDG() == particle->getPDG()) m_arich.scatter = 1;
            }
          }
        }
      }
 
 
      // get reconstructed photons associated with track
      const std::vector<ARICHPhoton>& photons = arichTrack.getPhotons();
      m_arich.nRec = photons.size();
      for (auto it = photons.begin(); it != photons.end(); ++it) {
        ARICHPhoton iph = *it;
        if (iph.getHitID() < arichHits.getEntries()) {
          auto ihit = *arichHits[iph.getHitID()];
          int module  = ihit.getModule();
          int channel = ihit.getChannel();
          int chid = 1 << 16 | module << 8 | channel;
          iph.setHitID(chid);
        }
        m_arich.photons.push_back(iph);
      }
 
      ROOT::Math::XYZVector recPos = arichTrack.getPosition();
      m_arich.recHit.x = recPos.X();
      m_arich.recHit.y = recPos.Y();
      m_arich.recHit.z = recPos.Z();
 
      ROOT::Math::XYZVector recMom = arichTrack.getDirection() * arichTrack.getMomentum();
      m_arich.recHit.p = recMom.R();
      m_arich.recHit.theta = recMom.Theta();
      m_arich.recHit.phi = recMom.Phi();
 
      if (aeroHit) {
        ROOT::Math::XYZVector truePos = aeroHit->getPosition();
        m_arich.mcHit.x = truePos.X();
        m_arich.mcHit.y = truePos.Y();
        m_arich.mcHit.z = truePos.Z();
 
        ROOT::Math::XYZVector trueMom = aeroHit->getMomentum();
        m_arich.mcHit.p = trueMom.R();
        m_arich.mcHit.theta = trueMom.Theta();
        m_arich.mcHit.phi = trueMom.Phi();
 
        m_arich.mcHit.PDG = aeroHit->getPDG();
        m_arich.status += 1000;
 
        if (!particle) {
          particle = aeroHit->getRelated<MCParticle>();
          if (particle) {
            m_arich.PDG = particle->getPDG();
            MCParticle* mother = particle->getMother();
            if (mother) m_arich.motherPDG = mother->getPDG();
            m_arich.primary = particle->getStatus(MCParticle::c_PrimaryParticle);
            m_arich.seen = particle->hasSeenInDetector(Const::ARICH);
            ROOT::Math::XYZVector prodVertex = particle->getProductionVertex();
            m_arich.rhoProd = prodVertex.Rho();
            m_arich.zProd = prodVertex.Z();
            m_arich.phiProd = prodVertex.Phi();
            ROOT::Math::XYZVector decVertex = particle->getDecayVertex();
            m_arich.rhoDec = decVertex.Rho();
            m_arich.zDec = decVertex.Z();
            m_arich.phiDec = decVertex.Phi();
 
            std::vector<Belle2::MCParticle*> daughs =  particle->getDaughters();
            for (const auto daugh : daughs) {
              if (daugh->getPDG() == particle->getPDG()) m_arich.scatter = 1;
            }
          }
        }
      }
      m_tree->Fill();
    }
  }

event() [7/12]

void event ( void  )

overridevirtual

Event processor.

Reimplemented from Module.

Definition at line 76 of file ARICHPackerModule.cc.

  {
 
    StoreObjPtr<EventMetaData> evtMetaData;
    StoreArray<ARICHDigit> digits(m_inputDigitsName);
    StoreArray<RawARICH> rawData(m_outputRawDataName);
 
    int nModules = N_MERGERS * N_FEB2MERGER;
 
    std::vector<const ARICHDigit*>* sortedDigits = new std::vector<const ARICHDigit*>[nModules];
    for (const auto& digit : digits) {
      int moduleID = digit.getModuleID();
      unsigned mergerID = m_mergerMap->getMergerID(moduleID);
      if (!mergerID) { B2WARNING("No module2merger mapping for module ID: " <<  moduleID << "; Digit will not be packed!"); continue;}
      if (!m_copperMap->getCopperID(mergerID)) { B2WARNING("No merger2copper mapping for merger ID: " <<  mergerID << "; Digit will not be packed!"); continue;}
      sortedDigits[moduleID - 1].push_back(&digit);
    }
 
    int buffer[4][ARICH_BUFFER_NWORDS];
 
    for (const auto& copperID : m_copperMap->getCopperIDs()) {
 
      int bufferSize[4]  = {0, 0, 0, 0};
      for (int finesse = 0; finesse < 4; finesse++) {
 
        unsigned ibyte = 0;
 
        for (int j = 0; j < ARICH_BUFFER_NWORDS; j++) {
          buffer[finesse][j] = 0;
        }
 
        auto* buf = buffer[finesse];
 
        // get corresponding merger ID
        unsigned mergerID = m_copperMap->getMergerID(copperID, finesse);
        unsigned mergerSN = m_mergerMap->getMergerSN(mergerID);
        if (!mergerID) continue;
 
        ARICHRawHeader mergerHead;
        unsigned dataFormat = m_nonSuppressed + 1; // production data -> TODO: use enum
        mergerHead.type = dataFormat;
        mergerHead.version = m_version;
        mergerHead.mergerID = mergerSN;
        mergerHead.FEBSlot = 0;
        mergerHead.trigger = evtMetaData->getEvent();
        // mergerHead.length dont forget
 
        ibyte += ARICHRAW_HEADER_SIZE;
 
        int nboards = N_FEB2MERGER;
 
        for (int k = 0; k < nboards; k++) {
 
          int moduleID = m_mergerMap->getModuleID(mergerID, k + 1);
          if (moduleID <= 0) continue;
 
          // FEB header
          ARICHRawHeader FEBHead;
          FEBHead.type = dataFormat;
          FEBHead.version = m_version;
          FEBHead.mergerID = mergerSN;
          FEBHead.FEBSlot = k; // board slots go from 0-5 for now, if firmware is updated to 1-6 add +1 !!
          FEBHead.trigger = evtMetaData->getEvent();
 
          if (m_nonSuppressed) {
            // data length in bytes
            FEBHead.length = 144 + ARICHFEB_HEADER_SIZE; // 144ch * 1 byte
            writeHeader(buf, ibyte, FEBHead);
            ibyte += ARICHFEB_HEADER_SIZE; // leave slot for FEB header (FEB header to be implemented!)
 
            // write data
            for (const auto& digit : sortedDigits[moduleID - 1]) {
              unsigned chn = digit->getChannelID();
              //std::cout << "pack: mod: " << boardID << " ch " << chn<< std::endl;
              unsigned shift = 143 - chn;
              unsigned bitmap = (unsigned)digit->getBitmap();
              buf[(ibyte + shift) / 4] += (bitmap << (3 - (ibyte + shift) % 4) * 8);
            }
            ibyte += 144;
          } else {
            FEBHead.length = ARICHFEB_HEADER_SIZE + sortedDigits[moduleID - 1].size() * 2; // each hit is 2 bytes! channel + bitmap
            writeHeader(buf, ibyte, FEBHead);
            ibyte += ARICHFEB_HEADER_SIZE; // leave slot for FEB header (FEB header to be implemented!)
 
            for (const auto& digit : sortedDigits[moduleID - 1]) {
              unsigned chn = digit->getChannelID();
              //std::cout << "pack: mod: " << boardID << " ch " << chn<< std::endl;
              unsigned shift = (3 - ibyte % 4) * 8;
              buf[ibyte / 4] += (chn << shift);
              ibyte++;
              shift = (3 - ibyte % 4) * 8;
              unsigned bitmap = (unsigned)digit->getBitmap();
              buf[ibyte / 4] += (bitmap << shift);
              ibyte++;
            }
          }
        }
 
        unsigned merg = 0;
        mergerHead.length = ibyte - ARICHRAW_HEADER_SIZE;
        writeHeader(buf, merg, mergerHead);
 
        bufferSize[finesse] = ceil(ibyte / 4.);
 
        if (m_debug) {
          std::cout << "Pack finesse: " << finesse << std::endl;
          for (int i = 0; i < bufferSize[finesse]; i++) {
            std::cout << i << "-th word bitset " << std::bitset<32>(buf[i]) << std::endl;
          }
        }
      }
 
      RawCOPPERPackerInfo info;
      info.exp_num = evtMetaData->getExperiment();
      info.run_subrun_num = (evtMetaData->getRun() << 8) +
                            (evtMetaData->getSubrun() & 0xFF); // run number : 14bits, subrun # : 8bits
      info.eve_num = evtMetaData->getEvent();
      info.node_id = ARICH_ID + copperID;
      info.tt_ctime = 0;
      info.tt_utime = 0;
      info.b2l_ctime = 0;
      info.hslb_crc16_error_bit = 0;
      info.truncation_mask = 0;
      info.type_of_data = 0;
 
      auto* raw = rawData.appendNew();
      raw->PackDetectorBuf(buffer[0], bufferSize[0],
                           buffer[1], bufferSize[1],
                           buffer[2], bufferSize[2],
                           buffer[3], bufferSize[3],
                           info);
 
    }
 
    delete [] sortedDigits;
 
  }

event() [8/12]

void event ( void  )

overridevirtual

Event processor.

Reimplemented from HistoModule.

Definition at line 113 of file ARICHRateCalModule.cc.

  {
    StoreObjPtr<EventMetaData> evtmetadata;
    StoreObjPtr<ARICHInfo> arichinfo;
    if (!arichinfo->getthscan_mode()) return;
    double vth_thscan = arichinfo->getvth_thscan();
 
    //int runno = m_internalmode ? m_run_count : evtmetadata->getRun();
    //int raw_evtno = m_internalmode ? m_evt_count : evtmetadata->getEvent();
    int runno = evtmetadata->getRun();
    //int raw_evtno = evtmetadata->getEvent();
 
    ARICHThParam param(runno, m_dth, m_th0, m_nrun);
    StoreArray<ARICHRawDigit> rawdigits;
    for (auto& rawdigit : rawdigits) {
      const int mrgid = rawdigit.getBoardId();
      //const int nfebs = rawdigit.getNFEBs();
      //B2INFO("MB="<<mrgid<<" nfeb="<<nfebs);
      std::vector<ARICHRawDigit::FEBDigit>& febs(rawdigit.getFEBs());
      for (auto& feb : febs) {
        const int febno = feb.febno;
        std::vector<ARICHRawDigit::FEBDigit::ChannelDigit>& channels(feb());
        for (auto& channel : channels) {
          if (channel.val > 0) {
            //B2INFO("MB="<<mrgid<<" ch="<< channel.chno);
            double vth = m_internalmode ? vth_thscan * 1000 : param.getVth() * 1000 ;
            h_rate2D[mrgid]->Fill(channel.chno + febno * 144, vth);
          }
        }
      }
    }
 
    // evt,run counter for internal thscan mode
    m_evt_count++;
    if (m_evt_count == m_nevents) {
      m_evt_count = 0;
      m_run_count++;
    }
 
  }

event() [9/12]

void event ( void  )

overridevirtual

Event processor.

Reimplemented from HistoModule.

Definition at line 62 of file ARICHRawUnpackerModule.cc.

  {
    StoreArray<RawARICH> rawdata;
    StoreArray<ARICHRawDigit> rawdigits;
    for (auto& raw : rawdata) {
      for (int finesse = 0; finesse < 4; finesse++) {
        const int* buf = (const int*)raw.GetDetectorBuffer(0, finesse);
        int bufSize = raw.GetDetectorNwords(0, finesse);
        if (bufSize < 1) continue;
        m_ibyte = 0;
        // read merger header
        int type = calbyte(buf);
        int ver = calbyte(buf);
        int boardid = calbyte(buf);
        int febno = calbyte(buf);
        unsigned int length_all = calword(buf);
        unsigned int mrg_evtno = calword(buf);
        ARICHRawDigit* rawdigit = rawdigits.appendNew(type, ver, boardid, febno, length_all, mrg_evtno);
        rawdigit->setCopperId(raw.GetNodeID(0));
        rawdigit->setHslbId(finesse);
        while (m_ibyte < length_all) {
          int type_feb = calbyte(buf);
          ver = calbyte(buf);
          boardid = calbyte(buf);
          febno = calbyte(buf);
          unsigned int length = calword(buf);
          int evtno = calword(buf);
          unsigned int ibyte = 0;
          std::stringstream ss;
          ss << "type=" << type_feb << ", ver=" << ver << " "
             << ", boardid=" << boardid << ", febno=" << febno
             << ", length=" << length << ", evtno=" << evtno << " ";
          bool hasHit = false;
          long long feb_trigno = 0;
          for (int i = 0; i < 10; i++) {
            int val = calbyte(buf);
            ibyte++;
            if (i > 1 && i < 6) {
              feb_trigno |= (0xff & val) << (5 - i) * 8;
            }
          }
          ARICHRawDigit::FEBDigit feb;
          if (type_feb == 0x02) {//Raw mode
            int ch = 143;
            //B2INFO("raw mode");
            while (ibyte < length) {
              int val = calbyte(buf);
              if (val != 0) {
                ibyte++;
                ss << "ch# " << ch << "(" << val << ") ";
                hasHit = true;
                if (febno < 0 || febno > 6) {
                  B2ERROR("FEB is bad:" << LogVar("FEB", std::to_string(febno) + " hslb-" + std::to_string(finesse))
                          << LogVar("type", type_feb) << LogVar("ver", ver)
                          << LogVar("boardid", boardid) << LogVar("febno", febno)
                          << LogVar("length", length) << LogVar("evtno", evtno));
                }
                feb.push_back(ch, val);
              }
              ch--;
              if (ch < 0) break;
            }
          } else if (type_feb == 0x01) { // Suppressed mode
            if (length > 144 * 2 + 10) B2FATAL("error " << length);
            //B2INFO("suppreed mode");
            while (ibyte < length) {
              int ch = calbyte(buf);
              ibyte++;
              int val = calbyte(buf);
              ibyte++;
              if (val != 0) {
                ss << "ch# " << ch << "(" << val << ") ";
                hasHit = true;
                if (febno < 0 || febno > 6) {
                  B2ERROR("FEB is bad:" << LogVar("FEB", std::to_string(febno) + " hslb-" + std::to_string(finesse))
                          << LogVar("type", type_feb) << LogVar("ver", ver)
                          << LogVar("boardid", boardid) << LogVar("febno", febno)
                          << LogVar("length", length) << LogVar("evtno", evtno));
                  return;
                }
                feb.push_back(ch, val);
              }
            }
          }
          rawdigit->addFEB(feb, type, ver, boardid, febno, length, evtno, feb_trigno);
          if (m_debug && hasHit) {
            B2INFO(ss.str());
          }
        }
      }
    }
  }

event() [10/12]

void event ( void  )

overridevirtual

Event processor.

Reimplemented from Module.

Definition at line 119 of file ARICHReconstructorModule.cc.

  {
    // using track information form tracking system (mdst Track)
    if (m_inputTrackType == 0) {
 
      Const::EDetector myDetID = Const::EDetector::ARICH; // arich
      Const::ChargedStable hypothesis = Const::pion;
      int pdgCode = abs(hypothesis.getPDGCode());
 
      for (const Track& track : m_Tracks) {
        const TrackFitResult* fitResult = track.getTrackFitResultWithClosestMass(hypothesis);
        if (!fitResult) {
          B2ERROR("No TrackFitResult for " << hypothesis.getPDGCode());
          continue;
        }
 
        const MCParticle* mcParticle = track.getRelated<MCParticle>();
 
        ARICHAeroHit* aeroHit = NULL;
        if (mcParticle) aeroHit = mcParticle->getRelated<ARICHAeroHit>();
 
        RelationVector<ExtHit> extHits = DataStore::getRelationsWithObj<ExtHit>(&track);
        ARICHTrack* arichTrack = NULL;
 
        //const ExtHit* arich1stHit = NULL;
        const ExtHit* arich2ndHit = NULL;
        const ExtHit* arichWinHit = NULL;
 
        for (unsigned i = 0; i < extHits.size(); i++) {
          const ExtHit* extHit = extHits[i];
          if (abs(extHit->getPdgCode()) != pdgCode or
              extHit->getDetectorID() != myDetID or
              extHit->getStatus() != EXT_EXIT or // particles registered at the EXIT of the Al plate
              extHit->getMomentum().Z() < 0.0 or // track passes in backward
              extHit->getCopyID() == 12345) { continue;}
          if (extHit->getCopyID() == 6789) {  arich2ndHit = extHit; continue;}
          arichWinHit = extHit;
        }
 
        if (arich2ndHit) {
          // if aeroHit cannot be found using MCParticle check if it was already related to the extHit (by ARICHRelate module)
          if (!aeroHit) aeroHit = arich2ndHit->getRelated<ARICHAeroHit>();
 
          // make new ARICHTrack
          arichTrack = m_ARICHTracks.appendNew(arich2ndHit);
          if (arichWinHit) arichTrack->setHapdWindowHit(arichWinHit);
        }
 
        // skip if track has no extHit in ARICH
        if (!arichTrack) continue;
        // transform track parameters to ARICH local frame
        m_ana->transformTrackToLocal(*arichTrack, m_align);
        // make new ARICHLikelihood
        ARICHLikelihood* like = m_ARICHLikelihoods.appendNew();
        // calculate and set likelihood values
        m_ana->likelihood2(*arichTrack, m_ARICHHits, *like);
        // make relations
        track.addRelationTo(like);
        arichTrack->addRelationTo(arich2ndHit);
        if (aeroHit) arichTrack->addRelationTo(aeroHit);
 
      } // Tracks loop
    } // input type if
 
 
    // using track information form MC (stored in ARICHAeroHit)
    else {
 
      // Loop over all ARICHAeroHits
      for (const ARICHAeroHit& aeroHit : m_aeroHits) {
 
        // make new ARICHTrack
        ARICHTrack* arichTrack = m_ARICHTracks.appendNew(&aeroHit);
        // smearing of track parameters (to mimic tracking system resolutions)
        m_ana->smearTrack(*arichTrack);
        // transform track parameters to ARICH local frame
        m_ana->transformTrackToLocal(*arichTrack, m_align);
        // make associated ARICHLikelihood
        ARICHLikelihood* like = m_ARICHLikelihoods.appendNew();
        // calculate and set likelihood values
        m_ana->likelihood2(*arichTrack, m_ARICHHits, *like);
        // make relation
        arichTrack->addRelationTo(like);
        arichTrack->addRelationTo(&aeroHit);
      } // for iTrack
 
    }
  }

event() [11/12]

void event ( void  )

overridevirtual

Event processor.

Reimplemented from Module.

Definition at line 60 of file ARICHRelateModule.cc.

  {
    int nHits = m_aeroHits.getEntries();
    B2DEBUG(50, "No. of hits " << nHits);
 
    for (const ARICHAeroHit& aeroHit : m_aeroHits) {
      const MCParticle* particle = DataStore::getRelated<MCParticle>(&aeroHit);
      if (!particle) {
        B2DEBUG(50, "No MCParticle for AeroHit!");
        continue;
      }
 
      // Find the track produced by MCParticle
      const Track* track = DataStore::getRelated<Track>(particle);
      if (!track) continue;
 
      const TrackFitResult* fitResult = track->getTrackFitResultWithClosestMass(Const::pion);
      if (!fitResult) continue;
 
      RelationVector<ExtHit> extHits = DataStore::getRelationsWithObj<ExtHit>(track);
      const Const::EDetector arich = Const::EDetector::ARICH;
      for (const ExtHit& extHit : extHits) {
        if (abs(extHit.getPdgCode()) != Const::pion.getPDGCode() or
            extHit.getDetectorID() != arich or
            extHit.getCopyID() != 6789 or // aerogel Al support plate
            extHit.getStatus() != EXT_EXIT) continue; // particles registered at the EXIT of the Al plate
        extHit.addRelationTo(&aeroHit);
        break;
      }
    }
  }

event() [12/12]

void event ( void  )

overridevirtual

Event processor.

temporary! FEB Slots on merger should be 1-6 (now 0-5). Remove when firmware is updated!

Reimplemented from Module.

Definition at line 90 of file ARICHUnpackerModule.cc.

  {
 
    StoreArray<RawARICH> rawData(m_inputRawDataName);
    StoreArray<ARICHDigit> digits(m_outputDigitsName);
    StoreArray<ARICHRawDigit> rawdigits(m_outputRawDigitsName);
    StoreObjPtr<ARICHInfo> arichinfo(m_outputarichinfoName);
    arichinfo.create();
    StoreObjPtr<EventMetaData> evtMetaData;
 
    digits.clear();
    bool m_pciedata = false;
    int trgtype = 16;
    double vth_thscan = 0.0;
 
    if (m_debug) {
      std::cout << std::endl << "------------------------" << std::endl;
      std::cout << "Run: " << evtMetaData->getRun() << " Event: " << evtMetaData->getEvent()  << std::endl;
      std::cout << "------------------------" << std::endl << std::endl;
    }
 
    unsigned thscan_mode = 0;
    // regular Unpacker mode, fill ARICHDigit
//    if (m_disable_unpacker == 0) {
 
    for (auto& raw : rawData) {
      // Check PCIe40 data or Copper data
      if (raw.GetMaxNumOfCh(0) == 48) { m_pciedata = true; } // Could be 36 or 48
      else if (raw.GetMaxNumOfCh(0) == 4) { m_pciedata = false; }
      else { B2FATAL("ARICHUnpackerModule: Invalid value of GetMaxNumOfCh from raw data: " << LogVar("Number of ch: ", raw.GetMaxNumOfCh(0))); }
 
      for (int finesse = 0; finesse < raw.GetMaxNumOfCh(0); finesse++) {
        const int* buffer = raw.GetDetectorBuffer(0, finesse);
        int bufferSize = raw.GetDetectorNwords(0, finesse);
 
        if (bufferSize < 1)
          continue;
 
        // record the trigger type from the B2L data
        trgtype = raw.GetTRGType(0);
 
        // read merger header
        unsigned ibyte = 0;
        ARICHRawHeader head;
 
        readHeader(buffer, ibyte, head);
 
        if (m_debug > 1) printBits(buffer, bufferSize);
 
        if (m_debug) {
          std::cout << "Merger header" << std::endl;
          head.print();
        }
        //-- RawDigit for Merger info
        int type = (int)head.type;
        int ver = (int)head.version;
        int boardid = (int)head.mergerID;
        int febno = (int)head.FEBSlot;
        unsigned int length_all = (unsigned int)head.length;
        unsigned int mrg_evtno = (unsigned int)head.trigger;
        ARICHRawDigit* rawdigit = rawdigits.appendNew(type, ver, boardid, febno, length_all, mrg_evtno);
 
        if (!m_pciedata) {
          rawdigit->setCopperId(raw.GetNodeID(0));
          rawdigit->setHslbId(finesse);
        } else {
          if (raw.GetNodeID(0) == 0x4000001) {
            rawdigit->setCopperId(raw.GetNodeID(0) + (int)(finesse / 4));
          } else if (raw.GetNodeID(0) == 0x4000002) {
            rawdigit->setCopperId(0x400000A + (int)(finesse / 4));
          } else {
            B2FATAL("ARICHUnpackerModule: Invalid Node ID from readout: " <<  LogVar("NodeID: ", raw.GetNodeID(0)));
          }
 
          rawdigit->setHslbId(finesse % 4);
          rawdigit->setPcieId(raw.GetNodeID(0));
          rawdigit->setPcieChId(finesse);
        }
 
        //-- end of RawDigit for Merger info
 
        // record the ibyte here
        unsigned begin = ibyte;
 
        while (ibyte < head.length) {
 
          // new feb
          ARICHRawHeader febHead;
          readFEHeader(buffer, ibyte, febHead);
          if (febHead.thscan_mode) {thscan_mode++;}
          if (m_debug) febHead.print();
 
          if (/*febHead.type != head.type ||*/ febHead.version != head.version || febHead.mergerID != head.mergerID
                                               || febHead.trigger != head.trigger) {
            B2ERROR("ARICHUnpackerModule: data in FEB header not consistent with data in merger HEADER " << LogVar("FEB ID",
                    (unsigned)febHead.FEBSlot) <<
                    LogVar("merger ID", (unsigned)head.mergerID)); break;
          }
 
          // feb header shift
          ibyte += ARICHFEB_HEADER_SIZE;
          int dataLen = febHead.length - ARICHFEB_HEADER_SIZE;
 
          febHead.FEBSlot += 1; 
 
          unsigned mergID = m_mergerMap->getMergerIDfromSN((unsigned)head.mergerID);
 
          if (mergID == 99) { B2ERROR("ARICHUnpackerModule: unknown merger number. Merger data will be skipped.  " << LogVar("merger ID", mergID) << LogVar("Serial Number", (unsigned)head.mergerID)); break;}
 
          unsigned moduleID = m_mergerMap->getModuleID(mergID, (unsigned)febHead.FEBSlot);
 
          if (!moduleID) { B2ERROR("ARICHUnpackerModule: no merger to FEB mapping. Merger data will be skipped. " << LogVar("merger ID", mergID) << LogVar("Serial Number", (unsigned)head.mergerID) << LogVar("FEB slot", (unsigned)febHead.FEBSlot)); break;}
 
          // read data
          if (m_debug) std::cout << "Hit channels: " << std::endl;
          if (febHead.type == 1) {
            for (int i = 0; i < dataLen / 2; i++) {
              int shift = (3 - ibyte % 4) * 8;
              uint8_t asicCh = buffer[ibyte / 4] >> shift;
              ibyte++;
              shift = (3 - ibyte % 4) * 8;
              uint8_t hitBitSet = buffer[ibyte / 4] >> shift;
              if (m_debug && hitBitSet) std::cout << "ch: " << (unsigned)asicCh << " " <<  std::bitset<8>(hitBitSet) << std::endl;
              // store digit
              digits.appendNew(moduleID, (unsigned)asicCh, hitBitSet);
              ibyte++;
            }
          } else if (febHead.type == 2) {
            unsigned asicCh = 143;
            for (int i = 0; i < dataLen; i++) {
              int shift = (3 - ibyte % 4) * 8;
              uint8_t hitBitSet = buffer[ibyte / 4] >> shift;
              // store digit if hit
              if (hitBitSet & m_bitMask) {
                digits.appendNew(moduleID, asicCh, hitBitSet);
              }
              asicCh--;
              ibyte++;
            }
          } else B2ERROR("ARICHUnpackerModule: Unknown data type" << LogVar("type", febHead.type));
 
        }
 
        if (ceil(ibyte / 4.) != (unsigned)bufferSize)
          B2WARNING("ARICHUnpackerModule: data buffer size mismatch " <<  LogVar("size from copper", bufferSize) << LogVar("size from merger",
                    ceil(
                      ibyte / 4.)));
 
 
        //-- If thscan_mode detected from header: proceed to fill ARICHRawDigit
        //-- If m_rawmode is set to 1: proceed to fill ARICHRawDigit
        if (thscan_mode == 0 && m_rawmode == 0) continue;
        //-- go back to beginning again for second loop in case of thscan
        m_ibyte = begin;
 
        while (m_ibyte < length_all) {
          ARICHRawHeader febHead;
          readFEHeader(buffer, m_ibyte, febHead);
          int type_feb = febHead.type;
          ver = febHead.version;
          boardid = febHead.mergerID;
          febno = febHead.FEBSlot;
 
          vth_thscan = (febHead.vth * 0.0024) - 1.27;
          unsigned int length = febHead.length;
          int evtno = febHead.trigger;
          unsigned int jbyte = 0;
          std::stringstream ss;
          ss << "type=" << type_feb << ", ver=" << ver << " "
             << ", boardid=" << boardid << ", febno=" << febno
             << ", length=" << length << ", evtno=" << evtno << " ";
          bool hasHit = false;
          long long feb_trigno = 0;
          for (int i = 0; i < 10; i++) {
            int val = calbyte(buffer);
            jbyte++;
            if (i > 1 && i < 6) {
              feb_trigno |= (0xff & val) << (5 - i) * 8;
            }
          }
          ARICHRawDigit::FEBDigit feb;
          if (type_feb == 0x02) {//Raw mode
            int ch = 143;
            //B2INFO("raw mode");
            while (jbyte < length) {
              int val = calbyte(buffer);
              if (val != 0) {
                jbyte++;
                ss << "ch# " << ch << "(" << val << ") ";
                hasHit = true;
                if (febno < 0 || febno > 6) {
                  B2ERROR("FEB is bad : " << LogVar("FEB no.", febno) << LogVar("hslb", finesse) << LogVar("type", type_feb) << LogVar("ver",
                          ver) << LogVar("boardid", boardid) << LogVar("febno", febno) << LogVar("length", length) << LogVar("evtno", evtno));
                }
                feb.push_back(ch, val);
              }
              ch--;
              if (ch < 0) break;
            }
          } else if (type_feb == 0x01) { // Suppressed mode
            // The below line is commented since it sometimes causes problem during processing threshold scan data.
            // No harm to comment this line since it is only utilized for threshold scan data.
            //if (length > 144 * 2 + 10) B2FATAL("error " << LogVar("length", length));
            //B2INFO("suppreed mode");
            while (jbyte < length) {
              int ch = calbyte(buffer);
              jbyte++;
              int val = calbyte(buffer);
              jbyte++;
              if (val != 0) {
                ss << "ch# " << ch << "(" << val << ") ";
                hasHit = true;
                if (febno < 0 || febno > 6) {
                  B2ERROR("FEB is bad : " << LogVar("FEB no.", febno) << LogVar("hslb", finesse) << LogVar("type", type_feb) << LogVar("ver",
                          ver) << LogVar("boardid", boardid) << LogVar("febno", febno) << LogVar("length", length) << LogVar("evtno", evtno));
                  return;
                }
                feb.push_back(ch, val);
              }
            }
          }
          rawdigit->addFEB(feb, type, ver, boardid, febno, length, evtno, feb_trigno);
          if (m_debug && hasHit) {
            B2INFO(ss.str());
          }
        }
 
 
 
      }
    } // end of rawData loop
 
//    } // end of regular unpacker
    /*
        // RawUnpacker mode, fill ARICHRawDigit
        if (m_rawmode == 1) {
          for (auto& raw : rawData) {
            for (int finesse = 0; finesse < 4; finesse++) {
              const int* buf = (const int*)raw.GetDetectorBuffer(0, finesse);
              int bufSize = raw.GetDetectorNwords(0, finesse);
              if (bufSize < 1) continue;
              m_ibyte = 0;
              // read merger header
              int type = calbyte(buf);
              int ver = calbyte(buf);
              int boardid = calbyte(buf);
              int febno = calbyte(buf);
              unsigned int length_all = calword(buf);
              unsigned int mrg_evtno = calword(buf);
              ARICHRawDigit* rawdigit = rawdigits.appendNew(type, ver, boardid, febno, length_all, mrg_evtno);
              rawdigit->setCopperId(raw.GetNodeID(0));
              rawdigit->setHslbId(finesse);
              int nfebs = 0;
    //--done
              while (m_ibyte < length_all) {
                int type_feb = calbyte(buf);
                ver = calbyte(buf);
                boardid = calbyte(buf);
                febno = calbyte(buf);
 
                // first line: vth value
                unsigned int vth_int = cal2byte(buf);
                if (vth_int > 0) { vth_thscan = (vth_int * 0.0024) - 1.27; }
                // second line: length
                unsigned int length = cal2byte(buf);
                int evtno = calword(buf);
                unsigned int ibyte = 0;
                std::stringstream ss;
                ss << "type=" << type_feb << ", ver=" << ver << " "
                   << ", boardid=" << boardid << ", febno=" << febno
                   << ", length=" << length << ", evtno=" << evtno << " ";
                bool hasHit = false;
                long long feb_trigno = 0;
                for (int i = 0; i < 10; i++) {
                  int val = calbyte(buf);
                  ibyte++;
                  if (i < 6) {
                    feb_trigno |= (0xff & val) << (5 - i) * 8;
                  }
                }
                ARICHRawDigit::FEBDigit feb;
                nfebs++;
                if (type_feb == 0x02) {//Raw mode
                  int ch = 143;
                  //B2INFO("raw mode");
                  while (ibyte < length) {
                    int val = calbyte(buf);
                    if (val != 0) {
                      ibyte++;
                      ss << "ch# " << ch << "(" << val << ") ";
                      hasHit = true;
                      if (febno < 0 || febno > 6) {
                        B2ERROR("FEB is bad : " << LogVar("FEB no.", febno) << LogVar("hslb", finesse) << LogVar("type", type_feb) << LogVar("ver",
                                ver) << LogVar("boardid", boardid) << LogVar("febno", febno) << LogVar("length", length) << LogVar("evtno", evtno));
                      }
                      feb.push_back(ch, val);
                    }
                    ch--;
                    if (ch < 0) break;
                  }
                } else if (type_feb == 0x01) { // Suppressed mode
                  // The below line is commented since it sometimes causes problem during processing threshold scan data.
                  // No harm to comment this line since it is only utilized for threshold scan data.
                  //if (length > 144 * 2 + 10) B2FATAL("error " << LogVar("length", length));
                  //B2INFO("suppreed mode");
                  while (ibyte < length) {
                    int ch = calbyte(buf);
                    ibyte++;
                    int val = calbyte(buf);
                    ibyte++;
                    if (val != 0) {
                      ss << "ch# " << ch << "(" << val << ") ";
                      hasHit = true;
                      if (febno < 0 || febno > 6) {
                        B2ERROR("FEB is bad : " << LogVar("FEB no.", febno) << LogVar("hslb", finesse) << LogVar("type", type_feb) << LogVar("ver",
                                ver) << LogVar("boardid", boardid) << LogVar("febno", febno) << LogVar("length", length) << LogVar("evtno", evtno));
                        return;
                      }
                      feb.push_back(ch, val);
                    }
                  }
                }
                rawdigit->addFEB(feb, type, ver, boardid, febno, length, evtno, feb_trigno);
                if (m_debug && hasHit) {
                  B2INFO(ss.str());
                }
              }
            }
          }
 
        } // end of raw unpacker
    */
    arichinfo->settrgtype(trgtype);
    arichinfo->setpciedata(m_pciedata);
    if (vth_thscan > -1.27) { arichinfo->setvth_thscan(vth_thscan); }
    arichinfo->setntrack(0);
    arichinfo->setnexthit(0);
    arichinfo->setnhit(0);
    if (thscan_mode > 0 || m_rawmode != 0)
    { arichinfo->setthscan_mode(true); }
    else
    { arichinfo->setthscan_mode(false); }
 
  }

FastTracking()

ROOT::Math::XYZVector FastTracking ( ROOT::Math::XYZVector  dirf, ROOT::Math::XYZVector  r, double *  refind, double *  z, int  n, int  opt  )

private

Calculates the intersection of the Cherenkov photon emitted from point "r" in "dirf" direction with the detector plane.

Parameters

[in]	r	Vector of photon emission point.
[in]	dirf	Direction of photon emission.
[in]	refind	Array of layers refractive indices.
[in]	z	Array of z coordinates of borders between layers.
[in]	n	Number of aerogel layers through which photon passes.
[in]	opt	Parameter can be set to 1 to return empty TVector3 in case of errors

Definition at line 151 of file ARICHReconstruction.cc.

  {
    //
    // Description:
    // The method calculates the intersection  of the cherenkov photon
    // with the detector plane
 
    //  z[n+1]
    //  z[0] .. 1st aerogel exit
    //  z[n-1] .. 2nd aerogel exit
 
    double angmir = 0; int section[2] = {0, 0};
 
    unsigned tileID = m_arichgp->getAerogelPlane().getAerogelTileID(r.X(), r.Y());
 
    if (tileID == 0 && opt == 1) return ROOT::Math::XYZVector();
 
    int nmir = m_arichgp->getMirrors().getNMirrors();
    if (nmir > 0) {
      double dangle = 2 * M_PI / nmir;
      angmir = m_arichgp->getMirrors().getStartAngle() - dangle / 2.;
 
      double trkangle = r.Phi() - angmir;
      if (trkangle < 0) trkangle += 2 * M_PI;
      if (trkangle > 2 * M_PI) trkangle -= 2 * M_PI;
 
      section[1]  = int(trkangle / dangle) + 1;
    }
 
    bool reflok = false; bool refl = false;
    double path = (z[0] - r.Z()) / dirf.Z();
    r   += dirf * path;
    for (int a = 1; a <= n + 1 ; a++) {
      double rind = refractiveInd[a] / refractiveInd[a - 1];
      dirf = Refraction(dirf, rind);
      if (dirf.R() == 0) return ROOT::Math::XYZVector();
      path = (z[a] - r.Z()) / dirf.Z();
      ROOT::Math::XYZVector r0 = r;
      if (a == n && opt == 1) {
        if (m_arichgp->getAerogelPlane().getAerogelTileID(r.X(), r.Y()) != tileID) return ROOT::Math::XYZVector();
      }
      r += dirf * path;
      ROOT::Math::XYVector rxy(r);
      // check for possible reflections
      if (a != n || nmir == 0) continue;
      double angle = rxy.Phi() - angmir;
      if (angle < 0) angle += 2 * M_PI;
      if (angle > 2 * M_PI) angle -= 2 * M_PI;
      double dangle = 2 * M_PI / nmir;
      section[0] = int(angle / dangle) + 1;
      if (r.R() > (r - 2 * m_mirrorPoints[section[0] - 1]).R()) {
        refl = true;
        int nrefl = 2;
        if (section[0] == section[1]) nrefl = 1;
        for (int k = 0; k < nrefl; k++) {
          if (!HitsMirror(r0, dirf, section[k])) continue;
          ROOT::Math::XYZVector mirpoint = m_mirrorPoints[section[k] - 1];
          ROOT::Math::XYZVector mirnorm = m_mirrorNorms[section[k] - 1];
          double s = dirf.Dot(mirnorm);
          double s1 = (mirpoint - r0).Dot(mirnorm);
          r = r0 + s1 / s * dirf;
          dirf = dirf - 2 * (dirf.Dot(mirnorm)) * mirnorm;
          path = (z[a] - r.Z()) / dirf.Z();
          r += dirf * path;
          reflok = true;
          break;
        }
      }
    }
 
    if (!reflok && refl) return ROOT::Math::XYZVector();
    return r;
  }

getMirrorNorm()

ROOT::Math::XYZVector getMirrorNorm ( int  mirrorID )

private

Returns normal vector of the mirror plate with id mirrorID.

Definition at line 751 of file ARICHReconstruction.cc.

  {
    if (m_alignMirrors && m_mirrAlign.isValid()) {
 
      ROOT::Math::XYZVector mirnorm = m_arichgp->getMirrors().getNormVector(mirrorID);
 
      VectorUtil::setMagThetaPhi(mirnorm,
                                 mirnorm.R(),
                                 mirnorm.Theta() + m_mirrAlign->getAlignmentElement(mirrorID).getAlpha(),
                                 mirnorm.Phi() + m_mirrAlign->getAlignmentElement(mirrorID).getBeta());
 
      return mirnorm;
 
    }
    return m_arichgp->getMirrors().getNormVector(mirrorID);
  }

getMirrorPoint()

ROOT::Math::XYZVector getMirrorPoint ( int  mirrorID )

private

Returns point on the mirror plate with id mirrorID.

Definition at line 742 of file ARICHReconstruction.cc.

  {
 
    ROOT::Math::XYZVector mirpoint = m_arichgp->getMirrors().getPoint(mirrorID);
    if (m_alignMirrors && m_mirrAlign.isValid()) mirpoint += m_mirrAlign->getAlignmentElement(mirrorID).getTranslation();
    return mirpoint;
 
  }

getTrack()

int getTrack ( int  mask, ROOT::Math::XYZVector &  r, ROOT::Math::XYZVector &  dir  )

protected

Beamtest Track reconstruction.

Definition at line 318 of file arichBtestModule.cc.

  {
    int retval = 0;
    //const int trgch = 13;
    const double  t0 = 0;// m_tdc[trgch];
    for (int i = 0; i < 4; i++) {
      ARICHTracking* w = &m_mwpc[i];
 
      for (int k = 0; k < 4; k++) mwpc_tdc[i][k]->Fill(m_tdc[w->tdc[k]] - t0);
      mwpc_tdc[i][4]->Fill(m_tdc[w->atdc] - t0);
 
 
      for (int k = 0; k < 2; k++) {  // axis x,y
        w->status[k] = 1;
        w->diff[k] = m_tdc[w->tdc[2 * k]]  - m_tdc[w->tdc[2 * k + 1]];
        w->sum[k] = m_tdc[w->tdc[2 * k + 1]] + m_tdc[w->tdc[2 * k]] - 2 * m_tdc[w->atdc];
        mwpc_sum[i][k]->Fill(w->sum[k]);
        if (w->sum[k] < w->cutll[k] || w->sum[k] > w->cutul[k]) continue;
        mwpc_sum_cut[i][k]->Fill(w->sum[k]);
        w->status[k] = 0;
        w->reco[k] = w->diff[k] * w->slp[k] + w->offset[k];
        w->reco[k] += w->pos[k];
 
        mwpc_diff[i][k]->Fill(w->diff[k]);
      }
 
      w->reco[2] = w->pos[2]; // update z axis
      if (!w->status[0] &&  !w->status[1])  mwpc_xy[i]->Fill(w->reco[0], w->reco[1]);
 
      if (mask & (1 << i)) {
        if (!w->status[0] &&  !w->status[1]) {
          // add point to a fitter
        } else {
          retval = 1;
        }
      }
    }
 
    // replace by fitter
    int ii[4] = {0, 1, 0, 0};
    int ncnt = 0;
    for (int i = 0; i < 4; i++) {
      if (mask & (1 << i)) {
        //B2INFO(ncnt << " " << i << " Mask " << mask);
        ii[ncnt++] = i;
      }
    }
    int i0 = ii[0];
    int i1 = ii[1];
 
    r.SetXYZ(m_mwpc[i0].reco[1], m_mwpc[i0].reco[0], m_mwpc[i0].reco[2]);
    dir.SetXYZ(m_mwpc[i1].reco[1] - m_mwpc[i0].reco[1],
               m_mwpc[i1].reco[0] - m_mwpc[i0].reco[0],
               m_mwpc[i1].reco[2] - m_mwpc[i0].reco[2]);
 
    dir = dir.Unit();
 
// end replace by fitter
    if (dir.Z() != 0) {
      for (int i = 0; i < 4; i++) {
        ARICHTracking* w = &m_mwpc[i];
        double l = (w->reco[2] - r.Z()) / dir.Z() ;
        ROOT::Math::XYZVector rext = r + dir * l;
        if (!w->status[0])  mwpc_residuals[i][0]->Fill(w->reco[0] - rext.Y());
        if (!w->status[1])  mwpc_residuals[i][1]->Fill(w->reco[1] - rext.X());
 
        TAxis* axis =  mwpc_residualsz[i][1]->GetYaxis();
        for (int k = 0; k < axis->GetNbins(); k++) {
          double ll = (w->reco[2] + axis->GetBinCenter(k + 1) - r.Z()) / dir.Z();
          ROOT::Math::XYZVector rextt = r + dir * ll;
          mwpc_residualsz[i][0]->Fill(w->reco[0] - rextt.Y(), axis->GetBinCenter(k + 1));
          mwpc_residualsz[i][1]->Fill(w->reco[1] - rextt.X(), axis->GetBinCenter(k + 1));
 
        }
      }
    }
 
    return retval;
  }

getTrackMeanEmissionPosition()

ROOT::Math::XYZVector getTrackMeanEmissionPosition ( const ARICHTrack &  track, int  iAero  )

private

Returns mean emission position of Cherenkov photons from i-th aerogel layer.

Definition at line 699 of file ARICHReconstruction.cc.

  {
    // Emission length measured from aerogel exit
 
    ROOT::Math::XYZVector dir = track.getDirection();
    if (dir.Z() == 0) return ROOT::Math::XYZVector();
    double d   = m_thickness[iAero] / dir.Z() / m_transmissionLen[iAero];
    double dmean = 1 - d / expm1(d);
    //double dmean = -log((1 + exp(-d)) / 2.);
    double mel = dmean * m_transmissionLen[iAero];
 
    return (getTrackPositionAtZ(track, m_zaero[iAero]) - mel * dir);
  }

getTrackPositionAtZ()

ROOT::Math::XYZVector getTrackPositionAtZ ( const ARICHTrack &  track, double  zout  )

private

Returns track direction at point with z coordinate "zout" (assumes straight track).

Definition at line 713 of file ARICHReconstruction.cc.

  {
    ROOT::Math::XYZVector dir = track.getDirection();
    ROOT::Math::XYZVector pos = track.getPosition();
    if (dir.Z() == 0) return ROOT::Math::XYZVector(0, 0, 0);
    double path = (zout - pos.Z()) / dir.Z();
    return pos + dir * path;
  }

HitsMirror()

bool HitsMirror ( const ROOT::Math::XYZVector &  pos, const ROOT::Math::XYZVector &  dir, int  mirrorID  )

private

Returns true if photon at position pos with direction dir hits mirror plate with ID mirrorID.

Parameters

[in]	pos	Photon position.
[in]	dir	Photon direction.
[in]	mirrorID	ID of mirror plate.

Definition at line 236 of file ARICHReconstruction.cc.

  {
 
    ROOT::Math::XYZVector mirnorm = m_mirrorNorms[mirrorID - 1];
    ROOT::Math::XYZVector mirpoint = m_mirrorPoints[mirrorID - 1];
    ROOT::Math::Rotation3D rot = TransformToFixed(mirnorm);
    ROOT::Math::XYZVector dirTr = rot * dir;
    if (dirTr.Z() < 0) return 0; // comes from outter side
    ROOT::Math::XYZVector posTr =  rot * (pos - mirpoint);
    ROOT::Math::XYZVector pointOnMirr = posTr - (posTr.Z() / dirTr.Z()) * dirTr;
    if (fabs(pointOnMirr.Y()) < m_arichgp->getMirrors().getPlateLength() / 2.
        && fabs(pointOnMirr.X()) < m_arichgp->getMirrors().getPlateWidth() / 2.) return 1;
 
    return 0;
  }

HitVirtualPosition()

ROOT::Math::XYZVector HitVirtualPosition ( const ROOT::Math::XYZVector &  hitpos, int  mirrorID  )

private

Returns the hit virtual position, assuming that it was reflected from mirror.

Parameters

[in]	hitpos	Vector of hit positions.
[in]	mirrorID	Id of mirror from which the photon was reflected.

Definition at line 226 of file ARICHReconstruction.cc.

  {
 
    if (mirrorID == 0) return hitpos;
    ROOT::Math::XYZVector mirpoint = m_mirrorPoints[mirrorID - 1];
    ROOT::Math::XYZVector mirnorm = m_mirrorNorms[mirrorID - 1];
    return hitpos - 2 * ((hitpos - mirpoint).Dot(mirnorm)) * mirnorm;
  }

initialize() [1/13]

void initialize

(

)

Read geometry parameters from xml and initialize class members.

Definition at line 59 of file ARICHReconstruction.cc.

  {
 
    for (const auto& part : Const::chargedStableSet) {
      p_mass[part.getIndex()] = part.getMass();
    }
 
    m_nAerogelLayers = m_arichgp->getAerogelPlane().getNLayers();
    m_thickness[m_nAerogelLayers] = 0;
    for (unsigned int i = 0; i < m_nAerogelLayers; i++) {
      m_refractiveInd[i] = m_arichgp->getAerogelPlane().getLayerRefIndex(i + 1);
      m_anorm[i] = ROOT::Math::XYZVector(0, 0, 1);
      m_thickness[i] = m_arichgp->getAerogelPlane().getLayerThickness(i + 1);
      if (i == 0) m_zaero[i] = m_arichgp->getAerogelPlane().getAerogelZPosition() + m_thickness[i];
      else  m_zaero[i] = m_zaero[i - 1] +  m_thickness[i];
 
      m_transmissionLen[i] = m_arichgp->getAerogelPlane().getLayerTrLength(i + 1) ; // aerogel transmission length;
 
      // measured FOM; aeroMerit is number of detected photons for beam of beta=1 and perpedicular incidence to aerogel tile
      // (corrected for geometrical acceptance). n0[i] is then calculated from
      // aeroMerit[i]=n0[i]*sin2(thc)*transmissionLen[i] * (1 - exp(-thickness[i] / transmissionLen[i])
      m_n0[i] = m_recPars->getAerogelFOM(i) / ((1. - 1. / m_refractiveInd[i] / m_refractiveInd[i]) * m_transmissionLen[i] * (1 - exp(
                                                 -m_thickness[i] / m_transmissionLen[i])));
      m_thickness[m_nAerogelLayers]   += m_thickness[i];
    }
    m_refractiveInd[m_nAerogelLayers  ]   = 1.0;
    m_refractiveInd[m_nAerogelLayers + 1]   = m_arichgp->getHAPDGeometry().getWinRefIndex();
    m_zaero[m_nAerogelLayers  ] = m_arichgp->getDetectorZPosition();
    m_zaero[m_nAerogelLayers + 1] = m_zaero[m_nAerogelLayers] + m_arichgp->getHAPDGeometry().getWinThickness();
 
    if (m_mirrAlign.hasChanged()) {
      m_mirrorNorms.clear();
      m_mirrorPoints.clear();
      for (unsigned i = 1; i < m_arichgp->getMirrors().getNMirrors() + 1; i++) {
        m_mirrorNorms.push_back(getMirrorNorm(i));
        m_mirrorPoints.push_back(getMirrorPoint(i));
      }
    }
 
    if (m_tileAlign) {
      if (m_tileAlign.hasChanged()) {
        for (int iTile = 1; iTile < 125; iTile++) {
          m_tilePars[iTile - 1][0] = m_tileAlign->getAlignmentElement(iTile).getAlpha();
          m_tilePars[iTile - 1][1] = m_tileAlign->getAlignmentElement(iTile).getBeta();
        }
      }
    }
  }

initialize() [2/13]

void initialize ( void  )

overridevirtual

Initialize the Module.

Function is called only once at the beginning of your job at the beginning of the corresponding module. Things that can be done here, should be done here, as it is relatively cheap in terms of CPU time.

Reimplemented from Module.

Definition at line 93 of file arichBtestModule.cc.

  {
    B2INFO("arichBtestModule::initialize()");
    StoreArray<ARICHAeroHit> aeroHits;
    StoreArray<ARICHDigit> digits;
    aeroHits.registerInDataStore();
    digits.registerInDataStore();
 
    m_file = new TFile(m_outFile.c_str(), "RECREATE");
    m_tuple = new TNtuple("nt", "Btest data hits", "x:y:xrec:yrec:m:c:gx:gy");
    char hapdName[256];
    for (int i = 0; i < 6; i++) {
      sprintf(hapdName, "hapd%d", i);
      hapd[i] = new TH1F(hapdName, hapdName, 145, -0.5, 144.5);
    }
    char name[256];
 
    for (int i = 0; i < 4; i++) {
      for (int k = 0; k < 2; k++) {
        sprintf(name, "mwpc%d_a%d_", i, k);
        mwpc_diff[i][k]     = new TH1F(strcat(name, "diff"), name, 300, -150, 150);
        sprintf(name, "mwpc%d_a%d_", i, k);
        mwpc_sum[i][k]      = new TH1F(strcat(name, "sum"), name, 200, -0.5, 199.5);
        sprintf(name, "mwpc%d_a%d_", i, k);
        mwpc_sum_cut[i][k]      = new TH1F(strcat(name, "sum_cut"), name, 200, -0.5, 199.5);
        sprintf(name, "mwpc%d_a%d_", i, k);
        mwpc_residuals[i][k] = new TH1F(strcat(name, "resd"), name, 200, -100, 100);
        sprintf(name, "mwpc%d_a%d_", i, k);
        mwpc_residualsz[i][k] = new TH2F(strcat(name, "resd_z"), name, 200, -25, 25, 400, -1000, 1000);
      }
      for (int k = 0; k < 5; k++) {
        sprintf(name, "mwpc%d_a%d_", i, k);
        mwpc_tdc[i][k]     = new TH1F(strcat(name, "tdc"), name, 1024, -1024, 1024);
      }
      sprintf(name, "mwpc%d_", i);
      mwpc_xy[i] = new TH2F(strcat(name, "xy"), name, 120, -30, 30, 120, -30, 30);
    }
    //m_hapdmap = new TGraph("arich/modules/arichBtest/geometry/hapd.map");
    //m_el2pos = new TGraph("arich/modules/arichBtest/geometry/hapdchmap_v0.dat");
 
    /*
     arich::ARICHGeometryPar* _arichgp = arich::ARICHGeometryPar::Instance();
 
     ofstream dout;
     dout.open ("ChannelCenterGlob.txt");
     for (int i=0;i<6;i++){
      for (int k=0;k<144;k++){
        ROOT::Math::XYZVector r = _arichgp->getChannelCenterGlob(i + 1, k);
        dout  << r.X() << " " << r.Y() << endl;
      }
     }
     dout.close();
     */
    time(&m_timestart);
  }

initialize() [3/13]

void initialize ( void  )

overridevirtual

Initialize the Module.

This method is called at the beginning of data processing.

Reimplemented from Module.

Definition at line 73 of file ARICHDigitizerModule.cc.

  {
    // QE at 400nm (3.1eV) applied in SensitiveDetector
    m_maxQE = m_simPar->getQE(3.1);
 
    m_ARICHSimHits.isRequired();
    m_ARICHDigits.registerInDataStore();
 
    m_bgOverlay = false;
    StoreObjPtr<BackgroundInfo> bgInfo("", DataStore::c_Persistent);
    if (bgInfo.isValid()) {
      if (bgInfo->getMethod() == BackgroundInfo::c_Overlay) m_bgOverlay = true;
    }
 
  }

initialize() [4/13]

void initialize ( void  )

overridevirtual

Initialize the Module.

This method is called at the beginning of data processing.

better use isRequired(), but RawFTSW is not in sim

Reimplemented from HistoModule.

Definition at line 160 of file ARICHDQMModule.cc.

  {
    REG_HISTOGRAM
    // Only need the Tracks via relations, but since DQM is only really possible if the
    // Tracks StoreArray is present, make it a requirement instead of optional.
    StoreArray<Track> tracks;
    tracks.isOptional();
 
    m_arichHits.isRequired();
    m_arichDigits.isOptional();
    m_arichTracks.isOptional();
    m_arichLikelihoods.isOptional();
    m_rawFTSW.isOptional(); 
  }

initialize() [5/13]

void initialize ( void  )

overridevirtual

Initialize the Module.

This method is called at the beginning of data processing.

Reimplemented from Module.

Definition at line 57 of file ARICHFillHitsModule.cc.

  {
 
    StoreArray<ARICHDigit> digits;
    digits.isRequired();
 
    StoreArray<ARICHHit> arichHits;
    arichHits.registerInDataStore();
 
    arichHits.registerRelationTo(digits);
  }

initialize() [6/13]

void initialize ( void  )

overridevirtual

Initialize the Module.

This method is called at the beginning of data processing.

Reimplemented from Module.

Definition at line 52 of file ARICHMCParticlesModule.cc.

  {
    // Dependecies check
    m_tracks.isRequired();
    m_extHits.isRequired();
 
    m_arichMCPs.registerInDataStore();
    m_tracks.registerRelationTo(m_arichMCPs);
 
  }

initialize() [7/13]

void initialize ( void  )

overridevirtual

Initialize the Module.

This method is called at the beginning of data processing.

Reimplemented from Module.

Definition at line 69 of file ARICHNtupleModule.cc.

  {
 
    if (m_file) delete m_file;
    m_file = new TFile(m_outputFile.c_str(), "RECREATE");
    if (m_file->IsZombie()) {
      B2FATAL("Couldn't open file '" << m_outputFile << "' for writing!");
      return;
    }
 
    m_tree = new TTree("arich", "ARICH validation ntuple");
 
    m_tree->Branch("evt", &m_arich.evt, "evt/I");
    m_tree->Branch("run", &m_arich.run, "run/I");
    m_tree->Branch("exp", &m_arich.exp, "exp/I");
 
    m_tree->Branch("charge", &m_arich.charge, "charge/S");
    m_tree->Branch("pValue", &m_arich.pValue, "pValue/F");
    m_tree->Branch("d0", &m_arich.z0, "d0/F");
    m_tree->Branch("z0", &m_arich.d0, "z0/F");
 
#ifdef ALIGNMENT_USING_BHABHA
    m_tree->Branch("eop", &m_arich.eop, "eop/F");
    m_tree->Branch("e9e21", &m_arich.e9e21, "e9e21/F");
    m_tree->Branch("etot", &m_arich.etot, "etot/F");
#endif
 
    m_tree->Branch("PDG", &m_arich.PDG, "PDG/I");
    m_tree->Branch("motherPDG", &m_arich.motherPDG, "motherPDG/I");
    m_tree->Branch("primary", &m_arich.primary, "primary/S");
    m_tree->Branch("seen", &m_arich.seen, "seen/S");
    m_tree->Branch("scatter", &m_arich.scatter, "scatter/I");
 
    m_tree->Branch("rhoProd", &m_arich.rhoProd, "rhoProd/F");
    m_tree->Branch("zProd",   &m_arich.zProd,   "zProd/F");
    m_tree->Branch("phiProd", &m_arich.phiProd, "phiProd/F");
    m_tree->Branch("rhoDec", &m_arich.rhoDec, "rhoDec/F");
    m_tree->Branch("zDec",   &m_arich.zDec,   "zDec/F");
    m_tree->Branch("phiDec", &m_arich.phiDec, "phiDec/F");
    m_tree->Branch("status", &m_arich.status, "status/I");
 
    m_tree->Branch("detPhot",  &m_arich.detPhot,  "detPhot/I");
    m_tree->Branch("numBkg",  &m_arich.numBkg,  "e/F:mu:pi:K:p:d");
    m_tree->Branch("expPhot",  &m_arich.expPhot,  "e/F:mu:pi:K:p:d");
    m_tree->Branch("logL",  &m_arich.logL,  "e/F:mu:pi:K:p:d");
 
    m_tree->Branch("recHit",  &m_arich.recHit,  "PDG/I:x/F:y:z:p:theta:phi");
    m_tree->Branch("mcHit",  &m_arich.mcHit,  "PDG/I:x/F:y:z:p:theta:phi");
    m_tree->Branch("winHit",  &m_arich.winHit,  "x/F:y");
    m_tree->Branch("nrec", &m_arich.nRec, "nRec/I");
    m_tree->Branch("nCDC", &m_arich.nCDC, "nCDC/I");
    m_tree->Branch("inAcceptance", &m_arich.inAcc, "inAcc/O");
    m_tree->Branch("photons", "std::vector<Belle2::ARICHPhoton>", &m_arich.photons);
 
    // required input
    m_arichTracks.isRequired();
    m_arichLikelihoods.isRequired();
 
    // optional input
    m_tracks.isOptional();
    m_arichMCPs.isOptional();
    m_arichAeroHits.isOptional();
    m_arichInfo.isOptional();
 
    StoreArray<ARICHHit> arichHits;
    arichHits.isRequired();
 
  }

initialize() [8/13]

void initialize ( void  )

overridevirtual

Initialize the Module.

This method is called at the beginning of data processing.

Reimplemented from Module.

Definition at line 66 of file ARICHPackerModule.cc.

  {
 
    StoreArray<ARICHDigit> digits(m_inputDigitsName);
    digits.isRequired();
 
    StoreArray<RawARICH> rawData(m_outputRawDataName);
    rawData.registerInDataStore();
  }

initialize() [9/13]

void initialize ( void  )

overridevirtual

Initialize the Module.

This method is called at the beginning of data processing.

Reimplemented from HistoModule.

Definition at line 93 of file ARICHRateCalModule.cc.

  {
 
    REG_HISTOGRAM
    StoreArray<ARICHRawDigit> rawdata;
    rawdata.isRequired();
    //StoreArray<RawARICH> rawdata;
    //rawdata.isRequired();
  }

initialize() [10/13]

void initialize ( void  )

overridevirtual

Initialize the Module.

This method is called at the beginning of data processing.

Reimplemented from HistoModule.

Definition at line 53 of file ARICHRawUnpackerModule.cc.

  {
    REG_HISTOGRAM;
    StoreArray<RawARICH> rawdata;
    rawdata.isRequired();
    StoreArray<ARICHRawDigit> rawdigit;
    rawdigit.registerInDataStore();
  }

initialize() [11/13]

void initialize ( void  )

overridevirtual

Initialize the Module.

Reimplemented from Module.

Definition at line 74 of file ARICHReconstructorModule.cc.

  {
    // Initialize variables
    if (m_ana) delete m_ana;
    m_ana = new ARICHReconstruction(m_storePhot);
    m_ana->setTrackPositionResolution(m_trackPositionResolution);
    m_ana->setTrackAngleResolution(m_trackAngleResolution);
    m_ana->useMirrorAlignment(m_alignMirrors);
 
    // Input: ARICHDigits
    m_ARICHHits.isRequired();
 
    m_Tracks.isOptional();
    m_ExtHits.isOptional();
    m_aeroHits.isOptional();
 
    if (!m_aeroHits.isOptional()) {
      m_Tracks.isRequired();
      m_ExtHits.isRequired();
    }
 
    StoreArray<MCParticle> MCParticles;
    MCParticles.isOptional();
 
    // Output - log likelihoods
    m_ARICHLikelihoods.registerInDataStore();
 
    m_ARICHTracks.registerInDataStore();
 
    if (m_inputTrackType) m_ARICHTracks.registerRelationTo(m_ARICHLikelihoods);
    else {
      m_ARICHTracks.registerRelationTo(m_ExtHits);
      m_Tracks.registerRelationTo(m_ARICHLikelihoods);
    }
    //m_ARICHTracks.registerInDataStore(DataStore::c_DontWriteOut);
    //m_ARICHTracks.registerRelationTo(m_ARICHLikelihoods, DataStore::c_Event, DataStore::c_DontWriteOut);
    m_ARICHTracks.registerRelationTo(m_aeroHits);
    printModuleParams();
  }

initialize() [12/13]

void initialize ( void  )

overridevirtual

Initialize the Module.

This method is called at the beginning of data processing.

Reimplemented from Module.

Definition at line 49 of file ARICHRelateModule.cc.

  {
    // Dependecies check
    m_mcParticles.isRequired();
    m_mdstTracks.isRequired();
    m_aeroHits.isRequired();
    m_extHits.isRequired();
    m_extHits.registerRelationTo(m_aeroHits);
  }

initialize() [13/13]

void initialize ( void  )

overridevirtual

Initialize the Module.

This method is called at the beginning of data processing.

Reimplemented from Module.

Definition at line 73 of file ARICHUnpackerModule.cc.

  {
 
    StoreArray<RawARICH> rawData(m_inputRawDataName);
    rawData.isRequired();
 
    StoreArray<ARICHDigit> digits(m_outputDigitsName);
    digits.registerInDataStore();
 
    StoreArray<ARICHRawDigit> rawdigits(m_outputRawDigitsName);
    rawdigits.registerInDataStore();
 
    StoreObjPtr<ARICHInfo> arichinfo(m_outputarichinfoName);
    arichinfo.registerInDataStore();
 
  }

InsideDetector()

int InsideDetector ( ROOT::Math::XYZVector  a, int  copyno  )

private

Returns 1 if vector "a" lies on "copyno"-th detector active surface of detector and 0 else.

Definition at line 109 of file ARICHReconstruction.cc.

  {
    if (copyno == -1) return 0;
    ROOT::Math::XYVector origin;
    origin.SetXY(m_arichgp->getDetectorPlane().getSlotR(copyno) * std::cos(m_arichgp->getDetectorPlane().getSlotPhi(copyno)),
                 m_arichgp->getDetectorPlane().getSlotR(copyno) * std::sin(m_arichgp->getDetectorPlane().getSlotPhi(copyno)));
    ROOT::Math::XYVector a2(a);
    double phi = m_arichgp->getDetectorPlane().getSlotPhi(copyno);
    ROOT::Math::XYVector diff = a2 - origin;
    diff.Rotate(-phi);
    const double size = m_arichgp->getHAPDGeometry().getAPDSizeX();
    if (fabs(diff.X()) < size / 2. && fabs(diff.Y()) < size / 2.) {
      int chX, chY;
      m_arichgp->getHAPDGeometry().getXYChannel(diff.X(), diff.Y(), chX, chY);
      if (chX < 0 || chY < 0) return 0;
      int asicChannel = m_chnMap->getAsicFromXY(chX, chY);
      // eliminate un-active channels
      if (asicChannel < 0 || !m_chnMask->isActive(copyno, asicChannel)) return 0;
      return 1;
    }
    return 0;
  }

likelihood2()

int likelihood2	(	ARICHTrack &	arichTrack,
		const StoreArray< ARICHHit > &	arichHits,
		ARICHLikelihood &	arichLikelihood
	)

Computes the value of identity likelihood function for different particle hypotheses.

Definition at line 349 of file ARICHReconstruction.cc.

  {
 
    const unsigned int nPhotonHits = arichHits.getEntries(); // number of detected photons
 
    if (m_nAerogelLayers + 1 > c_noOfAerogels) B2ERROR("ARICHReconstrucion: number of aerogel layers defined in the xml file exceeds "
                                                         << c_noOfAerogels);
 
    double  logL[c_noOfHypotheses] = {0.0};
    double  nSig_w_acc[c_noOfHypotheses][c_noOfAerogels] = { {0.0} }; // expected no. of signal photons, including geometrical acceptance
    double  nSig_wo_acc[c_noOfHypotheses][c_noOfAerogels][20] = { { {0.0} } }; // expected no. of signal photons, without geometrical acceptance, divided in 20 phi bins (used for PDF normalization)
    double  nSig_wo_accInt[c_noOfHypotheses][c_noOfAerogels] = { {0.0} }; // expected no. of signal photons, without geometrical acceptance, integrated over phi
    double  esigi[c_noOfHypotheses] = {0.0}; // expected number of signal photons in hit pixel
    double  thetaCh[c_noOfHypotheses][c_noOfAerogels] = { {0.0} }; // expected Cherenkov angle
 
    // read some geometry parameters
    double padSize = m_arichgp->getHAPDGeometry().getPadSize();
    int nMirSeg = m_arichgp->getMirrors().getNMirrors();
    double angmir  = m_arichgp->getMirrors().getStartAngle();
 
    // Detected photons in cherenkov ring (integrated over 0.1<theta<0.5)
    int nDetPhotons = 0;
 
    double ebgri[c_noOfHypotheses] = {0.0}; // number of expected background photons per pad
    double  nBgr[c_noOfHypotheses] = {0.0}; // total number of expected background photons (in 0.1-0.5 rad ring)
 
    // reconstructed track direction
    ROOT::Math::XYZVector edir = arichTrack.getDirection();
    if (edir.Z() < 0) return 0;
    double momentum = arichTrack.getMomentum();
 
    double thcResolution = m_recPars->getThcResolution(momentum);
    if (thcResolution < 0) thcResolution = 0.01; // happens for spurious tracks with 100 GeV momentum!
 
    double wideGaussFract = (m_recPars->getParameters())[0];
    double wideGaussSigma = (m_recPars->getParameters())[1];
 
    unsigned tileID = m_arichgp->getAerogelPlane().getAerogelTileID(arichTrack.getPosition().X(), arichTrack.getPosition().Y());
    double r = arichTrack.getPosition().Rho();
    if (tileID > 0) correctEmissionPoint(tileID, r);
 
    //------------------------------------------------------
    // Calculate number of expected detected photons (emmited x geometrical acceptance).
    // -----------------------------------------------------
 
    float nphot_scaling = 20.; // number of photons to be traced is (expected number of emitted photons * nphot_scaling)
    int nStep = 5;             // number of steps in one aerogel layer
 
    // loop over all particle hypotheses
    for (int iHyp = 0; iHyp < c_noOfHypotheses; iHyp++) {
 
      // absorbtion factor (scattering)
      double abs = 1;
 
      // loop over aerogel layers
      for (int iAerogel = m_nAerogelLayers - 1; iAerogel >= 0; iAerogel--) {
 
        thetaCh[iHyp][iAerogel] = ExpectedCherenkovAngle(momentum,
                                                         p_mass[iHyp],
                                                         m_refractiveInd[iAerogel]);
 
        // track length in the radiator
        double pathLengthRadiator = arichTrack.getDirection().Dot(m_anorm[iAerogel]);
        if (pathLengthRadiator)  pathLengthRadiator = m_thickness[iAerogel] / pathLengthRadiator;
 
        // step length
        double dxx = pathLengthRadiator / double(nStep);
        // number of photons to be emmited per step (number of expected photons * nphot_scaling)
        double nPhot = m_n0[iAerogel] * sin(thetaCh[iHyp][iAerogel]) * sin(thetaCh[iHyp][iAerogel]) * dxx * nphot_scaling;
        ROOT::Math::XYZVector exit_point = getTrackPositionAtZ(arichTrack, m_zaero[iAerogel]);
 
        // loop over emmision point steps
        for (int iepoint = 0; iepoint < nStep; iepoint++) {
 
          ROOT::Math::XYZVector epoint = exit_point - (0.5 + iepoint) * dxx * edir;
          abs *= exp(-dxx / m_transmissionLen[iAerogel]);
          unsigned int genPhot = nPhot * abs; // number of photons to emmit in current step, including scattering  correction
 
          // loop over emmited "photons"
          for (unsigned int iPhoton = 0; iPhoton < genPhot; iPhoton++) {
            double fi = 2 * M_PI * iPhoton / float(genPhot); // uniformly distributed in phi
            ROOT::Math::XYZVector adirf = setThetaPhi(thetaCh[iHyp][iAerogel], fi); // photon direction in track system
            adirf =  TransformFromFixed(edir) * adirf;  // photon direction in global system
            int ifi = int (fi * 20 / 2. / M_PI); // phi bin
            // track photon from emission point to the detector plane
            ROOT::Math::XYZVector dposition = FastTracking(adirf, epoint, &m_refractiveInd[iAerogel], &m_zaero[iAerogel],
                                                           m_nAerogelLayers - iAerogel, 1);
            if (dposition.R() > 1.0) {nSig_wo_acc[iHyp][iAerogel][ifi] += 1; nSig_wo_accInt[iHyp][iAerogel] += 1;}
            else continue;
            unsigned  copyno =  m_arichgp->getDetectorPlane().pointSlotID(dposition.X(), dposition.Y());
            if (!copyno) continue;
            // check if photon fell on photosensitive area
            if (InsideDetector(dposition, copyno)) nSig_w_acc[iHyp][iAerogel] += 1;
          }
        }
 
        // scale the obtained numbers
        for (int ik = 0; ik < 20; ik++) {
          nSig_wo_acc[iHyp][iAerogel][ik] /= nphot_scaling;
        }
        nSig_w_acc[iHyp][iAerogel] /= nphot_scaling;
        nSig_wo_accInt[iHyp][iAerogel] /= nphot_scaling;
      } // for (unsigned int iAerogel=0;iAerogel<m_nAerogelLayers;iAerogel++)
 
      // sum up contribution from all aerogel layers
      for (unsigned int iAerogel = 0; iAerogel < m_nAerogelLayers; iAerogel++) {
        nSig_w_acc[iHyp][m_nAerogelLayers] += nSig_w_acc[iHyp][iAerogel];
        nSig_wo_accInt[iHyp][m_nAerogelLayers] += nSig_wo_accInt[iHyp][iAerogel];
 
        for (int ik = 0; ik < 20; ik++) {
          nSig_wo_acc[iHyp][m_nAerogelLayers][ik] += nSig_wo_acc[iHyp][iAerogel][ik];
        }
      }
 
      // get number of expected background photons in ring (integrated from 0.1 to 0.5 rad)
      std::vector<double> bkgPars = {momentum / sqrt(p_mass[iHyp]*p_mass[iHyp] + momentum * momentum), double(arichTrack.hitsWindow())};
      nBgr[iHyp] = m_recPars->getExpectedBackgroundHits(bkgPars);
 
    }  // for (int iHyp=0;iHyp < c_noOfHypotheses; iHyp++ )
    //#####################################################
 
    ROOT::Math::XYZVector track_at_detector = getTrackPositionAtZ(arichTrack, m_zaero[m_nAerogelLayers + 1]);
 
    // the id numbers of mirrors from which the photons could possibly reflect are calculated
    int mirrors[3];
    mirrors[0] = 0; // for no reflection
    int refl = 1;
 
    // only if particle track on detector is at radius larger than 850mm (for now hardcoded)
    // possible reflections are taken into account.
    if (track_at_detector.Rho() > 85.0) {
      double trackang = track_at_detector.Phi() - angmir;
      if (trackang < 0) trackang += 2 * M_PI;
      if (trackang > 2 * M_PI) trackang -= 2 * M_PI;
      int section1 = int(trackang * nMirSeg / 2. / M_PI) + 1;
      int section2 = section1 + 1;
      if (section1 == nMirSeg)  section2 = 1;
      mirrors[1] = section1; mirrors[2] = section2;
      refl = 3;
    }
 
    // loop over all detected photon hits
 
    for (unsigned int iPhoton = 0; iPhoton < nPhotonHits; iPhoton++) {
 
      ARICHHit* h = arichHits[iPhoton];
      int modID = h->getModule();
      int channel = h->getChannel();
      ROOT::Math::XYZVector hitpos = m_arichgp->getMasterVolume().pointToLocal(h->getPosition());
      bool bkgAdded = false;
      int nfoo = nDetPhotons;
      for (int iHyp = 0; iHyp < c_noOfHypotheses; iHyp++) { esigi[iHyp] = 0; ebgri[iHyp] = 0;}
 
      bool reflOK = true; // remove window photons from reflected hypothesis
 
      // loop over possible mirror reflections
      for (int mirr = 0; mirr < refl; mirr++) {
 
        if (!reflOK) break; // photon from window so break
 
        // calculate fi_ch for a given track refl
        ROOT::Math::XYZVector virthitpos =  HitVirtualPosition(hitpos, mirrors[mirr]);
 
        // if hit is more than 25cm from the track position on the detector plane, skip it.
        // (not reconstructing hits with irrelevantly large Cherenkov angle)
        if ((track_at_detector - virthitpos).R() > 25.0) continue;
 
        double sigExpArr[c_noOfHypotheses] = {0.0}; // esigi for given mirror hypothesis only
        double th_cer_all[c_noOfAerogels] = {0.0};
        double fi_cer_all[c_noOfAerogels] = {0.0};
 
        double weight[c_noOfHypotheses][c_noOfAerogels] = { {0.0} };
        double weight_sum[c_noOfHypotheses] = {0.0};
        int proc = 0;
        double fi_cer_trk = 0.;
 
        // loop over all aerogel layers
        for (unsigned int iAerogel = 0; iAerogel < m_nAerogelLayers; iAerogel++) {
 
          ROOT::Math::XYZVector initialrf = getTrackPositionAtZ(arichTrack, m_zaero[iAerogel]);
          ROOT::Math::XYZVector epoint = getTrackMeanEmissionPosition(arichTrack, iAerogel);
          ROOT::Math::XYZVector edirr  = arichTrack.getDirection();
          ROOT::Math::XYZVector photonDirection; // calculated photon direction
 
          if (CherenkovPhoton(epoint, virthitpos, initialrf, photonDirection, &m_refractiveInd[iAerogel], &m_zaero[iAerogel],
                              m_nAerogelLayers - iAerogel, mirrors[mirr]) < 0) break;
 
          ROOT::Math::XYZVector dirch = TransformToFixed(edirr) * photonDirection;
          double fi_cer = dirch.Phi();
          double th_cer = dirch.Theta();
 
 
          th_cer_all[iAerogel] = th_cer;
          fi_cer_all[iAerogel] = fi_cer;
          auto deltaPhi = dirch.Phi() - edir.Phi();
          if (deltaPhi > M_PI)
            deltaPhi -= 2 * M_PI;
          if (deltaPhi < -M_PI)
            deltaPhi += 2 * M_PI;
          fi_cer_trk = deltaPhi;
 
          if (mirr == 0 && th_cer < 0.1) reflOK = false;
          // skip photons with irrelevantly large/small Cherenkov angle
          if (th_cer > 0.5 || th_cer < 0.1) continue;
 
          // count photons with 0.1<thc<0.5
          if (nfoo == nDetPhotons) nDetPhotons++;
 
          if (fi_cer < 0) fi_cer += 2 * M_PI;
          double fii = fi_cer;
          if (mirr > 0) {
            double fi_mir = m_mirrorNorms[mirrors[mirr] - 1].Phi();
            fii = 2 * fi_mir - fi_cer - M_PI;
          }
 
 
          // loop over all particle hypotheses
          for (int iHyp = 0; iHyp < c_noOfHypotheses; iHyp++) {
 
            // track a photon from the mean emission point to the detector surface
            ROOT::Math::XYZVector photonDirection1 = setThetaPhi(thetaCh[iHyp][iAerogel], fi_cer);  // particle system
            photonDirection1 = TransformFromFixed(edirr) * photonDirection1;  // global system
            int ifi = int (fi_cer * 20 / 2. / M_PI);
            ROOT::Math::XYZVector detector_position;
 
            detector_position = FastTracking(photonDirection1, epoint, &m_refractiveInd[iAerogel], &m_zaero[iAerogel],
                                             m_nAerogelLayers - iAerogel, 0);
 
            ROOT::Math::XYZVector meanr             = detector_position - epoint;
            double   path              = meanr.R();
            meanr                      = meanr.Unit();
 
            double meanpath = (m_recPars->getParameters())[2];
            if (iAerogel == 1) meanpath = meanpath - m_thickness[iAerogel];
 
            double detector_sigma    = thcResolution * meanpath / meanr.Z();
            double wide_sigma = wideGaussSigma * path / meanr.Z();
            // calculate pad orientation and distance relative to that photon
            double modphi =  m_arichgp->getDetectorPlane().getSlotPhi(modID);
 
            double      pad_fi = fii - modphi;
            double      dx     = (detector_position - hitpos).R();
            double  dr = (track_at_detector - detector_position).R();
 
            if (dr > 0.01) {
              double normalizacija = nSig_wo_acc[iHyp][iAerogel][ifi] * padSize / (0.1 * M_PI * dr * meanr.Z());
              weight[iHyp][iAerogel] = normalizacija;
              weight_sum[iHyp] += weight[iHyp][iAerogel];
              double integralMain = SquareInt(padSize, pad_fi, dx, detector_sigma) / sqrt(2.);
              double integralWide = SquareInt(padSize, pad_fi, dx, wide_sigma) / sqrt(2.);
              // expected number of signal photons in each pixel
              double exp = normalizacija * ((1 - wideGaussFract) * integralMain + wideGaussFract * integralWide);
              esigi[iHyp] += exp;
              sigExpArr[iHyp] += exp;
            } // if (dr>0 && thetaCh[iHyp][iAerogel])
 
          }// for (int iHyp=0;iHyp< c_noOfHypotheses; iHyp++)
          if (iAerogel == m_nAerogelLayers - 1) proc = 1; // successfully processed for all layers
        }// for (unsigned int iAerogel=0; iAerogel<m_nAerogelLayers;iAerogel++)
 
        if (proc) {
          // add background contribution if not yet (add only once)
          if (!bkgAdded) {
            for (int iHyp = 0; iHyp < c_noOfHypotheses; iHyp++) {
              std::vector<double> pars = {momentum / sqrt(p_mass[iHyp]*p_mass[iHyp] + momentum * momentum), double(arichTrack.hitsWindow())};
              ebgri[iHyp] += m_recPars->getBackgroundPerPad(th_cer_all[1], pars);
            }
            bkgAdded = true;
          }
        }
        // create ARICHPhoton if desired
        if (m_storePhot && th_cer_all[1] > 0 && th_cer_all[1] < 0.6) {
          double n_cos_theta_ch[c_noOfHypotheses] = {0.0};
          double phi_ch[c_noOfHypotheses] = {0.0};
          for (int iHyp = 0; iHyp < c_noOfHypotheses; iHyp++) {
            if (weight_sum[iHyp] > 0) {
              for (unsigned int iAerogel = 0; iAerogel < m_nAerogelLayers; iAerogel++) {
                double emission_prob = weight[iHyp][iAerogel] / weight_sum[iHyp];
                n_cos_theta_ch[iHyp] += emission_prob * m_refractiveInd[iAerogel] * cos(th_cer_all[iAerogel]);
                phi_ch[iHyp] += emission_prob * fi_cer_all[iAerogel];
              }
              //std::cout << iHyp << " " <<  n_cos_theta_ch[iHyp] << " " << phi_ch[iHyp] << std::endl;
            } else {
              n_cos_theta_ch[iHyp] = -99999.;
              phi_ch[iHyp] = -99999.;
            }
          }
          ARICHPhoton phot(iPhoton, th_cer_all[1], fi_cer_all[1], mirrors[mirr]); // th_cer of the first aerogel layer assumption is stored
          phot.setBkgExp(ebgri); // store expected number of background hits
          phot.setSigExp(sigExpArr); // store expected number of signal hits
          phot.setPhiCerTrk(fi_cer_trk); // store phi angle in track coordinates
          phot.setNCosThetaCh(n_cos_theta_ch); // store n cos(theta_th) for all particle hypotheses
          phot.setPhiCh(phi_ch); // store phi_ch for all particle hypotheses
          phot.setXY(hitpos.X(), hitpos.Y()); // store x-y hit position
          phot.setModuleID(modID); // store module id
          phot.setChannel(channel); // store channel
          arichTrack.addPhoton(phot);
        }
 
 
      }// for (int mirr = 0; mirr < refl; mirr++)
 
      //******************************************
      // LIKELIHOOD construction
      //*******************************************
 
      for (int iHyp = 0; iHyp < c_noOfHypotheses; iHyp++) {
        double expected = esigi[iHyp] + ebgri[iHyp];
        if (bkgAdded) logL[iHyp] += expected + log(1 - exp(-expected));
      }
 
    } // for (unsigned  int iPhoton=0; iPhoton< nPhotonHits; iPhoton++)
 
    //*********************************************
    // add constant term to the LIKELIHOOD function
    //*********************************************
    double exppho[c_noOfHypotheses] = {0.0};
    for (int iHyp = 0; iHyp < c_noOfHypotheses; iHyp++) {
      exppho[iHyp] = nSig_w_acc[iHyp][m_nAerogelLayers] * (1 - wideGaussFract) + wideGaussFract * 0.7 *
                     nSig_wo_accInt[iHyp][m_nAerogelLayers] + nBgr[iHyp];
      logL[iHyp] -= exppho[iHyp];
      if (isnan(logL[iHyp]) || isinf(logL[iHyp])) {
        B2WARNING("ARICHReconstruction: log likelihood value infinite! Flat background hit probability is " << ebgri[iHyp] << "!");
        logL[iHyp] = 0;
      }
    }
 
    //******************************************
    // store LikeliHOOD info
    //******************************************
 
    int flag = 1;
    if ((thetaCh[0][0] > 0 || thetaCh[0][1] > 0) &&  nSig_w_acc[0][m_nAerogelLayers] == 0) flag = 0;
 
    // set values of ARICHLikelihood
    arichLikelihood.setValues(flag, logL, nDetPhotons, exppho);
 
    return 1;
  }

magFieldCorrection()

void magFieldCorrection

(

ROOT::Math::XYZVector &

hitpos

)

Corrects hit position for distorsion due to non-perpendicular magnetic field component.

Definition at line 155 of file ARICHFillHitsModule.cc.

  {
    ROOT::Math::XYZVector Bfield = BFieldManager::getField(hitpos);
    ROOT::Math::XYZVector shift = m_geoPar->getHAPDGeometry().getPhotocathodeApdDistance() / abs(Bfield.Z()) * Bfield;
    hitpos.SetX(hitpos.X() - shift.X());
    hitpos.SetY(hitpos.Y() - shift.Y());
  }

magFieldDistorsion()

void magFieldDistorsion	(	ROOT::Math::XYVector &	hit,
		int	copyno
	)

Apply correction to hit position due to non-perpendicular component of magnetic field.

Parameters

hit	local position of simhit
copyno	copy number of hapd

Definition at line 190 of file ARICHDigitizerModule.cc.

  {
 
    ROOT::Math::XYVector hitGlob;
    double phi = m_geoPar->getDetectorPlane().getSlotPhi(copyno);
    double r = m_geoPar->getDetectorPlane().getSlotR(copyno);
    double z = m_geoPar->getDetectorZPosition() + m_geoPar->getHAPDGeometry().getWinThickness();
    hitGlob.SetXY(r * std::cos(phi), r * std::sin(phi));
    ROOT::Math::XYVector shift = hit;
    shift.Rotate(phi);
    hitGlob += shift;
    ROOT::Math::XYZVector Bfield = BFieldManager::getField(m_geoPar->getMasterVolume().pointToGlobal(ROOT::Math::XYZVector(hitGlob.X(),
                                                           hitGlob.Y(), z)));
    double cc = m_geoPar->getHAPDGeometry().getPhotocathodeApdDistance() / abs(Bfield.Z());
    shift.SetX(cc * Bfield.X());
    shift.SetY(cc * Bfield.Y());
    shift.Rotate(phi);
    hit.SetX(hit.X() + shift.X());
    hit.SetY(hit.Y() + shift.Y());
  }

mergerClusterHitMap1D()

TH1 * mergerClusterHitMap1D	(	TH1 *	hitMap,
		int	mergerID
	)

Make 1D hit map of specified Merger Board.

Parameters

hitMap	1D hit map of all channels.
mergerID	Merger board identifier [1:72].

Returns

1D hit map of the merger board.

Definition at line 95 of file hitMapMaker.cc.

  {
 
    int m_mergerID = mergerID;
 
    DBObjPtr<ARICHChannelMapping> arichChMap;
    DBObjPtr<ARICHMergerMapping> arichMergerMap;
 
    TH1* m_hitMap = hitMap;
 
    std::vector<int> moduleIDs;
    for (int i = 1; i < 7; i++) {
      moduleIDs.push_back(arichMergerMap->getModuleID(m_mergerID, i));
    }
 
    TH1D* m_mergerHitMap1D = new TH1D("MergerHitMap1D", Form("Hit map in Merger Board %d", m_mergerID), 144 * 6, -0.5, 144 * 6 - 0.5);
    for (int i = 1; i < 7; i++) {
      for (int j = 0; j < 144; j++) {
        int hitsNum = m_hitMap->GetBinContent((moduleIDs[i] - 1) * 144 + i);
        m_mergerHitMap1D->Fill(144 * (i - 1) + j, hitsNum);
      }
    }
    return m_mergerHitMap1D;
  }

mergerClusterHitMap2D()

TCanvas * mergerClusterHitMap2D	(	TH1 *	hitMap,
		int	mergerID
	)

Make display of 6 HAPDs' 2D hit map of the Merger Board.

Parameters

hitMap	1D hit map of all channels.
mergerID	Merger board identifier [1:72].

Returns

Display of 6 HAPDs' 2D hit map.

Definition at line 120 of file hitMapMaker.cc.

  {
    int m_mergerID = mergerID;
 
    DBObjPtr<ARICHChannelMapping> arichChMap;
    DBObjPtr<ARICHMergerMapping> arichMergerMap;
 
    TH1* m_hitMap = hitMap;
 
 
    std::vector<int> moduleIDs;
    for (int i = 1; i < 7; i++) {
      moduleIDs.push_back(arichMergerMap->getModuleID(m_mergerID, i));
    }
 
    TCanvas* m_mergerHitMap = new TCanvas("MergerHitMap", "Hit map in Merger Board", 600, 400);
    m_mergerHitMap->Divide(3, 2);
    for (int i = 1; i < 7; i++) {
      m_mergerHitMap->cd(i);
      moduleHitMap(m_hitMap, moduleIDs[i])->Draw("coloz");
    }
    return m_mergerHitMap;
  }

moduleDeadMap()

TH2 * moduleDeadMap	(	TH1 *	hitMap,
		int	moduleID
	)

Make chip dead/alive map in HAPD view (2*2 chips).

Parameters

hitMap	1D hit map of all channels.
moduleID	Module identifier [1:420].

Returns

2D dead/alive map on HAPD.

Definition at line 62 of file hitMapMaker.cc.

  {
    int m_moduleID = moduleID;
 
    DBObjPtr<ARICHChannelMapping> arichChMap;
 
    TH2* m_moduleDeadMap = new TH2D(Form("HAPDDeadMapMod%d", moduleID), Form("HAPD alive/dead module %d;nChX;nChY", moduleID), 2, 0.5,
                                    2.5,
                                    2, 0.5, 2.5);
    TH1* m_hitMap = hitMap;
 
    int deadCh[2][2] = {};
 
    for (int i = 0; i < 144; i++) {
      int hitsNum = m_hitMap->GetBinContent((m_moduleID - 1) * 144 + i);
      int xChn, yChn;
      arichChMap->getXYFromAsic(i, xChn, yChn);
      if (hitsNum == 0) deadCh[(int)xChn / 6][(int)yChn / 6]++;
    }
 
    for (int j = 0; j < 2; j++) {
      for (int k = 0; k < 2; k++) {
        if (deadCh[j][k] > 18) {
          m_moduleDeadMap->Fill(j + 1, k + 1, 1);
        } else {
          m_moduleDeadMap->Fill(j + 1, k + 1, 10);
        }
      }
    }
 
    return m_moduleDeadMap;
  }

moduleHitMap()

TH2 * moduleHitMap	(	TH1 *	hitMap,
		int	moduleID
	)

Make hit map in HAPD view (12*12 channels).

Parameters

hitMap	1D hit map of all channels.
moduleID	Module identifier [1:420].

Returns

2D hit map on HAPD.

Definition at line 38 of file hitMapMaker.cc.

  {
 
    int m_moduleID = moduleID;
 
    DBObjPtr<ARICHChannelMapping> arichChMap;
 
    TH2* m_moduleHitMap = new TH2D(Form("HAPDHitMapMod%d", moduleID), Form("HAPD hit map module %d;nChX;nChY", moduleID), 12, 0.5, 12.5,
                                   12,
                                   0.5, 12.5);
    TH1* m_hitMap = hitMap;
 
    for (int i = 0; i < 144; i++) {
      int hitsNum = m_hitMap->GetBinContent((m_moduleID - 1) * 144 + i);
      int xChn, yChn;
      arichChMap->getXYFromAsic(i, xChn, yChn);
      if (hitsNum != 0) {
        m_moduleHitMap->Fill(xChn + 1, yChn + 1, hitsNum);
      }
    }
 
    return m_moduleHitMap;
  }

printBits()

void printBits ( const int *  buffer, int  bufferSize  )

private

Unpack raw data given in production format.

Parameters

buffer	raw data buffer
bufferSize	buffer size

Definition at line 529 of file ARICHUnpackerModule.cc.

  {
    for (int i = 0; i < bufferSize; i++) {
      std::cout << i << "-th word bitset: " << std::bitset<32>(*(buffer + i)) << std::endl;
    }
  }

printModuleParams()

void printModuleParams ( )

protected

Print module parameters.

Definition at line 208 of file ARICHReconstructorModule.cc.

  {
    if (m_inputTrackType == 0) { B2DEBUG(100, "ARICHReconstructorModule: track infromation is taken from mdst Tracks.");}
    else  B2DEBUG(100, "ARICHReconstructorModule: track information is taken from MC (ARICHAeroHit).");
  }

readdata()

int readdata ( gzFile  fp, int  rec_id, int  print  )

protected

Read the data from the file (can be compressed)

Definition at line 398 of file arichBtestModule.cc.

  {
 
    unsigned int len, data[10000];
    gzread(fp, &len, sizeof(unsigned int));
    //if (print) printf( "[%3d]   %d: ", len, rec_id );
    gzread(fp, data, sizeof(unsigned int)*len);
 
    ROOT::Math::XYZVector r;
    ROOT::Math::XYZVector dir;
    if (rec_id == 1) {
      readmwpc(data, len);
      int retval = getTrack(*(m_MwpcTrackMask.begin()), r, dir);
      //dir = ROOT::Math::XYZVector(0,0,1);
 
      if (!retval) {
        // global transf, add track to datastore
        StoreArray<ARICHAeroHit> arichAeroHits;
 
        int particleId = 11;// geant4
        dir *= m_beamMomentum * Unit::GeV;
        r *= Unit::mm /*/ CLHEP::mm*/;
        static ARICHBtestGeometryPar* _arichbtgp = ARICHBtestGeometryPar::Instance();
        static ROOT::Math::XYZVector dr =  _arichbtgp->getTrackingShift();
 
        r += dr;
 
        //----------------------------------------
        // Track rotation
        //
        ROOT::Math::Rotation3D rot  =  _arichbtgp->getFrameRotation();
        ROOT::Math::XYZVector  rc   =  _arichbtgp->getRotationCenter();
 
        ROOT::Math::XYZVector rrel  =  rc - rot * rc;
        r = rot * r + rrel;
        dir = rot * dir;
        r.SetX(-r.X()); dir.SetX(-dir.X());
 
        //
        // end track rotation
        //----------------------------------------
        r.SetY(-r.Y());
        dir.SetY(-dir.Y());
        B2DEBUG(50, "-----------> " <<  rc.X() <<  " " << rc.Y() << " " <<   rc.Z() << "::::" << rrel.X() <<  " " << rrel.Y() << " " <<
                rrel.Z()  << " ----> R " <<   r.X() <<  " " << r.Y() << " " <<   r.Z() << " ----> S " <<   dir.X() <<  " " << dir.Y() << " " <<
                dir.Z());
 
        // Add new ARIHCAeroHit to datastore
        arichAeroHits.appendNew(particleId, r, dir);
 
      }
    }
 
    if (rec_id == 2) {
      readhapd(len, data);
    }
 
    for (unsigned int j = 0; j < len; j++) {
      //if( j%8==0 && j!= 0 ) printf( "       " );
      //printf( " %08x", data[j] );
      //if( j%8==7 || j==len-1 ) putchar('\n');
    }
    return len + 1;
  }

readFEHeader()

void readFEHeader ( const int *  buffer, unsigned &  ibyte, ARICHRawHeader &  head  )

private

Read FE header.

Definition at line 483 of file ARICHUnpackerModule.cc.

  {
 
    // read the first line of header
    char line1[4];
    int shift;
    for (int i = 0; i < 4; i++) {
      shift = (3 - ibyte % 4) * 8;
      line1[3 - i] = buffer[ibyte / 4] >> shift;
      ibyte++;
    }
 
    head.type = line1[3];
    head.version = line1[2];
    head.mergerID = line1[1];
    head.FEBSlot = line1[0];
 
    // data length
    unsigned char len[4];
    for (int i = 0; i < 4; i++) {
      shift = (3 - ibyte % 4) * 8;
      len[3 - i] = buffer[ibyte / 4] >> shift;
      ibyte++;
    }
 
    unsigned vth_info = len[3] * 256 + len[2];
    if (vth_info >= 32768) { head.thscan_mode = true; vth_info -= 32768; }
    head.vth = vth_info;
    // This line (16 bits) is actaully not used for data length.
    len[2] = 0;
    len[3] = 0;
    uint32_t* tmp = (uint32_t*)len;
    head.length = *tmp;
 
    // trigger number
    char trg[4];
    for (int i = 0; i < 4; i++) {
      shift = (3 - ibyte % 4) * 8;
      trg[3 - i] = buffer[ibyte / 4] >> shift;
      ibyte++;
    }
    tmp = (uint32_t*)trg;
    head.trigger = *tmp;
 
  }

readhapd()

int readhapd ( unsigned int  len, unsigned int *  data  )

protected

Read the HAPD hits from the data buffer.

Definition at line 268 of file arichBtestModule.cc.

  {
 
    ARICHGeometryPar* _arichgp = ARICHGeometryPar::Instance();
    ARICHBtestGeometryPar* _arichbtgp = ARICHBtestGeometryPar::Instance();
    //-----------------------------------------------------
 
    unsigned int bmask = 0xF;
 
    for (int module = 0; module < 6; module++) {
      int istart = 19 * module;
      int istop = 19 * module + 18;
      for (int i = istart; i < istop; i++) {
        for (int j = 0; j < 8; j++) {
          unsigned int kdata = data[i] & (bmask << (j * 4));
          if (kdata) {
            int channelID = j + 8 * (i - istart);
 
            if (_arichgp->isActive(module,  channelID)) {
 
              hapd[module]->Fill(channelID);
 
              StoreArray<ARICHDigit> arichDigits;
              double globalTime = 0;
 
              double rposx = 0, rposy = 0;
              std::pair<int, int> eposhapd(_arichbtgp->GetHapdElectronicMap(module * 144 + channelID));
              int channel = eposhapd.second;
 
              if ((channel < 108 && channel > 71) || channel < 36) channel = 108 - (int(channel / 6) * 2 + 1) * 6 +
                    channel;
              else channel = 144 - (int((channel - 36) / 6) * 2 + 1) * 6 + channel - 36;
              ROOT::Math::XYVector loc = _arichgp->getChannelCenterLoc(channel);
              if (abs(loc.X()) > 2.3 || abs(loc.Y()) > 2.3) continue;
 
              arichDigits.appendNew(module + 1, channel, globalTime);
 
              ROOT::Math::XYZVector rechit = _arichgp->getChannelCenterGlob(module + 1, channel);
              std::pair<double, double> poshapd(_arichbtgp->GetHapdChannelPosition(module * 144 + channelID));
              m_tuple ->Fill(-poshapd.first, poshapd.second, rechit.X(), rechit.Y(), module, channelID, rposx, rposy);
            }
          }
        }
      }
    }
 
    return len;
  }

readHeader()

void readHeader ( const int *  buffer, unsigned &  ibyte, ARICHRawHeader &  head  )

private

Read Merger header.

Definition at line 435 of file ARICHUnpackerModule.cc.

  {
 
    // read the first line of header
    char line1[4];
    int shift;
    for (int i = 0; i < 4; i++) {
      shift = (3 - ibyte % 4) * 8;
      line1[3 - i] = buffer[ibyte / 4] >> shift;
      ibyte++;
    }
 
    head.type = line1[3];
    head.version = line1[2];
    head.mergerID = line1[1];
    head.FEBSlot = line1[0];
 
    // data length
    unsigned char len[4];
    for (int i = 0; i < 4; i++) {
      shift = (3 - ibyte % 4) * 8;
      len[3 - i] = buffer[ibyte / 4] >> shift;
      ibyte++;
    }
 
    unsigned seu = len[2];
    // This line (16 bits) is actaully not used for data length.
    len[2] = 0;
    len[3] = 0;
    uint32_t* tmp = (uint32_t*)len;
    head.length = *tmp;
 
    for (int i = 0; i < 6; i ++) {
      head.SEU_FEB[i] = (seu & (1 << i)) != 0;
    }
 
    // trigger number
    char trg[4];
    for (int i = 0; i < 4; i++) {
      shift = (3 - ibyte % 4) * 8;
      trg[3 - i] = buffer[ibyte / 4] >> shift;
      ibyte++;
    }
    tmp = (uint32_t*)trg;
    head.trigger = *tmp;
 
  }

readmwpc()

void readmwpc ( unsigned int *  dbuf, unsigned int  len, int  print = 0  )

protected

Read the MWPC information from the data buffer.

Definition at line 198 of file arichBtestModule.cc.

  {
 
    const int unsigned MAXV1290 = 32;
 
    unsigned int edge;
    unsigned int tdcch;
    unsigned int tdcl;
    unsigned int nhits;
 
    for (int i = 0; i < 32; i++) m_tdc[i] = 0XFFFF;
    //TDC V1290
    //for (int i=0;i<len;i++)   printf("V1290 %d(%d) 0x%08x \n",i,len, dbuf[i],n);
    for (unsigned int i = 0; i < len; i++) {
      int rid = (dbuf[i] >> 27) & 0x1F;
      switch (rid) {
        case 0x18:   // Filler
          if (print1290) printf("Filler  0x%08x %u.data\n", dbuf[i], i);
          break;
        case 0x8:   // Global Header
          if (print1290) printf("Global Header  0x%08x %u.data\n", dbuf[i], i);
          break;
        case 0x10:  // Global Trailer --- Last word of data
          nhits = ((dbuf[i] >> 5) & 0xFFFF); //Word Count = bit 31...21< Word Count: 20 ... 5 > 4...0
          if (print1290) printf("Global Trailer  0x%08x %u.data  STATUS=0x%03x nhits=%u\n", dbuf[i], i, (dbuf[i] >> 24) & 0x7, nhits);
 
          if (nhits != len) {
            if (print1290) printf("V1290 nhits!=len %u %u\n", nhits, len);
          };
          break;
        case 0x11:  // Global Trigger TimeTag
          if (print1290) printf("Global Trigger TimeTag  0x%08x %u.data\n", dbuf[i], i);
          break;
        case 0x1:   // TDC header
          if (print1290) printf("TDC header  0x%08x %u.data evid=%u wc=%u\n", dbuf[i], i, (dbuf[i] >> 12) & 0xFFF, dbuf[i] & 0xFFF);
 
          break;
        case 0x3:   // TDC trailer
          if (print1290) printf("TDC trailer  0x%08x %u.data\n", dbuf[i], i);
 
          break;
        case 0x4:   // TDC Error
          if (print1290) printf("TDC Error  0x%08x %u.data ERR=0x%x\n", dbuf[i], i, dbuf[i] & 0x3FFF);
 
          break;
        case 0x0  :   // TDC data
 
          edge  = (dbuf[i] >> 26) & 0x1;
          tdcch = (dbuf[i] >> 21) & 0x1F;
          tdcl  =   dbuf[i] & 0x1FFFFF;
          if (tdcch < MAXV1290) {
            //if (tdcch>15) gV1290[tdcch]->Fill(tdcl/40);
            //if (tdcch>15) if (tdcl< trg[tdcch-16] && tdcl>16000 && tdcl<20000 ) trg[tdcch-16]=tdcl;
            if (print1290)
              printf("V1290 0x%08x %u. [ch=%2u] edge=%u data=%u \n", dbuf[i], i, tdcch, edge, tdcl);
            m_tdc[tdcch] = tdcl / 40;
          }
 
          break;
 
        default:
          if (print1290) printf("Unknown  0x%08x %u.data\n", dbuf[i], i);
 
          break;
 
      }
    }
 
  }

REG_MODULE() [1/4]

REG_MODULE

(

ARICHDQM

)

---

REG_MODULE() [2/4]

REG_MODULE

(

ARICHPacker

)

---

REG_MODULE() [3/4]

REG_MODULE

(

ARICHRateCal

)

---

REG_MODULE() [4/4]

REG_MODULE

(

ARICHRawUnpacker

)

---

sectorDeadMap()

TCanvas * sectorDeadMap	(	TH1 *	hitMap,
		int	sector
	)

Make display of 70 HAPDs' 2D dead/alive map of the sector.

Parameters

hitMap	1D hit map of all channels.
sector	Sector identifier [1:6].

Returns

Display of 70 HAPDs' 2D dead/alive map.

Definition at line 176 of file hitMapMaker.cc.

  {
    //TH1* m_hitMap = NULL;
    TH2* m_moduleDeadMap = NULL;
    TExec* ex1 = NULL;
 
    //m_hitMap = hitMap;
    TPad* p_hitMaps[421] = {};
    TCanvas* m_sectorDeadMap = new TCanvas(Form("c_deadMap%d", sector - 1), Form("Dead chip map of sector%d", sector - 1), 600, 400);
    for (int i = 1; i < 421; i++) {
      for (int iRing = 0; iRing < 7; iRing++) {
        if (((iRing + 13)*iRing / 2) * 6 + (iRing + 7) * (sector - 1) < i && i <= ((iRing + 13)*iRing / 2) * 6 + (iRing + 7) * (sector)) {
          m_sectorDeadMap->cd();
          p_hitMaps[i] = new TPad(Form("p_deadMap%d", i), "",
                                  (double)((double)(6 - iRing) / 2 + ((i - ((iRing + 13)*iRing / 2) * 6 + (iRing + 7) * (sector - 1) - 1) % (iRing + 7))) / 13,
                                  (double)iRing / 7, (double)((double)(8 - iRing) / 2 + ((i - ((iRing + 13)*iRing / 2) * 6 + (iRing + 7) * (sector - 1) - 1) %
                                                              (iRing + 7))) / 13, (double)(iRing + 1) / 7);
          p_hitMaps[i]->Draw();
          p_hitMaps[i]->SetNumber(i);
          m_sectorDeadMap->cd(i);
          m_moduleDeadMap = moduleDeadMap(hitMap, i);
          m_moduleDeadMap->SetTitleSize(0, "xyz");
          m_moduleDeadMap->SetTitle(0);
          m_moduleDeadMap->SetLineWidth(1);
          m_moduleDeadMap->SetStats(0);
          gStyle->SetLabelColor(0, "xyz");
          ex1 = new TExec("ex1", "deadPalette();");
          ex1->Draw();
          m_moduleDeadMap->Draw("colz");
        }
      }
    }
    return m_sectorDeadMap;
  }

sectorHitMap()

TCanvas * sectorHitMap	(	TH1 *	hitMap,
		int	sector
	)

Make display of 70 HAPDs' 2D hit map of the sector.

Parameters

hitMap	1D hit map of all channels.
sector	Sector identifier [1:6].

Returns

Display of 70 HAPDs' 2D hit map.

Definition at line 144 of file hitMapMaker.cc.

  {
    //    TH1* m_hitMap = NULL;
    TH2* m_moduleHitMap = NULL;
 
    //m_hitMap = hitMap;
    TPad* p_hitMaps[421] = {};
    TCanvas* m_sectorHitMap = new TCanvas(Form("c_hitMap%d", sector - 1), Form("Hit map of sector%d", sector - 1), 600, 400);
    for (int i = 1; i < 421; i++) {
      for (int iRing = 0; iRing < 7; iRing++) {
        if (((iRing + 13)*iRing / 2) * 6 + (iRing + 7) * (sector - 1) < i && i <= ((iRing + 13)*iRing / 2) * 6 + (iRing + 7) * (sector)) {
          m_sectorHitMap->cd();
          p_hitMaps[i] = new TPad(Form("p_hitMap%d", i), "",
                                  (double)((double)(6 - iRing) / 2 + ((i - ((iRing + 13)*iRing / 2) * 6 + (iRing + 7) * (sector - 1) - 1) % (iRing + 7))) / 13,
                                  (double)iRing / 7, (double)((double)(8 - iRing) / 2 + ((i - ((iRing + 13)*iRing / 2) * 6 + (iRing + 7) * (sector - 1) - 1) %
                                                              (iRing + 7))) / 13, (double)(iRing + 1) / 7);
          p_hitMaps[i]->Draw();
          p_hitMaps[i]->SetNumber(i);
          m_sectorHitMap->cd(i);
          m_moduleHitMap = moduleHitMap(hitMap, i);
          m_moduleHitMap->SetTitleSize(0, "xyz");
          m_moduleHitMap->SetTitle(0);
          m_moduleHitMap->SetLineWidth(1);
          gStyle->SetLabelColor(0, "xyz");
          m_moduleHitMap->SetStats(0);
          m_moduleHitMap->Draw("col");
        }
      }
    }
    return m_sectorHitMap;
  }

setTrackAngleResolution()

void setTrackAngleResolution

(

double

aRes

)

Sets track direction resolution (from tracking).

Definition at line 694 of file ARICHReconstruction.cc.

  {
    m_trackAngRes = aRes;
  }

setTrackPositionResolution()

void setTrackPositionResolution

(

double

pRes

)

Sets track position resolution (from tracking).

Definition at line 690 of file ARICHReconstruction.cc.

  {
    m_trackPosRes = pRes;
  }

skipdata()

int skipdata ( gzFile  fp )

protected

Skip the data part of the record.

Definition at line 190 of file arichBtestModule.cc.

  {
    unsigned int u;
    gzread(fp, &u, sizeof(unsigned int));
    gzseek(fp, u * sizeof(unsigned int), SEEK_CUR);
    return u + 1;
  }

smearTrack()

int smearTrack

(

ARICHTrack &

arichTrack

)

Smears track parameters ("simulate" the uncertainties of tracking).

Definition at line 133 of file ARICHReconstruction.cc.

  {
    double a = gRandom->Gaus(0, m_trackAngRes);
    double b = gRandom->Gaus(0, m_trackAngRes);
    ROOT::Math::XYZVector dirf(a, b, sqrt(1 - a * a - b * b));
    double dx = gRandom->Gaus(0, m_trackPosRes);
    double dy = gRandom->Gaus(0, m_trackPosRes);
    ROOT::Math::XYZVector mod(dx, dy, 0);
    ROOT::Math::XYZVector rpoint = arichTrack.getPosition() + mod;
    ROOT::Math::XYZVector odir  = arichTrack.getDirection();
    double omomentum  = arichTrack.getMomentum();
    ROOT::Math::XYZVector rdir = TransformFromFixed(odir) * dirf;  // global system
    double rmomentum = omomentum;
    arichTrack.setReconstructedValues(rpoint, ROOT::Math::XYZVector(rdir), rmomentum);
    return 1;
  }

terminate() [1/2]

void terminate ( void  )

overridevirtual

Is called at the end of your Module.

Function is called only once at the end of your job at the end of the corresponding module. This function is for cleaning up, closing files, etc.

Reimplemented from Module.

Definition at line 568 of file arichBtestModule.cc.

  {
    int j = 1;
    BOOST_FOREACH(const std::string & fname, m_runList) {
      B2INFO(m_eveList[j] << " events processed from file " << fname);
      j++;
    }
    for (int i = 0; i < 4; i++) {
      //ARICHTracking* w = &m_mwpc[i];
      B2INFO(i << " a1=" << m_mwpc[i].tdc[0] << " a2="  << m_mwpc[i].tdc[1] << " a3=" << m_mwpc[i].tdc[2] << " a2="  << m_mwpc[i].tdc[3]
             << " A=" << m_mwpc[i].atdc);
    }
 
  }

terminate() [2/2]

void terminate ( void  )

overridevirtual

Termination action.

Clean-up, close files, summarize statistics, etc.

Reimplemented from Module.

Definition at line 335 of file ARICHNtupleModule.cc.

  {
    m_file->cd();
    m_tree->Write();
    m_file->Close();
  }

transformTrackToLocal()

void transformTrackToLocal	(	ARICHTrack &	arichTrack,
		bool	align
	)

Transforms track parameters from global Belle2 to ARICH local frame.

Definition at line 722 of file ARICHReconstruction.cc.

  {
    // tranform track from Belle II to local ARICH frame
    ROOT::Math::XYZVector locPos = m_arichgp->getMasterVolume().pointToLocal(arichTrack.getPosition());
    ROOT::Math::XYZVector locDir = m_arichgp->getMasterVolume().momentumToLocal(arichTrack.getDirection());
 
    // apply the alignment correction
    if (align && m_alignp.isValid()) {
      // apply global alignment correction
      locPos = m_alignp->pointToLocal(locPos);
      locDir = m_alignp->momentumToLocal(locDir);
    }
 
    // set parameters and return
    // is it needed to extrapolate to z of aerogel in local frame?? tabun
    arichTrack.setReconstructedValues(locPos, locDir, arichTrack.getMomentum());
    return;
  }

writeHeader()

void writeHeader	(	int *	buffer,
		unsigned &	ibyte,
		const ARICHRawHeader &	head
	)

TODO!

Definition at line 214 of file ARICHPackerModule.cc.

  {
 
    unsigned char line1[4];
    int shift;
 
    line1[3] = head.type;
    line1[2] = head.version;
    line1[1] = head.mergerID;
    line1[0] = head.FEBSlot;
 
    for (int i = 0; i < 4; i++) {
      shift = (3 - ibyte % 4) * 8;
      buffer[ibyte / 4] |= line1[3 - i] << shift;
      ibyte++;
    }
 
    auto len = reinterpret_cast<const unsigned char*>(&head.length);
 
    for (int i = 0; i < 4; i++) {
      shift = (3 - ibyte % 4) * 8;
      buffer[ibyte / 4] |= len[3 - i] << shift;
      ibyte++;
    }
 
    // trigger number
    auto trg = reinterpret_cast<const unsigned char*>(&head.trigger);
    for (int i = 0; i < 4; i++) {
      shift = (3 - ibyte % 4) * 8;
      buffer[ibyte / 4] |= trg[3 - i] << shift;
      ibyte++;
    }
 
  }

~ARICHDigitizerModule()

~ARICHDigitizerModule ( )

virtual

Destructor.

Definition at line 68 of file ARICHDigitizerModule.cc.

  {
 
  }

~ARICHDQMModule()

~ARICHDQMModule ( )

virtual

Destructor.

Definition at line 57 of file ARICHDQMModule.cc.

  {
  }

~ARICHFillHitsModule()

~ARICHFillHitsModule ( )

virtual

Destructor.

Definition at line 53 of file ARICHFillHitsModule.cc.

  {
  }

~ARICHMCParticlesModule()

~ARICHMCParticlesModule ( )

virtual

Destructor.

Definition at line 48 of file ARICHMCParticlesModule.cc.

  {
  }

~ARICHNtupleModule()

~ARICHNtupleModule ( )

virtual

Destructor.

Definition at line 65 of file ARICHNtupleModule.cc.

  {
  }

~ARICHPackerModule()

~ARICHPackerModule ( )

virtual

Destructor.

Definition at line 62 of file ARICHPackerModule.cc.

  {
  }

~ARICHRateCalModule()

~ARICHRateCalModule ( )

virtual

Destructor.

Definition at line 60 of file ARICHRateCalModule.cc.

  {
  }

~ARICHRawUnpackerModule()

~ARICHRawUnpackerModule ( )

virtual

Destructor.

Definition at line 41 of file ARICHRawUnpackerModule.cc.

  {
  }

~ARICHReconstructorModule()

~ARICHReconstructorModule ( )

virtual

Destructor.

Definition at line 69 of file ARICHReconstructorModule.cc.

  {
    if (m_ana) delete m_ana;
  }

~ARICHUnpackerModule()

~ARICHUnpackerModule ( )

virtual

Destructor.

Definition at line 69 of file ARICHUnpackerModule.cc.

  {
  }

Variable Documentation

hapd

TH1F* hapd[6]

histogram of hits for each hapd

Definition at line 82 of file arichBtestModule.cc.

m_tuple

TNtuple* m_tuple

ntuple containing hapd hits

Definition at line 81 of file arichBtestModule.cc.

mwpc_diff

TH1F* mwpc_diff[4][2]

tdc difference from mwpcs

Definition at line 84 of file arichBtestModule.cc.

mwpc_residuals

TH1F* mwpc_residuals[4][2]

residuals from mwpcs

Definition at line 87 of file arichBtestModule.cc.

mwpc_residualsz

TH2F* mwpc_residualsz[4][2]

z-residuals from mwpcs

Definition at line 89 of file arichBtestModule.cc.

mwpc_sum

TH1F* mwpc_sum[4][2]

tdc sum from mwpcs

Definition at line 85 of file arichBtestModule.cc.

mwpc_sum_cut

TH1F* mwpc_sum_cut[4][2]

tdc sum from mwpcs, with sum cut applied

Definition at line 86 of file arichBtestModule.cc.

mwpc_tdc

TH1F* mwpc_tdc[4][5]

tdc information from mwpcs

Definition at line 83 of file arichBtestModule.cc.

mwpc_xy

TH2F* mwpc_xy[4]

calculated x-y track positions

Definition at line 88 of file arichBtestModule.cc.