Changeset 2599 for trunk


Ignore:
Timestamp:
12/04/03 12:17:22 (21 years ago)
Author:
gaug
Message:
*** empty log message ***
Location:
trunk/MagicSoft/Mars
Files:
14 edited

Legend:

Unmodified
Added
Removed

  • trunk/MagicSoft/Mars/Changelog

    r2596 r2599  
    44
    55                                                 -*-*- END OF LINE -*-*-
     6 2003/12/04: Markus Gaug
     7
     8  * manalysis/MCalibration*
     9     - implemented some of Thomas Bretz suggestions to make the code nicer
     10     - implemented the possibility to have cosmics in the calibration data
     11       and still fit it
     12     - conversion HiGain/LoGain now left out (it seems to be already done in
     13        the data)
     14
     15  * mhist/MHCalibration*
     16     - implemented some of Thomas Bretz suggestions to make the code nicer
     17     - implemented the possibility to have cosmics in the calibration data
     18       and still fit it
     19
     20  * macros/calibration.C
     21     - MStatusDisplay of calibration histograms a little bit more readable
     22
    623 2003/12/03: Abelardo Moralejo
    724

  • trunk/MagicSoft/Mars/macros/calibration.C

    r2581 r2599  
    6262    tlist.AddToList(&fill);
    6363
    64     MStatusDisplay *d1 = new MStatusDisplay;
     64    //    MStatusDisplay *d1 = new MStatusDisplay;
    6565 
    6666    // Set update time to 3s
    67     d1->SetUpdateTime(3000);
     67    //    d1->SetUpdateTime(3000);
    6868 
    6969    //
    … …  
    7272    MEvtLoop evtloop;
    7373    evtloop.SetParList(&plist);
    74     evtloop.SetDisplay(d1);     
     74    //    evtloop.SetDisplay(d1);       
    7575
    7676    //
    … …  
    8282    tlist.PrintStatistics();
    8383
    84     MPedestalCam *ped = plist.FindObject("MPedestalCam");
    85     ped.Print();
    8684
    8785    //
    … …  
    122120
    123121    MHCamEvent hist2;
    124     hist2.SetType(0);
     122    hist2.SetType(8);
    125123    plist2.AddToList(&hist2);
    126124    MFillH fill2("MHCamEvent", "MCalibrationCam");
    … …  
    145143    //
    146144    MCalibrationCam *cam = plist2.FindObject("MCalibrationCam");
    147     MCalibrationPix *pix = cam->GetCalibrationPix(523);
    148     pix->Draw();
    149 
     145    cam.Print();
     146
     147    Int_t pixnr;
     148
     149    while (1)
     150      {
     151       
     152        cout << "Which pixel number do you want to display? (Press 0 to exit)" << endl;
     153        cin >> pixnr;
     154        if (pixnr == 0)
     155          break;
     156       
     157        if (pixnr >= 577)
     158          break;
     159       
     160        MCalibrationPix *pix = cam->GetCalibrationPix(pixnr);
     161        pix->Draw();
     162   
     163      }
     164   
    150165    //
    151166    // Here we are confronted to a serious bug in ROOT:
    … …  
    170185    MHCamera disp11  (geomcam, "MCalibrationCam;rq", "Reduced Charges");
    171186    MHCamera disp12  (geomcam, "MCalibrationCam;errrq", "Error of Reduced Charges");
     187    MHCamera disp13  (geomcam, "MCalibrationCam;phe", "Nr. of Phe's (F-Factor Method)");
     188    MHCamera disp14  (geomcam, "MCalibrationCam;convphe", "Conversion Factor (F-Factor Method)");
    172189
    173190    disp1.SetCamContent(*cam, 0);
    … …  
    183200    disp11.SetCamContent(*cam, 10);
    184201    disp12.SetCamContent(*cam, 11);
     202    disp12.SetCamContent(*cam, 12);
     203    disp13.SetCamContent(*cam, 13);
    185204
    186205    disp1.SetYTitle("Q [FADC counts]");
    … …  
    196215    disp11.SetYTitle("Q [FADC counts]");
    197216    disp12.SetYTitle("\\Delta_{Q} [FADC counts]");
     217    disp13.SetYTitle("Nr Phe's");
     218    disp14.SetYTitle("Conversion Factor [Phe/FADC count]");
    198219
    199220    MStatusDisplay *d2 = new MStatusDisplay;
    … …  
    202223    d2->SetUpdateTime(1000);
    203224
    204     TCanvas *c1 = &d2->AddTab("Fitted Charges");
    205     c1->Divide(5, 2);
     225    TCanvas *c1 = &d2->AddTab("Charges Mean");
     226    c1->Divide(2, 2);
    206227
    207228    TObject *obj;
    … …  
    211232    obj=disp1.DrawCopy("hist");
    212233
    213     c1->cd(6);
     234    c1->cd(3);
    214235    gPad->SetBorderMode(0);
    215236    obj->Draw();
    … …  
    219240    obj=disp2.DrawCopy("hist");
    220241
    221     c1->cd(7);
    222     gPad->SetBorderMode(0);
    223     obj->Draw();
    224 
    225     c1->cd(3);
     242    c1->cd(4);
     243    gPad->SetBorderMode(0);
     244    obj->Draw();
     245
     246    TCanvas *c11 = &d2->AddTab("Charges Sigma");
     247    c11->Divide(2, 2);
     248
     249    c11->cd(1);
    226250    gStyle->SetOptStat(1101);
    227251    obj=disp3.DrawCopy("hist");
    228252
    229     c1->cd(8);
    230     gPad->SetBorderMode(0);
    231     obj->Draw();
    232 
    233     c1->cd(4);
     253    c11->cd(3);
     254    gPad->SetBorderMode(0);
     255    obj->Draw();
     256
     257    c11->cd(2);
    234258    gStyle->SetOptStat(1101);
    235259    obj=disp4.DrawCopy("hist");
    236260
    237     c1->cd(9);
    238     gPad->SetBorderMode(0);
    239     obj->Draw();
    240 
    241     c1->cd(5);
     261    c11->cd(4);
     262    gPad->SetBorderMode(0);
     263    obj->Draw();
     264
     265
     266    TCanvas *c12 = &d2->AddTab("Fit Prob.");
     267    c12->Divide(1, 2);
     268
     269    c12->cd(1);
    242270    gStyle->SetOptStat(1101);
    243271    obj=disp5.DrawCopy("hist");
    244272
    245     c1->cd(10);
     273    c12->cd(2);
    246274    gPad->SetBorderMode(0);
    247275    obj->Draw();
    … …  
    282310
    283311    c3->cd(2);
    284     gStyle->SetOptStat(1101);
     312    gStyle->SetOptStat(1111);
    285313    obj=disp10.DrawCopy("hist");
    286314
    … …  
    305333    obj->Draw();
    306334
     335    TCanvas *c5 = &d2->AddTab("F-Factor Method");
     336    c5->Divide(2, 2);
     337
     338    c5->cd(1);
     339    gStyle->SetOptStat(1111);
     340    obj=disp13.DrawCopy("hist");
     341
     342    c5->cd(3);
     343    obj->Draw();
     344
     345    c5->cd(2);
     346    gStyle->SetOptStat(1101);
     347    obj=disp14.DrawCopy("hist");
     348
     349    c5->cd(4);
     350    obj->Draw();
     351
     352
    307353#endif
    308354

  • trunk/MagicSoft/Mars/manalysis/MCalibrationBlindPix.h

    r2525 r2599  
    5151  Float_t GetErrT()      const    { return fErrT;      }
    5252 
    53   Bool_t FillQ(Int_t q)            { return fHist->FillBPQ(q); }
    54   Bool_t FillT(Int_t t)            { return fHist->FillBPT(t); } 
    55   Bool_t FillRQvsT(Float_t rq, Int_t t) { return fHist->FillBPQvsN(rq,t); }   
     53  Bool_t FillQ(Int_t q)            { return fHist->FillBlindPixelQ(q); }
     54  Bool_t FillT(Int_t t)            { return fHist->FillBlindPixelT(t); } 
     55  Bool_t FillRQvsT(Float_t rq, Int_t t) { return fHist->FillBlindPixelQvsN(rq,t); }   
    5656 
    5757  Bool_t IsValid()                 { return fLambda > 0. || fErrLambda > 0.; }

  • trunk/MagicSoft/Mars/manalysis/MCalibrationCalc.cc

    r2581 r2599  
    163163      }
    164164
    165     //
    166     // fRawEvt->GetNumPixels() does not work !!!!!!!!!!
    167     //
    168     fCalibrations->InitSize(577);
    169165
    170166    switch (fColor)
    … …  
    189185      }
    190186   
    191 
    192     // MTime does not work!!!!
    193 
    194     //    fEvtTime = (MTime*)pList->FindObject("MRawEvtTime");
    195     //    if (!fEvtTime)
    196     //      {
    197     //        *fLog << dbginf << "MTime could not be created ... ¡!¡!¡!¡!" << endl;       
    198     //      }
    199    
    200     // fEvtTime->SetTime(0.);
    201 
    202187   
    203188    return kTRUE;
    … …  
    222207    fNumHiGainSamples =  fRunHeader->GetNumSamplesHiGain();
    223208    fNumLoGainSamples =  fRunHeader->GetNumSamplesLoGain();
    224  
     209
     210    //
     211    // FIXME: The next statement simply does not work:
     212    //        fRawEvt->GetNumPixels() returns always 0
     213    // fCalibrations->InitSize(fRawEvt->GetNumPixels());
     214    //
     215
     216    fCalibrations->InitSize(577);   
     217
     218    for (Int_t i=0;i<577;i++)
     219      {
     220        MCalibrationPix &pix = (*fCalibrations)[i];
     221        pix.ChangePixId(i);
     222      }
     223   
    225224    return kTRUE;
    226225 
    … …  
    238237    fEvents++;
    239238
    240     Bool_t cosmic = kFALSE;
     239    Int_t cosmicpix = 0;
    241240
    242241    MCalibrationBlindPix &blindpixel = *(fCalibrations->GetBlindPixel());
    … …  
    251250        const Int_t pixid = pixel.GetPixelId();
    252251
    253         if (!fCalibrations->IsPixelUsed(pixid))
    254           if (!fCalibrations->AddPixel(pixid))
    255             *fLog << err << dbginf << "MCalibrationPix(" << pixid << ") could not be created !!" << endl;
    256        
    257252        Byte_t *ptr = pixel.GetHiGainSamples();
    258253        Byte_t mid  = pixel.GetIdxMaxHiGainSample();
    … …  
    270265            mid = pixel.GetIdxMaxLoGainSample();
    271266
    272             sum = (max == gkSaturationLimit              // overflow of LoGain ??? -> GimmeABreak!!!
    273                    ? fHistOverFlow++, gkLoGainOverFlow   // OUCH (Florian was maybe right)
    274                    : sum*gkConversionHiLo    );          // OUFF (Florian was wrong) !!
    275 
    276             //            *fLog << warn << "Warning: Saturation of HiGain reached in slice " << (int)mid << " !!! " << endl;
    277             //            *fLog << warn << "Max = " << max << endl;
     267            //
     268            //  FIXME: It seems the conversion HiGain LoGain is already
     269            //         performed in the data?!?
     270            //
     271            sum = (max > gkSaturationLimit              // overflow of LoGain ??? -> GimmeABreak!!!
     272                   ? fHistOverFlow++, gkLoGainOverFlow  // OUCH (Florian was maybe right)
     273                   : sum  );                            // OUFF (Florian was wrong) !!
     274            //                   : sum*gkConversionHiLo    );          // OUFF (Florian was wrong) !!
    278275
    279276            if (fHistOverFlow)
    … …  
    282279
    283280          }
    284        
    285281
    286282        MPedestalPix    &ped = (*fPedestals)[pixid];
    … …  
    298294
    299295        //
    300         // This is a very primitive check for the number of cosmics
     296        // This is a very primitive check for the number of cosmicpixs
    301297        // The cut will be applied in the fit, but for the blind pixel,
    302298        // we need to remove this event
    … …  
    305301        //
    306302
    307         if ((float)sum < pedes+5.*pedrms)
    308            cosmic = kTRUE;
     303        if ((float)sum < pedes+4.*pedrms)
     304           cosmicpix++;
    309305
    310306        Float_t rsum      = (float)sum - pedes;
    … …  
    314310           
    315311          case gkCalibrationBlindPixelId:
     312
    316313            //
    317314            // FIXME: This works only when the blind pixel ID is much larger than
    318315            //        the rest of the pixels (which is the case right now)
    319316            //
    320 //          if (!cosmic)
    321                if (!blindpixel.FillQ(sum))
     317            if (cosmicpix < 100.)
     318              {
     319                if (!blindpixel.FillQ(sum))
    322320                  *fLog << warn <<
    323                   "Overflow or Underflow occurred filling Blind Pixel sum = " << sum << endl;
    324 
    325             if (!blindpixel.FillT((int)mid))
    326               *fLog << warn <<
    327                 "Overflow or Underflow occurred filling Blind Pixel time = " << (int)mid << endl;
    328 
    329             if (!blindpixel.FillRQvsT(rsum,fEvents))
    330               *fLog << warn <<
    331                 "Overflow or Underflow occurred filling Blind Pixel eventnr = " << fEvents << endl;
    332 
     321                    "Overflow or Underflow occurred filling Blind Pixel sum = " << sum << endl;
     322
     323                if (!blindpixel.FillT((int)mid))
     324                  *fLog << warn <<
     325                    "Overflow or Underflow occurred filling Blind Pixel time = " << (int)mid << endl;
     326
     327                if (!blindpixel.FillRQvsT(rsum,fEvents))
     328                  *fLog << warn <<
     329                    "Overflow or Underflow occurred filling Blind Pixel eventnr = " << fEvents << endl;
     330              }
     331           
    333332          case gkCalibrationPINDiodeId:
    334333            if (!pindiode.FillQ(sum))
    … …  
    364363      } /* while (pixel.Next()) */
    365364
    366     if (cosmic)
     365    if (cosmicpix > 300.)
    367366        fCosmics++;
    368367
    369     fCalibrations->SetReadyToSave();
    370    
    371368    return kTRUE;
    372369}
    … …  
    405402
    406403  *fLog << GetDescriptor() << " Fitting the Normal Pixels" << endl;
     404
    407405  //
    408406  // loop over the pedestal events and check if we have calibration
    … …  
    411409    {
    412410
    413       if (fCalibrations->IsPixelUsed(pixid))
    414         {
    415 
    416           MCalibrationPix &pix = (*fCalibrations)[pixid];
    417 
    418           const Float_t ped    = (*fPedestals)[pixid].GetPedestal() * fNumHiGainSamples;
    419           const Float_t prms   = (*fPedestals)[pixid].GetPedestalRms() * fNumHiGainSamples;
    420 
    421           pix.SetPedestal(ped,prms);
    422 
    423           if (TESTBIT(fFlags,kUseTFits))
    424             pix.FitT();
    425 
    426           if (!pix.FitQ())
    427              continue;
    428 
    429         }
     411      MCalibrationPix &pix = (*fCalibrations)[pixid];
     412
     413      const Float_t ped    = (*fPedestals)[pixid].GetPedestal() * fNumHiGainSamples;
     414      const Float_t prms   = (*fPedestals)[pixid].GetPedestalRms() * fNumHiGainSamples;
     415
     416      pix.SetPedestal(ped,prms);
     417
     418      if (TESTBIT(fFlags,kUseTFits))
     419        pix.FitT();
     420     
     421      pix.FitQ();
    430422    }
    431423

  • trunk/MagicSoft/Mars/manalysis/MCalibrationCam.cc

    r2582 r2599  
    6161    fTitle = title ? title : "Storage container for the Calibration Information in the camera";
    6262
    63     fPixels     = new TClonesArray("MCalibrationPix",0);
     63    fPixels     = new TClonesArray("MCalibrationPix",1);
    6464    fBlindPixel = new MCalibrationBlindPix();
    6565    fPINDiode   = new MCalibrationPINDiode();
    … …  
    9292
    9393  //
    94   // Here it is important to use GetEntriesFast,
    95   // We get the array index of the last "filled" fPixel entry
    96   // (Not the number of "filled" fPixels!!)
     94  // Here we get the number of "filled" fPixels!!
    9795  //
    9896  return fPixels->GetEntriesFast();
    … …  
    116114    return *static_cast<MCalibrationPix*>(fPixels->UncheckedAt(i));
    117115}
    118 
    119 Bool_t MCalibrationCam::AddPixel(Int_t idx)
    120 {
    121 
    122   //
    123   // Check bounds first
    124   //
    125   if (!CheckBounds(idx))
    126     InitSize(idx);
    127 
    128   //
    129   // in case, new cannot allocate memory,
    130   // return FALSE
    131   //
    132   return (new ((*fPixels)[idx]) MCalibrationPix(idx));
    133 
    134 }
    135 
    136116
    137117// --------------------------------------------------------------------------
    … …  
    173153Bool_t MCalibrationCam::IsPixelFitted(Int_t idx) const
    174154{
    175   if (!IsPixelUsed(idx))
    176     return kFALSE;
    177 
    178155  return ((*this)[idx].GetRQ() > 0. && (*this)[idx].GetErrRQ() > 0.);
    179156}
    … …  
    208185    {
    209186
    210       Int_t id = 0;
    211      
    212187      TIter Next(fPixels);
    213188      MCalibrationPix *pix;
    214189      while ((pix=(MCalibrationPix*)Next()))
    215190        {
    216           if (!IsPixelUsed(id++))
    217             continue;
    218191          if (pix->FitQ())
    219192            nsuccess++;
    … …  
    222195  else                  // fit only the pixel with index i
    223196    {
    224       if (!IsPixelUsed(i))
    225         return 0;
    226197      if((*this)[i].FitQ())
    227198        nsuccess++;
    … …  
    247218  UShort_t nsuccess = 0;
    248219
    249   Int_t id = 0;
    250 
    251220  TIter Next(fPixels);
    252221  MCalibrationPix *pix;
    253222  while ((pix=(MCalibrationPix*)Next()))
    254223    {
    255       if (!IsPixelUsed(id++))
    256         continue;
    257224      if (pix->FitQ())
    258225        nsuccess++;
    … …  
    293260    {
    294261
    295       Int_t id = 0;
    296      
    297262      TIter Next(fPixels);
    298263      MCalibrationPix *pix;
    299264      while ((pix=(MCalibrationPix*)Next()))
    300265        {
    301           if (!IsPixelUsed(id++))
    302             continue;
    303           if (pix->FitT())
     266         if (pix->FitT())
    304267            nsuccess++;
    305268        }
    … …  
    331294  UShort_t nsuccess = 0;
    332295
    333   Int_t id = 0;
    334 
    335296  TIter Next(fPixels);
    336297  MCalibrationPix *pix;
    337298  while ((pix=(MCalibrationPix*)Next()))
    338299    {
    339       if (!IsPixelUsed(id++))
    340         continue;
    341300      if (pix->FitT())
    342301        nsuccess++;
    … …  
    391350  // we want to keep all pixel not used with a NULL pointer.
    392351  //
    393   fPixels->Expand(size);
    394   //
    395   // Set all entries to the null pointer. 
    396   // Later we fill the array per pixId with the contruction: new (fPixels[i]) MCalibrationPix(pixid)
    397   // To check, if pixels is already filled, we test the NULL pointer (see IsPixelUsed)
    398   //
    399   for (Int_t i=0; i< size; i++)
    400     {
    401       MCalibrationPix *pix = &(*this)[i];
    402       pix = NULL;
    403     }
     352  fPixels->ExpandCreate(size);
     353
    404354}
    405355 
    … …  
    413363    while ((pix=(MCalibrationPix*)Next()))
    414364    {
    415         if (!IsPixelUsed(id))
    416           continue;
    417 
    418         *fLog << id << ": ";
     365
     366        *fLog << id << ": " << pix->GetPed() << " " << pix->GetPedRms() << " Charges: " ;
    419367        *fLog << pix->GetQ() << " " << pix->GetRQ() << endl;
    420368
    … …  
    426374Bool_t MCalibrationCam::GetPixelContent(Double_t &val, Int_t idx, const MGeomCam &cam, Int_t type) const
    427375{
    428 
    429   if (!IsPixelFitted(idx))
    430     return kFALSE;
    431376
    432377  switch (type)
    … …  
    468413      val = (*this)[idx].GetErrRQ();
    469414      break;
     415    case 12:
     416      val = (*this)[idx].GetPheFFactorMethod();
     417      break;
     418    case 13:
     419      val = (*this)[idx].GetConversionFFactorMethod();
     420      break;
    470421    default:
    471422      return kFALSE;
    … …  
    474425}
    475426
    476 void MCalibrationCam::DrawPixelContent(Int_t num) const
    477 {
    478     *fLog << warn << "MCalibrationCam::DrawPixelContent - not available." << endl;
     427void MCalibrationCam::DrawPixelContent(Int_t idx) const
     428{
     429
     430  (*this)[idx].Draw();
     431
    479432}
    480433

  • trunk/MagicSoft/Mars/manalysis/MCalibrationCam.h

    r2525 r2599  
    7777  Bool_t CheckBounds(Int_t i) const;
    7878
    79   Bool_t AddPixel(Int_t idx);
    80 
    8179  void Print(Option_t *o="") const;
    8280 

  • trunk/MagicSoft/Mars/manalysis/MCalibrationPix.cc

    r2581 r2599  
    4040using namespace std;
    4141
    42 static const TString gsDefHTitle = "Calibration Histograms Pixel ";
    4342// --------------------------------------------------------------------------
    4443//
    4544// Default Constructor.
    4645//
    47 MCalibrationPix::MCalibrationPix(Int_t pix, const char *name, const char *title)
    48     : fPixId(pix)
     46MCalibrationPix::MCalibrationPix(const char *name, const char *title)
     47    : fPixId(-1),
     48      fQ(-1.),
     49      fErrQ(-1.),
     50      fSigmaQ(-1.),
     51      fErrSigmaQ(-1.),
     52      fQProb(-1.),
     53      fPed(-1.),
     54      fPedRms(-1.),
     55      fT(-1.),
     56      fSigmaT(-1.),
     57      fTProb(-1.),
     58      fRQ(-1.),
     59      fErrRQ(-1.),
     60      fFactor(1.3),
     61      fPheFFactorMethod(-1.),
     62      fConversionFFactorMethod(-1.)
    4963{
    5064
    … …  
    5266  fTitle = title ? title : "Container of the MHCalibrationPixels and the fit results";
    5367
    54   fHist = new MHCalibrationPixel(fPixId,"MHCalibrationPixel",gsDefHTitle.Data()+fPixId);
    55 
    56   fQ   = fErrQ     = 0.;
    57   fPed = fPedRms   = 0.;
    58   fT   = fSigmaT  = 0.;
    59   fRQ  = fErrRQ = 0.;
    60   fSigmaQ = fErrSigmaQ = 0.;
    61   fQProb  = 0.;
    62 
     68  fHist = new MHCalibrationPixel("MHCalibrationPixel","Calibration Histograms Pixel");
    6369}
    6470
    … …  
    6773  delete fHist;
    6874}
     75
     76
     77void MCalibrationPix::ChangePixId(Int_t i)
     78{
     79 
     80  fPixId = i;
     81  fHist->ChangeHistId(i);
     82 
     83}
     84
    6985
    7086// ------------------------------------------------------------------------
    … …  
    84100
    85101  if (fPed && fPedRms)
    86     fHist->SetLowerFitRange(fPed + 3.5*fPedRms);
     102    fHist->SetLowerFitRange(fPed + 2.0*fPedRms);
    87103  else
    88104    *fLog << warn << "Cannot set lower fit range to suppress cosmics: Pedestals not available" << endl;
    … …  
    101117  fQProb     = fHist->GetQProb();
    102118
    103   if (fPed)
     119  if ((fPed > 0.)  && (fPedRms > 0.))
     120    {
     121     
    104122    fRQ      = fQ - fPed;
    105   if (fPedRms)
    106123    fErrRQ   = TMath::Sqrt(fErrQ*fErrQ + fPedRms*fPedRms);
    107  
     124
     125    fPheFFactorMethod =
     126      fFactor
     127      * fRQ * fRQ
     128      / (fSigmaQ * fSigmaQ - fPedRms*fPedRms) ;
     129
     130    fConversionFFactorMethod = fPheFFactorMethod / fRQ ;
     131   
     132    }
    108133
    109134  return kTRUE;
    … …  
    142167}
    143168
    144 void MCalibrationPix::Test()
    145 {
    146   *fLog << "TEST: pixid: " << fPixId << endl; 
    147 }

  • trunk/MagicSoft/Mars/manalysis/MCalibrationPix.h

    r2544 r2599  
    1212private:
    1313
    14   //
    15   // FIXME: Derive class from MCerphotPix ??
    16   //
    1714  Int_t   fPixId;           // the pixel Id
    1815 
    … …  
    3229  Float_t fRQ;               // The reduced mean charge after the fit
    3330  Float_t fErrRQ;            // The error of the reduced mean charge after the fit 
    34  
     31
     32  Float_t fFactor;            // The F-factor
     33  Float_t fPheFFactorMethod;  // The number of Phe's calculated after the F-factor method
     34  Float_t fConversionFFactorMethod; // The conversion factor to Phe's calculated after the F-factor method
     35   
    3536  MHCalibrationPixel *fHist; // Pointer to the histograms performing the fits, etc. 
    3637 
    3738public:
    3839
    39   MCalibrationPix(Int_t pix=-1, const char *name=NULL, const char *title=NULL);
     40  MCalibrationPix(const char *name=NULL, const char *title=NULL);
    4041  ~MCalibrationPix();
    4142 
    … …  
    6667  Bool_t IsValid()      const           { return fRQ >=0 || fErrRQ >= 0; }
    6768  Int_t  GetPixId()     const           { return fPixId;   }
    68 
     69  void ChangePixId(Int_t i);
     70 
    6971  Bool_t FitQ();
    7072  Bool_t FitT();
    … …  
    7375  virtual void Draw(Option_t *opt="")     { fHist->Draw(opt); }
    7476 
    75   void Test();
    76    
     77  Float_t GetPheFFactorMethod() const        { return fPheFFactorMethod;  } 
     78  Float_t GetConversionFFactorMethod() const { return fConversionFFactorMethod;  }
     79 
    7780  ClassDef(MCalibrationPix, 1)  // Storage Container for Calibration information of one pixel
    7881};

  • trunk/MagicSoft/Mars/manalysis/MPedCalcPedRun.cc

    r2536 r2599  
    6060
    6161    AddToBranchList("fHiGainPixId");
    62     //AddToBranchList("fLoGainPixId");
    6362    AddToBranchList("fHiGainFadcSamples");
    64     //AddToBranchList("fLoGainFadcSamples");
    6563}
    6664
    … …  
    9896    fNumHiGainSamples =  fRunheader->GetNumSamplesHiGain();
    9997
     98    fPedestals->InitSize(fRunheader->GetNumPixel());
     99
    100100    return kTRUE;
    101101
    … …  
    105105Int_t MPedCalcPedRun::Process()
    106106{
     107
    107108    MRawEvtPixelIter pixel(fRawEvt);
    108109
    109110    while (pixel.Next())
    110111    {
    111               Byte_t *ptr = pixel.GetHiGainSamples();
     112             Byte_t *ptr = pixel.GetHiGainSamples();
    112113        const Byte_t *end = ptr + fRawEvt->GetNumHiGainSamples();
    113114
    … …  
    115116        const Float_t higainrms = CalcHiGainRms(ptr, end, higainped);
    116117
    117         //const Float_t higainpederr = CalcHiGainMeanErr(higainrms);
    118         //const Float_t higainrmserr = CalcHiGainRmsErr(higainrms);
    119 
    120118        const UInt_t pixid = pixel.GetPixelId();
    121119        MPedestalPix &pix = (*fPedestals)[pixid];
    122120
    123121        pix.Set(higainped, higainrms);
    124        //pix.SetPedestalRms(higainpederr, higainrmserr);
     122
    125123    }
    126124

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

    r2525 r2599  
    6969
    7070    // Create a large number of bins, later we will rebin
    71     fBPQfirst = 0;
    72     fBPQlast  = gkStartBlindPixelBinNr;
    73     fBPQnbins = gkStartBlindPixelBinNr;
    74 
    75     fHBPQ = new TH1I("HBPQ","Distribution of Summed FADC Slices",fBPQnbins,fBPQfirst,fBPQlast);
    76     fHBPQ->SetXTitle("Sum FADC Slices");
    77     fHBPQ->SetYTitle("Nr. of events");
    78     fHBPQ->Sumw2();
    79 
    80     fErrBPQfirst = 0.;
    81     fErrBPQlast  = gkStartBlindPixelBinNr;
    82     fErrBPQnbins = gkStartBlindPixelBinNr;
    83 
    84     fHBPErrQ = new TH1F("HBPErrQ","Distribution of Variances of Summed FADC Slices",
    85                         fErrBPQnbins,fErrBPQfirst,fErrBPQlast);
    86     fHBPErrQ->SetXTitle("Variance Summed FADC Slices");
    87     fHBPErrQ->SetYTitle("Nr. of events");
    88     fHBPErrQ->Sumw2();
     71    fBlindPixelQfirst = 0;
     72    fBlindPixelQlast  = gkStartBlindPixelBinNr;
     73    fBlindPixelQnbins = gkStartBlindPixelBinNr;
     74
     75    fHBlindPixelQ = new TH1I("HBlindPixelQ","Distribution of Summed FADC Slices",fBlindPixelQnbins,fBlindPixelQfirst,fBlindPixelQlast);
     76    fHBlindPixelQ->SetXTitle("Sum FADC Slices");
     77    fHBlindPixelQ->SetYTitle("Nr. of events");
     78    fHBlindPixelQ->Sumw2();
     79
     80    fErrBlindPixelQfirst = 0.;
     81    fErrBlindPixelQlast  = gkStartBlindPixelBinNr;
     82    fErrBlindPixelQnbins = gkStartBlindPixelBinNr;
     83
     84    fHBlindPixelErrQ = new TH1F("HBlindPixelErrQ","Distribution of Variances of Summed FADC Slices",
     85                        fErrBlindPixelQnbins,fErrBlindPixelQfirst,fErrBlindPixelQlast);
     86    fHBlindPixelErrQ->SetXTitle("Variance Summed FADC Slices");
     87    fHBlindPixelErrQ->SetYTitle("Nr. of events");
     88    fHBlindPixelErrQ->Sumw2();
    8989
    9090    Axis_t tfirst = -0.5;
    … …  
    9292    Int_t nbins   = 16;
    9393
    94     fHBPT = new TH1I("HBPT","Distribution of Mean Arrival Times",nbins,tfirst,tlast);
    95     fHBPT->SetXTitle("Mean Arrival Times [FADC slice nr]");
    96     fHBPT->SetYTitle("Nr. of events");
    97     fHBPT->Sumw2();
     94    fHBlindPixelT = new TH1I("HBlindPixelT","Distribution of Mean Arrival Times",nbins,tfirst,tlast);
     95    fHBlindPixelT->SetXTitle("Mean Arrival Times [FADC slice nr]");
     96    fHBlindPixelT->SetYTitle("Nr. of events");
     97    fHBlindPixelT->Sumw2();
    9898
    9999    // We define a reasonable number and later enlarge it if necessary
    … …  
    102102    Axis_t nlast  = (Axis_t)nbins - 0.5;
    103103
    104     fHBPQvsN = new TH1I("HBPQvsN","Sum of Charges vs. Event Number",nbins,nfirst,nlast);
    105     fHBPQvsN->SetXTitle("Event Nr.");
    106     fHBPQvsN->SetYTitle("Sum of FADC slices");
     104    fHBlindPixelQvsN = new TH1I("HBlindPixelQvsN","Sum of Charges vs. Event Number",nbins,nfirst,nlast);
     105    fHBlindPixelQvsN->SetXTitle("Event Nr.");
     106    fHBlindPixelQvsN->SetYTitle("Sum of FADC slices");
    107107
    108108    fgSinglePheFitFunc = &gfKto8;
    … …  
    113113{
    114114
    115   delete fHBPQ;
    116   delete fHBPT;
    117   delete fHBPErrQ;
     115  delete fHBlindPixelQ;
     116  delete fHBlindPixelT;
     117  delete fHBlindPixelErrQ;
    118118 
    119119  if (fSinglePheFit)
    … …  
    127127void MHCalibrationBlindPixel::ResetBin(Int_t i)
    128128{
    129     fHBPQ->SetBinContent (i, 1.e-20);
    130     fHBPErrQ->SetBinContent  (i, 1.e-20);
    131     fHBPT->SetBinContent(i, 1.e-20);
     129    fHBlindPixelQ->SetBinContent (i, 1.e-20);
     130    fHBlindPixelErrQ->SetBinContent  (i, 1.e-20);
     131    fHBlindPixelT->SetBinContent(i, 1.e-20);
    132132}
    133133
    … …  
    209209    gPad->SetTicks();
    210210
    211     fHBPQ->DrawCopy(opt);
     211    fHBlindPixelQ->DrawCopy(opt);
    212212   
    213213    if (fSinglePheFit)
    … …  
    231231    gPad->SetLogy(1);
    232232    gPad->SetBorderMode(0);
    233     fHBPT->DrawCopy(opt);
    234 
    235     if (fHBPT->GetFunction("GausTime"))
     233    fHBlindPixelT->DrawCopy(opt);
     234
     235    if (fHBlindPixelT->GetFunction("GausTime"))
    236236      {
    237         TF1 *tfit = fHBPT->GetFunction("GausTime");
     237        TF1 *tfit = fHBlindPixelT->GetFunction("GausTime");
    238238        if (tfit->GetProb() < 0.01)
    239239          tfit->SetLineColor(kRed);
    … …  
    248248    c->cd(4);
    249249
    250     fHBPQvsN->DrawCopy(opt);
     250    fHBlindPixelQvsN->DrawCopy(opt);
    251251
    252252    c->Modified();
    … …  
    260260  gRandom->SetSeed();
    261261
    262   if (fHBPQ->GetEntries() != 0)
     262  if (fHBlindPixelQ->GetEntries() != 0)
    263263    {
    264       *fLog << err << "Histogram " << fHBPQ->GetTitle() << " is already filled. " << endl;
     264      *fLog << err << "Histogram " << fHBlindPixelQ->GetTitle() << " is already filled. " << endl;
    265265      *fLog << err << "Create new class MHCalibrationBlindPixel for simulation! " << endl;
    266266      return kFALSE;
    … …  
    268268 
    269269  TF1 *simulateSinglePhe = new TF1("simulateSinglePhe",fgSinglePheFitFunc,
    270                                    fBPQfirst,fBPQlast,fgSinglePheFitNPar);
     270                                   fBlindPixelQfirst,fBlindPixelQlast,fgSinglePheFitNPar);
    271271 
    272272  simulateSinglePhe->SetParameters(lambda,mu0,mu1,sigma0,sigma1);
    273273  simulateSinglePhe->SetParNames("#lambda","#mu_0","#mu_1","#sigma_0","#sigma_1");
    274   simulateSinglePhe->SetNpx(fBPQnbins); 
     274  simulateSinglePhe->SetNpx(fBlindPixelQnbins); 
    275275
    276276  for (Int_t i=0;i<10000; i++)
    277277    {
    278       fHBPQ->Fill(simulateSinglePhe->GetRandom());
     278      fHBlindPixelQ->Fill(simulateSinglePhe->GetRandom());
    279279    }
    280280 
    … …  
    283283
    284284
    285 void MHCalibrationBlindPixel::ChangeFitFunc(BPFitFunc fitfunc, Int_t par)
     285void MHCalibrationBlindPixel::ChangeFitFunc(BlindPixelFitFunc fitfunc, Int_t par)
    286286{
    287287 
    … …  
    303303  // Get the fitting ranges
    304304  //
    305   rmin = (rmin != 0.) ? rmin : fBPQfirst;
    306   rmax = (rmax != 0.) ? rmax : fBPQlast;
     305  rmin = (rmin != 0.) ? rmin : fBlindPixelQfirst;
     306  rmax = (rmax != 0.) ? rmax : fBlindPixelQlast;
    307307
    308308  //
    … …  
    310310  // otherwise the fit goes gaga because of high number of dimensions ...
    311311  //
    312   const Stat_t   entries      = fHBPQ->GetSumOfWeights();
     312  const Stat_t   entries      = fHBlindPixelQ->GetSumOfWeights();
    313313  const Double_t lambda_guess = 0.2;
    314   const Double_t mu_0_guess = fHBPQ->GetBinCenter(fHBPQ->GetMaximumBin());
    315   const Double_t si_0_guess = mu_0_guess/500.;
    316   const Double_t mu_1_guess = mu_0_guess + 2.;
     314  const Double_t mu_0_guess = fHBlindPixelQ->GetBinCenter(fHBlindPixelQ->GetMaximumBin());
     315  const Double_t si_0_guess = mu_0_guess/10.;
     316  const Double_t mu_1_guess = mu_0_guess + 50.;
    317317  const Double_t si_1_guess = si_0_guess + si_0_guess;
    318318
    … …  
    325325  fSinglePheFit->SetParLimits(3,1.0,rmax-rmin);
    326326  fSinglePheFit->SetParLimits(4,1.7,rmax-rmin);
    327   fSinglePheFit->SetParLimits(5,0.,2.*entries);
     327  fSinglePheFit->SetParLimits(5,0.,2.5*entries);
    328328  //
    329329  // Normalize the histogram to facilitate faster fitting of the area
    … …  
    333333  // ROOT gives us another nice example of user-unfriendly behavior:
    334334  // Although the normalization of the function fSinglePhe and the
    335   // Histogram fHBPQ agree (!!), the fit does not normalize correctly INTERNALLY
     335  // Histogram fHBlindPixelQ agree (!!), the fit does not normalize correctly INTERNALLY
    336336  // in the fitting procedure !!!
    337337  //
    … …  
    342342  // So, WE have to adapt to that internal flaw of ROOT:
    343343  //
    344   const Int_t  npx     = fSinglePheFit->GetNpx();
    345   const Int_t  bins    = fHBPQ->GetXaxis()->GetLast()-fHBPQ->GetXaxis()->GetFirst();
    346   //  fHBPQ->Scale(gkSq2Pi*(float)bins/npx/entries);
     344  //  const Int_t  npx     = fSinglePheFit->GetNpx();
     345  //  const Int_t  bins    = fHBlindPixelQ->GetXaxis()->GetLast()-fHBlindPixelQ->GetXaxis()->GetFirst();
     346  //  fHBlindPixelQ->Scale(gkSq2Pi*(float)bins/npx/entries);
    347347
    348348  //
    … …  
    353353  //  fSinglePheFit->SetNpx(fQnbins); 
    354354
    355   fHBPQ->Fit("SinglePheFit",opt);
     355  fHBlindPixelQ->Fit(fSinglePheFit,opt);
    356356
    357357  fLambda = fSinglePheFit->GetParameter(0);
    … …  
    403403  // If you find another solution which WORKS!!, please tell me!!
    404404  //
    405   Int_t nbins = 100;
    406 
    407   *fLog << "New number of bins in HSinQ: " << CutEdges(fHBPQ,nbins) << endl;
    408 
    409   fBPQfirst = fHBPQ->GetBinLowEdge(fHBPQ->GetXaxis()->GetFirst());
    410   fBPQlast  = fHBPQ->GetBinLowEdge(fHBPQ->GetXaxis()->GetLast())+fHBPQ->GetBinWidth(0);
    411   fBPQnbins = nbins;
    412 
    413   *fLog << "New number of bins in HErrQ: " << CutEdges(fHBPErrQ,30) << endl;
    414   fErrBPQfirst = fHBPErrQ->GetBinLowEdge(fHBPErrQ->GetXaxis()->GetFirst());
    415   fErrBPQlast  = fHBPErrQ->GetBinLowEdge(fHBPErrQ->GetXaxis()->GetLast())+fHBPErrQ->GetBinWidth(0);
    416   fErrBPQnbins = nbins;
    417 
    418   CutEdges(fHBPQvsN,0);
     405  Int_t nbins = 50;
     406
     407  *fLog << "New number of bins in HSinQ: " << CutEdges(fHBlindPixelQ,nbins) << endl;
     408
     409  fBlindPixelQfirst = fHBlindPixelQ->GetBinLowEdge(fHBlindPixelQ->GetXaxis()->GetFirst());
     410  fBlindPixelQlast  = fHBlindPixelQ->GetBinLowEdge(fHBlindPixelQ->GetXaxis()->GetLast())+fHBlindPixelQ->GetBinWidth(0);
     411  fBlindPixelQnbins = nbins;
     412
     413  *fLog << "New number of bins in HErrQ: " << CutEdges(fHBlindPixelErrQ,30) << endl;
     414  fErrBlindPixelQfirst = fHBlindPixelErrQ->GetBinLowEdge(fHBlindPixelErrQ->GetXaxis()->GetFirst());
     415  fErrBlindPixelQlast  = fHBlindPixelErrQ->GetBinLowEdge(fHBlindPixelErrQ->GetXaxis()->GetLast())+fHBlindPixelErrQ->GetBinWidth(0);
     416  fErrBlindPixelQnbins = nbins;
     417
     418  CutEdges(fHBlindPixelQvsN,0);
    419419
    420420}
    … …  
    426426    return kFALSE;
    427427
    428   rmin = (rmin != 0.) ? rmin : 0.;
    429   rmax = (rmax != 0.) ? rmax : 16.;
    430 
    431   const Stat_t   entries     = fHBPT->GetEntries();
    432   const Double_t mu_guess    = fHBPT->GetBinCenter(fHBPT->GetMaximumBin());
     428  rmin = (rmin != 0.) ? rmin : 4.;
     429  rmax = (rmax != 0.) ? rmax : 9.;
     430
     431  const Stat_t   entries     = fHBlindPixelT->GetEntries();
     432  const Double_t mu_guess    = fHBlindPixelT->GetBinCenter(fHBlindPixelT->GetMaximumBin());
    433433  const Double_t sigma_guess = (rmax - rmin)/2.;
    434434  const Double_t area_guess  = entries/gkSq2Pi;
    … …  
    441441  fTimeGausFit->SetParLimits(2,0.,rmax-rmin);
    442442
    443   fHBPT->Fit("GausTime",opt);
     443  fHBlindPixelT->Fit(fTimeGausFit,opt);
     444
     445  rmin = fTimeGausFit->GetParameter(1) - 2.*fTimeGausFit->GetParameter(2);
     446  rmax = fTimeGausFit->GetParameter(1) + 2.*fTimeGausFit->GetParameter(2);
     447  fTimeGausFit->SetRange(rmin,rmax); 
     448
     449  fHBlindPixelT->Fit(fTimeGausFit,opt);
     450
    444451
    445452  fMeanT     = fTimeGausFit->GetParameter(2);

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

    r2525 r2599  
    3232private:
    3333
    34   TH1I* fHBPQ;        //-> Histogram with the single Phe spectrum
    35   TH1F* fHBPErrQ;     //-> Variance of summed FADC slices
    36   TH1I* fHBPT;        //-> Variance of summed FADC slices
    37   TH1I* fHBPQvsN;     //-> Summed Charge vs. Event Nr.
     34  TH1I* fHBlindPixelQ;        //-> Histogram with the single Phe spectrum
     35  TH1F* fHBlindPixelErrQ;     //-> Variance of summed FADC slices
     36  TH1I* fHBlindPixelT;        //-> Variance of summed FADC slices
     37  TH1I* fHBlindPixelQvsN;     //-> Summed Charge vs. Event Nr.
    3838 
    3939  TF1 *fSinglePheFit;
    4040  TF1 *fTimeGausFit; 
    4141
    42   Axis_t  fBPQfirst;
    43   Axis_t  fBPQlast;
    44   Int_t   fBPQnbins;
     42  Axis_t  fBlindPixelQfirst;
     43  Axis_t  fBlindPixelQlast;
     44  Int_t   fBlindPixelQnbins;
    4545 
    46   Axis_t fErrBPQfirst;
    47   Axis_t fErrBPQlast;
    48   Int_t  fErrBPQnbins;
     46  Axis_t fErrBlindPixelQfirst;
     47  Axis_t fErrBlindPixelQlast;
     48  Int_t  fErrBlindPixelQnbins;
    4949
    5050  void ResetBin(Int_t i);
    … …  
    5454  Bool_t fFitOK; 
    5555 
    56   BPFitFunc fgSinglePheFitFunc;     // In the beginning,
     56  BlindPixelFitFunc fgSinglePheFitFunc;     // In the beginning,
    5757  Int_t     fgSinglePheFitNPar;     // we want to be flexible using different functions
    5858
    … …  
    8383  ~MHCalibrationBlindPixel();
    8484
    85   Bool_t FillBPQ(Int_t q)         { return fHBPQ->Fill(q) > -1;  } 
    86   Bool_t FillErrBPQ(Float_t errq) { return fHBPErrQ->Fill(errq) > -1; }
    87   Bool_t FillBPT(Int_t t)         { return fHBPT->Fill(t) > -1;  }
    88   Bool_t FillBPQvsN(Stat_t rq, Int_t t) { return fHBPQvsN->Fill(t,rq) > -1;  } 
     85  Bool_t FillBlindPixelQ(Int_t q)         { return fHBlindPixelQ->Fill(q) > -1;  } 
     86  Bool_t FillErrBlindPixelQ(Float_t errq) { return fHBlindPixelErrQ->Fill(errq) > -1; }
     87  Bool_t FillBlindPixelT(Int_t t)         { return fHBlindPixelT->Fill(t) > -1;  }
     88  Bool_t FillBlindPixelQvsN(Stat_t rq, Int_t t) { return fHBlindPixelQvsN->Fill(t,rq) > -1;  } 
    8989 
    9090  const Double_t GetLambda()   const { return fLambda; }
    … …  
    109109  const Double_t GetSigmaTErr()    const { return fSigmaTErr; }
    110110
    111   const TH1F *GetHErrQ() { return fHBPErrQ; }
    112   const TH1F *GetHErrQ() const { return fHBPErrQ; }
     111  const TH1F *GetHErrQ() { return fHBlindPixelErrQ; }
     112  const TH1F *GetHErrQ() const { return fHBlindPixelErrQ; }
    113113 
    114114  Bool_t SimulateSinglePhe(Double_t lambda,
    … …  
    121121  Bool_t FitT(Axis_t rmin=0., Axis_t rmax=0.,Option_t *opt="R0+");
    122122
    123   void ChangeFitFunc(BPFitFunc fitfunc, Int_t par=5);
     123  void ChangeFitFunc(BlindPixelFitFunc fitfunc, Int_t par=5);
    124124
    125125

  • trunk/MagicSoft/Mars/mhist/MHCalibrationConfig.h

    r2533 r2599  
    1212
    1313// Global rejection criteria for the acceptance of a fit: Prob=0.01 == 99% Probability
    14 const Float_t gkProbLimit = 0.01;
     14const Float_t gkProbLimit = 0.001;
    1515
    1616// Starting number of bins for the histo:
    … …  
    3030
    3131// typedef to the fitting functions for the blind pixel
    32 typedef Double_t (*BPFitFunc)(Double_t *, Double_t *);
     32typedef Double_t (*BlindPixelFitFunc)(Double_t *, Double_t *);
    3333
    3434#endif /* MARS_MHCalibrationBlindPixelConfig */

  • trunk/MagicSoft/Mars/mhist/MHCalibrationPixel.cc

    r2581 r2599  
    6060// Default Constructor.
    6161//
    62 MHCalibrationPixel::MHCalibrationPixel(Int_t pix, const char *name, const char *title)
    63       : fPixId(pix),
     62MHCalibrationPixel::MHCalibrationPixel(const char *name, const char *title)
     63      : fPixId(-1),
    6464        fQGausFit(NULL),
    6565        fTGausFit(NULL),
    6666        fFitLegend(NULL),
    6767        fLowerFitRange(0.),
    68         fFitOK(kFALSE)
     68        fFitOK(kFALSE),
     69        fQChisquare(-1.),
     70        fQProb(-1.),
     71        fQNdf(-1),
     72        fTChisquare(-1.),
     73        fTProb(-1.),
     74        fTNdf(-1)
    6975{
    7076
    … …  
    7278    fTitle = title ? title : "Fill the accumulated charges and times of all events and perform fits";
    7379
    74     TString qname  = "HQ";
    7580    TString qtitle = "Distribution of Summed FADC Slices Pixel ";
    76     qname  += pix;
    77     qtitle += pix;
    7881
    7982    // Create a large number of bins, later we will rebin
    … …  
    8285    fQnbins = gkStartPixelBinNr;
    8386
    84     fHQ = new TH1I( qname.Data(),qtitle.Data(),
     87    fHQ = new TH1I("HQ",qtitle.Data(),
    8588                   fQnbins,fQfirst,fQlast);
    8689    fHQ->SetXTitle("Sum FADC Slices");
    … …  
    8891    fHQ->Sumw2();
    8992
    90     TString tname  = "HT";
     93    fHQ->SetDirectory(NULL);
     94
    9195    TString ttitle = "Distribution of Mean Arrival Times Pixel ";
    92     tname  += pix;
    93     ttitle += pix;
    9496
    9597    Axis_t tfirst = -0.5;
    … …  
    9799    Int_t nbins   = 16;
    98100
    99     fHT = new TH1I(tname.Data(),ttitle.Data(),
     101    fHT = new TH1I("HT",ttitle.Data(),
    100102                  nbins,tfirst,tlast);
    101103    fHT->SetXTitle("Mean Arrival Times [FADC slice nr]");
    … …  
    103105    fHT->Sumw2();
    104106
    105     TString qvsnname  = "HQvsN";
     107    fHT->SetDirectory(NULL);
     108
    106109    TString qvsntitle = "Sum of Charges vs. Event Number Pixel ";
    107     qvsnname  += pix;
    108     qvsntitle += pix;
    109110
    110111    // We define a reasonable number and later enlarge it if necessary
    … …  
    113114    Axis_t nlast  = (Axis_t)nbins - 0.5;
    114115
    115     fHQvsN = new TH1I(qvsnname.Data(),qvsntitle.Data(),
     116    fHQvsN = new TH1I("HQvsN",qvsntitle.Data(),
    116117                     nbins,nfirst,nlast);
    117118    fHQvsN->SetXTitle("Event Nr.");
    118119    fHQvsN->SetYTitle("Sum of FADC slices");
    119120
    120     fQChisquare = -1.;
    121     fQProb      = -1.;
    122     fQNdf       = -1;
    123    
    124     fTChisquare = -1.;
    125     fTProb      = -1.;
    126     fTNdf       = -1;
     121    fHQvsN->SetDirectory(NULL);
    127122
    128123}
    … …  
    144139}
    145140
     141
     142void MHCalibrationPixel::ChangeHistId(Int_t id)
     143{
     144
     145  fPixId = id;
     146 
     147  TString nameQ = TString(fHQ->GetName());
     148  nameQ += id;
     149  fHQ->SetName(nameQ.Data());
     150
     151  TString nameT = TString(fHT->GetName());
     152  nameT += id;
     153  fHT->SetName(nameT.Data());
     154
     155  TString nameQvsN  = TString(fHQvsN->GetName());
     156  nameQvsN += id;
     157  fHQvsN->SetName(nameQvsN.Data());
     158}
     159
     160
    146161void MHCalibrationPixel::Reset()
    147162{
    … …  
    206221
    207222  char line8[32];
    208   sprintf(line8,"Probability: %4.2f ",fQProb);
     223  sprintf(line8,"Probability: %4.3f ",fQProb);
    209224  fFitLegend->AddText(line8);
    210225
    … …  
    266281    fHT->DrawCopy(opt);
    267282
    268     if (fHT->GetFunction("GausTime"))
     283    if (fTGausFit)
    269284      {
    270         TF1 *tfit = fHT->GetFunction("GausTime");
    271         if (tfit->GetProb() < 0.01)
    272           tfit->SetLineColor(kRed);
     285        if (fTChisquare > 1.)
     286          fTGausFit->SetLineColor(kRed);
    273287        else
    274           tfit->SetLineColor(kGreen);
    275 
    276         tfit->DrawCopy("same");
     288          fTGausFit->SetLineColor(kGreen);
     289
     290        fTGausFit->DrawCopy("same");
    277291        c->Modified();
    278292        c->Update();
    … …  
    290304Bool_t MHCalibrationPixel::FitT(Axis_t rmin, Axis_t rmax, Option_t *option)
    291305{
    292  
     306
    293307  if (fTGausFit)
    294308    return kFALSE;
    295309
    296   rmin = (rmin != 0.) ? rmin : -0.5;
    297   rmax = (rmax != 0.) ? rmax : 15.5;
     310  rmin = (rmin != 0.) ? rmin : 4.;
     311  rmax = (rmax != 0.) ? rmax : 9.;
    298312
    299313  const Stat_t   entries     = fHT->GetEntries();
    … …  
    302316  const Double_t area_guess  = entries/gkSq2Pi;
    303317
    304   fTGausFit = new TF1("GausTime","gaus",rmin,rmax); 
     318  TString name = TString("GausTime");
     319  name += fPixId;
     320  fTGausFit = new TF1(name.Data(),"gaus",rmin,rmax); 
    305321
    306322  if (!fTGausFit)
    … …  
    314330  fTGausFit->SetParLimits(0,0.,entries);
    315331  fTGausFit->SetParLimits(1,rmin,rmax);
    316   fTGausFit->SetParLimits(2,0.,rmax-rmin);
    317 
    318   fHT->Fit("GausTime",option);
     332  fTGausFit->SetParLimits(2,0.,(rmax-rmin)/2.);
     333  fTGausFit->SetRange(rmin,rmax);
     334
     335  fHT->Fit(fTGausFit,option);
     336
     337  rmin = fTGausFit->GetParameter(1) - 3.*fTGausFit->GetParameter(2);
     338  rmax = fTGausFit->GetParameter(1) + 3.*fTGausFit->GetParameter(2);
     339  fTGausFit->SetRange(rmin,rmax); 
     340
     341  fHT->Fit(fTGausFit,option);
    319342
    320343  fTChisquare = fTGausFit->GetChisquare();
    … …  
    324347  fTSigma     = fTGausFit->GetParameter(2);
    325348
    326   if (fTProb < gkProbLimit)
     349  if (fTChisquare > 1.)
    327350    {
    328351      *fLog << warn << "Fit of the Arrival times failed ! " << endl;
    … …  
    344367  //
    345368  Axis_t rmin = (fLowerFitRange != 0.) ? fLowerFitRange : fQfirst;
     369  Axis_t rmax = 0.;
    346370
    347371  //
    … …  
    352376  const Double_t ar_guess = entries/gkSq2Pi;
    353377  const Double_t mu_guess = fHQ->GetBinCenter(fHQ->GetMaximumBin());
    354   const Double_t si_guess = mu_guess/500.;
    355 
    356   fQGausFit = new TF1("QGausFit","gaus",rmin,fQlast);
     378  const Double_t si_guess = mu_guess/50.;
     379
     380  TString name = TString("QGausFit");
     381  name += fPixId;
     382
     383  fQGausFit = new TF1(name.Data(),"gaus",rmin,fQlast);
    357384
    358385  if (!fQGausFit)
    … …  
    367394  fQGausFit->SetParLimits(1,rmin,fQlast);
    368395  fQGausFit->SetParLimits(2,0.,fQlast-rmin);
    369 
    370   fHQ->Fit("QGausFit",option);
     396  fQGausFit->SetRange(rmin,fQlast);
     397  fQGausFit->Update();
     398
     399  fHQ->Fit(fQGausFit,option);
     400 
     401  rmin = fQGausFit->GetParameter(1) - 3.*fQGausFit->GetParameter(2);
     402  rmax = fQGausFit->GetParameter(1) + 3.*fQGausFit->GetParameter(2);
     403
     404  fQGausFit->SetRange(rmin,rmax); 
     405  fHQ->Fit(fQGausFit,option);
     406
     407  rmin = fQGausFit->GetParameter(1) - 3.5*fQGausFit->GetParameter(2);
     408  rmax = fQGausFit->GetParameter(1) + 3.5*fQGausFit->GetParameter(2);
     409
     410  fQGausFit->SetRange(rmin,rmax); 
     411  fHQ->Fit(fQGausFit,option);
    371412
    372413  fQChisquare = fQGausFit->GetChisquare();
    … …  
    399440{
    400441
    401   //
    402   // The number 100 is necessary because it is the internal binning
    403   // of ROOT functions. A call to SetNpx() does NOT help
    404   // If you find another solution which WORKS!!, please tell me!!
    405   //
    406   Int_t nbins = 100;
     442  Int_t nbins = 50;
    407443
    408444  CutEdges(fHQ,nbins);

  • trunk/MagicSoft/Mars/mhist/MHCalibrationPixel.h

    r2581 r2599  
    6868public:
    6969
    70   MHCalibrationPixel(Int_t pix=-1, const char *name=NULL, const char *title=NULL);
     70  MHCalibrationPixel(const char *name=NULL, const char *title=NULL);
    7171  ~MHCalibrationPixel();
     72
     73  void ChangeHistId(Int_t i);
    7274 
    7375  Bool_t SetupFill(const MParList *pList);
    … …  
    105107  const TH1I *GetHQvsN() const { return fHQvsN; }
    106108 
    107   Bool_t FitQ(Option_t *option="RQ0+"); 
    108   Bool_t FitT(Axis_t rmin=0, Axis_t rmax=0, Option_t *option="RQ0+");   
     109  Bool_t FitQ(Option_t *option="RQ0"); 
     110  Bool_t FitT(Axis_t rmin=0, Axis_t rmax=0, Option_t *option="RQ0");   
    109111
    110112  virtual void Draw(Option_t *option="");
Note: See TracChangeset for help on using the changeset viewer.