New script to do common-mode-noise analysis of raw data.

The script can also simulate the effect of common mode noise.

This required a few changes to AliFMDInput.

diff --git a/FMD/AliFMDInput.cxx b/FMD/AliFMDInput.cxx
index 20b759b61c9095a826e1895f31e9a3d7e5622597..d1c8daa5257d469d42401b50305bf02183fec6f0 100644 (file)
--- a/FMD/AliFMDInput.cxx
+++ b/FMD/AliFMDInput.cxx
@@ -442,11 +442,12 @@ AliFMDInput::Begin(Int_t event)
 
   // Possibly load FMD Digit information 
   if (TESTBIT(fTreeMask, kRaw) || TESTBIT(fTreeMask, kRawCalib)) {
 
   // Possibly load FMD Digit information 
   if (TESTBIT(fTreeMask, kRaw) || TESTBIT(fTreeMask, kRawCalib)) {

+    Bool_t mon = fRawFile.Contains("mem://");

     // AliInfo("Getting FMD raw data digits");
     // AliInfo("Getting FMD raw data digits");

-    std::cout << "Waiting for event ..." << std::endl;
+    if (mon) std::cout << "Waiting for event ..." << std::flush;

     do { 
       if (!fReader->NextEvent()) { 
     do { 
       if (!fReader->NextEvent()) { 

-       if (fRawFile.Contains("mem://")) { 
+       if (mon) { 

          gSystem->Sleep(3);
          continue;
        }
          gSystem->Sleep(3);
          continue;
        }


@@ -457,6 +458,7 @@ AliFMDInput::Begin(Int_t event)
         eventType == AliRawEventHeaderBase::kCalibrationEvent) 
        break;
     } while (true);
         eventType == AliRawEventHeaderBase::kCalibrationEvent) 
        break;
     } while (true);

+    if (mon) std::cout << "got it" << std::endl;

     // AliFMDRawReader r(fReader, 0);
     fArrayA->Clear();
     fFMDReader->ReadAdcs(fArrayA);
     // AliFMDRawReader r(fReader, 0);
     fArrayA->Clear();
     fFMDReader->ReadAdcs(fArrayA);


@@ -497,6 +499,8 @@ AliFMDInput::Event()
     if (!ProcessRecPoints()) return kFALSE;
   if (TESTBIT(fTreeMask, kESD))
     if (!ProcessESDs()) return kFALSE;
     if (!ProcessRecPoints()) return kFALSE;
   if (TESTBIT(fTreeMask, kESD))
     if (!ProcessESDs()) return kFALSE;

+  if (TESTBIT(fTreeMask, kUser))
+    if (!ProcessUsers()) return kFALSE;

   
   return kTRUE;
 }
   
   return kTRUE;
 }


@@ -782,6 +786,27 @@ AliFMDInput::ProcessESDs()
   return kTRUE;
 }
 
   return kTRUE;
 }
 

+//____________________________________________________________________
+Bool_t 
+AliFMDInput::ProcessUsers()
+{
+  // Process event summary data
+  for (UShort_t det = 1; det <= 3; det++) {
+    Char_t rings[] = { 'I', (det == 1 ? '\0' : 'O'), '\0' };
+    for (Char_t* rng = rings; *rng != '\0'; rng++) {
+      UShort_t nsec = (*rng == 'I' ?  20 :  40);
+      UShort_t nstr = (*rng == 'I' ? 512 : 256);
+      for (UShort_t sec = 0; sec < nsec; sec++) {
+       for (UShort_t str = 0; str < nstr; str++) {
+         Float_t v  = GetSignal(det,*rng,sec,str);
+         if (!ProcessUser(det, *rng, sec, str, v)) continue;
+       }
+      }
+    }
+  }
+  return kTRUE;
+}
+

 //____________________________________________________________________
 Bool_t
 AliFMDInput::End()
 //____________________________________________________________________
 Bool_t
 AliFMDInput::End()

diff --git a/FMD/AliFMDInput.h b/FMD/AliFMDInput.h
index 2a99a1415ccd8abddc812e1b308ecc3a8267feb0..8a931bcdd5091c5b326b4aff1aec6266156975c1 100644 (file)
--- a/FMD/AliFMDInput.h
+++ b/FMD/AliFMDInput.h
@@ -116,7 +116,8 @@ public:
     kGeometry,        // Not really a tree 
     kTracks,         // Hits and tracs - for BG study  
     kTrackRefs,       // Track references - also for BG study
     kGeometry,        // Not really a tree 
     kTracks,         // Hits and tracs - for BG study  
     kTrackRefs,       // Track references - also for BG study

-    kRawCalib         // Read raws and calibrate them
+    kRawCalib,        // Read raws and calibrate them
+    kUser

   };
   /** CTOR  */
   AliFMDInput();
   };
   /** CTOR  */
   AliFMDInput();


@@ -194,6 +195,9 @@ public:
   /** Loop over all ESD data, and call ProcessESD for each entry.
       @return  @c false on error  */
   virtual Bool_t ProcessESDs();
   /** Loop over all ESD data, and call ProcessESD for each entry.
       @return  @c false on error  */
   virtual Bool_t ProcessESDs();

+  /** Loop over all strips and ask user routine to supply the data.
+      @return  @c false on error  */
+  virtual Bool_t ProcessUsers();

 
   /** Process one hit, and optionally it's corresponding kinematics
       track.  Users should over this to process each hit. 
 
   /** Process one hit, and optionally it's corresponding kinematics
       track.  Users should over this to process each hit. 


@@ -248,8 +252,18 @@ public:
       @param eta  Psuedo-rapidity 
       @param mult Psuedo-multiplicity 
       @return  @c false on error  */
       @param eta  Psuedo-rapidity 
       @param mult Psuedo-multiplicity 
       @return  @c false on error  */

-  virtual Bool_t ProcessESD(UShort_t, Char_t, UShort_t, UShort_t, 
-                           Float_t, Float_t);
+  virtual Bool_t ProcessESD(UShort_t d, Char_t r, UShort_t s, UShort_t t, 
+                           Float_t eta, Float_t mult);
+  /** Process User data for the FMD.  Users should overload this to
+      deal with ESD data. 
+      @param d    Detector number (1-3)
+      @param r    Ring identifier ('I' or 'O')
+      @param s    Sector number (0-19, or 0-39)
+      @param t    Strip number (0-511, or 0-255)
+      @param v    Value
+      @return  @c false on error  */
+  virtual Bool_t ProcessUser(UShort_t d, Char_t r, UShort_t s, UShort_t t, 
+                            Float_t v);

   /** Service function to make a logarithmic axis. 
       @param n    Number of bins 
       @param min  Minimum of axis 
   /** Service function to make a logarithmic axis. 
       @param n    Number of bins 
       @param min  Minimum of axis 


@@ -301,6 +315,17 @@ protected:
   /** Assignement operator 
       @return  REference to this */
   AliFMDInput& operator=(const AliFMDInput&) { return *this; }
   /** Assignement operator 
       @return  REference to this */
   AliFMDInput& operator=(const AliFMDInput&) { return *this; }

+  /** 
+   * Get user supplued data
+   * 
+   * @param d Detector
+   * @param r Ring
+   * @param s Sector
+   * @param t Strip
+   * 
+   * @return Value
+   */
+  virtual Float_t GetSignal(UShort_t d, Char_t r, UShort_t s, UShort_t t);

 
   TString          fGAliceFile; // File name of gAlice file
   AliRunLoader*    fLoader;     // Loader of FMD data 
 
   TString          fGAliceFile; // File name of gAlice file
   AliRunLoader*    fLoader;     // Loader of FMD data 


@@ -348,6 +373,10 @@ inline Bool_t AliFMDInput::ProcessRawCalibDigit(AliFMDDigit*) { return kTRUE; }
 inline Bool_t AliFMDInput::ProcessRecPoint(AliFMDRecPoint*) { return kTRUE; }
 inline Bool_t AliFMDInput::ProcessESD(UShort_t,Char_t,UShort_t,UShort_t,
                                      Float_t,Float_t) { return kTRUE; }
 inline Bool_t AliFMDInput::ProcessRecPoint(AliFMDRecPoint*) { return kTRUE; }
 inline Bool_t AliFMDInput::ProcessESD(UShort_t,Char_t,UShort_t,UShort_t,
                                      Float_t,Float_t) { return kTRUE; }

+inline Bool_t AliFMDInput::ProcessUser(UShort_t,Char_t,UShort_t,UShort_t,
+                                      Float_t) { return kTRUE; }
+inline Float_t AliFMDInput::GetSignal(UShort_t, Char_t, UShort_t, UShort_t) { 
+  return 0.; }

 
 
 #endif
 
 
 #endif

diff --git a/FMD/scripts/FindCommonModeNoise.C b/FMD/scripts/FindCommonModeNoise.C
new file mode 100644 (file)
index 0000000..71a7d62
--- /dev/null
+++ b/FMD/scripts/FindCommonModeNoise.C
@@ -0,0 +1,368 @@
+//____________________________________________________________________
+//
+// $Id: DrawDigits.C 30718 2009-01-22 16:07:40Z cholm $
+//
+// Script that contains a class to draw eloss from hits, versus ADC
+// counts from digits, using the AliFMDInputHits class in the util library. 
+//
+// It draws the energy loss versus the p/(mq^2).  It can be overlayed
+// with the Bethe-Bloc curve to show how the simulation behaves
+// relative to the expected. 
+//
+// Use the script `Compile.C' to compile this class using ACLic. 
+//
+#include <TH1F.h>
+#include <TMath.h>
+#include <TStyle.h>
+#include <TArrayF.h>
+#include <TSystem.h>
+#include <TFile.h>
+#include <TCanvas.h>
+#include <TSystem.h>
+#include <TRandom.h>
+#include <TPad.h>
+#include <TLegend.h>
+
+#include <AliCDBManager.h>
+#include <AliFMDDigit.h>
+#include <AliFMDInput.h>
+#include <AliFMDParameters.h>
+#include <AliRawReader.h>
+#include <AliLog.h>
+
+#include <iostream>
+
+/** @class DrawDigits
+    @brief Draw pedestal, noise, and both common mode noise corrected
+    @code 
+    Root> .L Compile.C
+    Root> Compile("FindCommonModeNoise.C")
+    Root> FindCommonModeNoise fcmn("file");
+    Root> fcmn.Run();
+    @endcode
+
+    If null is passed instead of a file, then a simulation of common
+    mode noise is run instead.
+
+    @ingroup FMD_script
+ */
+class FindCommonModeNoise : public AliFMDInput
+{
+private:
+  TString fCalibDir;
+  Int_t   fCount;
+  TFile*  fOut;
+  Float_t fShift;
+  Float_t fCmn;                // Common mode noise component
+  TH1F*   fPed;                // 
+  TH1F*   fNoise;              //
+  TH1F*   fCenteredPed;        //
+  TH1F*   fCorrectedNoise;     //
+  TH1F*   fSummedAdc;          //
+  TH1F*   fSummedAdc2;         //
+  TH1F*   fSummedCenteredAdc;  //
+  TH1F*   fSummedCenteredAdc2; //
+  TH1F*   fAdc;                //
+
+public:
+  //__________________________________________________________________
+  FindCommonModeNoise(const char*  file, 
+                     const char*  calibDir=0) 
+    : AliFMDInput(), 
+      fCalibDir(""), 
+      fOut(0)
+  { 
+    fCmn = 2;
+    if (calibDir) fCalibDir = calibDir;
+    if (file) { 
+      AddLoad(kRaw);
+      SetRawFile(file);
+    }
+    else 
+      AddLoad(kUser);
+    
+    TString outName(Form("histo_%s", gSystem->BaseName(fRawFile.Data())));
+    if (!outName.EndsWith(".root")) outName.Append(".root");
+    fOut  = TFile::Open(outName.Data(),"RECREATE");
+
+    fPed                = MakeHist("ped_dist","Pedestal Distribution");
+    fNoise              = MakeHist("noise_dist",
+                                  "Noise Distribution (Not CMN corrected)");
+    fCenteredPed        = MakeHist("cmn_ped_dist",
+                                  "Pedestal - Average Pedestal Distribution");
+    fCorrectedNoise     = MakeHist("noise_corrected_dist",
+                                  "Noise Distribution (CMN corrected)");
+    fSummedAdc          = MakeHist("sum_adc_dist","Sum ADC Distribution");
+    fSummedAdc2         = MakeHist("sum_adc2_dist",
+                                  "Sum ADC Squared Distribution");
+    fSummedCenteredAdc  = MakeHist("sum_cmn_adc_dist",
+                                  "Sum CMN Adjusted ADC Distribution");
+    fSummedCenteredAdc2 = MakeHist("sum_cmn_adc2_dist",
+                                  "Sum CMN Adjusted ADC Squared Distribution");
+    fAdc                = MakeHist("adc_dist","Current ADC Distribution");
+  }
+  //__________________________________________________________________
+  Int_t NEvents() const 
+  { 
+    Int_t nEv = (fReader ? fReader->GetNumberOfEvents() : -1);
+    if (nEv > 0) return TMath::Min(nEv, 1000);
+    return 1000;
+  }
+  //__________________________________________________________________
+  TH1F* MakeHist(const char* name, const char* title)
+  {
+    TH1F* h = new TH1F(name,title,51200,-0.5,51199.5);
+    h->SetXTitle("Strip number");
+    h->SetFillColor(kRed);
+    h->SetFillStyle(3001);
+    h->SetDirectory(0);
+    h->SetStats(kFALSE);
+    return h;
+  }
+  //__________________________________________________________________
+  Bool_t Init()
+  {
+    AliCDBManager* cdb = AliCDBManager::Instance();
+    cdb->SetRun(0);
+    cdb->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
+
+    AliFMDParameters* pars = AliFMDParameters::Instance();
+    if (fCalibDir.IsNull())
+      pars->Init(kFALSE);
+    else 
+      pars->Init(fCalibDir.Data(), kFALSE);
+
+    fCount = 0;
+
+    return AliFMDInput::Init();
+  }
+
+  //__________________________________________________________________
+  Bool_t Begin(Int_t e)
+  {
+    fCount++;
+    fAdc->Reset();
+    // Common mode noise generation
+    fShift = gRandom->Gaus(0, fCmn);
+    return AliFMDInput::Begin(e);
+  }
+
+  //__________________________________________________________________
+  Int_t BinCenter(UShort_t d, Char_t r, UShort_t s, UShort_t t) const
+  {
+    if (d == 1 && !(r == 'I' || r == 'i')) return -1;
+
+    Int_t idx = -1;
+    switch (d) { 
+    case 1: idx = 0; break;
+    case 2: idx = 10240; break;
+    case 3: idx = 3 * 10240; break;
+    default: 
+      return -1;
+    }
+    
+    UShort_t q    = (r == 'I' || r == 'i' ? 0 : 1);
+    // UShort_t nSec = (q == 0 ?  20 :  40);
+    UShort_t nStr = (q == 0 ? 512 : 256);
+
+    idx += (q == 0 ? 0 : 10240);
+    idx += s * nStr + t;
+
+    if (idx >= 51200) return -1;
+    return idx;
+  }
+  
+  //__________________________________________________________________
+  Int_t BinCenter(UShort_t chip, UShort_t strip) const 
+  {
+    return chip * 128 + strip;
+  }
+  //__________________________________________________________________
+  Int_t BinNumber(UShort_t chip, UShort_t strip) const 
+  {
+    return 1 + BinCenter(chip, strip);
+  }
+  //__________________________________________________________________
+  Float_t GetSignal(UShort_t d, Char_t, UShort_t s, UShort_t t)
+  {
+    // Calculate signal value.  Only the first 256 strips of FMD are
+    // filled. 
+    if (d > 1 || s > 0 || t > 256) return 0;
+    Float_t c = 0.005;
+    Float_t x = (t % 128);
+    Float_t p = c * (64 * 64 - 128 * x + x * x) + 100;
+    Float_t n = (p-100)/40 + 1;
+
+    // In case we simulate a common mode noise component. 
+    if (t < 128) 
+      return fShift + gRandom->Gaus(p, n);
+   
+    // In case all noise is uncorrelated. 
+    return gRandom->Gaus(p, TMath::Sqrt(n*n + fCmn*fCmn));
+  }    
+  //__________________________________________________________________
+  Bool_t ProcessRawDigit(AliFMDDigit* digit)
+  {
+    // From data 
+    if (!digit) return kTRUE;
+
+    UShort_t d             =  digit->Detector();
+    Char_t   r             =  digit->Ring();
+    UShort_t s             =  digit->Sector();
+    UShort_t t             =  digit->Strip();
+    return ProcessOne(d, r, s, t, digit->Counts());
+  }
+  //__________________________________________________________________
+  Bool_t ProcessUser(UShort_t d, Char_t r, UShort_t s, UShort_t t, Float_t v)
+  {
+    // From our simulator in GetValue
+    return ProcessOne(d, r, s, t, v);
+  }
+  //__________________________________________________________________
+  Bool_t ProcessOne(UShort_t d, Char_t r, UShort_t s, UShort_t t, Float_t v)
+  {
+    Int_t    i             =  BinCenter(d, r, s, t);
+    if (i < 0) { 
+      std::cout << "Invalid index " << i << " returned for FMD" 
+               << d << r << '[' << s << ',' << t << ']' << std::endl;
+      return kFALSE;
+    }
+
+    fAdc->Fill(i, v);
+
+    return kTRUE;
+  }
+  //__________________________________________________________________
+  Bool_t End()
+  {
+    // At the end of each event, fill the summed histograms. 
+    for (UShort_t chip = 0; chip < 400; chip++) { 
+      Double_t chipAv = 0;
+      for (UShort_t strip = 0; strip < 128; strip++) 
+       chipAv += fAdc->GetBinContent(1+BinCenter(chip, strip));
+      chipAv /= 128;
+      
+      for (UShort_t strip = 0; strip < 128; strip++) {
+       Int_t    i    = BinCenter(chip, strip);
+       Double_t adc  = fAdc->GetBinContent(i+1);
+       Double_t cAdc = chipAv - adc;
+       fSummedAdc->Fill(i, adc);
+       fSummedAdc2->Fill(i, adc * adc);
+       fSummedCenteredAdc->Fill(i, cAdc);
+       fSummedCenteredAdc2->Fill(i, cAdc*cAdc);
+      }
+    }
+
+    return AliFMDInput::End();
+  }
+  //__________________________________________________________________
+  Bool_t Finish()
+  {
+    if (fEventCount == 0) { 
+      std::cerr << "No events!" << std::endl;
+      return kFALSE;
+    }
+
+    for (UInt_t bin = 1; bin <= 51200; bin++) { 
+      if (bin % (51200 / 10) == 0) 
+       std::cout << "Bin # " << bin << std::endl;
+      Double_t sumAdc  = fSummedAdc->GetBinContent(bin);
+      Double_t sumAdc2 = fSummedAdc2->GetBinContent(bin);
+      Double_t avAdc   = sumAdc / fEventCount;
+      Double_t avAdc2  = sumAdc2 / fEventCount;
+      Double_t noise2  = avAdc2 - avAdc * avAdc;
+    
+      fPed->SetBinContent(bin, avAdc);
+      fNoise->SetBinContent(bin, noise2 >= 0 ? TMath::Sqrt(noise2) : 0);
+      
+      Double_t sumCAdc  = fSummedCenteredAdc->GetBinContent(bin);
+      Double_t sumCAdc2 = fSummedCenteredAdc2->GetBinContent(bin);
+      Double_t avCAdc   = sumCAdc / fEventCount;
+      Double_t avCAdc2  = sumCAdc2 / fEventCount;
+      Double_t cNoise2  = avCAdc2 - avCAdc * avCAdc;
+      fCenteredPed->SetBinContent(bin, avCAdc);
+      fCorrectedNoise->SetBinContent(bin, TMath::Sqrt(cNoise2));
+    }
+    
+    gStyle->SetPalette(1);
+    gStyle->SetOptTitle(0);
+    gStyle->SetCanvasColor(0);
+    gStyle->SetCanvasBorderSize(0);
+    gStyle->SetPadColor(0);
+    gStyle->SetPadBorderSize(0);
+  
+    TCanvas* c = new TCanvas("c", "C", 800, 600);
+    c->SetFillColor(0);
+    c->SetBorderMode(0);
+    c->SetBorderSize(0);
+    c->SetTopMargin(0);
+    c->SetBottomMargin(0);
+    c->Divide(2, 1);
+    
+    TPad* p1 = static_cast<TPad*>(c->cd(1));
+    p1->SetTopMargin(0.05);
+    p1->SetRightMargin(0.05);
+    if (!fReader) fNoise->GetXaxis()->SetRangeUser(-0.5, 255.5);
+    fNoise->SetMinimum(0);
+    fNoise->SetMaximum(1.5*fNoise->GetMaximum());
+    fNoise->DrawCopy();
+    fCorrectedNoise->SetMinimum(0);
+    fCorrectedNoise->SetFillColor(kBlue);
+    fCorrectedNoise->DrawCopy("same");
+    TLegend* l1 = new TLegend(0.2, 0.7, 0.945, 0.945);
+    l1->SetFillColor(0);
+    l1->SetFillStyle(0);
+    l1->SetBorderSize(0);
+    l1->AddEntry(fNoise, "Noise", "f");
+    l1->AddEntry(fCorrectedNoise, "Corrected noise", "f");
+    l1->Draw();
+      
+
+    TPad* p2 = static_cast<TPad*>(c->cd(2));
+    p2->SetTopMargin(0.05);
+    p2->SetRightMargin(0.05);
+    if (!fReader) fPed->GetXaxis()->SetRangeUser(-0.5, 255.5);
+    fPed->SetMinimum(fCenteredPed->GetMinimum());
+    fPed->SetMaximum(1.5*fPed->GetMaximum());
+    fPed->DrawCopy();
+    fCenteredPed->SetFillColor(kBlue);
+    fCenteredPed->DrawCopy("same");
+    TLegend* l2 = new TLegend(0.2, 0.7, 0.945, 0.945);
+    l2->SetFillColor(0);
+    l2->SetFillStyle(0);
+    l2->SetBorderSize(0);
+    l2->AddEntry(fPed, "Pedestal", "f");
+    l2->AddEntry(fCenteredPed, "Pedestal - #bar{Pedestal}", "f");
+    l2->Draw();
+
+    c->Modified();
+    c->Update();
+    c->cd();
+    c->Print("cmn.png");
+
+    if (fOut) { 
+      std::cout << "Flusing to disk ... " << std::flush;
+      fOut->cd();
+      fPed->Write();                //
+      fNoise->Write();              //
+      fCenteredPed->Write();        //
+      fCorrectedNoise->Write();     //
+      fSummedAdc->Write();          //
+      fSummedAdc2->Write();         //
+      fSummedCenteredAdc->Write();  //
+      fSummedCenteredAdc2->Write(); //
+      fAdc->Write();                //
+      fOut->Write();
+      fOut->Close();
+      std::cout << "done" << std::endl;
+    }
+    return kTRUE;
+  }
+
+  ClassDef(FindCommonModeNoise,0);
+};
+
+//____________________________________________________________________
+//
+// EOF
+//