FMD/scripts/FindCommonModeNoise.C

   1 //____________________________________________________________________
   2 //
   3 // $Id: DrawDigits.C 30718 2009-01-22 16:07:40Z cholm $
   4 //
   5 // Script that contains a class to draw eloss from hits, versus ADC
   6 // counts from digits, using the AliFMDInputHits class in the util library.
   7 //
   8 // It draws the energy loss versus the p/(mq^2).  It can be overlayed
   9 // with the Bethe-Bloc curve to show how the simulation behaves
  10 // relative to the expected.
  11 //
  12 // Use the script `Compile.C' to compile this class using ACLic.
  13 //
  14 #include <TH1F.h>
  15 #include <TMath.h>
  16 #include <TStyle.h>
  17 #include <TArrayF.h>
  18 #include <TSystem.h>
  19 #include <TFile.h>
  20 #include <TCanvas.h>
  21 #include <TSystem.h>
  22 #include <TRandom.h>
  23 #include <TPad.h>
  24 #include <TLegend.h>
  25
  26 #include <AliCDBManager.h>
  27 #include <AliFMDDigit.h>
  28 #include <AliFMDInput.h>
  29 #include <AliFMDParameters.h>
  30 #include <AliRawReader.h>
  31 #include <AliLog.h>
  32
  33 #include <iostream>
  34
  35 /** @class DrawDigits
  36     @brief Draw pedestal, noise, and both common mode noise corrected
  37     @code
  38     Root> .L Compile.C
  39     Root> Compile("FindCommonModeNoise.C")
  40     Root> FindCommonModeNoise fcmn("file");
  41     Root> fcmn.Run();
  42     @endcode
  43
  44     If null is passed instead of a file, then a simulation of common
  45     mode noise is run instead.
  46
  47     @ingroup FMD_script
  48  */
  49 class FindCommonModeNoise : public AliFMDInput
  50 {
  51 private:
  52   TString fCalibDir;
  53   Int_t   fCount;
  54   TFile*  fOut;
  55   Float_t fShift;
  56   Float_t fCmn;                // Common mode noise component
  57   TH1F*   fPed;                //
  58   TH1F*   fNoise;              //
  59   TH1F*   fCenteredPed;        //
  60   TH1F*   fCorrectedNoise;     //
  61   TH1F*   fSummedAdc;          //
  62   TH1F*   fSummedAdc2;         //
  63   TH1F*   fSummedCenteredAdc;  //
  64   TH1F*   fSummedCenteredAdc2; //
  65   TH1F*   fAdc;                //
  66
  67 public:
  68   //__________________________________________________________________
  69   FindCommonModeNoise(const char*  file,
  70                       const char*  calibDir=0)
  71     : AliFMDInput(),
  72       fCalibDir(""),
  73       fOut(0)
  74   {
  75     fCmn = 2;
  76     if (calibDir) fCalibDir = calibDir;
  77     if (file) {
  78       AddLoad(kRaw);
  79       SetRawFile(file);
  80     }
  81     else
  82       AddLoad(kUser);
  83
  84     TString outName(Form("histo_%s", gSystem->BaseName(fRawFile.Data())));
  85     if (!outName.EndsWith(".root")) outName.Append(".root");
  86     fOut  = TFile::Open(outName.Data(),"RECREATE");
  87
  88     fPed                = MakeHist("ped_dist","Pedestal Distribution");
  89     fNoise              = MakeHist("noise_dist",
  90                                    "Noise Distribution (Not CMN corrected)");
  91     fCenteredPed        = MakeHist("cmn_ped_dist",
  92                                    "Pedestal - Average Pedestal Distribution");
  93     fCorrectedNoise     = MakeHist("noise_corrected_dist",
  94                                    "Noise Distribution (CMN corrected)");
  95     fSummedAdc          = MakeHist("sum_adc_dist","Sum ADC Distribution");
  96     fSummedAdc2         = MakeHist("sum_adc2_dist",
  97                                    "Sum ADC Squared Distribution");
  98     fSummedCenteredAdc  = MakeHist("sum_cmn_adc_dist",
  99                                    "Sum CMN Adjusted ADC Distribution");
 100     fSummedCenteredAdc2 = MakeHist("sum_cmn_adc2_dist",
 101                                    "Sum CMN Adjusted ADC Squared Distribution");
 102     fAdc                = MakeHist("adc_dist","Current ADC Distribution");
 103   }
 104   //__________________________________________________________________
 105   Int_t NEvents() const
 106   {
 107     Int_t nEv = (fReader ? fReader->GetNumberOfEvents() : -1);
 108     if (nEv > 0) return TMath::Min(nEv, 1000);
 109     return 1000;
 110   }
 111   //__________________________________________________________________
 112   TH1F* MakeHist(const char* name, const char* title)
 113   {
 114     TH1F* h = new TH1F(name,title,51200,-0.5,51199.5);
 115     h->SetXTitle("Strip number");
 116     h->SetFillColor(kRed);
 117     h->SetFillStyle(3001);
 118     h->SetDirectory(0);
 119     h->SetStats(kFALSE);
 120     return h;
 121   }
 122   //__________________________________________________________________
 123   Bool_t Init()
 124   {
 125     AliCDBManager* cdb = AliCDBManager::Instance();
 126     cdb->SetRun(0);
 127     cdb->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
 128
 129     AliFMDParameters* pars = AliFMDParameters::Instance();
 130     if (fCalibDir.IsNull())
 131       pars->Init(kFALSE);
 132     else
 133       pars->Init(fCalibDir.Data(), kFALSE);
 134
 135     fCount = 0;
 136
 137     return AliFMDInput::Init();
 138   }
 139
 140   //__________________________________________________________________
 141   Bool_t Begin(Int_t e)
 142   {
 143     fCount++;
 144     fAdc->Reset();
 145     // Common mode noise generation
 146     fShift = gRandom->Gaus(0, fCmn);
 147     return AliFMDInput::Begin(e);
 148   }
 149
 150   //__________________________________________________________________
 151   Int_t BinCenter(UShort_t d, Char_t r, UShort_t s, UShort_t t) const
 152   {
 153     if (d == 1 && !(r == 'I' || r == 'i')) return -1;
 154
 155     Int_t idx = -1;
 156     switch (d) {
 157     case 1: idx = 0; break;
 158     case 2: idx = 10240; break;
 159     case 3: idx = 3 * 10240; break;
 160     default:
 161       return -1;
 162     }
 163
 164     UShort_t q    = (r == 'I' || r == 'i' ? 0 : 1);
 165     // UShort_t nSec = (q == 0 ?  20 :  40);
 166     UShort_t nStr = (q == 0 ? 512 : 256);
 167
 168     idx += (q == 0 ? 0 : 10240);
 169     idx += s * nStr + t;
 170
 171     if (idx >= 51200) return -1;
 172     return idx;
 173   }
 174
 175   //__________________________________________________________________
 176   Int_t BinCenter(UShort_t chip, UShort_t strip) const
 177   {
 178     return chip * 128 + strip;
 179   }
 180   //__________________________________________________________________
 181   Int_t BinNumber(UShort_t chip, UShort_t strip) const
 182   {
 183     return 1 + BinCenter(chip, strip);
 184   }
 185   //__________________________________________________________________
 186   Float_t GetSignal(UShort_t d, Char_t, UShort_t s, UShort_t t)
 187   {
 188     // Calculate signal value.  Only the first 256 strips of FMD are
 189     // filled.
 190     if (d > 1 || s > 0 || t > 256) return 0;
 191     Float_t c = 0.005;
 192     Float_t x = (t % 128);
 193     Float_t p = c * (64 * 64 - 128 * x + x * x) + 100;
 194     Float_t n = (p-100)/40 + 1;
 195
 196     // In case we simulate a common mode noise component.
 197     if (t < 128)
 198       return fShift + gRandom->Gaus(p, n);
 199
 200     // In case all noise is uncorrelated.
 201     return gRandom->Gaus(p, TMath::Sqrt(n*n + fCmn*fCmn));
 202   }
 203   //__________________________________________________________________
 204   Bool_t ProcessRawDigit(AliFMDDigit* digit)
 205   {
 206     // From data
 207     if (!digit) return kTRUE;
 208
 209     UShort_t d             =  digit->Detector();
 210     Char_t   r             =  digit->Ring();
 211     UShort_t s             =  digit->Sector();
 212     UShort_t t             =  digit->Strip();
 213     return ProcessOne(d, r, s, t, digit->Counts());
 214   }
 215   //__________________________________________________________________
 216   Bool_t ProcessUser(UShort_t d, Char_t r, UShort_t s, UShort_t t, Float_t v)
 217   {
 218     // From our simulator in GetValue
 219     return ProcessOne(d, r, s, t, v);
 220   }
 221   //__________________________________________________________________
 222   Bool_t ProcessOne(UShort_t d, Char_t r, UShort_t s, UShort_t t, Float_t v)
 223   {
 224     Int_t    i             =  BinCenter(d, r, s, t);
 225     if (i < 0) {
 226       std::cout << "Invalid index " << i << " returned for FMD"
 227                 << d << r << '[' << s << ',' << t << ']' << std::endl;
 228       return kFALSE;
 229     }
 230
 231     fAdc->Fill(i, v);
 232
 233     return kTRUE;
 234   }
 235   //__________________________________________________________________
 236   Bool_t End()
 237   {
 238     // At the end of each event, fill the summed histograms.
 239     for (UShort_t chip = 0; chip < 400; chip++) {
 240       Double_t chipAv = 0;
 241       for (UShort_t strip = 0; strip < 128; strip++)
 242         chipAv += fAdc->GetBinContent(1+BinCenter(chip, strip));
 243       chipAv /= 128;
 244
 245       for (UShort_t strip = 0; strip < 128; strip++) {
 246         Int_t    i    = BinCenter(chip, strip);
 247         Double_t adc  = fAdc->GetBinContent(i+1);
 248         Double_t cAdc = chipAv - adc;
 249         fSummedAdc->Fill(i, adc);
 250         fSummedAdc2->Fill(i, adc * adc);
 251         fSummedCenteredAdc->Fill(i, cAdc);
 252         fSummedCenteredAdc2->Fill(i, cAdc*cAdc);
 253       }
 254     }
 255
 256     return AliFMDInput::End();
 257   }
 258   //__________________________________________________________________
 259   Bool_t Finish()
 260   {
 261     if (fEventCount == 0) {
 262       std::cerr << "No events!" << std::endl;
 263       return kFALSE;
 264     }
 265
 266     for (UInt_t bin = 1; bin <= 51200; bin++) {
 267       if (bin % (51200 / 10) == 0)
 268         std::cout << "Bin # " << bin << std::endl;
 269       Double_t sumAdc  = fSummedAdc->GetBinContent(bin);
 270       Double_t sumAdc2 = fSummedAdc2->GetBinContent(bin);
 271       Double_t avAdc   = sumAdc / fEventCount;
 272       Double_t avAdc2  = sumAdc2 / fEventCount;
 273       Double_t noise2  = avAdc2 - avAdc * avAdc;
 274
 275       fPed->SetBinContent(bin, avAdc);
 276       fNoise->SetBinContent(bin, noise2 >= 0 ? TMath::Sqrt(noise2) : 0);
 277
 278       Double_t sumCAdc  = fSummedCenteredAdc->GetBinContent(bin);
 279       Double_t sumCAdc2 = fSummedCenteredAdc2->GetBinContent(bin);
 280       Double_t avCAdc   = sumCAdc / fEventCount;
 281       Double_t avCAdc2  = sumCAdc2 / fEventCount;
 282       Double_t cNoise2  = avCAdc2 - avCAdc * avCAdc;
 283       fCenteredPed->SetBinContent(bin, avCAdc);
 284       fCorrectedNoise->SetBinContent(bin, TMath::Sqrt(cNoise2));
 285     }
 286
 287     gStyle->SetPalette(1);
 288     gStyle->SetOptTitle(0);
 289     gStyle->SetCanvasColor(0);
 290     gStyle->SetCanvasBorderSize(0);
 291     gStyle->SetPadColor(0);
 292     gStyle->SetPadBorderSize(0);
 293
 294     TCanvas* c = new TCanvas("c", "C", 800, 600);
 295     c->SetFillColor(0);
 296     c->SetBorderMode(0);
 297     c->SetBorderSize(0);
 298     c->SetTopMargin(0);
 299     c->SetBottomMargin(0);
 300     c->Divide(2, 1);
 301
 302     TPad* p1 = static_cast<TPad*>(c->cd(1));
 303     p1->SetTopMargin(0.05);
 304     p1->SetRightMargin(0.05);
 305     if (!fReader) fNoise->GetXaxis()->SetRangeUser(-0.5, 255.5);
 306     fNoise->SetMinimum(0);
 307     fNoise->SetMaximum(1.5*fNoise->GetMaximum());
 308     fNoise->DrawCopy();
 309     fCorrectedNoise->SetMinimum(0);
 310     fCorrectedNoise->SetFillColor(kBlue);
 311     fCorrectedNoise->DrawCopy("same");
 312     TLegend* l1 = new TLegend(0.2, 0.7, 0.945, 0.945);
 313     l1->SetFillColor(0);
 314     l1->SetFillStyle(0);
 315     l1->SetBorderSize(0);
 316     l1->AddEntry(fNoise, "Noise", "f");
 317     l1->AddEntry(fCorrectedNoise, "Corrected noise", "f");
 318     l1->Draw();
 319
 320
 321     TPad* p2 = static_cast<TPad*>(c->cd(2));
 322     p2->SetTopMargin(0.05);
 323     p2->SetRightMargin(0.05);
 324     if (!fReader) fPed->GetXaxis()->SetRangeUser(-0.5, 255.5);
 325     fPed->SetMinimum(fCenteredPed->GetMinimum());
 326     fPed->SetMaximum(1.5*fPed->GetMaximum());
 327     fPed->DrawCopy();
 328     fCenteredPed->SetFillColor(kBlue);
 329     fCenteredPed->DrawCopy("same");
 330     TLegend* l2 = new TLegend(0.2, 0.7, 0.945, 0.945);
 331     l2->SetFillColor(0);
 332     l2->SetFillStyle(0);
 333     l2->SetBorderSize(0);
 334     l2->AddEntry(fPed, "Pedestal", "f");
 335     l2->AddEntry(fCenteredPed, "Pedestal - #bar{Pedestal}", "f");
 336     l2->Draw();
 337
 338     c->Modified();
 339     c->Update();
 340     c->cd();
 341     c->Print("cmn.png");
 342
 343     if (fOut) {
 344       std::cout << "Flusing to disk ... " << std::flush;
 345       fOut->cd();
 346       fPed->Write();                //
 347       fNoise->Write();              //
 348       fCenteredPed->Write();        //
 349       fCorrectedNoise->Write();     //
 350       fSummedAdc->Write();          //
 351       fSummedAdc2->Write();         //
 352       fSummedCenteredAdc->Write();  //
 353       fSummedCenteredAdc2->Write(); //
 354       fAdc->Write();                //
 355       fOut->Write();
 356       fOut->Close();
 357       std::cout << "done" << std::endl;
 358     }
 359     return kTRUE;
 360   }
 361
 362   ClassDef(FindCommonModeNoise,0);
 363 };
 364
 365 //____________________________________________________________________
 366 //
 367 // EOF
 368 //