Changeset 14245


Ignore:
Timestamp:
06/28/12 12:12:04 (12 years ago)
Author:
Jens Buss
Message:
whitespace changes
File:
1 edited

Legend:

Unmodified
Added
Removed

  • fact/tools/marsmacros/singlepe.C

    r14242 r14245  
    4444struct Single
    4545{
    46     float    fSignal;
    47     float    fTime;
     46  float    fSignal;
     47  float    fTime;
    4848};
    4949
    … …  
    5454    vector<vector<Single>> fData;
    5555
    56 public:
     56  public:
    5757    MSingles(const char *name=NULL, const char *title=NULL) : fIntegrationWindow(0)
    5858    {
    59         fName = name ? name : "MSingles";
     59      fName = name ? name : "MSingles";
    6060    }
    6161
    6262    void InitSize(const UInt_t i)
    6363    {
    64         fData.resize(i);
    65     }
    66 
    67     vector<Single> &operator[](UInt_t i) { return fData[i]; }
    68     vector<Single> &GetVector(UInt_t i)  { return fData[i]; }
    69 
    70     UInt_t GetNumPixels() const { return fData.size(); }
    71 
    72     void SetIntegrationWindow(Int_t w) { fIntegrationWindow = w; }
    73     Int_t GetIntegrationWindow() const { return fIntegrationWindow; }
     64      fData.resize(i);
     65    }
     66
     67    vector<Single> &operator[](UInt_t i)
     68    {
     69      return fData[i];
     70    }
     71    vector<Single> &GetVector(UInt_t i)
     72    {
     73      return fData[i];
     74    }
     75
     76    UInt_t GetNumPixels() const
     77    {
     78      return fData.size();
     79    }
     80
     81    void SetIntegrationWindow(Int_t w)
     82    {
     83      fIntegrationWindow = w;
     84    }
     85    Int_t GetIntegrationWindow() const
     86    {
     87      return fIntegrationWindow;
     88    }
    7489
    7590    Bool_t GetPixelContent(Double_t &val, Int_t idx, const MGeomCam &cam, Int_t type=0) const
    7691    {
    77         return kTRUE;
     92      return kTRUE;
    7893    }
    7994    void   DrawPixelContent(Int_t num) const { }
    … …  
    8499class MH2F : public TH2F
    85100{
    86 public:
     101  public:
    87102    Int_t Fill(Double_t x,Double_t y)
    88103    {
    89         Int_t binx, biny, bin;
    90         fEntries++;
    91         binx = TMath::Nint(x+1);
    92         biny = fYaxis.FindBin(y);
    93         bin  = biny*(fXaxis.GetNbins()+2) + binx;
    94         AddBinContent(bin);
    95         if (fSumw2.fN) ++fSumw2.fArray[bin];
    96         if (biny == 0 || biny > fYaxis.GetNbins()) {
    97             if (!fgStatOverflows) return -1;
    98         }
    99         ++fTsumw;
    100         ++fTsumw2;
    101         fTsumwy  += y;
    102         fTsumwy2 += y*y;
    103         return bin;
     104      Int_t binx, biny, bin;
     105      fEntries++;
     106      binx = TMath::Nint(x+1);
     107      biny = fYaxis.FindBin(y);
     108      bin  = biny*(fXaxis.GetNbins()+2) + binx;
     109      AddBinContent(bin);
     110
     111      if (fSumw2.fN) ++fSumw2.fArray[bin];
     112
     113      if (biny == 0 || biny > fYaxis.GetNbins())
     114        {
     115          if (!fgStatOverflows) return -1;
     116        }
     117
     118      ++fTsumw;
     119      ++fTsumw2;
     120      fTsumwy  += y;
     121      fTsumwy2 += y*y;
     122      return bin;
    104123    }
    105124    ClassDef(MH2F, 1);
    … …  
    108127class MProfile2D : public TProfile2D
    109128{
    110 public:
     129  public:
    111130    Int_t Fill(Double_t x, Double_t y, Double_t z)
    112131    {
    113         Int_t bin,binx,biny;
    114         fEntries++;
    115         binx =TMath::Nint(x+1);
    116         biny =fYaxis.FindBin(y);
    117         bin = GetBin(binx, biny);
    118         fArray[bin] += z;
    119         fSumw2.fArray[bin] += z*z;
    120         fBinEntries.fArray[bin] += 1;
    121         if (fBinSumw2.fN)  fBinSumw2.fArray[bin] += 1;
    122         if (biny == 0 || biny > fYaxis.GetNbins()) {
    123             if (!fgStatOverflows) return -1;
    124         }
    125         ++fTsumw;
    126         ++fTsumw2;
    127         fTsumwy  += y;
    128         fTsumwy2 += y*y;
    129         fTsumwz  += z;
    130         fTsumwz2 += z*z;
    131         return bin;
     132      Int_t bin,binx,biny;
     133      fEntries++;
     134      binx =TMath::Nint(x+1);
     135      biny =fYaxis.FindBin(y);
     136      bin = GetBin(binx, biny);
     137      fArray[bin] += z;
     138      fSumw2.fArray[bin] += z*z;
     139      fBinEntries.fArray[bin] += 1;
     140
     141      if (fBinSumw2.fN)  fBinSumw2.fArray[bin] += 1;
     142
     143      if (biny == 0 || biny > fYaxis.GetNbins())
     144        {
     145          if (!fgStatOverflows) return -1;
     146        }
     147
     148      ++fTsumw;
     149      ++fTsumw2;
     150      fTsumwy  += y;
     151      fTsumwy2 += y*y;
     152      fTsumwz  += z;
     153      fTsumwz2 += z*z;
     154      return bin;
    132155    }
    133156    ClassDef(MProfile2D, 1);
    … …  
    147170    MPedestalSubtractedEvt *fEvent;     //!
    148171
    149 public:
     172  public:
    150173    MHSingles() : fNumEntries(0), fSingles(0), fEvent(0)
    151174    {
    152         fName = "MHSingles";
    153 
    154         fSignal.SetName("Signal");
    155         fSignal.SetTitle("pulse integral spectrum");
    156         fSignal.SetXTitle("pixel [SoftID]");
    157         fSignal.SetYTitle("time [a.u.]");
    158         fSignal.SetDirectory(NULL);
    159 
    160         fTime.SetName("Time");
    161         fTime.SetTitle("arival time spectrum");
    162         fTime.SetXTitle("pixel [SoftID]");
    163         fTime.SetYTitle("time [a.u.]");
    164         fTime.SetDirectory(NULL);
    165 
    166         fPulse.SetName("Pulse");
    167         fPulse.SetTitle("average pulse shape spectrum");
    168         fPulse.SetXTitle("pixel [SoftID]");
    169         fPulse.SetYTitle("time [a.u.]");
    170         fPulse.SetDirectory(NULL);
     175      fName = "MHSingles";
     176
     177      fSignal.SetName("Signal");
     178      fSignal.SetTitle("pulse integral spectrum");
     179      fSignal.SetXTitle("pixel [SoftID]");
     180      fSignal.SetYTitle("time [a.u.]");
     181      fSignal.SetDirectory(NULL);
     182
     183      fTime.SetName("Time");
     184      fTime.SetTitle("arival time spectrum");
     185      fTime.SetXTitle("pixel [SoftID]");
     186      fTime.SetYTitle("time [a.u.]");
     187      fTime.SetDirectory(NULL);
     188
     189      fPulse.SetName("Pulse");
     190      fPulse.SetTitle("average pulse shape spectrum");
     191      fPulse.SetXTitle("pixel [SoftID]");
     192      fPulse.SetYTitle("time [a.u.]");
     193      fPulse.SetDirectory(NULL);
    171194    }
    172195
    173196    Bool_t SetupFill(const MParList *plist)
    174197    {
    175         fSingles = (MSingles*)plist->FindObject("MSingles");
    176         if (!fSingles)
    177         {
    178             *fLog << /* err << */ "MSingles not found... abort." << endl;
    179             return kFALSE;
    180         }
    181 
    182         fEvent = (MPedestalSubtractedEvt*)plist->FindObject("MPedestalSubtractedEvt");
    183         if (!fEvent)
    184         {
    185             *fLog << /* err << */ "MPedestalSubtractedEvt not found... abort." << endl;
    186             return kFALSE;
    187         }
    188 
    189         fNumEntries = 0;
    190 
    191         return kTRUE;
     198      fSingles = (MSingles*)plist->FindObject("MSingles");
     199
     200      if (!fSingles)
     201        {
     202          *fLog << /* err << */ "MSingles not found... abort." << endl;
     203          return kFALSE;
     204        }
     205
     206      fEvent = (MPedestalSubtractedEvt*)plist->FindObject("MPedestalSubtractedEvt");
     207
     208      if (!fEvent)
     209        {
     210          *fLog << /* err << */ "MPedestalSubtractedEvt not found... abort." << endl;
     211          return kFALSE;
     212        }
     213
     214      fNumEntries = 0;
     215
     216      return kTRUE;
    192217    }
    193218    Bool_t ReInit(MParList *plist)
    194219    {
    195         const MRawRunHeader *header = (MRawRunHeader*)plist->FindObject("MRawRunHeader");
    196         if (!header)
    197         {
    198             *fLog << /* err << */ "MRawRunHeader not found... abort." << endl;
    199             return kFALSE;
    200         }
    201 
    202         const Int_t w = fSingles->GetIntegrationWindow();
    203 
    204         MBinning binsx, binsy;
    205         binsx.SetEdges(fSingles->GetNumPixels(), -0.5, fSingles->GetNumPixels()-0.5);
    206         binsy.SetEdges(2*150, -2*7.5*w, 2*(60-7.5)*w);
    207 
    208         MH::SetBinning(fSignal, binsx, binsy);
    209 
    210         const UShort_t roi = header->GetNumSamples();
    211 
    212         // Correct for DRS timing!!
    213         MBinning binst(roi, -0.5, roi-0.5);
    214         MH::SetBinning(fTime, binsx, binst);
    215 
    216         MBinning binspy(2*roi, -roi-0.5, roi-0.5);
    217         MH::SetBinning(fPulse, binsx, binspy);
    218 
    219         return kTRUE;
     220      const MRawRunHeader *header = (MRawRunHeader*)plist->FindObject("MRawRunHeader");
     221
     222      if (!header)
     223        {
     224          *fLog << /* err << */ "MRawRunHeader not found... abort." << endl;
     225          return kFALSE;
     226        }
     227
     228      const Int_t w = fSingles->GetIntegrationWindow();
     229
     230      MBinning binsx, binsy;
     231
     232      binsx.SetEdges(fSingles->GetNumPixels(), -0.5, fSingles->GetNumPixels()-0.5);
     233
     234      binsy.SetEdges(2*150, -2*7.5*w, 2*(60-7.5)*w);
     235
     236      MH::SetBinning(fSignal, binsx, binsy);
     237
     238      const UShort_t roi = header->GetNumSamples();
     239
     240      // Correct for DRS timing!!
     241      MBinning binst(roi, -0.5, roi-0.5);
     242
     243      MH::SetBinning(fTime, binsx, binst);
     244
     245      MBinning binspy(2*roi, -roi-0.5, roi-0.5);
     246
     247      MH::SetBinning(fPulse, binsx, binspy);
     248
     249      return kTRUE;
    220250    }
    221251
    222252    Int_t Fill(const MParContainer *par, const Stat_t weight=1)
    223253    {
    224         const Float_t *ptr = fEvent->GetSamples(0);
    225         const Short_t  roi = fEvent->GetNumSamples();
    226 
    227         for (unsigned int i=0; i<fSingles->GetNumPixels(); i++)
    228         {
    229             const vector<Single> &result = fSingles->GetVector(i);
    230 
    231             for (unsigned int j=0; j<result.size(); j++)
     254      const Float_t *ptr = fEvent->GetSamples(0);
     255      const Short_t  roi = fEvent->GetNumSamples();
     256
     257      for (unsigned int i=0; i<fSingles->GetNumPixels(); i++)
     258        {
     259          const vector<Single> &result = fSingles->GetVector(i);
     260
     261          for (unsigned int j=0; j<result.size(); j++)
    232262            {
    233                 fSignal.Fill(i, result[j].fSignal);
    234                 fTime.Fill  (i, result[j].fTime);
    235 
    236                 if (!ptr)
    237                     continue;
    238 
    239                 const Short_t tm = floor(result[j].fTime);
    240 
    241                 for (int s=0; s<roi; s++)
    242                     fPulse.Fill(i, s-tm, ptr[s]);
     263              fSignal.Fill(i, result[j].fSignal);
     264              fTime.Fill  (i, result[j].fTime);
     265
     266              if (!ptr)
     267                continue;
     268
     269              const Short_t tm = floor(result[j].fTime);
     270
     271              for (int s=0; s<roi; s++)
     272                fPulse.Fill(i, s-tm, ptr[s]);
    243273            }
    244274
    245             ptr += roi;
    246         }
    247 
    248         fNumEntries++;
    249 
    250         return kTRUE;
    251     }
    252 
    253     TH2 *GetSignal() { return &fSignal; }
    254     TH2 *GetTime()   { return &fTime; }
    255     TH2 *GetPulse()  { return &fPulse; }
    256 
    257     UInt_t GetNumEntries() const { return fNumEntries; }
     275          ptr += roi;
     276        }
     277
     278      fNumEntries++;
     279
     280      return kTRUE;
     281    }
     282
     283    TH2 *GetSignal()
     284    {
     285      return &fSignal;
     286    }
     287    TH2 *GetTime()
     288    {
     289      return &fTime;
     290    }
     291    TH2 *GetPulse()
     292    {
     293      return &fPulse;
     294    }
     295
     296    UInt_t GetNumEntries() const
     297    {
     298      return fNumEntries;
     299    }
    258300
    259301    void Draw(Option_t *opt)
    260302    {
    261         TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this);
    262 
    263         AppendPad("");
    264 
    265         pad->Divide(2,2);
    266 
    267         pad->cd(1);
    268         gPad->SetBorderMode(0);
    269         gPad->SetFrameBorderMode(0);
    270         fSignal.Draw("colz");
    271 
    272         pad->cd(2);
    273         gPad->SetBorderMode(0);
    274         gPad->SetFrameBorderMode(0);
    275         fTime.Draw("colz");
    276 
    277         pad->cd(3);
    278         gPad->SetBorderMode(0);
    279         gPad->SetFrameBorderMode(0);
    280         fPulse.Draw("colz");
     303      TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this);
     304
     305      AppendPad("");
     306
     307      pad->Divide(2,2);
     308
     309      pad->cd(1);
     310      gPad->SetBorderMode(0);
     311      gPad->SetFrameBorderMode(0);
     312      fSignal.Draw("colz");
     313
     314      pad->cd(2);
     315      gPad->SetBorderMode(0);
     316      gPad->SetFrameBorderMode(0);
     317      fTime.Draw("colz");
     318
     319      pad->cd(3);
     320      gPad->SetBorderMode(0);
     321      gPad->SetFrameBorderMode(0);
     322      fPulse.Draw("colz");
    281323    }
    282324
    … …  
    295337Int_t PreProcess(MParList *plist)
    296338{
    297     fSingles = (MSingles*)plist->FindCreateObj("MSingles");
    298     if (!fSingles)
    299         return kFALSE;
    300 
    301     fEvent = (MPedestalSubtractedEvt*)plist->FindObject("MPedestalSubtractedEvt");
    302     if (!fEvent)
    303     {
    304         *fLog << /* err << */ "MPedestalSubtractedEvt not found... abort." << endl;
    305         return kFALSE;
    306     }
    307 
    308     d->AddTab("Fluct");
    309     fluct2->Draw();
    310     fluct1->Draw("same");
    311 
    312     fSingles->SetIntegrationWindow(20);
    313 
    314     //if (weights.GetN()==0)
    315     //    return kFALSE;
    316 
    317     return kTRUE;
     339  fSingles = (MSingles*)plist->FindCreateObj("MSingles");
     340
     341  if (!fSingles)
     342    return kFALSE;
     343
     344  fEvent = (MPedestalSubtractedEvt*)plist->FindObject("MPedestalSubtractedEvt");
     345
     346  if (!fEvent)
     347    {
     348      *fLog << /* err << */ "MPedestalSubtractedEvt not found... abort." << endl;
     349      return kFALSE;
     350    }
     351
     352  d->AddTab("Fluct");
     353  fluct2->Draw();
     354  fluct1->Draw("same");
     355
     356  fSingles->SetIntegrationWindow(20);
     357
     358  //if (weights.GetN()==0)
     359  //    return kFALSE;
     360
     361  return kTRUE;
    318362}
    319363
    320364Int_t Process()
    321365{
    322     const UInt_t roi = fEvent->GetNumSamples();
    323 
    324     const size_t navg             = 10;
    325     const float  threshold        = 5;
    326     const UInt_t start_skip       = 20;
    327     const UInt_t end_skip         = 10;
    328     const UInt_t integration_size = 10*3;
    329     const UInt_t max_search_window = 30;
    330 
    331     vector<float> val(roi-navg);//result of first sliding average
    332     for (int pix=0; pix<fEvent->GetNumPixels(); pix++)
    333     {
    334         vector<Single> &result = fSingles->GetVector(pix);
    335         result.clear();
    336 
    337         const Float_t *ptr = fEvent->GetSamples(pix);
    338 
    339         //Sliding window
    340         for (size_t i=0; i<roi-navg; i++)
    341         {
    342             val[i] = 0;
    343             for (size_t j=i; j<i+navg; j++)
    344                 val[i] += ptr[j];
    345             val[i] /= navg;
    346         }
    347 
    348         //peak finding via threshold
    349         for (UInt_t i=start_skip; i<roi-navg-end_skip-30; i++)
    350         {
    351             //search for threshold crossings
    352             if (val[i+0]>threshold ||
    353                 val[i+4]>threshold ||
    354 
    355                 val[i+5]<threshold ||
    356                 val[i+9]<threshold)
    357                 continue;
    358 
    359             //search for maximum after threshold crossing
    360             UInt_t k_max = i+5;
    361             for (UInt_t k=i+5; k<i+max_search_window; k++)
     366  const UInt_t roi = fEvent->GetNumSamples();
     367
     368  const size_t navg             = 10;
     369  const float  threshold        = 5;
     370  const UInt_t start_skip       = 20;
     371  const UInt_t end_skip         = 10;
     372  const UInt_t integration_size = 10*3;
     373  const UInt_t max_search_window = 30;
     374
     375  vector<float> val(roi-navg);//result of first sliding average
     376
     377  for (int pix=0; pix<fEvent->GetNumPixels(); pix++)
     378    {
     379      vector<Single> &result = fSingles->GetVector(pix);
     380      result.clear();
     381
     382      const Float_t *ptr = fEvent->GetSamples(pix);
     383
     384      //Sliding window
     385      for (size_t i=0; i<roi-navg; i++)
     386        {
     387          val[i] = 0;
     388
     389          for (size_t j=i; j<i+navg; j++)
     390            val[i] += ptr[j];
     391
     392          val[i] /= navg;
     393        }
     394
     395      //peak finding via threshold
     396      for (UInt_t i=start_skip; i<roi-navg-end_skip-30; i++)
     397        {
     398          //search for threshold crossings
     399          if (val[i+0]>threshold ||
     400              val[i+4]>threshold ||
     401
     402              val[i+5]<threshold ||
     403              val[i+9]<threshold)
     404            continue;
     405
     406          //search for maximum after threshold crossing
     407          UInt_t k_max = i+5;
     408
     409          for (UInt_t k=i+5; k<i+max_search_window; k++)
    362410            {
    363                 if (val[k] > val[k_max])
    364                     k_max = k;
     411              if (val[k] > val[k_max])
     412                k_max = k;
    365413            }
    366414
    367             if (k_max == i+5 || k_max == i + max_search_window-1) continue;
    368 
    369             //search for half maximum before maximum
    370             UInt_t k_half_max = 0;
    371             for (UInt_t k=k_max; k>k_max-25; k--)
     415          if (k_max == i+5 || k_max == i + max_search_window-1) continue;
     416
     417          //search for half maximum before maximum
     418          UInt_t k_half_max = 0;
     419
     420          for (UInt_t k=k_max; k>k_max-25; k--)
    372421            {
    373                 if (val[k-1] < val[k_max]/2 &&
    374                     val[k] >= val[k_max]/2 )
     422              if (val[k-1] < val[k_max]/2 &&
     423                  val[k] >= val[k_max]/2 )
    375424                {
    376                     k_half_max = k;
    377                     break;
     425                  k_half_max = k;
     426                  break;
    378427                }
    379428            }
    380             if (k_half_max > roi-navg-end_skip-max_search_window - integration_size)
    381                 continue;
    382             if (k_half_max == 0)
    383                 continue;
    384             if (k_max - k_half_max > 14)
    385                 continue;
    386 
    387             fluct2->Fill(k_max - k_half_max);
    388 
    389             // Evaluate arrival time more precisely!!!
    390             // Get a better integration window
    391 
    392             Single single;
    393             single.fSignal = 0;
    394             single.fTime   = k_half_max;
    395 
    396 
    397             // Crossing of the threshold is at 5
    398             for (UInt_t j=k_half_max; j<k_half_max+integration_size; j++)
    399                 single.fSignal += ptr[j];
    400             result.push_back(single);
    401 
    402             i += 5+integration_size;
    403         }
    404     }
    405     return kTRUE;
     429
     430          if (k_half_max > roi-navg-end_skip-max_search_window - integration_size)
     431            continue;
     432
     433          if (k_half_max == 0)
     434            continue;
     435
     436          if (k_max - k_half_max > 14)
     437            continue;
     438
     439          fluct2->Fill(k_max - k_half_max);
     440
     441          // Evaluate arrival time more precisely!!!
     442          // Get a better integration window
     443
     444          Single single;
     445          single.fSignal = 0;
     446          single.fTime   = k_half_max;
     447
     448
     449          // Crossing of the threshold is at 5
     450          for (UInt_t j=k_half_max; j<k_half_max+integration_size; j++)
     451            single.fSignal += ptr[j];
     452
     453          result.push_back(single);
     454
     455          i += 5+integration_size;
     456        }
     457    }
     458
     459  return kTRUE;
    406460}
    407461
    408462Int_t PostProcess()
    409463{
    410     return kTRUE;
     464  return kTRUE;
    411465}
    412466
    … …  
    418472
    419473Double_t MedianOfH1 (
    420         TH1*            inputHisto
    421         );
     474  TH1*            inputHisto
     475);
    422476
    423477int singlepe(
    424478//        const char *runfile, const char *drsfile, const char *outfile
    425         )
     479)
    426480{
    427481
    428     // ======================================================
    429 
    430     const char *drsfile = "/fact/raw/2012/05/18/20120518_003.drs.fits.gz";
     482  // ======================================================
     483
     484  const char *drsfile = "/fact/raw/2012/05/18/20120518_003.drs.fits.gz";
    431485//    const char *drsfile = "/fact/raw/2012/06/25/20120625_112.drs.fits.gz";
    432486
    433     MDirIter iter;
    434     iter.AddDirectory("/fact/raw/2012/05/18", "20120518_005.fits.gz");
     487  MDirIter iter;
     488  iter.AddDirectory("/fact/raw/2012/05/18", "20120518_005.fits.gz");
    435489//    iter.AddDirectory(gSystem->DirName(runfile), gSystem->BaseName(runfile));
    436490
    … …  
    449503//     iter.AddDirectory("/fact/raw/2012/06/25", "20120625_119.fits.gz");
    450504
    451     // ======================================================
    452 
    453     // true:  Display correctly mapped pixels in the camera displays
    454     //        but the value-vs-index plot is in software/spiral indices
    455     // false: Display pixels in hardware/linear indices,
    456     //        but the order is the camera display is distorted
    457     bool usemap = true;
    458 
    459     // map file to use (get that from La Palma!)
    460     const char *map = usemap ? "/scratch_nfs/FACTmap111030.txt" : NULL;
    461 
    462     // ------------------------------------------------------
    463 
    464     Long_t max1 = 0;
    465 
    466     // ======================================================
    467 
    468     if (map && gSystem->AccessPathName(map, kFileExists))
    469     {
    470         gLog << "ERROR - Cannot access mapping file '" << map << "'" << endl;
    471         return 1;
    472     }
    473     if (gSystem->AccessPathName(drsfile, kFileExists))
    474     {
    475         gLog << "ERROR - Cannot access drsfile file '" << drsfile << "'" << endl;
    476         return 2;
    477     }
    478 
    479     // --------------------------------------------------------------------------------
    480 
    481     gLog.Separator("Singles");
    482     gLog << "Extract singles with '" << drsfile << "'" << endl;
    483     gLog << endl;
    484 
    485     // ------------------------------------------------------
    486 
    487     MDrsCalibration drscalib300;
    488     if (!drscalib300.ReadFits(drsfile))
    489         return 4;
    490 
    491     // ------------------------------------------------------
    492 
    493     iter.Sort();
    494     iter.Print();
    495 
    496     // ======================================================
    497 
    498     //MStatusDisplay *d = new MStatusDisplay;
    499 
    500     MBadPixelsCam badpixels;
    501     badpixels.InitSize(1440);
    502     badpixels[ 424].SetUnsuitable(MBadPixelsPix::kUnsuitable);
    503     badpixels[ 583].SetUnsuitable(MBadPixelsPix::kUnsuitable);
    504     badpixels[ 830].SetUnsuitable(MBadPixelsPix::kUnsuitable);
    505     badpixels[ 923].SetUnsuitable(MBadPixelsPix::kUnsuitable);
    506     badpixels[1208].SetUnsuitable(MBadPixelsPix::kUnsuitable);
    507     badpixels[1399].SetUnsuitable(MBadPixelsPix::kUnsuitable);
    508 
    509     //  Twin pixel
    510     //     113
    511     //     115
    512     //     354
    513     //     423
    514     //    1195
    515     //    1393
    516 
    517     //MDrsCalibrationTime timecam;
    518 
    519     MH::SetPalette();
    520 
    521     // ======================================================
    522 
    523     gLog << endl;
    524     gLog.Separator("Processing external light pulser run");
    525 
    526     MTaskList tlist1;
    527 
    528     MSingles singles;
    529     MRawRunHeader header;
    530 
    531     MParList plist1;
    532     plist1.AddToList(&tlist1);
    533     plist1.AddToList(&drscalib300);
    534     plist1.AddToList(&badpixels);
    535     plist1.AddToList(&singles);
    536     plist1.AddToList(&header);
    537 
    538     TString Title;
    539     Title =  iter.Next();
    540     iter.Reset();
    541     Title += "; ";
    542     Title += max1;
    543 
    544     MEvtLoop loop1(Title);
    545     loop1.SetDisplay(d);
    546     loop1.SetParList(&plist1);
    547 
    548     // ------------------ Setup the tasks ---------------
    549 
    550     MRawFitsRead read1;
    551     read1.LoadMap(map);
    552     read1.AddFiles(iter);
    553 
    554     MContinue cont1("(MRawEvtHeader.GetTriggerID&0xff00)!=0x100", "SelectCal");
    555 
    556     MGeomApply apply1;
    557 
    558     MDrsCalibApply drsapply1;
    559 
    560     MTaskInteractive mytask;
    561     mytask.SetPreProcess(PreProcess);
    562     mytask.SetProcess(Process);
    563     mytask.SetPostProcess(PostProcess);
    564 
    565     MFillH fill("MHSingles");
    566 
    567     // ------------------ Setup eventloop and run analysis ---------------
    568 
    569     tlist1.AddToList(&read1);
     505  // ======================================================
     506
     507  // true:  Display correctly mapped pixels in the camera displays
     508  //        but the value-vs-index plot is in software/spiral indices
     509  // false: Display pixels in hardware/linear indices,
     510  //        but the order is the camera display is distorted
     511  bool usemap = true;
     512
     513  // map file to use (get that from La Palma!)
     514  const char *map = usemap ? "/scratch_nfs/FACTmap111030.txt" : NULL;
     515
     516  // ------------------------------------------------------
     517
     518  Long_t max1 = 0;
     519
     520  // ======================================================
     521
     522  if (map && gSystem->AccessPathName(map, kFileExists))
     523    {
     524      gLog << "ERROR - Cannot access mapping file '" << map << "'" << endl;
     525      return 1;
     526    }
     527
     528  if (gSystem->AccessPathName(drsfile, kFileExists))
     529    {
     530      gLog << "ERROR - Cannot access drsfile file '" << drsfile << "'" << endl;
     531      return 2;
     532    }
     533
     534  // --------------------------------------------------------------------------------
     535
     536  gLog.Separator("Singles");
     537  gLog << "Extract singles with '" << drsfile << "'" << endl;
     538  gLog << endl;
     539
     540  // ------------------------------------------------------
     541
     542  MDrsCalibration drscalib300;
     543
     544  if (!drscalib300.ReadFits(drsfile))
     545    return 4;
     546
     547  // ------------------------------------------------------
     548
     549  iter.Sort();
     550  iter.Print();
     551
     552  // ======================================================
     553
     554  //MStatusDisplay *d = new MStatusDisplay;
     555
     556  MBadPixelsCam badpixels;
     557  badpixels.InitSize(1440);
     558  badpixels[ 424].SetUnsuitable(MBadPixelsPix::kUnsuitable);
     559  badpixels[ 583].SetUnsuitable(MBadPixelsPix::kUnsuitable);
     560  badpixels[ 830].SetUnsuitable(MBadPixelsPix::kUnsuitable);
     561  badpixels[ 923].SetUnsuitable(MBadPixelsPix::kUnsuitable);
     562  badpixels[1208].SetUnsuitable(MBadPixelsPix::kUnsuitable);
     563  badpixels[1399].SetUnsuitable(MBadPixelsPix::kUnsuitable);
     564
     565  //  Twin pixel
     566  //     113
     567  //     115
     568  //     354
     569  //     423
     570  //    1195
     571  //    1393
     572
     573  //MDrsCalibrationTime timecam;
     574
     575  MH::SetPalette();
     576
     577  // ======================================================
     578
     579  gLog << endl;
     580  gLog.Separator("Processing external light pulser run");
     581
     582  MTaskList tlist1;
     583
     584  MSingles singles;
     585  MRawRunHeader header;
     586
     587  MParList plist1;
     588  plist1.AddToList(&tlist1);
     589  plist1.AddToList(&drscalib300);
     590  plist1.AddToList(&badpixels);
     591  plist1.AddToList(&singles);
     592  plist1.AddToList(&header);
     593
     594  TString Title;
     595  Title =  iter.Next();
     596  iter.Reset();
     597  Title += "; ";
     598  Title += max1;
     599
     600  MEvtLoop loop1(Title);
     601  loop1.SetDisplay(d);
     602  loop1.SetParList(&plist1);
     603
     604  // ------------------ Setup the tasks ---------------
     605
     606  MRawFitsRead read1;
     607  read1.LoadMap(map);
     608  read1.AddFiles(iter);
     609
     610  MContinue cont1("(MRawEvtHeader.GetTriggerID&0xff00)!=0x100", "SelectCal");
     611
     612  MGeomApply apply1;
     613
     614  MDrsCalibApply drsapply1;
     615
     616  MTaskInteractive mytask;
     617  mytask.SetPreProcess(PreProcess);
     618  mytask.SetProcess(Process);
     619  mytask.SetPostProcess(PostProcess);
     620
     621  MFillH fill("MHSingles");
     622
     623  // ------------------ Setup eventloop and run analysis ---------------
     624
     625  tlist1.AddToList(&read1);
    570626//    tlist1.AddToList(&cont1);
    571     tlist1.AddToList(&apply1);
    572     tlist1.AddToList(&drsapply1);
    573     tlist1.AddToList(&mytask);
    574     tlist1.AddToList(&fill);
    575 
    576     if (!loop1.Eventloop(max1))
    577         return 10;
    578 
    579     if (!loop1.GetDisplay())
    580         return 11;
    581 
    582     gStyle->SetOptFit(kTRUE);
    583 
    584     MHSingles* h = (MHSingles*) plist1.FindObject("MHSingles");
    585     if (!h)
    586         return 0;
    587 
    588     TH2 *hsignal = h->GetSignal();
    589     TH2 *htime   = h->GetTime();
    590     TH2 *hpulse  = h->GetPulse();
    591 
    592     d->AddTab("Time");
    593     TH1 *ht = htime->ProjectionY();
    594     ht->SetYTitle("Counts [a.u.]");
    595     ht->Scale(1./1440);
    596     ht->Draw();
    597 
    598     d->AddTab("Pulse");
    599     TH1 *hp = hpulse->ProjectionY();
    600     hp->SetYTitle("Counts [a.u.]");
    601     hp->Scale(1./1440);
    602     hp->Draw();
    603 
    604 
    605     // ------------------ Spectrum Fit ---------------
    606     // AFTER this the Code for the spektrum fit follows
    607     // access the spektrumhistogram by the pointer *hist
    608     UInt_t maxOrder     = 8;
    609 
    610 
    611     MGeomCamFACT fact;
    612     MHCamera hcam(fact);
    613     MHCamera hcam2(fact);
    614 
    615     //built an array of Amplitude spectra
    616     TH1F hAmplitude ("Rate",      "Average number of dark counts per event",
    617                      200,  0,  20);
    618     hAmplitude.SetXTitle("Amplitude [a.u.]");
    619     hAmplitude.SetYTitle("Rate");
    620 
    621     TH1F hGain      ("Gain",      "Gain distribution",
    622                      50,  100,   300);
    623     hGain.SetXTitle("gain [a.u.]");
    624     hGain.SetYTitle("Rate");
    625 
    626     TH1F hGainRMS   ("RelSigma",  "Distribution of rel. Sigma",
    627                      100,   0,   30);
    628     hGainRMS.SetXTitle("gainRMS [a.u.]");
    629     hGainRMS.SetYTitle("Rate");
    630     hGainRMS.SetStats(false);
    631 
    632     TH1F hCrosstlk  ("Crosstalk", "Distribution of crosstalk probability",
    633                      100,   0,    25);
    634     hCrosstlk.SetXTitle("crosstalk [a.u.]");
    635     hCrosstlk.SetYTitle("Rate");
    636 
    637     TH1F hOffset    ("Offset",    "Baseline offset distribution",
    638                      50, -7.5,  2.5);
    639     hOffset.SetXTitle("baseline offset [a.u.]");
    640     hOffset.SetYTitle("Rate");
    641 
    642     TH1F hFitProb   ("FitProb",   "FitProb distribution",
    643                      50,   0,  100);
    644     hFitProb.SetXTitle("FitProb p-value [a.u.]");
    645     hFitProb.SetYTitle("Rate");
    646     hFitProb.SetStats(false);
    647 
    648 
    649     TH1F hSum       ("Sum1",      "Sum of all pixel's' integral spectra",
    650                      100,    0,  2000);
    651     hSum.SetXTitle("pulse integral [a.u.]");
    652     hSum.SetYTitle("Rate");
    653     hSum.SetStats(false);
    654     hSum.Sumw2();
    655 
    656 
    657     TH1F hSumScale  ("Sum2",
    658                      "Sum of all pixel's' integral spectra (scaled with gain)",
    659                      100,    0.05,   10.05);
    660     hSumScale.SetXTitle("pulse integral [pe]");
    661     hSumScale.SetYTitle("Rate");
    662     hSumScale.SetStats(false);
    663     hSumScale.Sumw2();
    664 
    665     // define fit function for Amplitudespectrum
    666     TF1 func("spektrum", fcn, 0, 1000, 5);
    667     func.SetParNames("Maximum", "Gain", "Sigma", "XtalkProb", "Offset");
    668     func.SetLineColor(kBlue);
    669 
    670     // define fit function for Amplitudespectrum
    671     TF1 funcSum("sum_spektrum", fcn, 0, 2000, 5);
    672     funcSum.SetParNames("Maximum", "Gain", "Sigma", "XtalkProb", "Offset", "PowOfXtalk");
    673     funcSum.SetLineColor(kBlue);
    674 
    675     // define fit function for Amplitudespectrum
    676     TF1 funcScaled("gain_sum_spektrum", fcn, 0, 10, 5);
    677     funcScaled.SetParNames("Maximum", "Gain", "Sigma", "XtalkProb", "Offset");
    678     funcScaled.SetLineColor(kBlue);
    679 
    680 
    681 
    682     // Map for which pixel shall be plotted and which not
    683     TArrayC usePixel(1440);
    684     memset(usePixel.GetArray(), 1, 1440);
    685 
    686     // List of Pixel that should be ignored in camera view
    687     usePixel[424]  = 0;
    688     usePixel[923]  = 0;
    689     usePixel[1208] = 0;
    690     usePixel[583]  = 0;
    691     usePixel[830]  = 0;
    692     usePixel[1399] = 0;
    693     usePixel[113]  = 0;
    694     usePixel[115]  = 0;
    695     usePixel[354]  = 0;
    696     usePixel[423]  = 0;
    697     usePixel[1195] = 0;
    698     usePixel[1393] = 0;
     627  tlist1.AddToList(&apply1);
     628  tlist1.AddToList(&drsapply1);
     629  tlist1.AddToList(&mytask);
     630  tlist1.AddToList(&fill);
     631
     632  if (!loop1.Eventloop(max1))
     633    return 10;
     634
     635  if (!loop1.GetDisplay())
     636    return 11;
     637
     638  gStyle->SetOptFit(kTRUE);
     639
     640  MHSingles* h = (MHSingles*) plist1.FindObject("MHSingles");
     641
     642  if (!h)
     643    return 0;
     644
     645  TH2 *hsignal = h->GetSignal();
     646  TH2 *htime   = h->GetTime();
     647  TH2 *hpulse  = h->GetPulse();
     648
     649  d->AddTab("Time");
     650  TH1 *ht = htime->ProjectionY();
     651  ht->SetYTitle("Counts [a.u.]");
     652  ht->Scale(1./1440);
     653  ht->Draw();
     654
     655  d->AddTab("Pulse");
     656  TH1 *hp = hpulse->ProjectionY();
     657  hp->SetYTitle("Counts [a.u.]");
     658  hp->Scale(1./1440);
     659  hp->Draw();
     660
     661
     662  // ------------------ Spectrum Fit ---------------
     663  // AFTER this the Code for the spektrum fit follows
     664  // access the spektrumhistogram by the pointer *hist
     665  UInt_t maxOrder     = 8;
     666
     667
     668  MGeomCamFACT fact;
     669  MHCamera hcam(fact);
     670  MHCamera hcam2(fact);
     671
     672  //built an array of Amplitude spectra
     673  TH1F hAmplitude ("Rate",      "Average number of dark counts per event",
     674                   200,  0,  20);
     675  hAmplitude.SetXTitle("Amplitude [a.u.]");
     676  hAmplitude.SetYTitle("Rate");
     677
     678  TH1F hGain      ("Gain",      "Gain distribution",
     679                   400,  100,   300);
     680  hGain.SetXTitle("gain [a.u.]");
     681  hGain.SetYTitle("Rate");
     682
     683  TH1F hGainRMS   ("RelSigma",  "Distribution of rel. Sigma",
     684                   100,   0,   30);
     685  hGainRMS.SetXTitle("gainRMS [a.u.]");
     686  hGainRMS.SetYTitle("Rate");
     687  hGainRMS.SetStats(false);
     688
     689  TH1F hCrosstlk  ("Crosstalk", "Distribution of crosstalk probability",
     690                   100,   0,    25);
     691  hCrosstlk.SetXTitle("crosstalk [a.u.]");
     692  hCrosstlk.SetYTitle("Rate");
     693
     694  TH1F hOffset    ("Offset",    "Baseline offset distribution",
     695                   50, -7.5,  2.5);
     696  hOffset.SetXTitle("baseline offset [a.u.]");
     697  hOffset.SetYTitle("Rate");
     698
     699  TH1F hFitProb   ("FitProb",   "FitProb distribution",
     700                   50,   0,  100);
     701  hFitProb.SetXTitle("FitProb p-value [a.u.]");
     702  hFitProb.SetYTitle("Rate");
     703  hFitProb.SetStats(false);
     704
     705
     706  TH1F hSum       ("Sum1",      "Sum of all pixel's' integral spectra",
     707                   100,    0,  2000);
     708  hSum.SetXTitle("pulse integral [a.u.]");
     709  hSum.SetYTitle("Rate");
     710  hSum.SetStats(false);
     711  hSum.Sumw2();
     712
     713
     714  TH1F hSumScale  ("Sum2",
     715                   "Sum of all pixel's' integral spectra (scaled with gain)",
     716                   100,    0.05,   10.05);
     717  hSumScale.SetXTitle("pulse integral [pe]");
     718  hSumScale.SetYTitle("Rate");
     719  hSumScale.SetStats(false);
     720  hSumScale.Sumw2();
     721
     722  // define fit function for Amplitudespectrum
     723  TF1 func("spektrum", fcn, 0, 1000, 5);
     724  func.SetParNames("Maximum", "Gain", "Sigma", "XtalkProb", "Offset");
     725  func.SetLineColor(kBlue);
     726
     727  // define fit function for Amplitudespectrum
     728  TF1 funcSum("sum_spektrum", fcn, 0, 2000, 5);
     729  funcSum.SetParNames("Maximum", "Gain", "Sigma", "XtalkProb", "Offset", "PowOfXtalk");
     730  funcSum.SetLineColor(kBlue);
     731
     732  // define fit function for Amplitudespectrum
     733  TF1 funcScaled("gain_sum_spektrum", fcn, 0, 10, 5);
     734  funcScaled.SetParNames("Maximum", "Gain", "Sigma", "XtalkProb", "Offset");
     735  funcScaled.SetLineColor(kBlue);
     736
     737  TF1 funcScaledBL("gain_sum_spektrum", fcn, 0, 10, 5);
     738  funcScaledBL.SetParNames("Maximum", "Gain", "Sigma", "XtalkProb", "Offset");
     739  funcScaledBL.SetLineColor(kRed);
     740
     741
     742
     743  // Map for which pixel shall be plotted and which not
     744  TArrayC usePixel(1440);
     745  memset(usePixel.GetArray(), 1, 1440);
     746
     747  // List of Pixel that should be ignored in camera view
     748  usePixel[424]  = 0;
     749  usePixel[923]  = 0;
     750  usePixel[1208] = 0;
     751  usePixel[583]  = 0;
     752  usePixel[830]  = 0;
     753  usePixel[1399] = 0;
     754  usePixel[113]  = 0;
     755  usePixel[115]  = 0;
     756  usePixel[354]  = 0;
     757  usePixel[423]  = 0;
     758  usePixel[1195] = 0;
     759  usePixel[1393] = 0;
    699760
    700761//    usePixel[990] = 0;
    701762
    702     // Map for which pixel which were suspicous
    703     bool suspicous[1440];
    704     for (int i=0; i<1440; i++) suspicous[i]=false;
    705 
    706     // List of Pixel that should be ignored in camera view
    707     //    suspicous[802]       = true;
    708 
    709     //------------------------------------------------------------------------
    710     //  Fill and calculate sum spectrum
    711 
    712     // Begin of Loop over Pixels
    713     for (UInt_t pixel = 0; pixel < singles.GetNumPixels(); pixel++)
    714     {
    715         //jump to next pixel if the current is marked as faulty
    716         if (usePixel[pixel]==0)
    717             continue;
    718 
    719         TH1D *hist = hsignal->ProjectionY("proj", pixel+1, pixel+1);
    720         hist->SetDirectory(0);
    721         //loop over slices
    722         for (int b=1; b<=hist->GetNbinsX(); b++)
    723         {
    724             //Fill sum spectrum
    725             hSum.Fill(hist->GetBinCenter(b), hist->GetBinContent(b));
    726         }
    727     }
    728     for(UInt_t bin = 1; bin < (UInt_t)hSum.GetNbinsX()+1; bin++)
    729     {
    730         hSum.SetBinError(bin, hSum.GetBinContent(bin)*0.1);
    731     }
    732 
    733     gROOT->SetSelectedPad(0);
    734     TCanvas &c11 = d->AddTab("SumHist");
    735     c11.cd();
    736     gPad->SetLogy();
    737     gPad->SetGridx();
    738     gPad->SetGridy();
    739     hSum.DrawCopy("hist");
    740     //------------------------fit sum spectrum-------------------------------
    741     const Int_t    maxbin   = hSum.GetMaximumBin();
    742     const Double_t binwidth = hSum.GetBinWidth(maxbin);
    743     const Double_t ampl     = hSum.GetBinContent(maxbin);
    744 
    745     double fwhmSum = 0;
    746 
    747     //Calculate full width half Maximum
    748     for (int i=1; i<maxbin; i++)
    749         if (hSum.GetBinContent(maxbin-i)+hSum.GetBinContent(maxbin+i)<ampl*2)
    750         {
    751             fwhmSum = (2*i+1)*binwidth;
    752             break;
    753         }
    754     if (fwhmSum==0)
    755     {
    756         cout << "Could not determine start value for sigma." << endl;
    757     }
    758 
    759     //Set fit parameters
    760     Double_t sum_par[6] =
    761     {
    762         ampl, hSum.GetBinCenter(maxbin), fwhmSum*1.4, 0.1, -10, 1
    763     };
    764     funcSum.SetParameters(sum_par);
    765 
    766     //calculate fit range
    767     const double min = sum_par[1]-fwhmSum*2;
    768     const double max = sum_par[1]*(maxOrder);
    769     funcSum.SetRange(min, max);
    770     funcSum.FixParameter(5,0.1);
    771 
    772     //Fit and draw spectrum
    773     hSum.Fit(&funcSum, "NOQS", "", min, max);
    774     funcSum.DrawCopy("SAME");
    775     funcSum.GetParameters(sum_par);
    776 
    777     //Set Parameters for following fits
     763  // Map for which pixel which were suspicous
     764  bool suspicous[1440];
     765
     766  for (int i=0; i<1440; i++) suspicous[i]=false;
     767
     768  // List of Pixel that should be ignored in camera view
     769  //    suspicous[802]       = true;
     770
     771  //------------------------------------------------------------------------
     772  //  Fill and calculate sum spectrum
     773
     774  // Begin of Loop over Pixels
     775  for (UInt_t pixel = 0; pixel < singles.GetNumPixels(); pixel++)
     776    {
     777      //jump to next pixel if the current is marked as faulty
     778      if (usePixel[pixel]==0)
     779        continue;
     780
     781      TH1D *hist = hsignal->ProjectionY("proj", pixel+1, pixel+1);
     782      hist->SetDirectory(0);
     783
     784      //loop over slices
     785      for (int b=1; b<=hist->GetNbinsX(); b++)
     786        {
     787          //Fill sum spectrum
     788          hSum.Fill(hist->GetBinCenter(b), hist->GetBinContent(b));
     789        }
     790    }
     791
     792  for (UInt_t bin = 1; bin < (UInt_t)hSum.GetNbinsX()+1; bin++)
     793    {
     794      hSum.SetBinError(bin, hSum.GetBinContent(bin)*0.1);
     795    }
     796
     797  gROOT->SetSelectedPad(0);
     798  TCanvas &c11 = d->AddTab("SumHist");
     799  c11.cd();
     800  gPad->SetLogy();
     801  gPad->SetGridx();
     802  gPad->SetGridy();
     803  hSum.DrawCopy("hist");
     804  //------------------------fit sum spectrum-------------------------------
     805  const Int_t    maxbin   = hSum.GetMaximumBin();
     806  const Double_t binwidth = hSum.GetBinWidth(maxbin);
     807  const Double_t ampl     = hSum.GetBinContent(maxbin);
     808
     809  double fwhmSum = 0;
     810
     811  //Calculate full width half Maximum
     812  for (int i=1; i<maxbin; i++)
     813    if (hSum.GetBinContent(maxbin-i)+hSum.GetBinContent(maxbin+i)<ampl*2)
     814      {
     815        fwhmSum = (2*i+1)*binwidth;
     816        break;
     817      }
     818
     819  if (fwhmSum==0)
     820    {
     821      cout << "Could not determine start value for sigma." << endl;
     822    }
     823
     824  //Set fit parameters
     825  Double_t sum_par[6] =
     826  {
     827    ampl, hSum.GetBinCenter(maxbin), fwhmSum*1.4, 0.1, -10, 1
     828  };
     829  funcSum.SetParameters(sum_par);
     830
     831  //calculate fit range
     832  const double min = sum_par[1]-fwhmSum*2;
     833  const double max = sum_par[1]*(maxOrder);
     834  funcSum.SetRange(min, max);
     835  funcSum.FixParameter(5,0.1);
     836
     837  //Fit and draw spectrum
     838  hSum.Fit(&funcSum, "NOQS", "", min, max);
     839  funcSum.DrawCopy("SAME");
     840  funcSum.GetParameters(sum_par);
     841
     842  //Set Parameters for following fits
    778843//    const float  Amplitude      = sum_par[0];
    779     const float  gain           = sum_par[1];
     844  const float  gain           = sum_par[1];
    780845//    const float  GainRMS        = sum_par[2];
    781     const float  Crosstlk       = sum_par[3];
    782     const float  Offset         = sum_par[4];
     846  const float  Crosstlk       = sum_par[3];
     847  const float  Offset         = sum_par[4];
    783848//    const float  PowCrosstlk    = sum_par[5];
    784849
    785850
    786     //--------fit gausses to peaks of sum specetrum -----------
    787 
    788     // create histo for height of Maximum amplitudes vs. pe
    789     double min_bin  =   hSum.GetBinCenter(maxbin);
    790     double max_bin  =   hSum.GetBinCenter((maxOrder-2)*(maxbin));
    791     double bin_width=   (max_bin - min_bin)/10;
    792 
    793     TH1D hMaxHeightsSum("MaxHeightSum",      "peak heights",
    794                         10,  min_bin-bin_width, max_bin*2-bin_width/2 );
    795 
    796     FitGaussiansToSpectrum
    797             (
    798                 maxOrder-2,
    799                 hSum,
    800                 hMaxHeightsSum,
    801                 maxbin,
    802                 fwhmSum,
    803                 gain
    804                 );
    805 
    806     //EOF: fit sum spectrum-----------------------------
    807 
    808     hMaxHeightsSum.SetLineColor(kRed);
     851  //--------fit gausses to peaks of sum specetrum -----------
     852
     853  // create histo for height of Maximum amplitudes vs. pe
     854  double min_bin  =   hSum.GetBinCenter(maxbin);
     855  double max_bin  =   hSum.GetBinCenter((maxOrder-2)*(maxbin));
     856  double bin_width=   (max_bin - min_bin)/10;
     857
     858  TH1D hMaxHeightsSum("MaxHeightSum",      "peak heights",
     859                      10,  min_bin-bin_width, max_bin*2-bin_width/2 );
     860
     861  FitGaussiansToSpectrum
     862  (
     863    maxOrder-2,
     864    hSum,
     865    hMaxHeightsSum,
     866    maxbin,
     867    fwhmSum,
     868    gain
     869  );
     870
     871  //EOF: fit sum spectrum-----------------------------
     872
     873  hMaxHeightsSum.SetLineColor(kRed);
    809874//    hMaxHeightsSum.DrawCopy("SAME");
    810875
    811     //Fit Peak heights
    812     TF1 fExpo1( "fExpo1","expo", min, max);
    813     fExpo1.SetLineColor(kRed);
    814     hMaxHeightsSum.Fit(&fExpo1, "N0QS" );
     876  //Fit Peak heights
     877  TF1 fExpo1( "fExpo1","expo", min, max);
     878  fExpo1.SetLineColor(kRed);
     879  hMaxHeightsSum.Fit(&fExpo1, "N0QS" );
    815880//    hMaxHeightsSum.DrawCopy("SAME");
    816881//    fExpo1.DrawCopy("SAME");
    817882
    818     //  EOF:    Fill and calculate sum spectrum
    819     //------------------------------------------------------------------------
    820 
    821 
    822 
    823     //------------------------------------------------------------------------
    824     //  Gain Calculation for each Pixel
    825 
    826     int count_ok = 0;
    827     // Begin of Loop over Pixels
    828     for (UInt_t pixel = 0; pixel < singles.GetNumPixels(); pixel++)
    829     {
    830 
    831         if (usePixel[pixel]==0)
    832             continue;
    833 
    834         //Projectipon of a certain Pixel out of Ampl.Spectrum
    835         TH1D *hist = hsignal->ProjectionY("proj", pixel+1, pixel+1);
    836         hist->SetDirectory(0);
    837 
    838         for(UInt_t bin = 1; bin < (UInt_t)hSum.GetNbinsX()+1; bin++)
    839         {
    840             hist->SetBinError(bin, hSum.GetBinContent(bin)*0.1);
    841         }
    842 
    843         const Int_t    maxbin   = hist->GetMaximumBin();
    844         const Double_t binwidth = hist->GetBinWidth(maxbin);
    845         const Double_t ampl     = hist->GetBinContent(maxbin);
    846         const Double_t maxPos   = hist->GetBinCenter(maxbin);
    847 
    848         double fwhm = 0;
    849         for (int i=1; i<maxbin; i++)
    850             if (hist->GetBinContent(maxbin-i)+hist->GetBinContent(maxbin+i)<ampl*2)
    851             {
    852                 fwhm = (2*i+1)*binwidth;
    853                 break;
    854             }
    855 
    856         if (fwhm==0)
    857         {
    858             cout << "Could not determine start value for sigma." << endl;
    859             continue;
    860         }
    861 
    862         // Amplitude, Gain, Sigma, Crosstalk Prob, Shift
    863         // par[1] + par[4] is the position of the first maximum
    864         Double_t par[5] =
    865         {
    866             ampl,
    867             hist->GetBinCenter(maxbin),
     883  //  EOF:    Fill and calculate sum spectrum
     884  //------------------------------------------------------------------------
     885
     886
     887
     888  //------------------------------------------------------------------------
     889  //  Gain Calculation for each Pixel
     890
     891  int count_ok = 0;
     892
     893  // Begin of Loop over Pixels
     894  for (UInt_t pixel = 0; pixel < singles.GetNumPixels(); pixel++)
     895    {
     896
     897      if (usePixel[pixel]==0)
     898        continue;
     899
     900      //Projectipon of a certain Pixel out of Ampl.Spectrum
     901      TH1D *hist = hsignal->ProjectionY("proj", pixel+1, pixel+1);
     902      hist->SetDirectory(0);
     903
     904      for (UInt_t bin = 1; bin < (UInt_t)hSum.GetNbinsX()+1; bin++)
     905        {
     906          hist->SetBinError(bin, hSum.GetBinContent(bin)*0.1);
     907        }
     908
     909      const Int_t    maxbin   = hist->GetMaximumBin();
     910
     911      const Double_t binwidth = hist->GetBinWidth(maxbin);
     912
     913      const Double_t ampl     = hist->GetBinContent(maxbin);
     914
     915      const Double_t maxPos   = hist->GetBinCenter(maxbin);
     916
     917      double fwhm = 0;
     918
     919      for (int i=1; i<maxbin; i++)
     920        if (hist->GetBinContent(maxbin-i)+hist->GetBinContent(maxbin+i)<ampl*2)
     921          {
     922            fwhm = (2*i+1)*binwidth;
     923            break;
     924          }
     925
     926      if (fwhm==0)
     927        {
     928          cout << "Could not determine start value for sigma." << endl;
     929          continue;
     930        }
     931
     932      // Amplitude, Gain, Sigma, Crosstalk Prob, Shift
     933      // par[1] + par[4] is the position of the first maximum
     934      Double_t par[5] =
     935      {
     936        ampl,
     937        hist->GetBinCenter(maxbin),
    868938//            gain,
    869             fwhm*2,
     939        fwhm*2,
    870940//            GainRMS,
    871             Crosstlk,
    872             Offset
    873         };
    874 
    875 
    876         func.SetParameters(par);
    877         const double fit_min = par[1]-fwhm*1.4;
    878         const double fit_max = par[1]*maxOrder;
    879         func.SetRange(fit_min-fwhm*2, fit_max);
    880         func.SetParLimits(1, maxPos - fwhm, maxPos + fwhm);
    881 
    882         if (suspicous[pixel] == true)
    883         {
    884             cout << "Maxbin\t" << maxbin << endl;
    885             cout << "min\t" << fit_min << endl;
    886             cout << "max\t" << fit_max << endl;
    887             cout << "Amplitude, 1.Max, Sigma, fwhm:\t"
    888                  << par[0] << "\t"
    889                  << par[1] << "\t"
    890                  << par[2] << "\t"
    891                  << fwhm << "\t" << endl;
    892         }
    893 
    894         //Rebin Projection
    895         hist->Rebin(2);
    896 
    897         //Fit Pixels spectrum
    898         const TFitResultPtr rc = hist->Fit(&func, "N0QS", "", fit_min, fit_max);
    899 
    900         const bool ok = int(rc)==0;
    901 
    902         const double fit_prob   = rc->Prob();
    903         const float  fGain      = func.GetParameter(1);
    904         const float  fAmplitude = func.Integral(0, 7*fGain)/h->GetNumEntries(); //func.GetParameter(0);
    905         const float  f2GainRMS  = func.GetParameter(2);
    906         const float  fCrosstlk  = func.GetParameter(3);
    907         const float  fOffset    = func.GetParameter(4)/singles.GetIntegrationWindow();
     941        Crosstlk,
     942        Offset
     943      };
     944
     945
     946      func.SetParameters(par);
     947      const double fit_min = par[1]-fwhm*1.4;
     948      const double fit_max = par[1]*maxOrder;
     949      func.SetRange(fit_min-fwhm*2, fit_max);
     950      func.SetParLimits(1, maxPos - fwhm, maxPos + fwhm);
     951
     952      if (suspicous[pixel] == true)
     953        {
     954          cout << "Maxbin\t" << maxbin << endl;
     955          cout << "min\t" << fit_min << endl;
     956          cout << "max\t" << fit_max << endl;
     957          cout << "Amplitude, 1.Max, Sigma, fwhm:\t"
     958               << par[0] << "\t"
     959               << par[1] << "\t"
     960               << par[2] << "\t"
     961               << fwhm << "\t" << endl;
     962        }
     963
     964      //Rebin Projection
     965      hist->Rebin(2);
     966
     967      //Fit Pixels spectrum
     968      const TFitResultPtr rc = hist->Fit(&func, "N0QS", "", fit_min, fit_max);
     969
     970      const bool ok = int(rc)==0;
     971
     972      const double fit_prob   = rc->Prob();
     973      const float  fGain      = func.GetParameter(1);
     974      const float  fAmplitude = func.Integral(0, 7*fGain)/h->GetNumEntries(); //func.GetParameter(0);
     975      const float  f2GainRMS  = func.GetParameter(2);
     976      const float  fCrosstlk  = func.GetParameter(3);
     977      const float  fOffset    = func.GetParameter(4)/singles.GetIntegrationWindow();
    908978//        const float  fCrossOrd  = func.GetParameter(5);
    909979
    910         hAmplitude.Fill(fAmplitude);
    911         hGain.Fill(fGain);
    912         //hGainRMS.Fill(f2GainRMS);
    913         hCrosstlk.Fill(fCrosstlk*100);
    914         hOffset.Fill(fOffset);
    915         hFitProb.Fill(100*fit_prob);
    916         hGainRMS.Fill(100*f2GainRMS/fGain);
    917 
    918         hcam.SetBinContent(pixel+1, fGain);
    919         hcam2.SetBinContent(pixel+1, 100*f2GainRMS/fGain);
    920 
    921         usePixel[pixel] = ok;
    922 
    923         if (!ok)
    924         {
    925             cout << "Pixel " << setw(4) << pixel << ": failed!" << endl;
    926             continue;
    927         }
    928 
    929         //Fill sum spectrum
    930         for (int b=1; b<=hist->GetNbinsX(); b++)
     980      hAmplitude.Fill(fAmplitude);
     981      hGain.Fill(fGain);
     982      //hGainRMS.Fill(f2GainRMS);
     983      hCrosstlk.Fill(fCrosstlk*100);
     984      hOffset.Fill(fOffset);
     985      hFitProb.Fill(100*fit_prob);
     986      hGainRMS.Fill(100*f2GainRMS/fGain);
     987
     988      hcam.SetBinContent(pixel+1, fGain);
     989      hcam2.SetBinContent(pixel+1, 100*f2GainRMS/fGain);
     990
     991      usePixel[pixel] = ok;
     992
     993      if (!ok)
     994        {
     995          cout << "Pixel " << setw(4) << pixel << ": failed!" << endl;
     996          continue;
     997        }
     998
     999      //Fill sum spectrum
     1000      for (int b=1; b<=hist->GetNbinsX(); b++)
    9311001        {
    9321002//            hSum.Fill(hist->GetBinCenter(b), hist->GetBinContent(b));
    933             hSumScale.Fill((hist->GetBinCenter(b)-fOffset)/fGain, (hist->GetBinContent(b)));
    934         }
    935 
    936         //plott special pixel
    937         if (
    938                 count_ok==0 ||
    939                 suspicous[pixel]
    940                 )
    941         {
    942             TCanvas &c = d->AddTab(Form("Pix%d", pixel));
    943             c.cd();
    944             gPad->SetLogy();
    945             gPad->SetGridx();
    946             gPad->SetGridy();
    947 
    948             hist->SetStats(false);
    949             hist->SetYTitle("Counts [a.u.]");
    950             hist->SetName(Form("Pix%d", pixel));
    951             hist->DrawCopy("hist")->SetDirectory(0);
     1003          hSumScale.Fill((hist->GetBinCenter(b)-fOffset)/fGain, (hist->GetBinContent(b)));
     1004        }
     1005
     1006      //plott special pixel
     1007      if (
     1008        count_ok==0 ||
     1009        suspicous[pixel]
     1010      )
     1011        {
     1012          TCanvas &c = d->AddTab(Form("Pix%d", pixel));
     1013          c.cd();
     1014          gPad->SetLogy();
     1015          gPad->SetGridx();
     1016          gPad->SetGridy();
     1017
     1018          hist->SetStats(false);
     1019          hist->SetYTitle("Counts [a.u.]");
     1020          hist->SetName(Form("Pix%d", pixel));
     1021          hist->DrawCopy("hist")->SetDirectory(0);
    9521022//            hist->Draw();
    953             func.DrawCopy("SAME");
    954         }
    955 
    956         count_ok++;
    957 
    958         delete hist;
     1023          func.DrawCopy("SAME");
     1024        }
     1025
     1026      count_ok++;
     1027
     1028      delete hist;
    9591029//        if (suspicous[pixel] == true) usePixel[pixel]=0;
    9601030    }
    961     // End of Loop over Pixel
    962 
    963     //------------------fill histograms-----------------------
    964 
    965     hcam.SetUsed(usePixel);
    966     hcam2.SetUsed(usePixel);
    967 
    968     TCanvas &canv = d->AddTab("Gain");
    969     canv.cd();
    970     canv.Divide(2,2);
    971 
    972     canv.cd(1);
    973     hcam.SetNameTitle( "gain","Gain distribution over whole camera");
    974     hcam.DrawCopy();
    975 
    976     canv.cd(2);
    977     hcam2.SetNameTitle( "gainRMS","GainRMS distribution over whole camera");
    978     hcam2.DrawCopy();
    979 
    980     TCanvas &cam_canv = d->AddTab("Camera_Gain");
    981 
    982     //------------------ Draw gain corrected sum specetrum -------------------
    983     gROOT->SetSelectedPad(0);
    984     TCanvas &c12 = d->AddTab("GainHist");
    985     c12.cd();
    986     gPad->SetLogy();
    987     gPad->SetGridx();
    988     gPad->SetGridy();
    989 
    990     for(UInt_t bin = 1; bin < (UInt_t)hSum.GetNbinsX()+1; bin++)
    991     {
    992         hSumScale.SetBinError(bin, hSum.GetBinContent(bin)*0.1);
    993     }
    994     hSumScale.DrawCopy("hist");
    995 
    996 
    997     //-------------------- fit gain corrected sum spectrum -------------------
    998     const Int_t    maxbin2   = hSumScale.GetMaximumBin();
    999     const Double_t binwidth2 = hSumScale.GetBinWidth(maxbin2);
    1000     const Double_t ampl2     = hSumScale.GetBinContent(maxbin2);
    1001 
    1002     //Calculate full width half Maximum
    1003     double fwhmScaled = 0;
    1004     for (int i=1; i<maxbin; i++)
    1005         if (hSumScale.GetBinContent(maxbin2-i)+hSumScale.GetBinContent(maxbin2+i)<ampl2*2)
    1006         {
    1007             fwhmScaled = (2*i+1)*binwidth2;
    1008             break;
    1009         }
    1010     if (fwhmScaled==0)
    1011     {
    1012         cout << "Could not determine start value for sigma." << endl;
    1013     }
    1014 
    1015     //Set fit parameters
    1016     Double_t par2[6] =
    1017     {
    1018         ampl2, 0.9, fwhmScaled*1.4, Crosstlk, 0, 1.1
    1019     };
    1020 
    1021     funcScaled.SetParameters(par2);
     1031
     1032  // End of Loop over Pixel
     1033
     1034  //------------------fill histograms-----------------------
     1035
     1036  hcam.SetUsed(usePixel);
     1037  hcam2.SetUsed(usePixel);
     1038
     1039  TCanvas &canv = d->AddTab("Gain");
     1040  canv.cd();
     1041  canv.Divide(2,2);
     1042
     1043  canv.cd(1);
     1044  hcam.SetNameTitle( "gain","Gain distribution over whole camera");
     1045  hcam.DrawCopy();
     1046
     1047  canv.cd(2);
     1048  hcam2.SetNameTitle( "gainRMS","GainRMS distribution over whole camera");
     1049  hcam2.DrawCopy();
     1050
     1051  TCanvas &cam_canv = d->AddTab("Camera_Gain");
     1052
     1053  //------------------ Draw gain corrected sum specetrum -------------------
     1054  gROOT->SetSelectedPad(0);
     1055  TCanvas &c12 = d->AddTab("GainHist");
     1056  c12.cd();
     1057  gPad->SetLogy();
     1058  gPad->SetGridx();
     1059  gPad->SetGridy();
     1060
     1061  for (UInt_t bin = 1; bin < (UInt_t)hSum.GetNbinsX()+1; bin++)
     1062    {
     1063      hSumScale.SetBinError(bin, hSum.GetBinContent(bin)*0.1);
     1064    }
     1065
     1066  hSumScale.DrawCopy("hist");
     1067
     1068
     1069  //-------------------- fit gain corrected sum spectrum -------------------
     1070  const Int_t    maxbin2   = hSumScale.GetMaximumBin();
     1071  const Double_t binwidth2 = hSumScale.GetBinWidth(maxbin2);
     1072  const Double_t ampl2     = hSumScale.GetBinContent(maxbin2);
     1073
     1074  //Calculate full width half Maximum
     1075  double fwhmScaled = 0;
     1076
     1077  for (int i=1; i<maxbin; i++)
     1078    if (hSumScale.GetBinContent(maxbin2-i)+hSumScale.GetBinContent(maxbin2+i)<ampl2*2)
     1079      {
     1080        fwhmScaled = (2*i+1)*binwidth2;
     1081        break;
     1082      }
     1083
     1084  if (fwhmScaled==0)
     1085    {
     1086      cout << "Could not determine start value for sigma." << endl;
     1087    }
     1088
     1089  //Set fit parameters
     1090  Double_t par2[6] =
     1091  {
     1092    ampl2, 0.9, fwhmScaled*1.4, Crosstlk, 0, 1.1
     1093  };
     1094
     1095  funcScaled.SetParameters(par2);
     1096  funcScaledBL.SetParameters(par2);
    10221097//    funcScaled.FixParameter(1,0.9);
    1023     funcScaled.FixParameter(4,0);
    1024     funcScaled.FixParameter(5,Crosstlk);
    1025 
    1026     const double minScaled = par2[1]-fwhmScaled;
    1027     const double maxScaled = par2[1]*maxOrder;
    1028     cout << "par2[1]" << par2[1] << endl;
    1029     cout << "maxScaled\t"           << maxScaled
    1030          << " \t\t minScaled\t"     << minScaled << endl;
    1031 
    1032     funcScaled.SetRange(minScaled, maxScaled);
    1033     hSumScale.Fit(&funcScaled, "N0QS", "", minScaled, maxScaled);
    1034     funcScaled.DrawCopy("SAME");
    1035     //--------fit gausses to peaks of gain corrected sum specetrum -----------
    1036 
    1037     // create histo for height of Maximum amplitudes vs. pe
    1038     TH1D hMaxHeights("MaxHeights",      "peak heights",
    1039                       10,  hSumScale.GetBinCenter(maxbin-1)-0.5,   10+hSumScale.GetBinCenter(maxbin-1)-0.5);
    1040 
    1041     FitGaussiansToSpectrum
    1042             (
    1043                 maxOrder-2,
    1044                 hSumScale,
    1045                 hMaxHeights,
    1046                 maxbin2,
    1047                 fwhmScaled,
    1048                 1
    1049                 );
     1098  funcScaled.FixParameter(4,0);
     1099  funcScaledBL.FixParameter(4,0);
     1100  funcScaled.FixParameter(5,Crosstlk);
     1101  funcScaledBL.FixParameter(5,Crosstlk);
     1102
     1103  const double minScaled = par2[1]-fwhmScaled;
     1104  const double maxScaled = par2[1]*maxOrder;
     1105  cout << "par2[1]" << par2[1] << endl;
     1106  cout << "maxScaled\t"           << maxScaled
     1107       << " \t\t minScaled\t"     << minScaled << endl;
     1108
     1109  funcScaled.SetRange(minScaled, maxScaled);
     1110  funcScaledBL.SetRange(minScaled, maxScaled);
     1111  hSumScale.Fit(&funcScaled, "N0QS", "", minScaled, maxScaled);
     1112  hSumScale.Fit(&funcScaledBL, "WLN0QS", "", minScaled, maxScaled);
     1113  funcScaled.DrawCopy("SAME");
     1114  funcScaledBL.DrawCopy("SAME");
     1115  //--------fit gausses to peaks of gain corrected sum specetrum -----------
     1116
     1117  // create histo for height of Maximum amplitudes vs. pe
     1118  TH1D hMaxHeights("MaxHeights",      "peak heights",
     1119                   10,  hSumScale.GetBinCenter(maxbin-1)-0.5,   10+hSumScale.GetBinCenter(maxbin-1)-0.5);
     1120
     1121  FitGaussiansToSpectrum
     1122  (
     1123    maxOrder-2,
     1124    hSumScale,
     1125    hMaxHeights,
     1126    maxbin2,
     1127    fwhmScaled,
     1128    1
     1129  );
    10501130
    10511131
    … …  
    10531133//    TCanvas &c16 = d->AddTab("PeakHeights");
    10541134//    gPad->SetLogy();
    1055     hMaxHeights.SetLineColor(kRed);
     1135  hMaxHeights.SetLineColor(kRed);
    10561136//    hMaxHeights.DrawCopy("SAME");
    1057     TF1 fExpo( "fExpo","expo", 0,  maxOrder-2);
    1058     fExpo.SetLineColor(kRed);
    1059     hMaxHeights.Fit(&fExpo, "N0QS" );
     1137  TF1 fExpo( "fExpo","expo", 0,  maxOrder-2);
     1138  fExpo.SetLineColor(kRed);
     1139  hMaxHeights.Fit(&fExpo, "N0QS" );
    10601140//    fExpo.DrawCopy("SAME");
    10611141//    c12.cd();
    10621142//    fExpo.DrawCopy("SAME");
    10631143
    1064     TCanvas &c2 = d->AddTab("Result");
    1065     c2.Divide(3,2);
    1066     c2.cd(1);
    1067     gPad->SetLogy();
    1068     hAmplitude.DrawCopy();
    1069 
    1070     c2.cd(2);
    1071     gPad->SetLogy();
     1144  TCanvas &c2 = d->AddTab("Result");
     1145  c2.Divide(3,2);
     1146  c2.cd(1);
     1147  gPad->SetLogy();
     1148  hAmplitude.DrawCopy();
     1149
     1150  c2.cd(2);
     1151  gPad->SetLogy();
    10721152//    hGain.DrawCopy();
    10731153
    10741154
    1075     TF1 GainGaus( "GainGaus", "gaus", 0, hGain.GetXaxis()->GetXmax());
    1076     hGain.Rebin(2);
    1077     hGain.Fit(&GainGaus, "N0QS");
     1155  TF1 GainGaus( "GainGaus", "gaus", 0, hGain.GetXaxis()->GetXmax());
     1156  hGain.Rebin(2);
     1157  hGain.Fit(&GainGaus, "N0QS");
    10781158
    10791159//    float gain_mean     = GainGaus.GetParameter(1);
    1080     float gain_median   = MedianOfH1(&hGain);
    1081     TH1F hNormGain          ("NormGain",      "median normalized Gain distribution",
    1082                             hGain.GetNbinsX(),
    1083                             (hGain.GetXaxis()->GetXmin())/gain_median,
    1084                             (hGain.GetXaxis()->GetXmax())/gain_median
    1085                             );
    1086     hNormGain.SetXTitle("gain [fraction of median gain value]");
    1087     hNormGain.SetYTitle("Counts");
    1088 
    1089 
    1090     hNormGain.Add(&hGain);
     1160  float gain_median   = MedianOfH1(&hGain);
     1161  TH1F hNormGain          ("NormGain",      "median normalized Gain distribution",
     1162                           hGain.GetNbinsX(),
     1163                           (hGain.GetXaxis()->GetXmin())/gain_median,
     1164                           (hGain.GetXaxis()->GetXmax())/gain_median
     1165                          );
     1166  hNormGain.SetXTitle("gain [fraction of median gain value]");
     1167  hNormGain.SetYTitle("Counts");
     1168
     1169
     1170  hNormGain.Add(&hGain);
    10911171//    hNormGain.SetStats(false);
    1092     hNormGain.GetXaxis()->SetRangeUser(
    1093                             1-(hNormGain.GetXaxis()->GetXmax()-1),
    1094                             hNormGain.GetXaxis()->GetXmax()
    1095                             );
    1096     hNormGain.DrawCopy();
    1097     hNormGain.Fit(&GainGaus, "N0QS");
    1098     GainGaus.DrawCopy("SAME");
    1099 
    1100     //plot normailzed camera view
    1101     cam_canv.cd();
    1102     MHCamera hNormCam(fact);
    1103     hNormCam.SetStats(false);
    1104     for (UInt_t pixel =1; pixel <= 1440; pixel++)
    1105     {
    1106         hNormCam.SetBinContent(pixel, hcam.GetBinContent(pixel)/gain_median);
    1107     }
    1108 
    1109     hNormCam.SetNameTitle( "gain","normalized [Median] Gain distribution over whole camera");
    1110     hNormCam.SetZTitle("Horst");
    1111     hNormCam.SetUsed(usePixel);
    1112     hNormCam.DrawCopy();
     1172  hNormGain.GetXaxis()->SetRangeUser(
     1173    1-(hNormGain.GetXaxis()->GetXmax()-1),
     1174    hNormGain.GetXaxis()->GetXmax()
     1175  );
     1176  hNormGain.DrawCopy();
     1177  hNormGain.Fit(&GainGaus, "N0QS");
     1178  GainGaus.DrawCopy("SAME");
     1179
     1180  //plot normailzed camera view
     1181  cam_canv.cd();
     1182  MHCamera hNormCam(fact);
     1183  hNormCam.SetStats(false);
     1184
     1185  for (UInt_t pixel =1; pixel <= 1440; pixel++)
     1186    {
     1187      hNormCam.SetBinContent(pixel, hcam.GetBinContent(pixel)/gain_median);
     1188    }
     1189
     1190  hNormCam.SetNameTitle( "gain","normalized [Median] Gain distribution over whole camera");
     1191  hNormCam.SetZTitle("Horst");
     1192  hNormCam.SetUsed(usePixel);
     1193  hNormCam.DrawCopy();
    11131194//    hGain.Scale(1/GainGaus.GetParameter(1));
    11141195//    GainGaus.DrawCopy("SAME");
    11151196
    1116     c2.cd(3);
    1117     gPad->SetLogy();
    1118     hGainRMS.DrawCopy();
    1119 
    1120     c2.cd(4);
    1121     gPad->SetLogy();
    1122     hCrosstlk.DrawCopy();
    1123 
    1124     c2.cd(5);
    1125     gPad->SetLogy();
    1126     hOffset.DrawCopy();
    1127 
    1128     c2.cd(6);
    1129     gPad->SetLogy();
    1130     hFitProb.DrawCopy();
    1131 
    1132 
    1133     TCanvas &c25 = d->AddTab("Spectra");
    1134     c25.Divide(2,1);
    1135     c25.cd(1);
    1136     gPad->SetLogy();
    1137     gPad->SetGrid();
    1138     hSum.DrawCopy("hist");
    1139     funcSum.DrawCopy("SAME");
    1140 
    1141     c25.cd(2);
    1142     gPad->SetLogy();
    1143     gPad->SetGrid();
    1144     hSumScale.DrawCopy("hist");
    1145     funcScaled.DrawCopy("SAME");
    1146 
    1147 
    1148     /*
    1149     TCanvas &c3 = d->AddTab("Separation");
    1150     gPad->SetLogy();
    1151     hSep.DrawCopy();
    1152     */
    1153 
    1154 //    d->SaveAs(outfile);
    1155     return 0;
     1197  c2.cd(3);
     1198  gPad->SetLogy();
     1199  hGainRMS.DrawCopy();
     1200
     1201  c2.cd(4);
     1202  gPad->SetLogy();
     1203  hCrosstlk.DrawCopy();
     1204
     1205  c2.cd(5);
     1206  gPad->SetLogy();
     1207  hOffset.DrawCopy();
     1208
     1209  c2.cd(6);
     1210  gPad->SetLogy();
     1211//    hFitProb.DrawCopy();
     1212  hGain.DrawCopy();
     1213
     1214
     1215  TCanvas &c25 = d->AddTab("Spectra");
     1216  c25.Divide(2,1);
     1217  c25.cd(1);
     1218  gPad->SetLogy();
     1219  gPad->SetGrid();
     1220  hSum.DrawCopy("hist");
     1221  funcSum.DrawCopy("SAME");
     1222
     1223  c25.cd(2);
     1224  gPad->SetLogy();
     1225  gPad->SetGrid();
     1226  hSumScale.DrawCopy("hist");
     1227  funcScaled.DrawCopy("SAME");
     1228  funcScaledBL.DrawCopy("SAME");
     1229
     1230
     1231  /*
     1232  TCanvas &c3 = d->AddTab("Separation");
     1233  gPad->SetLogy();
     1234  hSep.DrawCopy();
     1235  */
     1236
     1237  d->SaveAs("/home_nfs/isdc/jbbuss/singlepe.root");
     1238  return 0;
    11561239}
    11571240
    11581241Double_t fcn(Double_t *xx, Double_t *par)
    11591242{
    1160     Double_t x     = xx[0];
    1161 
    1162     Double_t ampl   = par[0];
    1163     Double_t gain   = par[1];
    1164     Double_t sigma  = par[2];
    1165     Double_t cross  = par[3];
    1166     Double_t shift  = par[4];
     1243  Double_t x     = xx[0];
     1244
     1245  Double_t ampl   = par[0];
     1246  Double_t gain   = par[1];
     1247  Double_t sigma  = par[2];
     1248  Double_t cross  = par[3];
     1249  Double_t shift  = par[4];
    11671250//    Double_t k      = par[5];
    11681251
    1169     Double_t xTalk = 1;
    1170 
    1171     Double_t y = 0;
    1172     for (int order = 1; order < 14; order++)
    1173     {
    1174         Double_t arg   = (x - order*gain - shift)/sigma;
    1175 
    1176         Double_t gauss = TMath::Exp(-0.5 * arg * arg/order);
    1177 
    1178         y += xTalk*gauss;
     1252  Double_t xTalk = 1;
     1253
     1254  Double_t y = 0;
     1255
     1256  for (int order = 1; order < 14; order++)
     1257    {
     1258      Double_t arg   = (x - order*gain - shift)/sigma;
     1259
     1260      Double_t gauss = TMath::Exp(-0.5 * arg * arg/order);
     1261
     1262      y += xTalk*gauss;
    11791263
    11801264//        for (int j = 1; j < k; j++)
    11811265//        {
    1182             xTalk *= cross;
     1266      xTalk *= cross;
    11831267//        }
    11841268    }
    11851269
    1186     return ampl*y;
     1270  return ampl*y;
    11871271}
    11881272
    11891273Double_t fcn_cross(Double_t *xx, Double_t *par)
    11901274{
    1191     Double_t x     = xx[0];
    1192 
    1193     Double_t ampl   = par[0];
    1194     Double_t gain   = par[1];
    1195     Double_t sigma  = par[2];
    1196     Double_t cross  = par[3];
    1197     Double_t shift  = par[4];
     1275  Double_t x     = xx[0];
     1276
     1277  Double_t ampl   = par[0];
     1278  Double_t gain   = par[1];
     1279  Double_t sigma  = par[2];
     1280  Double_t cross  = par[3];
     1281  Double_t shift  = par[4];
    11981282//    Double_t powOfX = par[5];
    11991283
    1200     Double_t xTalk = 1;
    1201 
    1202     Double_t y = 0;
    1203     for (int order = 1; order < 14; order++)
    1204     {
    1205         Double_t arg   = (x - order*gain - shift)/sigma;
    1206 
    1207         Double_t gauss = TMath::Exp(-0.5 * arg * arg/order);
    1208 
    1209         y += xTalk*gauss;
     1284  Double_t xTalk = 1;
     1285
     1286  Double_t y = 0;
     1287
     1288  for (int order = 1; order < 14; order++)
     1289    {
     1290      Double_t arg   = (x - order*gain - shift)/sigma;
     1291
     1292      Double_t gauss = TMath::Exp(-0.5 * arg * arg/order);
     1293
     1294      y += xTalk*gauss;
    12101295
    12111296//        xTalk = pow(cross, order*powOfX);
    … …  
    12151300//        }
    12161301//        cross *= TMath::Exp(-powOfX*cross);
    1217         xTalk *= cross;
    1218     }
    1219 
    1220     return ampl*y;
     1302      xTalk *= cross;
     1303    }
     1304
     1305  return ampl*y;
    12211306}
    12221307
    … …  
    12241309void
    12251310FitGaussiansToSpectrum(
    1226         UInt_t      maxOrder,
    1227         TH1        &hSource,
    1228         TH1        &hDest,
    1229         Int_t       maxbin,
    1230         double      fwhm,
    1231         Double_t    gain
    1232         )
     1311  UInt_t      maxOrder,
     1312  TH1        &hSource,
     1313  TH1        &hDest,
     1314  Int_t       maxbin,
     1315  double      fwhm,
     1316  Double_t    gain
     1317)
    12331318{
    12341319
    1235     Double_t    peakposition = hSource.GetBinCenter(maxbin);
    1236     char        fname[20];
    1237 
    1238     // fit gauss functions to spectrum
    1239     for (UInt_t nr = 1; nr<maxOrder; nr++)
    1240     {
    1241         sprintf(fname,"gaus%d",nr);
    1242 
    1243         //create gauss functions
    1244         TF1 gaus( fname,"gaus", peakposition-fwhm,     peakposition+fwhm);
    1245         gaus.SetParLimits(1, peakposition-fwhm, peakposition+fwhm);
    1246         gaus.SetParLimits(2, -fwhm, fwhm );
    1247         //fit and draw gauss functions
    1248         hSource.Fit(  &gaus,   "N0QS", "", peakposition-fwhm, peakposition+fwhm);
     1320  Double_t    peakposition = hSource.GetBinCenter(maxbin);
     1321  char        fname[20];
     1322
     1323  // fit gauss functions to spectrum
     1324  for (UInt_t nr = 1; nr<maxOrder; nr++)
     1325    {
     1326      sprintf(fname,"gaus%d",nr);
     1327
     1328      //create gauss functions
     1329      TF1 gaus( fname,"gaus", peakposition-fwhm,     peakposition+fwhm);
     1330      gaus.SetParLimits(1, peakposition-fwhm, peakposition+fwhm);
     1331      gaus.SetParLimits(2, -fwhm, fwhm );
     1332      //fit and draw gauss functions
     1333      hSource.Fit(  &gaus,   "N0QS", "", peakposition-fwhm, peakposition+fwhm);
    12491334//        gaus.DrawCopy("SAME");
    12501335
    1251         peakposition = gaus.GetParameter(1);
    1252 
    1253         //fill histo for height of Maximum amplitudes vs. pe
    1254         hDest.SetBinContent(nr, gaus.GetParameter(0));
    1255 
    1256         peakposition += gain;
    1257     }
    1258     return;
     1336      peakposition = gaus.GetParameter(1);
     1337
     1338      //fill histo for height of Maximum amplitudes vs. pe
     1339      hDest.SetBinContent(nr, gaus.GetParameter(0));
     1340
     1341      peakposition += gain;
     1342    }
     1343
     1344  return;
    12591345}
    12601346
    12611347Double_t MedianOfH1 (
    1262         TH1*            inputHisto
    1263         )
     1348  TH1*            inputHisto
     1349)
    12641350{
    1265    //compute the median for 1-d histogram h1
    1266    Int_t nbins = inputHisto->GetXaxis()->GetNbins();
    1267    Double_t *x = new Double_t[nbins];
    1268    Double_t *y = new Double_t[nbins];
    1269    for (Int_t i=0;i<nbins;i++) {
     1351  //compute the median for 1-d histogram h1
     1352  Int_t nbins = inputHisto->GetXaxis()->GetNbins();
     1353  Double_t *x = new Double_t[nbins];
     1354  Double_t *y = new Double_t[nbins];
     1355
     1356  for (Int_t i=0; i<nbins; i++)
     1357    {
    12701358      x[i] = inputHisto->GetXaxis()->GetBinCenter(i+1);
    12711359      y[i] = inputHisto->GetBinContent(i+1);
    1272    }
    1273    Double_t median = TMath::Median(nbins,x,y);
    1274    delete [] x;
    1275    delete [] y;
    1276    return median;
     1360    }
     1361
     1362  Double_t median = TMath::Median(nbins,x,y);
     1363  delete [] x;
     1364  delete [] y;
     1365  return median;
    12771366}
    12781367// end of PlotMedianEachSliceOfPulse
Note: See TracChangeset for help on using the changeset viewer.