Changeset 2920 for trunk/MagicSoft/Mars


Ignore:
Timestamp:
01/26/04 19:32:58 (21 years ago)
Author:
gaug
Message:
*** empty log message ***
Location:
trunk/MagicSoft/Mars
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • trunk/MagicSoft/Mars/Changelog

    r2919 r2920  
    1111       in arrays: gsNCells, gsNTrigPixels,
    1212       gsNPixInCell, gsNLutInCell, gsNPixInLut, fNumPixCell.
    13      - Added method MMcTriggerLvl2::GetCellCompactPixel(int cell, MGeomCam *fCam)
     13     - Added method
     14       MMcTriggerLvl2::GetCellCompactPixel(int cell, MGeomCam *fCam)
    1415       which computes compact pixels into a given L2T macrocell.
    15      - Added method MMcTriggerLvl2::CalcBiggerCellPseudoSize() which
    16        computes fCellPseudoSize, the maximum Pseudo Size into L2T macrocells
    17      - Added method MMcTriggerLvl2::GetCellPseudoSize() const which
    18        returns fCellPseudoSize
    19      - Added method MMcTriggerLvl2::IsPixelInCell(Int_t pixel, Int_t cell),
     16     - Added method
     17       MMcTriggerLvl2::CalcBiggerCellPseudoSize()
     18       which computes fCellPseudoSize, the maximum Pseudo Size into L2T
     19       macrocells
     20     - Added method
     21       MMcTriggerLvl2::GetCellPseudoSize() const
     22       which returns fCellPseudoSize
     23     - Added method
     24       MMcTriggerLvl2::IsPixelInCell(Int_t pixel, Int_t cell),
    2025       which controls whether a pixel belongs to a given L2T cell.
    21      - Added method MMcTriggerLvl2::GetMaxCell() const which returns fMaxCell, the cell with
    22        the maximum fCellPseudoSize.
    23      
    24 
    25 
    26  2004/01/26: Markus Gaug
     26     - Added method
     27       MMcTriggerLvl2::GetMaxCell() const
     28       which returns fMaxCell, the cell with the maximum
     29       fCellPseudoSize.
     30     
     31
     32
     33 2004/01/26: Markus Gaug, Michele Doro
    2734
    2835   * manalysis/MArrivalTime.[h,cc], manalysis/MArrivalTimeCalc.[h,cc]:
    … …  
    3340
    3441   * mcalib/MHCalibrationBlindPixel.[h,cc]:
    35      - force mu_0 in Blind Pixel Fit to be around 0 in fKPoisson4
     42     - force mu_{0} in Blind Pixel Fit to be around 0 in fKPoisson4
     43     - implement combined Polya fit and Michele's back-scattered electron
     44       fit
    3645
    3746

  • trunk/MagicSoft/Mars/mcalib/MHCalibrationBlindPixel.cc

    r2913 r2920  
    260260{
    261261
    262     gStyle->SetOptFit(0);
     262    gStyle->SetOptFit(1);
    263263    gStyle->SetOptStat(111111);
    264264
    … …  
    366366      break;
    367367    case kEPolya:
    368       break;
    369 
     368      fSinglePheFit = new TF1("SinglePheFit",&fPolya,min,max,8);
     369      break;
    370370    case kEMichele:
    371371      break;
    … …  
    381381  const Double_t mu_1_guess = mu_0_guess + 50.;
    382382  const Double_t si_1_guess = si_0_guess + si_0_guess;
     383  // Michele
     384  const Double_t lambda_1cat_guess = 0.5;
     385  const Double_t lambda_1dyn_guess = 0.5;
     386  const Double_t mu_1cat_guess = mu_0_guess + 50.;
     387  const Double_t mu_1dyn_guess = mu_0_guess + 20.;
     388  const Double_t si_1cat_guess = si_0_guess + si_0_guess;
     389  const Double_t si_1dyn_guess = si_0_guess;
     390  // Polya
     391  const Double_t excessPoisson_guess = 0.5;
     392  const Double_t delta1_guess     = 8.;
     393  const Double_t delta2_guess     = 5.;
     394  const Double_t electronicAmp_guess  = 0.0025;
     395
    383396  //
    384397  // Initialize boundaries and start parameters
    … …  
    389402    case kEPoisson4:
    390403      fSinglePheFit->SetParameters(lambda_guess,mu_0_guess,mu_1_guess,si_0_guess,si_1_guess,norm);
    391       fSinglePheFit->SetParNames("#lambda","#mu_0","#mu_1","#sigma_0","#sigma_1","Area");
     404      fSinglePheFit->SetParNames("#lambda","#mu_{0}","#mu_{1}","#sigma_{0}","#sigma_{1}","Area");
    392405      fSinglePheFit->SetParLimits(0,0.,1.);
    393406      fSinglePheFit->SetParLimits(1,-2.,2.);
    … …  
    400413    case kEPoisson6:
    401414      fSinglePheFit->SetParameters(lambda_guess,mu_0_guess,mu_1_guess,si_0_guess,si_1_guess,norm);
    402       fSinglePheFit->SetParNames("#lambda","#mu_0","#mu_1","#sigma_0","#sigma_1","Area");
     415      fSinglePheFit->SetParNames("#lambda","#mu_{0}","#mu_{1}","#sigma_{0}","#sigma_{1}","Area");
    403416      fSinglePheFit->SetParLimits(0,0.,1.);
    404417      fSinglePheFit->SetParLimits(1,min,(max-min)/1.5);
    … …  
    410423
    411424    case kEPolya:
     425      fSinglePheFit->SetParameters(lambda_guess, excessPoisson_guess,
     426                                   delta1_guess,delta2_guess,
     427                                   electronicAmp_guess,
     428                                   si_0_guess,
     429                                   norm, mu_0_guess);
     430      fSinglePheFit->SetParNames("#lambda","b_{tot}",
     431                                 "#delta_{1}","#delta_{2}",
     432                                 "amp_{e}","#sigma_{0}",
     433                                 "Area", "#mu_{0}");
     434      fSinglePheFit->SetParLimits(0,0.,1.);
     435      fSinglePheFit->SetParLimits(1,0.,1.);
     436      fSinglePheFit->SetParLimits(2,6.,12.);   
     437      fSinglePheFit->SetParLimits(3,3.,8.);   
     438      fSinglePheFit->SetParLimits(4,0.,0.005);   
     439      fSinglePheFit->SetParLimits(5,min,(max-min)/1.5);
     440      fSinglePheFit->SetParLimits(6,norm-0.1,norm+0.1);
     441      fSinglePheFit->SetParLimits(7,-35.,15.);
    412442      break;
    413443
    414444    case kEMichele:
     445
     446
    415447      break;
    416448
    … …  
    450482      fSigma1Err = fSinglePheFit->GetParError(4);
    451483      break;
     484    case kEPolya:
     485      fLambda =  fSinglePheFit->GetParameter(0);
     486      fMu0    =  fSinglePheFit->GetParameter(7);
     487      fMu1    = 0.;
     488      fSigma0 =  fSinglePheFit->GetParameter(5);
     489      fSigma1 = 0.;
     490
     491      fLambdaErr = fSinglePheFit->GetParError(0);
     492      fMu0Err    = fSinglePheFit->GetParError(7);
     493      fMu1Err    = 0.;
     494      fSigma0Err = fSinglePheFit->GetParError(5);
     495      fSigma1Err = 0.;
    452496    default:
    453497      break;
    … …  
    479523
    480524  // Perform the cross-check fitting only the pedestal:
    481   fSinglePhePedFit = new TF1("GausPed","gaus",rmin,fMu0);
     525  fSinglePhePedFit = new TF1("GausPed","gaus",rmin,0.);
    482526  fHBlindPixelCharge->Fit(fSinglePhePedFit,opt);
    483527

  • trunk/MagicSoft/Mars/mcalib/MHCalibrationBlindPixel.h

    r2913 r2920  
    7171
    7272  //Getters
    73   const Double_t GetLambda()         const { return fLambda; }
    74   const Double_t GetLambdaCheck()     const { return fLambdaCheck; }
    75   const Double_t GetMu0()            const { return fMu0; }
    76   const Double_t GetMu1()            const { return fMu1; }
    77   const Double_t GetSigma0()         const { return fSigma0; }
    78   const Double_t GetSigma1()         const { return fSigma1; }
    79 
    80   const Double_t GetLambdaErr()      const { return fLambdaErr; }
    81   const Double_t GetLambdaCheckErr()  const { return fLambdaCheckErr; }
    82   const Double_t GetMu0Err()         const { return fMu0Err; }
    83   const Double_t GetMu1Err()         const { return fMu1Err; }
    84   const Double_t GetSigma0Err()      const { return fSigma0Err; }
    85   const Double_t GetSigma1Err()      const { return fSigma1Err; }
    86 
    87   const Double_t GetChiSquare()      const { return fChisquare; }
    88   const Double_t GetProb()         const { return fProb;      } 
    89   const Int_t    GetNdf()          const { return fNdf;       }   
    90 
    91   const Double_t GetMeanTime()      const { return fMeanTime; }
    92   const Double_t GetMeanTimeErr()    const { return fMeanTimeErr; }
    93   const Double_t GetSigmaTime()      const { return fSigmaTime; }
    94   const Double_t GetSigmaTimeErr()    const { return fSigmaTimeErr; }
    95 
    96   const Bool_t IsFitOK()                { return fFitOK; }
    97 
    98   const TH1I *GetHBlindPixelChargevsN()  const  { return fHBlindPixelChargevsN; }
    99   const TH1F *GetHBlindPixelPSD()      const  { return fHBlindPixelPSD;  }
     73  const Double_t GetLambda()         const { return fLambda;         }
     74  const Double_t GetLambdaCheck()    const { return fLambdaCheck;    }
     75  const Double_t GetMu0()            const { return fMu0;            }
     76  const Double_t GetMu1()            const { return fMu1;            }
     77  const Double_t GetSigma0()         const { return fSigma0;         }
     78  const Double_t GetSigma1()         const { return fSigma1;         }
     79
     80  const Double_t GetLambdaErr()      const { return fLambdaErr;      }
     81  const Double_t GetLambdaCheckErr() const { return fLambdaCheckErr; }
     82  const Double_t GetMu0Err()         const { return fMu0Err;         }
     83  const Double_t GetMu1Err()         const { return fMu1Err;         }
     84  const Double_t GetSigma0Err()      const { return fSigma0Err;      }
     85  const Double_t GetSigma1Err()      const { return fSigma1Err;      }
     86
     87  const Double_t GetChiSquare()      const { return fChisquare;      }
     88  const Double_t GetProb()           const { return fProb;           } 
     89  const Int_t    GetNdf()            const { return fNdf;            }   
     90
     91  const Double_t GetMeanTime()       const { return fMeanTime;      }
     92  const Double_t GetMeanTimeErr()    const { return fMeanTimeErr;    }
     93  const Double_t GetSigmaTime()      const { return fSigmaTime;      }
     94  const Double_t GetSigmaTimeErr()   const { return fSigmaTimeErr;  }
     95
     96  const Bool_t IsFitOK()             const { return fFitOK;          }
     97
     98  const TH1I *GetHBlindPixelChargevsN() const  { return fHBlindPixelChargevsN; }
     99  const TH1F *GetHBlindPixelPSD()       const  { return fHBlindPixelPSD;       }
    100100 
    101101  // Draws
    … …  
    131131  inline static Double_t fFitFuncMichele(Double_t *x, Double_t *par)
    132132    {
    133 
    134       const Double_t qe = 0.247;
    135 
    136       Double_t lambda = par[0]; 
    137 
    138       Double_t b1 = par[1];
    139       Double_t b2 = par[2];
    140      
    141       Double_t sum = 0.;
     133     
     134      Double_t lambda1cat = par[0]; 
     135      Double_t lambda1dyn = par[1];
     136      Double_t mu0        = par[2];
     137      Double_t mu1cat     = par[3];
     138      Double_t mu1dyn     = par[4];
     139      Double_t sigma0     = par[5];
     140      Double_t sigma1cat  = par[6];
     141      Double_t sigma1dyn  = par[7];
     142       
     143      Double_t sumcat = 0.;
     144      Double_t sumdyn = 0.;
    142145      Double_t arg = 0.;
    143146     
    144       return TMath::Exp(-1.*lambda)*par[5]*sum;
    145      
     147      if (mu1cat    < mu0)
     148        return fNoWay;
     149
     150      if (sigma1cat < sigma0)
     151        return fNoWay;
     152
     153      // if (sigma1cat < sigma1dyn)
     154      // return NoWay;
     155
     156      //if (mu1cat < mu1dyn)
     157      // return NoWay;
     158
     159      //      if (lambda1cat < lambda1dyn)
     160      // return NoWay;
     161
     162      Double_t mu2cat = (2.*mu1cat)-mu0; 
     163      Double_t mu2dyn = (2.*mu1dyn)-mu0; 
     164      Double_t mu3cat = (3.*mu1cat)-(2.*mu0);
     165      Double_t mu3dyn = (3.*mu1dyn)-(2.*mu0);
     166     
     167      Double_t sigma2cat = TMath::Sqrt((2.*sigma1cat*sigma1cat) - (sigma0*sigma0)); 
     168      Double_t sigma2dyn = TMath::Sqrt((2.*sigma1dyn*sigma1dyn) - (sigma0*sigma0)); 
     169      Double_t sigma3cat = TMath::Sqrt((3.*sigma1cat*sigma1cat) - (2.*sigma0*sigma0));
     170      Double_t sigma3dyn = TMath::Sqrt((3.*sigma1dyn*sigma1dyn) - (2.*sigma0*sigma0));
     171     
     172      Double_t lambda2cat = lambda1cat*lambda1cat;
     173      Double_t lambda2dyn = lambda1dyn*lambda1dyn;
     174      Double_t lambda3cat = lambda2cat*lambda1cat;
     175      Double_t lambda3dyn = lambda2dyn*lambda1dyn;
     176
     177     // k=0:
     178      arg = (x[0] - mu0)/sigma0;
     179      sumcat = TMath::Exp(-0.5*arg*arg)/sigma0;
     180      sumdyn =sumcat;
     181
     182      // k=1cat:
     183      arg = (x[0] - mu1cat)/sigma1cat;
     184      sumcat += lambda1cat*TMath::Exp(-0.5*arg*arg)/sigma1cat;
     185      // k=1dyn:
     186      arg = (x[0] - mu1dyn)/sigma1dyn;
     187      sumdyn += lambda1dyn*TMath::Exp(-0.5*arg*arg)/sigma1dyn;
     188     
     189      // k=2cat:
     190      arg = (x[0] - mu2cat)/sigma2cat;
     191      sumcat += 0.5*lambda2cat*TMath::Exp(-0.5*arg*arg)/sigma2cat;
     192      // k=2dyn:
     193      arg = (x[0] - mu2dyn)/sigma2dyn;
     194      sumdyn += 0.5*lambda2dyn*TMath::Exp(-0.5*arg*arg)/sigma2dyn;
     195 
     196     
     197      // k=3cat:
     198      arg = (x[0] - mu3cat)/sigma3cat;
     199      sumcat += 0.1666666667*lambda3cat*TMath::Exp(-0.5*arg*arg)/sigma3cat;
     200      // k=3dyn:
     201      arg = (x[0] - mu3dyn)/sigma3dyn;
     202      sumdyn += 0.1666666667*lambda3dyn*TMath::Exp(-0.5*arg*arg)/sigma3dyn; 
     203   
     204      sumcat = TMath::Exp(-1.*lambda1cat)*sumcat;
     205      sumdyn = TMath::Exp(-1.*lambda1dyn)*sumdyn;
     206     
     207      return par[8]*(sumcat+sumdyn)/2.;
     208
    146209    }
    147210   
    … …  
    339402     
    340403    }
     404
     405  inline static Double_t fPolya(Double_t *x, Double_t *par)
     406    {
     407
     408      const Double_t QEcat = 0.247;            // mean quantum efficiency
     409      const Double_t sqrt2 = 1.4142135623731;
     410      const Double_t sqrt3 = 1.7320508075689;
     411      const Double_t sqrt4 = 2.;
     412     
     413      const Double_t lambda = par[0];           // mean number of photons
     414     
     415      const Double_t excessPoisson = par[1];    // non-Poissonic noise contribution
     416      const Double_t delta1 = par[2];           // amplification first dynode
     417      const Double_t delta2 = par[3];           // amplification subsequent dynodes
     418     
     419      const Double_t electronicAmpl = par[4];   // electronic amplification and conversion to FADC charges
     420
     421      const Double_t pmtAmpl = delta1*delta2*delta2*delta2*delta2*delta2;  // total PMT gain
     422      const Double_t A = 1. + excessPoisson - QEcat                       
     423        + 1./delta1
     424                + 1./delta1/delta2
     425        + 1./delta1/delta2/delta2;                                  // variance contributions from PMT and QE
     426     
     427      const Double_t totAmpl = QEcat*pmtAmpl*electronicAmpl;        // Total gain and conversion
     428     
     429      const Double_t mu0 = par[7];                                      // pedestal
     430      const Double_t mu1 = totAmpl;                                 // single phe position
     431      const Double_t mu2 = 2*totAmpl;                               // double phe position
     432      const Double_t mu3 = 3*totAmpl;                               // triple phe position
     433      const Double_t mu4 = 4*totAmpl;                               // quadruple phe position
     434     
     435      const Double_t sigma0 = par[5];
     436      const Double_t sigma1 = electronicAmpl*pmtAmpl*TMath::Sqrt(QEcat*A);
     437      const Double_t sigma2 = sqrt2*sigma1;
     438      const Double_t sigma3 = sqrt3*sigma1;
     439      const Double_t sigma4 = sqrt4*sigma1;
     440     
     441      const Double_t lambda2 = lambda*lambda;
     442      const Double_t lambda3 = lambda2*lambda;
     443      const Double_t lambda4 = lambda3*lambda;
     444     
     445      //-- calculate the area----
     446      Double_t arg = (x[0] - mu0)/sigma0;
     447      Double_t sum = TMath::Exp(-0.5*arg*arg)/sigma0;
     448     
     449     // k=1:
     450      arg = (x[0] - mu1)/sigma1;
     451      sum += lambda*TMath::Exp(-0.5*arg*arg)/sigma1;
     452     
     453      // k=2:
     454      arg = (x[0] - mu2)/sigma2;
     455      sum += 0.5*lambda2*TMath::Exp(-0.5*arg*arg)/sigma2;
     456     
     457      // k=3:
     458      arg = (x[0] - mu3)/sigma3;
     459      sum += 0.1666666667*lambda3*TMath::Exp(-0.5*arg*arg)/sigma3;
     460     
     461      // k=4:
     462      arg = (x[0] - mu4)/sigma4;
     463      sum += 0.041666666666667*lambda4*TMath::Exp(-0.5*arg*arg)/sigma4;     
     464     
     465      return TMath::Exp(-1.*lambda)*par[6]*sum;
     466    }
     467 
     468
    341469 
    342470  ClassDef(MHCalibrationBlindPixel, 1)  // Histograms from the Calibration Blind Pixel
Note: See TracChangeset for help on using the changeset viewer.