00001 #include "SimPmtM64.h"
00002 #include "SimPmtM64CrosstalkTable.h"
00003 #include "SimPmtMaker.h"
00004 #include "MessageService/MsgService.h"
00005 #include "TMath.h"
00006
00007 CVSID("$Id: SimPmtM64.cxx,v 1.26 2008/11/18 17:31:02 rhatcher Exp $");
00008
00009 SimPmtMakerProxy<SimPmtM64> gSimPmtM64Proxy("SimPmtM64");
00010
00011 ClassImp(SimPmtM64)
00012
00013
00014 SimPmtM64CrosstalkTable* SimPmtM64::fsCrosstalkTable = NULL;
00015
00016 SimPmtM64::SimPmtM64(PlexPixelSpotId tube,
00017 VldContext context,
00018 TRandom* random ) :
00019 SimPmt(tube,context,random,64,1),
00020 fQieClock(context)
00021 {
00022
00023
00024
00025
00026
00027
00028
00029
00030 if(fsCrosstalkTable == NULL) {
00031 fsCrosstalkTable = new SimPmtM64CrosstalkTable(context);
00032 };
00033
00034
00035
00036 Double_t g, w;
00037 for(Int_t pix = 1; pix <= fNPixels; pix++) {
00038 GetGainAndWidth(pix, 1, g, w);
00039
00040 if(w>0.8) w = 0.8;
00041
00042
00043
00044
00045
00046 Float_t secemm = 1.0/(w*w); if(w>0.8) w = 0.8;
00047
00048 fPixelGains[pix] = g;
00049 fPixelSecEmm[pix] = secemm;
00050 };
00051 }
00052
00053 void SimPmtM64::Reset( const VldContext& context )
00054 {
00055
00056
00057 fQieClock.Reset(context);
00058
00059
00060 if(fsCrosstalkTable) fsCrosstalkTable->Reset(context);
00061 else fsCrosstalkTable = new SimPmtM64CrosstalkTable(context);
00062
00063 SimPmt::Reset(context);
00064
00065 }
00066
00067 Int_t SimPmtM64::TimeToBucket(Double_t time)
00068 {
00069
00070 return fQieClock.GetBucketForTime(time);
00071 }
00072
00073 Double_t SimPmtM64::BucketToStartTime(Int_t bucket)
00074 {
00075
00076 return fQieClock.GetTimeOfBucket(bucket);
00077 }
00078
00079 Double_t SimPmtM64::BucketToStopTime(Int_t bucket)
00080 {
00081
00082 return fQieClock.GetTimeOfBucket(bucket+1);
00083 }
00084
00085
00086 Float_t SimPmtM64::GetBucketDuration( Int_t )
00087 {
00088 return fQieClock.fNDTick;
00089 }
00090
00091 void SimPmtM64::Print(Option_t* ) const
00092 {
00093 printf(" SimPmtM64::Print() -- Tube: %s\n",fTube.AsString("t"));
00094
00095 SimPmtBucketIterator it(*(SimPmt*)this);
00096 printf(" Photoelectrons(Npe) Charge (fC) \n");
00097 for( ; !it.End(); it.Next() ) {
00098 int ibucket = it.BucketId();
00099 printf(" ---------------------------Time bucket %d----------------------\n",ibucket);
00100 for(int row=0; row < 8; row++) {
00101 printf(" %3.0f %3.0f %3.0f %3.0f %3.0f %3.0f %3.0f %3.0f",
00102 GetPe(row*8+1,1,ibucket),
00103 GetPe(row*8+2,1,ibucket),
00104 GetPe(row*8+3,1,ibucket),
00105 GetPe(row*8+4,1,ibucket),
00106 GetPe(row*8+5,1,ibucket),
00107 GetPe(row*8+6,1,ibucket),
00108 GetPe(row*8+7,1,ibucket),
00109 GetPe(row*8+8,1,ibucket));
00110
00111 printf(" | %5.0f %5.0f %5.0f %5.0f %5.0f %5.0f %5.0f %5.0f",
00112 GetCharge(row*8+1,ibucket)/Munits::fC,
00113 GetCharge(row*8+2,ibucket)/Munits::fC,
00114 GetCharge(row*8+3,ibucket)/Munits::fC,
00115 GetCharge(row*8+4,ibucket)/Munits::fC,
00116 GetCharge(row*8+5,ibucket)/Munits::fC,
00117 GetCharge(row*8+6,ibucket)/Munits::fC,
00118 GetCharge(row*8+7,ibucket)/Munits::fC,
00119 GetCharge(row*8+8,ibucket)/Munits::fC);
00120
00121 printf("\n");
00122 }
00123 }
00124 printf("\n");
00125 }
00126
00127
00128 Float_t SimPmtM64::GenChargeFromPE( Int_t pixel, Int_t , Float_t pe )
00129 {
00130
00131
00132
00133
00134
00135
00136
00137
00138
00139
00140 if(pe<=0) return 0;
00141
00142 Float_t gain = fPixelGains[pixel];
00143 Float_t secemm = fPixelSecEmm[pixel];
00144
00145
00146 Float_t mean2ndaryPe = pe * secemm;
00147 Float_t secondary_pe = RandomGenPoisson(mean2ndaryPe);
00148
00149 MSG("DetSim",Msg::kVerbose) << "Rolled: " << pe << " pe -> " << secondary_pe / secemm << " charge" << " with " << gain << " gain" << endl;
00150
00151
00152 Float_t q = (secondary_pe / secemm)
00153 * gain
00154 * Munits::e_SI;
00155
00156 return q;
00157 }
00158
00159 Float_t SimPmtM64::GenDarkNoiseCharge( Int_t , Float_t )
00160 {
00161
00162
00163
00164
00165
00166
00167
00168
00169 return 0;
00170 }
00171
00172
00173 Float_t SimPmtM64::GenNonlinearCharge( Int_t
00174 , Float_t inCharge )
00175 {
00176
00177
00178
00179
00180
00181
00182
00183
00184 return inCharge;
00185 }
00186
00187
00188
00189 void SimPmtM64::SimulateOpticalXtalk()
00190 {
00191
00192
00193
00194
00195
00196
00197
00198
00199
00200 for(SimPmtBucketIterator it(*this); !it.End(); it.Next()) {
00201 SimPmtTimeBucket& pmttb = it.Bucket();
00202
00203
00204 for(Int_t injPix = 1; injPix <= fNPixels; injPix++) {
00205 SimPixelTimeBucket& pixtb_inj = pmttb.GetPixelBucket(injPix);
00206
00207
00208
00209 float injPE = pixtb_inj.GetPE(1);
00210 if( injPE > 0) {
00211
00212
00213
00214 float prob = fsCrosstalkTable->fOptTotProb[injPix];
00215 prob *= fsPmtScaleOpticalCrosstalk;
00216 float ttalk = fRandom->Poisson(prob*injPE);
00217 int ntalk = TMath::Nint(ttalk);
00218
00219 for(int ipe=0;ipe<ntalk;ipe++) {
00220
00221
00222
00223 int whichPix = injPix;
00224 float r = fRandom->Uniform();
00225 for(Int_t xPix = 1; xPix <= fNPixels; xPix++) {
00226 if(r<fsCrosstalkTable->fOptDist[injPix][xPix]) {
00227 whichPix=xPix;
00228 break;
00229 }
00230 }
00231 if(whichPix==injPix){
00232 MSG("DetSim",Msg::kWarning) << "Random number error choosing xtalk pixel.";
00233 continue;
00234 }
00235
00236
00237
00238
00239
00240 SimPixelTimeBucket& pixtb = pmttb.GetPixelBucket(whichPix);
00241
00242
00243 MoveOpticalPE(1,
00244 pixtb_inj,1,
00245 pixtb, 1 );
00246
00247 pixtb.SetTruthBit(DigiSignal::kCrosstalkOptical);
00248
00249 }
00250
00251 }
00252 }
00253
00254
00255 }
00256
00257
00258
00259 CopyPEtoPEXtalk();
00260 }
00261
00262
00263 Float_t SimPmtM64::GenElectricalCrosstalk( Int_t injPixel,
00264 Int_t xPixel, Float_t injCharge )
00265 {
00266
00267
00268
00269
00270
00271
00272
00273
00274 if(injCharge <= 0) return 0;
00275
00276 float frac = fsCrosstalkTable->fElecFrac[injPixel][xPixel];
00277 float k = fsCrosstalkTable->fElecK[injPixel][xPixel];
00278
00279 frac *= fsPmtScaleChargeCrosstalk;
00280
00281
00282 float charge = frac*injCharge;
00283
00284
00285
00286
00287
00288
00289 const float kMinVariance2 = (2.0*Munits::fC)*(2.0*Munits::fC);
00290
00291 float width2 = fabs(k*charge);
00292 if(width2 > kMinVariance2 ){
00293 float width = sqrt(width2);
00294 return fRandom->Gaus(charge,width);
00295 };
00296
00297 return charge;
00298 }
00299