Ignore:
Timestamp:
04/15/04 11:17:54 (21 years ago)
Author:
gaug
Message:
*** empty log message ***
File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r3740 r3742  
    2222!
    2323\* ======================================================================== */
     24//const TString inpath = "/home/rootdata/Calib/2004_04_08/";
     25const TString inpath = "/home/rootdata/BlindPixel/";
     26const Int_t pedrun = 20491;
     27const Int_t calrun = 20494;
    2428
    25 void calibration()
     29void calibration(const Int_t prun=pedrun, const Int_t crun=calrun)
    2630{
    2731 
    28   const MCalibrationCam::PulserColor_t color = MCalibrationCam::kGREEN;
    29 
    30   const TString inpath = "/home/rootdata/Calib/2004_04_08/";
    31 
    3232  MRunIter pruns;
    3333  MRunIter cruns;
    3434
    35   pruns.AddRun(22254,inpath);
    36   cruns.AddRun(22253,inpath);
     35  pruns.AddRun(prun,inpath);
     36  cruns.AddRun(crun,inpath);
    3737
    38   gStyle->SetOptStat(1111);
    39   gStyle->SetOptFit();
    40  
     38  MCalibrationCam::PulserColor_t color;
     39
     40  if (crun < 20000)
     41    color = MCalibrationCam::kCT1;
     42  else
     43    color = FindColor((MDirIter*)&cruns);
     44
     45  if (color == MCalibrationCam::kNONE)
     46    return;
     47
    4148  MCalibrationQECam qecam;
    4249  MBadPixelsCam     badcam;
    … …  
    98105  if (!calloop.Process(pedloop.GetPedestalCam()))
    99106    return;
     107 
     108}
    100109
     110MCalibrationCam::PulserColor_t FindColor(MDirIter* run)
     111{
     112 
     113  MCalibrationCam::PulserColor_t col = MCalibrationCam::kNONE;
    101114
     115  TString filenames;
    102116
    103 #if 0
    104   //
    105   // The longer version:
    106   //
    107  
    108   //
    109   // Create a empty Parameter List and an empty Task List
    110   //
    111   MParList  plist;
    112   MTaskList tlist;
    113   plist.AddToList(&tlist);
    114   plist.AddToList(&pedloop.GetPedestalCam());
    115   plist.AddToList(&pedloop.GetBadPixels());   
     117  while (!(filenames=run->Next()).IsNull())
     118    {
    116119
    117   gLog << endl;;
    118   gLog << "Calculate MCalibrationCam from Runs " << cruns.GetRunsAsString() << endl;
    119   gLog << endl;
    120  
    121   MReadMarsFile read("Events");
    122   read.DisableAutoScheme();
    123   static_cast<MRead&>(read).AddFiles(cruns);
    124  
    125   MGeomCamMagic              geomcam;
    126   MExtractedSignalCam        sigcam;
    127   MArrivalTimeCam            timecam;
    128   MCalibrationRelTimeCam     relcam;
    129   MCalibrationQECam          qecam;
    130   MCalibrationChargeCam      calcam;
    131   MCalibrationChargePINDiode pindiode;   
    132   MCalibrationChargeBlindPix blindpix;   
    133  
    134   MHCalibrationRelTimeCam     histtime;
    135   MHCalibrationChargeCam      histcharge;
    136   MHCalibrationChargePINDiode histpin;
    137   MHCalibrationChargeBlindPix histblind;
    138   histcharge.SetPulserFrequency(500);
    139   histblind.SetSinglePheCut(600);
    140   //
    141   // Get the previously created MPedestalCam into the new Parameter List
    142   //
    143   plist.AddToList(&geomcam);
    144   plist.AddToList(&sigcam);
    145   plist.AddToList(&timecam);
    146   plist.AddToList(&relcam);
    147   plist.AddToList(&qecam);
    148   plist.AddToList(&calcam);
    149   plist.AddToList(&histtime);
    150   plist.AddToList(&histcharge);
    151   //    plist.AddToList(&histpin);
    152   plist.AddToList(&histblind);
    153  
    154   //
    155   // We saw that the signal jumps between slices,
    156   // thus take the sliding window
    157   //           
    158   MExtractSignal2        sigcalc2;
    159   MExtractPINDiode       pincalc;
    160   MExtractBlindPixel     blindcalc;
    161   sigcalc2.SetRange(2,15,6,5,14,6);
    162   blindcalc.SetRange(12,17);
     120      filenames.ToLower();
    163121
    164   MArrivalTimeCalc2      timecalc;
    165   MCalibrationChargeCalc calcalc;
    166   calcalc.SetPulserColor(color);
    167   MGeomApply             geomapl;
    168  
    169   MFillH filltime( "MHCalibrationRelTimeCam"    , "MArrivalTimeCam");
    170   //   MFillH fillpin  ("MHCalibrationChargePINDiode", "MExtractedSignalPINDiode");
    171   MFillH fillblind("MHCalibrationChargeBlindPix", "MExtractedSignalBlindPixel");
    172   MFillH fillcam  ("MHCalibrationChargeCam"     , "MExtractedSignalCam");
    173  
    174   //
    175   // Skip the HiGain vs. LoGain calibration
    176   //
    177   calcalc.SkipHiLoGainCalibration();
    178  
    179   //
    180   // Apply a filter against cosmics
    181   // (was directly in MCalibrationCalc in earlier versions)
    182   //
    183   MFCosmics            cosmics;
    184   MContinue            cont(&cosmics);
    185  
    186   tlist.AddToList(&read);
    187   tlist.AddToList(&geomapl);
    188   tlist.AddToList(&sigcalc2);
    189   tlist.AddToList(&blindcalc);
    190   //    tlist.AddToList(&pincalc);
    191   //
    192   // In case, you want to skip the cosmics rejection,
    193   // uncomment the next line
    194   //
    195   tlist.AddToList(&cont);
    196   //
    197   // In case, you want to skip the somewhat lengthy calculation
    198   // of the arrival times using a spline, uncomment the next two lines
    199   //
    200   tlist.AddToList(&timecalc);
    201   tlist.AddToList(&filltime);
    202   //    tlist.AddToList(&fillpin);
    203   tlist.AddToList(&fillblind);
    204   tlist.AddToList(&fillcam);
    205   //
    206   tlist.AddToList(&calcalc);
    207   //
    208   // Create and setup the eventloop
    209   //
    210   MEvtLoop evtloop;
    211   evtloop.SetParList(&plist);
    212   evtloop.SetDisplay(display);
    213  
    214   //
    215   // Execute second analysis
    216   //
    217   if (!evtloop.Eventloop())
    218     return;
    219  
    220   tlist.PrintStatistics();
    221  
    222   MBadPixelsCam *badpixels = (MBadPixelsCam*)plist->FindObject("MBadPixelsCam");
     122      if (filenames.Contains("green"))
     123        if (col == MCalibrationCam::kNONE)
     124          {
     125            cout << "Found colour: Green  in " << filenames << endl;
     126            col = MCalibrationCam::kGREEN;
     127          }
     128        else if (col != MCalibrationCam::kGREEN)
     129          {
     130            cout << "Different colour found in " << filenames << "... abort" << endl;
     131            return MCalibrationCam::kNONE;
     132          }
    223133
    224   //
    225   // print the most important results of all pixels to a file
    226   //
    227   /*
    228   MLog gauglog;
    229   gauglog.SetOutputFile(Form("%s%s",calcam.GetName(),".txt"),1);
    230   calcam.SetLogStream(&gauglog);
    231   badpixels->Print();
    232   calcam.SetLogStream(&gLog);
    233   */
    234   //
    235   // just one example how to get the plots of individual pixels
    236   //
    237   //  histblind.DrawClone("all");
    238   //  histcharge[400].DrawClone("all");
    239   //  histcharge(5).DrawClone("all");
    240   //  histtime[5].DrawClone("fourierevents");
    241   for (Int_t aidx=0;aidx<2;aidx++)
    242     {
    243       histcharge.GetAverageHiGainArea(aidx).DrawClone("all");
    244       histcharge.GetAverageLoGainArea(aidx).DrawClone("all");
    245     }
    246  
    247   for (Int_t sector=1;sector<7;sector++)
    248     {
    249       histcharge.GetAverageHiGainSector(sector).DrawClone("all");
    250       histcharge.GetAverageLoGainSector(sector).DrawClone("all");
     134      if (filenames.Contains("blue"))
     135        if (col == MCalibrationCam::kNONE)
     136          {
     137            cout << "Found colour: Blue  in " << filenames << endl;
     138            col = MCalibrationCam::kBLUE;
     139          }
     140        else if (col != MCalibrationCam::kBLUE)
     141          {
     142            cout << "Different colour found in " << filenames << "... abort" << endl;
     143            return MCalibrationCam::kNONE;
     144          }
     145
     146      if (filenames.Contains("uv"))
     147        if (col == MCalibrationCam::kNONE)
     148          {
     149            cout << "Found colour: Uv  in " << filenames << endl;
     150            col = MCalibrationCam::kUV;
     151          }
     152        else if (col != MCalibrationCam::kUV)
     153          {
     154            cout << "Different colour found in " << filenames << "... abort" << endl;
     155            return MCalibrationCam::kNONE;
     156          }
     157
     158      if (filenames.Contains("ct1"))
     159        if (col == MCalibrationCam::kNONE)
     160          {
     161            cout << "Found colour: Ct1  in " << filenames << endl;
     162            col = MCalibrationCam::kCT1;
     163          }
     164        else if (col != MCalibrationCam::kCT1)
     165          {
     166            cout << "Different colour found in " << filenames << "... abort" << endl;
     167            return MCalibrationCam::kNONE;
     168          }
     169     
    251170    }
    252171 
    253172
    254   // Create histograms to display
    255   MHCamera disp1  (geomcam, "Cal;Charge",               "Fitted Mean Charges");
    256   MHCamera disp2  (geomcam, "Cal;SigmaCharge",          "Sigma of Fitted Charges");
    257   MHCamera disp3  (geomcam, "Cal;FitProb",              "Probability of Fit");
    258   MHCamera disp4  (geomcam, "Cal;RSigma",               "Reduced Sigmas");
    259   MHCamera disp5  (geomcam, "Cal;RSigma/Charge",        "Reduced Sigma per Charge");
    260   MHCamera disp6  (geomcam, "Cal;FFactorPhe",           "Nr. of Photo-electrons (F-Factor Method)");
    261   MHCamera disp7  (geomcam, "Cal;FFactorConv",          "Conversion Factor to photons (F-Factor Method)");
    262   MHCamera disp8  (geomcam, "Cal;FFactorFFactor",       "Total F-Factor (F-Factor Method)");
    263   MHCamera disp9  (geomcam, "Cal;CascadesQEFFactor",    "Av. Quantum Efficiency (F-Factor Method)");
    264   MHCamera disp10 (geomcam, "Cal;QEFFactor",            "Measured QE (F-Factor Method)");
    265   MHCamera disp11 (geomcam, "Cal;PINDiodeConv",         "Conversion Factor tp photons (PIN Diode Method)");
    266   MHCamera disp12 (geomcam, "Cal;PINDiodeFFactor",      "Total F-Factor (PIN Diode Method)");
    267   MHCamera disp13 (geomcam, "Cal;Excluded",             "Pixels previously excluded");
    268   MHCamera disp14 (geomcam, "Cal;Unsuited",             "Unsuited Pixels  ");
    269   MHCamera disp15 (geomcam, "Cal;Unreliable",           "Unreliable Pixels");
    270   MHCamera disp16 (geomcam, "Cal;HiGainOscillating",    "Oscillating Pixels High Gain");
    271   MHCamera disp17 (geomcam, "Cal;LoGainOscillating",    "Oscillating Pixels Low Gain");
    272   MHCamera disp18 (geomcam, "Cal;HiGainPickup",         "Number Pickup events Hi Gain");
    273   MHCamera disp19 (geomcam, "Cal;LoGainPickup",         "Number Pickup events Lo Gain");
    274   MHCamera disp20 (geomcam, "Cal;Saturation",           "Pixels with saturated Hi Gain");
    275   MHCamera disp21 (geomcam, "Cal;FFactorValid",         "Pixels with valid F-Factor calibration");
    276   MHCamera disp22 (geomcam, "Cal;BlindPixelValid",      "Pixels with valid BlindPixel calibration");
    277   MHCamera disp23 (geomcam, "Cal;PINdiodeFFactorValid", "Pixels with valid PINDiode calibration");
     173     
     174  if (col == MCalibrationCam::kNONE)
     175    cout <<  "No colour found in filenames of runs: " << ((MRunIter*)run)->GetRunsAsString()
     176         << "... abort" << endl;
    278177 
    279   MHCamera disp24 (geomcam, "Cal;Ped",         "Pedestals");
    280   MHCamera disp25 (geomcam, "Cal;PedRms",      "Pedestal RMS");
    281  
    282   MHCamera disp26 (geomcam, "time;Time",        "Rel. Arrival Times");
    283   MHCamera disp27 (geomcam, "time;SigmaTime",   "Sigma of Rel. Arrival Times");
    284   MHCamera disp28 (geomcam, "time;TimeProb",    "Probability of Time Fit");
    285   MHCamera disp29 (geomcam, "time;NotFitValid", "Pixels with not valid fit results");
    286   MHCamera disp30 (geomcam, "time;Oscillating", "Oscillating Pixels");
    287  
    288   MHCamera disp31 (geomcam, "Cal;AbsTimeMean",  "Abs. Arrival Times");
    289   MHCamera disp32 (geomcam, "Cal;AbsTimeRms",   "RMS of Arrival Times");
    290  
    291   // Fitted charge means and sigmas
    292   disp1.SetCamContent(calcam,  0);
    293   disp1.SetCamError(  calcam,  1);
    294   disp2.SetCamContent(calcam,  2);
    295   disp2.SetCamError(  calcam,  3);
    296  
    297   // Fit probabilities
    298   disp3.SetCamContent(calcam,  4);
    299  
    300   // Reduced Sigmas and reduced sigmas per charge
    301   disp4.SetCamContent(calcam,  5);
    302   disp4.SetCamError(  calcam,  6);
    303   disp5.SetCamContent(calcam,  7);
    304   disp5.SetCamError(  calcam,  8);
    305  
    306   // F-Factor Method
    307   disp6.SetCamContent(calcam,  9);
    308   disp6.SetCamError(  calcam, 10);
    309   disp7.SetCamContent(calcam, 11);
    310   disp7.SetCamError(  calcam, 12);
    311   disp8.SetCamContent(calcam, 13);
    312   disp8.SetCamError(  calcam, 14);
    313  
    314   // Quantum efficiency
    315   disp9.SetCamContent( qecam, 0 );
    316   disp9.SetCamError(   qecam, 1 );
    317   disp10.SetCamContent(qecam, 8);
    318   disp10.SetCamError(  qecam, 9);
    319  
    320   // PIN Diode Method
    321   disp11.SetCamContent(calcam,21);
    322   disp11.SetCamError(  calcam,22);
    323   disp12.SetCamContent(calcam,23);
    324   disp12.SetCamError(  calcam,24);
    325  
    326   // Pixels with defects
    327   disp13.SetCamContent(calcam,26);
    328   disp14.SetCamContent(*badpixels,1);
    329   disp15.SetCamContent(*badpixels,3);
    330   disp16.SetCamContent(*badpixels,10);
    331   disp17.SetCamContent(*badpixels,11);
    332   disp18.SetCamContent(calcam,27);
    333   disp19.SetCamContent(calcam,28);
    334  
    335   // Lo Gain calibration
    336   disp20.SetCamContent(calcam,29);
    337  
    338   // Valid flags
    339   disp21.SetCamContent(calcam,15);
    340   disp22.SetCamContent(calcam,20);
    341   disp23.SetCamContent(calcam,25);
    342  
    343   // Pedestals
    344   disp24.SetCamContent(calcam,30);
    345   disp24.SetCamError(  calcam,31);
    346   disp25.SetCamContent(calcam,32);
    347   disp25.SetCamError(  calcam,33);
    348  
    349   // Relative Times
    350   disp26.SetCamContent(histtime,0);
    351   disp26.SetCamError(  histtime,1);
    352   disp27.SetCamContent(histtime,2);
    353   disp27.SetCamError(  histtime,3);
    354   disp28.SetCamContent(histtime,4);
    355   disp29.SetCamContent(histtime,5);
    356   disp30.SetCamContent(histtime,6);
    357  
    358   // Absolute Times
    359   disp31.SetCamContent(calcam,34);
    360   disp31.SetCamError(  calcam,35);
    361   disp32.SetCamContent(calcam,35);
    362  
    363   disp1.SetYTitle("Mean Charge [FADC Counts]");
    364   disp2.SetYTitle("\\sigma_{Charge} [FADC Counts]");
    365   disp3.SetYTitle("P_{Sum} [1]");
    366  
    367   disp4.SetYTitle("\\sqrt{\\sigma^{2}_{Charge} - RMS^{2}_{Ped}} [FADC Counts]");
    368   disp5.SetYTitle("Reduced Sigma / Mean Charge [1]");
    369  
    370   disp6.SetYTitle("Nr. Photo-electrons [1]");
    371   disp7.SetYTitle("Conversion Factor [Ph/FADC Count]");
    372   disp8.SetYTitle("\\sqrt{N_{Ph}}*\\sigma_{Charge}/\\mu_{Charge} [1] ");
    373  
    374   disp9.SetYTitle("Average QE Cascades [1]");
    375   disp10.SetYTitle("Measured QE (F-Factor Method) [1]");
    376  
    377   disp11.SetYTitle("Conversion Factor [Phot/FADC Count]");
    378   disp12.SetYTitle("\\sqrt{N_{Ph}}*\\sigma_{Charge}/\\mu_{Charge} [1]");
    379  
    380   disp13.SetYTitle("[1]");
    381   disp14.SetYTitle("[1]");
    382   disp15.SetYTitle("[1]");
    383   disp16.SetYTitle("[1]");
    384   disp17.SetYTitle("[1]");
    385   disp18.SetYTitle("[1]");
    386   disp19.SetYTitle("[1]");
    387   disp20.SetYTitle("[1]");
    388   disp21.SetYTitle("[1]");
    389   disp22.SetYTitle("[1]");
    390   disp23.SetYTitle("[1]");
    391  
    392   disp24.SetYTitle("Ped [FADC Counts ]");
    393   disp25.SetYTitle("RMS_{Ped} [FADC Counts ]");
    394  
    395   disp26.SetYTitle("Time Offset [ns]");
    396   disp27.SetYTitle("Timing resolution [ns]");
    397   disp28.SetYTitle("P_{Time} [1]");
    398  
    399   disp29.SetYTitle("[1]");
    400   disp30.SetYTitle("[1]");
    401  
    402   disp31.SetYTitle("Mean Abs. Time [FADC slice]");
    403   disp32.SetYTitle("RMS Abs. Time [FADC slices]");
    404  
    405   gStyle->SetOptStat(1111);
    406   gStyle->SetOptFit();
    407 
    408   // Charges
    409   TCanvas &c1 = display->AddTab("Fit.Charge");
    410   c1.Divide(2, 4);
    411  
    412   CamDraw(c1, disp1,calcam,1, 2 , 2);
    413   CamDraw(c1, disp2,calcam,2, 2 , 2);
    414  
    415   // Fit Probability
    416   TCanvas &c2 = display->AddTab("Fit.Prob");
    417   c2.Divide(1,4);
    418  
    419   CamDraw(c2, disp3,calcam,1,1,4);
    420  
    421   // Reduced Sigmas
    422   TCanvas &c3 = display->AddTab("Red.Sigma");
    423   c3.Divide(2,4);
    424  
    425   CamDraw(c3, disp4,calcam,1, 2 , 2);
    426   CamDraw(c3, disp5,calcam,2, 2 , 2);
    427 
    428  
    429   // F-Factor Method
    430   TCanvas &c4 = display->AddTab("F-Factor");
    431   c4.Divide(3,4);
    432  
    433   CamDraw(c4, disp6,calcam,1, 3 , 2);
    434   CamDraw(c4, disp7,calcam,2, 3 , 2);
    435   CamDraw(c4, disp8,calcam,3, 3 , 2);
    436  
    437 
    438   // Quantum Efficiencies
    439   TCanvas &c5 = display->AddTab("QE");
    440   c5.Divide(2, 4);
    441  
    442   CamDraw(c5, disp9 ,calcam,1,2, 2);
    443   CamDraw(c5, disp10,calcam,2,2, 2);
    444  
    445   // PIN Diode Method
    446   TCanvas &c6 = display->AddTab("PINDiode");
    447   c6.Divide(2,4);
    448  
    449   CamDraw(c6, disp11,calcam,1,2, 2);
    450   CamDraw(c6, disp12,calcam,2,2, 2);
    451  
    452   // Defects
    453   TCanvas &c7 = display->AddTab("Defects");
    454   c7.Divide(4,2);
    455  
    456   CamDraw(c7, disp13,calcam,1,4, 0);
    457   CamDraw(c7, disp14,calcam,2,4, 0);
    458   CamDraw(c7, disp18,calcam,3,4, 0);
    459   CamDraw(c7, disp19,calcam,4,4, 0);
    460  
    461   // BadCam
    462   TCanvas &c8 = display->AddTab("Defects");
    463   c8.Divide(3,2);
    464  
    465   CamDraw(c8, disp15,badcam,1,3, 0);
    466   CamDraw(c8, disp16,badcam,2,3, 0);
    467   CamDraw(c8, disp17,badcam,3,3, 0);
    468  
    469   // Valid flags
    470   TCanvas &c9 = display->AddTab("Validity");
    471   c9.Divide(4,2);
    472  
    473   CamDraw(c9, disp20,calcam,1,4,0);
    474   CamDraw(c9, disp21,calcam,2,4,0);
    475   CamDraw(c9, disp22,calcam,3,4,0);
    476   CamDraw(c9, disp23,calcam,4,4,0);
    477  
    478   // Pedestals
    479   TCanvas &c10 = display->AddTab("Pedestals");
    480   c10.Divide(2,4);
    481  
    482   CamDraw(c10,disp24,calcam,1,2,1);
    483   CamDraw(c10,disp25,calcam,2,2,2);
    484  
    485   // Rel. Times
    486   TCanvas &c11 = display->AddTab("Fitted Rel. Times");
    487   c11.Divide(3,4);
    488  
    489   CamDraw(c11,disp26,calcam,1,3,2);
    490   CamDraw(c11,disp27,calcam,2,3,2);
    491   CamDraw(c11,disp28,calcam,3,3,4);
    492  
    493   // Time Defects
    494   TCanvas &c12 = display->AddTab("Time Def.");
    495   c12.Divide(2,2);
    496  
    497   CamDraw(c12, disp29,calcam,1,2, 0);
    498   CamDraw(c12, disp30,calcam,2,2, 0);
    499  
    500   // Abs. Times
    501   TCanvas &c13 = display->AddTab("Abs. Times");
    502   c13.Divide(2,4);
    503  
    504   CamDraw(c13,disp31,calcam,1,2,2);
    505   CamDraw(c13,disp32,calcam,2,2,2);
    506 #endif
    507 
     178  return col;     
    508179}
    509 
    510 
    511 void CamDraw(TCanvas &c, MHCamera &cam, TObject &evt, Int_t i, Int_t j, Int_t fit)
    512 {
    513  
    514   TArrayI s0(6);
    515   s0[0] = 1;
    516   s0[1] = 2;
    517   s0[2] = 3;
    518   s0[3] = 4;
    519   s0[4] = 5;
    520   s0[5] = 6;
    521 
    522   TArrayI s1(3);
    523   s1[0] = 6;
    524   s1[1] = 1;
    525   s1[2] = 2;
    526  
    527   TArrayI s2(3);
    528   s2[0] = 3;
    529   s2[1] = 4;
    530   s2[2] = 5;
    531  
    532   TArrayI inner(1);
    533   inner[0] = 0;
    534  
    535   TArrayI outer(1);
    536   outer[0] = 1;
    537 
    538   c.cd(i);
    539   gPad->SetBorderMode(0);
    540   gPad->SetTicks();
    541   cam.GetXaxis()->SetLabelOffset(0.005);
    542   cam.GetXaxis()->SetLabelSize(0.06); 
    543   cam.GetYaxis()->SetLabelOffset(0.005);
    544   cam.GetYaxis()->SetLabelSize(0.06); 
    545   cam.GetXaxis()->SetTitleOffset(0.85);
    546   cam.GetXaxis()->SetTitleSize(0.06); 
    547   cam.GetYaxis()->SetTitleOffset(0.7);
    548   cam.GetYaxis()->SetTitleSize(0.06); 
    549   MHCamera *obj1 = (MHCamera*)cam.DrawCopy("hist");
    550   obj1->SetDirectory(NULL);
    551 
    552 
    553   c.cd(i+j);
    554   //  obj1->AddNotify(&evt);
    555   obj1->SetPrettyPalette();
    556   obj1->Draw();
    557 
    558   if (fit != 0)
    559     {
    560       c.cd(i+2*j);
    561       gPad->SetBorderMode(0);
    562       gPad->SetTicks();
    563       TProfile *obj2 = obj1->RadialProfile(Form("%s%s",obj1->GetName(),"_rad"));
    564       obj2->SetDirectory(NULL);
    565       obj2->GetXaxis()->SetLabelOffset(0.005);
    566       obj2->GetXaxis()->SetLabelSize(0.06); 
    567       obj2->GetYaxis()->SetLabelOffset(0.005);
    568       obj2->GetYaxis()->SetLabelSize(0.06); 
    569       obj2->GetXaxis()->SetTitleOffset(0.85);
    570       obj2->GetXaxis()->SetTitleSize(0.06); 
    571       obj2->GetYaxis()->SetTitleOffset(0.7);
    572       obj2->GetYaxis()->SetTitleSize(0.06); 
    573       obj2->Draw();
    574       obj2->SetBit(kCanDelete);
    575      
    576       TProfile *hprof[2];
    577       hprof[0] = obj1->RadialProfileS(s0, inner,Form("%s%s",obj1->GetName(), "Inner"));
    578       hprof[1] = obj1->RadialProfileS(s0, outer,Form("%s%s",obj1->GetName(), "Outer"));
    579 
    580      
    581       for (Int_t k=0; k<2; k++)     
    582         {
    583           Double_t min = cam.GetGeomCam().GetMinRadius(k);
    584           Double_t max = cam.GetGeomCam().GetMaxRadius(k);
    585 
    586           hprof[k]->SetLineColor(kRed+k);
    587           hprof[k]->SetDirectory(0);
    588           hprof[k]->SetBit(kCanDelete);
    589           hprof[k]->Draw("same");
    590           hprof[k]->Fit("pol1","Q","",min,max);
    591           hprof[k]->GetFunction("pol1")->SetLineColor(kRed+k);
    592           hprof[k]->GetFunction("pol1")->SetLineWidth(1);
    593         }
    594      
    595       gPad->Modified();
    596       gPad->Update();
    597 
    598       c.cd(i+3*j);
    599       gPad->SetBorderMode(0);
    600       gPad->SetTicks();
    601       TH1D *obj3 = (TH1D*)obj1->Projection(Form("%s%s",obj1->GetName(),"_py"));
    602       obj3->SetDirectory(NULL);
    603 //      obj3->Sumw2();
    604       obj3->GetXaxis()->SetLabelOffset(0.005);
    605       obj3->GetXaxis()->SetLabelSize(0.06); 
    606       obj3->GetYaxis()->SetLabelOffset(0.005);
    607       obj3->GetYaxis()->SetLabelSize(0.06); 
    608       obj3->GetXaxis()->SetTitleOffset(0.85);
    609       obj3->GetXaxis()->SetTitleSize(0.06); 
    610       obj3->GetYaxis()->SetTitleOffset(0.7);
    611       obj3->GetYaxis()->SetTitleSize(0.06); 
    612       obj3->Draw();
    613       obj3->SetBit(kCanDelete);
    614 
    615       gPad->Modified();
    616       gPad->Update();
    617 
    618       const Double_t min   = obj3->GetBinCenter(obj3->GetXaxis()->GetFirst());
    619       const Double_t max   = obj3->GetBinCenter(obj3->GetXaxis()->GetLast());
    620       const Double_t integ = obj3->Integral("width")/2.5066283;
    621       const Double_t mean  = obj3->GetMean();
    622       const Double_t rms   = obj3->GetRMS();
    623       const Double_t width = max-min;
    624 
    625       if (rms == 0. || width == 0. )
    626         return;
    627      
    628       switch (fit)
    629         {
    630         case 1:
    631           TF1 *sgaus = new TF1("sgaus","gaus(0)",min,max);
    632           sgaus->SetBit(kCanDelete);
    633           sgaus->SetParNames("Area","#mu","#sigma");
    634           sgaus->SetParameters(integ/rms,mean,rms);
    635           sgaus->SetParLimits(0,0.,integ);
    636           sgaus->SetParLimits(1,min,max);
    637           sgaus->SetParLimits(2,0,width/1.5);
    638           obj3->Fit("sgaus","QLR");
    639           obj3->GetFunction("sgaus")->SetLineColor(kYellow);
    640           break;
    641 
    642         case 2:
    643           TString dgausform = "([0]-[3])/[2]*exp(-0.5*(x-[1])*(x-[1])/[2]/[2])";
    644           dgausform += "+[3]/[5]*exp(-0.5*(x-[4])*(x-[4])/[5]/[5])";
    645           TF1 *dgaus = new TF1("dgaus",dgausform.Data(),min,max);
    646           dgaus->SetBit(kCanDelete);
    647           dgaus->SetParNames("A_{tot}","#mu_{1}","#sigma_{1}","A_{2}","#mu_{2}","#sigma_{2}");
    648           dgaus->SetParameters(integ,(min+mean)/2.,width/4.,
    649                                integ/width/2.,(max+mean)/2.,width/4.);
    650           // The left-sided Gauss
    651           dgaus->SetParLimits(0,integ-1.5,integ+1.5);
    652           dgaus->SetParLimits(1,min+(width/10.),mean);
    653           dgaus->SetParLimits(2,0,width/2.);
    654           // The right-sided Gauss
    655           dgaus->SetParLimits(3,0,integ);
    656           dgaus->SetParLimits(4,mean,max-(width/10.));
    657           dgaus->SetParLimits(5,0,width/2.);
    658           obj3->Fit("dgaus","QLRM");
    659           obj3->GetFunction("dgaus")->SetLineColor(kYellow);
    660           break;
    661          
    662         case 3:
    663           TString tgausform = "([0]-[3]-[6])/[2]*exp(-0.5*(x-[1])*(x-[1])/[2]/[2])";
    664           tgausform += "+[3]/[5]*exp(-0.5*(x-[4])*(x-[4])/[5]/[5])";
    665           tgausform += "+[6]/[8]*exp(-0.5*(x-[7])*(x-[7])/[8]/[8])";
    666           TF1 *tgaus = new TF1("tgaus",tgausform.Data(),min,max);
    667           tgaus->SetBit(kCanDelete);
    668           tgaus->SetParNames("A_{tot}","#mu_{1}","#sigma_{1}",
    669                              "A_{2}","#mu_{2}","#sigma_{2}",
    670                              "A_{3}","#mu_{3}","#sigma_{3}");
    671           tgaus->SetParameters(integ,(min+mean)/2,width/4.,
    672                                integ/width/3.,(max+mean)/2.,width/4.,
    673                                integ/width/3.,mean,width/2.);
    674           // The left-sided Gauss
    675           tgaus->SetParLimits(0,integ-1.5,integ+1.5);
    676           tgaus->SetParLimits(1,min+(width/10.),mean);
    677           tgaus->SetParLimits(2,width/15.,width/2.);
    678           // The right-sided Gauss
    679           tgaus->SetParLimits(3,0.,integ);
    680           tgaus->SetParLimits(4,mean,max-(width/10.));
    681           tgaus->SetParLimits(5,width/15.,width/2.);
    682           // The Gauss describing the outliers
    683           tgaus->SetParLimits(6,0.,integ);
    684           tgaus->SetParLimits(7,min,max);
    685           tgaus->SetParLimits(8,width/4.,width/1.5);
    686           obj3->Fit("tgaus","QLRM");
    687           obj3->GetFunction("tgaus")->SetLineColor(kYellow);
    688           break;
    689         case 4:
    690           obj3->Fit("pol0","Q");
    691           obj3->GetFunction("pol0")->SetLineColor(kYellow);
    692           break;
    693         case 9:
    694           break;
    695         default:
    696           obj3->Fit("gaus","Q");
    697           obj3->GetFunction("gaus")->SetLineColor(kYellow);
    698           break;
    699         }
    700      
    701 
    702 
    703         // Just to get the right (maximum) binning
    704         TH1D *half[4];
    705         half[0] = (TH1D*)obj1->ProjectionS(s1, inner, "Sector 6-1-2 Inner");
    706         half[1] = (TH1D*)obj1->ProjectionS(s2, inner, "Sector 3-4-5 Inner");
    707         half[2] = (TH1D*)obj1->ProjectionS(s1, outer, "Sector 6-1-2 Outer");
    708         half[3] = (TH1D*)obj1->ProjectionS(s2, outer, "Sector 3-4-5 Outer");
    709 
    710         for (Int_t k=0; k<4; k++)     
    711         {
    712             half[k]->SetLineColor(kRed+k);
    713             half[k]->SetDirectory(0);
    714             half[k]->SetBit(kCanDelete);
    715             half[k]->Draw("same");
    716         }
    717 
    718       gPad->Modified();
    719       gPad->Update();
    720      
    721     }
    722 }
    723 
Note: See TracChangeset for help on using the changeset viewer.