groot/macros/NoiseSimulator.C

0001 TH1D* BlobPosition;
0002 TH1D* BlobSigma;
0003 TH1D* BlobSigmaGen;
0004 TH1D* BlobPos;
0005 TH1D* BlobCharge;
0006 TF1* BlobFit;
0007 TF1* NORMGAUSS;
0008
0009 TH2D* ChargeQuality;
0010 TH2D* MeanQuality;
0011 TH2D* SigmaQuality;
0012
0013 TH1D* ChargeDiffBeforeNoise;
0014 TH1D* MeanDiffBeforeNoise;
0015 TH1D* SigmaDiffBeforeNoise;
0016
0017 TH1D* ChargeDiffAfterNoise;
0018 TH1D* MeanDiffAfterNoise;
0019 TH1D* SigmaDiffAfterNoise;
0020
0021 void NoiseSimulator(int Nevent = 10000, double NOISE = 9)
0022 {
0023   gSystem->Load("libgroot");
0024
0025   groot* Tree = groot::instance();
0026
0027   double Minphi =  9999;
0028   double Maxphi = -9999;
0029   for (int i=0; i<128; i++)
0030     {
0031       if (Tree->ZigzagMap2[0][i])
0032     {
0033       double phi = Tree->ZigzagMap2[0][i]->PCenter();
0034       if (phi < Minphi) Minphi = phi;
0035       if (phi > Maxphi) Maxphi = phi;
0036     }
0037     }
0038   double dPhi = fabs( Tree->ZigzagMap2[0][90]->PCenter() - Tree->ZigzagMap2[0][91]->PCenter() );
0039   cout << dPhi << endl;
0040   int NBINS = int( (Maxphi - Minphi)/dPhi + 1.5 );
0041   cout << NBINS << endl;
0042   double left  = Minphi - 0.5*dPhi;
0043   double right = Maxphi + 0.5*dPhi;
0044   BlobFit      = new TF1("BlobFit", "[0]*exp(-(x-[1])*(x-[1])/(2*[2]*[2]))", left, right);
0045   BlobPos      = new TH1D("BlobPos", "BlobPos", NBINS, left, right);
0046   BlobSigma    = new TH1D("BlobSigma", "BlobSigma", 250, 0.0, 8.0E-3);
0047   BlobSigmaGen = new TH1D("BlobSigmaGen", "BlobSigmaGen", 250, 0.0, 8.0E-3);
0048
0049   BlobPosition = new TH1D("BlobPosition", "BlobPosition", NBINS*10, left, right);
0050   BlobCharge   = new TH1D("BlobCharge", "BlobCharge", 256, 0, 2000);
0051
0052   ChargeQuality = new TH2D("ChargeQuality", "ChargeQuality", 256, 0, 2000, 256, 0, 2000);
0053   MeanQuality   = new TH2D("MeanQuality", "MeanQuality", NBINS*10, left, right, NBINS*10, left, right);
0054   SigmaQuality  = new TH2D("SigmaQuality", "SigmaQuality", 250, 0.0, 8.0E-3, 250, 0.0, 8.0E-3);
0055
0056   ChargeDiffBeforeNoise = new TH1D("ChargeDiffBeforeNoise", "ChargeDiffBeforeNoise", 300, -10, 10);
0057   MeanDiffBeforeNoise   = new TH1D("MeanDiffBeforeNoise", "MeanDiffBeforeNoise", 300, -dPhi, dPhi);
0058   SigmaDiffBeforeNoise  = new TH1D("SigmaDiffBeforeNoise", "SigmaDiffBeforeNoise", 300, -1e-4, 1e-4);
0059
0060   ChargeDiffAfterNoise = new TH1D("ChargeDiffAfterNoise", "ChargeDiffAfterNoise", 300, -100, 100);
0061   MeanDiffAfterNoise   = new TH1D("MeanDiffAfterNoise", "MeanDiffAfterNoise", 300, -3e-4, 3e-4);
0062   SigmaDiffAfterNoise  = new TH1D("SigmaDiffAfterNoise", "SigmaDiffAfterNoise", 300, -3e-4, 3e-4);
0063
0064   double sqrt2pi = 2.506628275;
0065   NORMGAUSS    = new TF1("NORMGAUSS", "[0]*0.00396825/([2]*2.506628275)*exp(-(x-[1])*(x-[1])/(2*[2]*[2]))", left, right);
0066
0067   TRandom randy;
0068
0069   double R5 = 467900.0;             // microns
0070
0071   //  These are values at 17.5 inches...
0072   //double BlobSigmaMean  = 33.5e-4;  // radians
0073   //double BlobSigmaSigma =  4.1e-4*4.1/4.66;  // radians
0074
0075   //  These are values at 14.0 inches...
0076   //double BlobSigmaMean  = 31.7e-4;  // radians
0077   //double BlobSigmaSigma =  3.5e-4*4.1/4.66;  // radians
0078
0079   //  These are values at 10.0 inches...
0080   //double BlobSigmaMean  = 29.4e-4;  // radians
0081   //double BlobSigmaSigma =  3.0e-4*4.1/4.66;  // radians
0082
0083   //  These are values at 4.0 inches...
0084   double BlobSigmaMean  = 25.1e-4;  // radians
0085   double BlobSigmaSigma =  2.0e-4*4.1/4.66;  // radians
0086
0087   double BlobLandauMPV   = 326.8;
0088   double BlobLandauSigma = 135.1;
0089
0090   for (int i=0; i<Nevent; i++)
0091     {
0092       double MEAN = (left + 0.1*(right-left)) + 0.8*(right-left)*randy.Rndm();
0093       BlobPosition->Fill(MEAN);
0094
0095       double SIGMA = randy.Gaus(BlobSigmaMean, BlobSigmaSigma);
0096       BlobSigmaGen->Fill(SIGMA); // radians
0097
0098       double CHARGE = randy.Landau(BlobLandauMPV, BlobLandauSigma);
0099       BlobCharge->Fill(CHARGE);
0100
0101       NORMGAUSS->SetParameter(0,CHARGE);
0102       NORMGAUSS->SetParameter(1,MEAN);
0103       NORMGAUSS->SetParameter(2,SIGMA);
0104
0105       for (int j=1; j<BlobPos->GetNbinsX()+1; j++)
0106     {
0107       double x = BlobPos->GetBinCenter(j);
0108       double height = NORMGAUSS->Eval(x);
0109       BlobPos->SetBinContent(j,height);
0110       BlobPos->SetBinError(j,9.0);
0111     }
0112
0113       double max   = BlobPos->GetMaximum();
0114       double mean  = BlobPos->GetBinCenter(BlobPos->GetMaximumBin());
0115       double sigma = 0.003;
0116       if (sigma < dPhi/3.0) sigma = dPhi/3.0;
0117       if (sigma > 5.0*dPhi) sigma = 5.0*dPhi;
0118
0119       BlobFit->SetParameter(0,max);
0120       BlobFit->SetParameter(1,mean);
0121       BlobFit->SetParameter(2,sigma);
0122
0123       BlobFit->SetParLimits(0, max/2.0, 3.0*max);
0124       BlobFit->SetParLimits(1, mean-1.0*dPhi, mean+1.0*dPhi);
0125       BlobFit->SetParLimits(2, 0.0007, 0.0055);
0126
0127       BlobPos->Fit(BlobFit,"Q0");
0128
0129       double FITTEDCHARGE = BlobFit->GetParameter(0)*BlobFit->GetParameter(2)*sqrt2pi/dPhi;
0130       double FITTEDMEAN   = BlobFit->GetParameter(1);
0131       double FITTEDSIGMA  = BlobFit->GetParameter(2);
0132
0133       ChargeQuality->Fill(CHARGE,FITTEDCHARGE);
0134       MeanQuality->Fill(MEAN,FITTEDMEAN);
0135       SigmaQuality->Fill(SIGMA,FITTEDSIGMA);
0136
0137       ChargeDiffBeforeNoise->Fill(CHARGE-FITTEDCHARGE);
0138       MeanDiffBeforeNoise->Fill(MEAN-FITTEDMEAN);
0139       SigmaDiffBeforeNoise->Fill(SIGMA-FITTEDSIGMA);
0140
0141       for (int j=1; j<BlobPos->GetNbinsX()+1; j++)
0142     {
0143       double oldheight = BlobPos->GetBinContent(j);
0144       double newheight = oldheight + randy.Gaus(0,NOISE);
0145       BlobPos->SetBinContent(j,newheight);
0146     }
0147
0148       max   = BlobPos->GetMaximum();
0149       mean  = BlobPos->GetBinCenter(BlobPos->GetMaximumBin());
0150       sigma = 0.003;
0151       if (sigma < dPhi/3.0) sigma = dPhi/3.0;
0152       if (sigma > 5.0*dPhi) sigma = 5.0*dPhi;
0153
0154       BlobFit->SetParameter(0,max);
0155       BlobFit->SetParameter(1,mean);
0156       BlobFit->SetParameter(2,sigma);
0157
0158       BlobFit->SetParLimits(0, max/2.0, 3.0*max);
0159       BlobFit->SetParLimits(1, mean-1.0*dPhi, mean+1.0*dPhi);
0160       BlobFit->SetParLimits(2, 0.0007, 0.0055);
0161
0162       BlobPos->Fit(BlobFit,"Q0");
0163
0164       double REFITTEDCHARGE = BlobFit->GetParameter(0)*BlobFit->GetParameter(2)*sqrt2pi;
0165       double REFITTEDMEAN   = BlobFit->GetParameter(1);
0166       double REFITTEDSIGMA  = BlobFit->GetParameter(2);
0167
0168       ChargeDiffAfterNoise->Fill(CHARGE-REFITTEDCHARGE);
0169       MeanDiffAfterNoise->Fill(MEAN-REFITTEDMEAN);
0170       SigmaDiffAfterNoise->Fill(SIGMA-REFITTEDSIGMA);
0171
0172       BlobSigma->Fill(REFITTEDSIGMA);
0173
0174     }
0175
0176 }