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