Example of low-level data processing

diff --git a/PHOS/macros/RawProcessing/PedAmp.C b/PHOS/macros/RawProcessing/PedAmp.C
new file mode 100644 (file)
index 0000000..cd38454
--- /dev/null
+++ b/PHOS/macros/RawProcessing/PedAmp.C
@@ -0,0 +1,250 @@
+#if !defined(__CINT__) || defined(__MAKECINT__)
+#include "AliRawReaderRoot.h"
+#include "AliCaloRawStream.h"
+#include "TString.h"
+#include "TH2F.h"
+#include "TH1S.h"
+#include "TFile.h"
+#include "TStopwatch.h"
+#include "iostream.h"
+
+#endif
+
+/* $Id$ */
+
+void PedAmp(const char *rawFile, const Int_t debug=0)
+{
+  // Read raw data, decode it to samples,
+  // calculate pedestals from presamples, 
+  // evaluate the signal amplitude as a maximum sample, 
+  // and fill histograms with pedestals and amplitudes
+  // This script should be compiled to speed up the data processing:
+  // .L PedAmp.C++
+  //___
+  // Yuri Kharlov. 6 September 2007
+
+  TStopwatch stopwatch;
+  stopwatch.Start();
+  
+  AliRawReaderRoot* rf = new AliRawReaderRoot(rawFile);
+  AliCaloRawStream in(rf,"PHOS");
+  in.SetOldRCUFormat(kTRUE);
+  const Int_t nPresample = 10;
+
+  TString baseNamePed ="hPed";
+  TString baseTitlePed="Ped in cell (";
+  TString baseNameLed ="hLED";
+  TString baseTitleLed="LED in cell (";
+  const char* sgain[2]={"high","low"};
+
+  const Int_t gainMax=2,modMax=5,colMax=56,rowMax=64;
+  TH1F *hPed[2][56][64];
+  TH1F *hLed[2][56][64];
+  for (Int_t gain=0; gain<gainMax; gain++) {
+    for (Int_t mod=0; mod<modMax; mod++) {
+      for (Int_t col=0; col<colMax; col++) {
+       for (Int_t row=0; row<rowMax; row++) {
+         hPed[gain][col][row] = 0;
+         hLed[gain][col][row] = 0;
+       }
+      }
+    }
+  }
+  TH1F *hPedHiMean1 = new TH1F("hPedHiMean1","Mean pedestals, high gain" ,100,0.,100.);
+  TH1F *hPedHiRMS1  = new TH1F("hPedHiRMS1" ,"RMS pedestals, high gain"  ,100,0.,50.);
+  TH1F *hPedLoMean1 = new TH1F("hPedLoMean1","Mean pedestals, low gain"  ,100,0.,100.);
+  TH1F *hPedLoRMS1  = new TH1F("hPedLoRMS1" ,"RMS pedestals, low gain"   ,100,0.,50.);
+  hPedHiMean1->Sumw2();
+  hPedHiRMS1 ->Sumw2();
+  hPedLoMean1->Sumw2();
+  hPedLoRMS1 ->Sumw2();
+
+  TH2F *hPedHiMean  = new TH2F("hPedHiMean","Mean pedestals, high gain",
+                              rowMax,0.,rowMax, colMax,0.,colMax);
+  TH2F *hPedHiRMS   = new TH2F("hPedHiRMS" ,"R.M.S. of pedestals, high gain",
+                              rowMax,0.,rowMax, colMax,0.,colMax);
+  TH2F *hPedHiNum   = new TH2F("hPedHiNum" ,"Number of pedestals, high gain",
+                              rowMax,0.,rowMax, colMax,0.,colMax);
+  TH2F *hPedLoMean  = new TH2F("hPedLoMean","Mean pedestals, low gain",
+                              rowMax,0.,rowMax, colMax,0.,colMax);
+  TH2F *hPedLoRMS   = new TH2F("hPedLoRMS" ,"R.M.S. of pedestals, low gain",
+                              rowMax,0.,rowMax, colMax,0.,colMax);
+  TH2F *hPedLoNum   = new TH2F("hPedLoNum" ,"Number of pedestals, low gain",
+                              rowMax,0.,rowMax, colMax,0.,colMax);
+  
+  TH2F *hLedHiMean  = new TH2F("hLedHiMean","Mean LED, high gain",
+                              rowMax,0.,rowMax, colMax,0.,colMax);
+  TH2F *hLedHiRMS   = new TH2F("hLedHiRMS" ,"R.M.S. of LED, high gain",
+                              rowMax,0.,rowMax, colMax,0.,colMax);
+  TH2F *hLedHiNum   = new TH2F("hLedHiNum" ,"Number of LED, high gain",
+                              rowMax,0.,rowMax, colMax,0.,colMax);
+  TH2F *hLedLoMean  = new TH2F("hLedLoMean","Mean LED, low gain",
+                              rowMax,0.,rowMax, colMax,0.,colMax);
+  TH2F *hLedLoRMS   = new TH2F("hLedLoRMS" ,"R.M.S. of LED, low gain",
+                              rowMax,0.,rowMax, colMax,0.,colMax);
+  TH2F *hLedLoNum   = new TH2F("hLedLoNum" ,"Number of LED, low gain",
+                              rowMax,0.,rowMax, colMax,0.,colMax);
+
+  Int_t iEvent=0;
+  Int_t iBin=0;
+  Int_t gain=-1,mod=-1,col=-1,row=-1;
+  Int_t runNum = 0;
+  Double_t maxAmp = 0;
+//   Double_t meanPed,rmsPed;
+  Float_t pedMean = 0;
+  Int_t   nPed = 0;
+  while (rf->NextEvent()) {
+    cout << "\n\n\tEvent "<< iEvent << " of type " << rf->GetType() << endl;
+    runNum = rf->GetRunNumber();
+    while ( in.Next() ) {
+      if (!in.IsNewHWAddress() && debug > 1) {
+       cout << " time=" << in.GetTime() << "/"<< in.GetTimeLength()
+            << " amp=" <<  in.GetSignal() <<endl;
+      }
+      if (in.IsNewHWAddress()) {
+       iBin   = 0;
+       if (gain!=-1 && col!=-1 && row!=-1) {
+         if (nPed < 1) continue;
+         maxAmp -= (Double_t)(pedMean/nPed); // "pedestal subtraction"
+         if (maxAmp > 5) hLed[gain][col][row]->Fill(maxAmp);
+       }
+       maxAmp  = 0;
+       pedMean = 0;
+       nPed    = 0;
+       gain = in.IsLowGain(); // 0-high, 1-low
+       mod  = in.GetModule();
+       col  = in.GetColumn();
+       row  = in.GetRow();
+       if (debug > 0) {
+         cout << "\n\tNew HW address: "<<in.GetHWAddress()
+              << ", gain="<< gain
+              << ", module="<<mod
+              << ", xz=("<<col<<","<<row<<")"<< endl;
+       }
+
+       // Create a new histo for the new cell, if it does not exist yet
+       if (!hPed[gain][col][row]) {
+         TString name  = baseNamePed;
+         TString title = baseTitlePed;
+         name +="_g"; name +=gain;
+         name +="_m"; name +=mod;
+         name +="_z"; name +=col;
+         name +="_x"; name +=row;
+
+         title +=mod; title +=",";
+         title +=col; title +=",";
+         title +=row; title +="), ";
+         title +=sgain[gain]; title +=" gain";
+
+         hPed[gain][col][row] = new TH1F(name,title,100,0.,100.);
+         hPed[gain][col][row]->Sumw2();
+         hPed[gain][col][row]->SetMarkerStyle(20);
+         hPed[gain][col][row]->SetOption("eph");
+       }
+       if (!hLed[gain][col][row]) {
+         TString name  = baseNameLed;
+         TString title = baseTitleLed;
+         name +="_g"; name +=gain;
+         name +="_m"; name +=mod;
+         name +="_z"; name +=col;
+         name +="_x"; name +=row;
+
+         title +=mod; title +=",";
+         title +=col; title +=",";
+         title +=row; title +="), ";
+         title +=sgain[gain]; title +=" gain";
+
+         hLed[gain][col][row] = new TH1F(name,title,1023,0.,1023.);
+         hLed[gain][col][row]->Sumw2();
+         hLed[gain][col][row]->SetMarkerStyle(20);
+         hLed[gain][col][row]->SetOption("eph");
+       }
+      }
+
+      iBin++;
+      // Assume that first nPresample are pedestals
+      if (in.GetTimeLength()-iBin < nPresample) {
+       hPed[gain][col][row]->Fill(in.GetSignal());
+       pedMean += in.GetSignal();
+       nPed++;
+      }
+      if (in.GetSignal() > maxAmp) maxAmp = in.GetSignal();
+    } // end of next signal
+    iEvent++;
+  } // end of next event
+
+  // Fill 2-dim histograms for mean, rms and n pedestals
+
+  gain = 0;
+  for (Int_t col=0; col<colMax; col++) {
+    for (Int_t row=0; row<rowMax; row++) {
+      if (hPed[gain][col][row] != 0) {
+       hPedHiMean1->Fill( hPed[gain][col][row]->GetMean());
+       hPedHiRMS1 ->Fill( hPed[gain][col][row]->GetRMS() );
+       hPedHiMean ->Fill( (Double_t)row, (Double_t)col, hPed[gain][col][row]->GetMean()    );
+       hPedHiRMS  ->Fill( (Double_t)row, (Double_t)col, hPed[gain][col][row]->GetRMS()     );
+       hPedHiNum  ->Fill( (Double_t)row, (Double_t)col, hPed[gain][col][row]->GetEntries() );
+      }
+      if (hLed[gain][col][row] != 0) {
+       hLedHiMean ->Fill( (Double_t)row, (Double_t)col, hLed[gain][col][row]->GetMean()    );
+       hLedHiRMS  ->Fill( (Double_t)row, (Double_t)col, hLed[gain][col][row]->GetRMS()     );
+       hLedHiNum  ->Fill( (Double_t)row, (Double_t)col, hLed[gain][col][row]->GetEntries() );
+      }
+    }
+  }
+  gain = 1;
+  for (Int_t col=0; col<colMax; col++) {
+    for (Int_t row=0; row<rowMax; row++) {
+      if (hPed[gain][col][row] != 0) {
+       hPedLoMean1->Fill( hPed[gain][col][row]->GetMean());
+       hPedLoRMS1 ->Fill( hPed[gain][col][row]->GetRMS() );
+       hPedLoMean ->Fill( (Double_t)row, (Double_t)col, hPed[gain][col][row]->GetMean()    );
+       hPedLoRMS  ->Fill( (Double_t)row, (Double_t)col, hPed[gain][col][row]->GetRMS()     );
+       hPedLoNum  ->Fill( (Double_t)row, (Double_t)col, hPed[gain][col][row]->GetEntries() );
+      }
+      if (hLed[gain][col][row] != 0) {
+       hLedLoMean ->Fill( (Double_t)row, (Double_t)col, hLed[gain][col][row]->GetMean()    );
+       hLedLoRMS  ->Fill( (Double_t)row, (Double_t)col, hLed[gain][col][row]->GetRMS()     );
+       hLedLoNum  ->Fill( (Double_t)row, (Double_t)col, hLed[gain][col][row]->GetEntries() );
+      }
+    }
+  }
+
+  // Write existing histograms to a root file
+
+  TString fileName = "ped";
+  fileName += runNum;
+  fileName += ".root";
+  TFile *file = new TFile(fileName,"RECREATE");
+
+  for (Int_t gain=0; gain<gainMax; gain++) {
+    for (Int_t col=0; col<colMax; col++) {
+      for (Int_t row=0; row<rowMax; row++) {
+       if (hPed[gain][col][row] != 0)
+         hPed[gain][col][row]->Write();
+       if (hLed[gain][col][row] != 0)
+         hLed[gain][col][row]->Write();
+      }
+    }
+  }
+  hPedHiMean1->Write();
+  hPedHiRMS1 ->Write();
+  hPedLoMean1->Write();
+  hPedLoRMS1 ->Write();
+  hPedHiMean ->Write();
+  hPedHiRMS  ->Write();
+  hPedHiNum  ->Write();
+  hPedLoMean ->Write();
+  hPedLoRMS  ->Write();
+  hPedLoNum  ->Write();
+  
+  hLedHiMean ->Write();
+  hLedHiRMS  ->Write();
+  hLedHiNum  ->Write();
+  hLedLoMean ->Write();
+  hLedLoRMS  ->Write();
+  hLedLoNum  ->Write();
+  
+  file->Close();
+  stopwatch.Print();
+}