1.The QA data created on demand according to the event species at filling time. 2...
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALQADataMakerRec.cxx
index 665a601..f7f2d7e 100644 (file)
@@ -26,9 +26,11 @@ Authors:  J.Klay (Cal Poly) May 2008
 #include <TH1F.h> 
 #include <TH1I.h> 
 #include <TH2F.h> 
+#include <TProfile.h> 
 
 // --- Standard library ---
 
+
 // --- AliRoot header files ---
 #include "AliESDCaloCluster.h"
 #include "AliESDCaloCells.h"
@@ -36,28 +38,32 @@ Authors:  J.Klay (Cal Poly) May 2008
 #include "AliLog.h"
 #include "AliEMCALQADataMakerRec.h"
 #include "AliQAChecker.h"
+#include "AliEMCALDigit.h" 
 #include "AliEMCALRecPoint.h" 
 #include "AliEMCALRawUtils.h"
 #include "AliEMCALReconstructor.h"
 #include "AliEMCALRecParam.h"
 #include "AliRawReader.h"
+#include "AliCaloRawStream.h"
 
 ClassImp(AliEMCALQADataMakerRec)
            
 //____________________________________________________________________________ 
   AliEMCALQADataMakerRec::AliEMCALQADataMakerRec() : 
-  AliQADataMakerRec(AliQA::GetDetName(AliQA::kEMCAL), "EMCAL Quality Assurance Data Maker")
+    AliQADataMakerRec(AliQAv1::GetDetName(AliQAv1::kEMCAL), "EMCAL Quality Assurance Data Maker"),
+    fSuperModules(4) // FIXME!!! number of SuperModules; 4 for 2009; update default to 12 for later runs..
 {
   // ctor
 }
 
 //____________________________________________________________________________ 
 AliEMCALQADataMakerRec::AliEMCALQADataMakerRec(const AliEMCALQADataMakerRec& qadm) :
-  AliQADataMakerRec()
+  AliQADataMakerRec(), fSuperModules()
 {
   //copy ctor 
   SetName((const char*)qadm.GetName()) ; 
   SetTitle((const char*)qadm.GetTitle()); 
+  fSuperModules = qadm.GetSuperModules();
 }
 
 //__________________________________________________________________
@@ -70,50 +76,71 @@ AliEMCALQADataMakerRec& AliEMCALQADataMakerRec::operator = (const AliEMCALQAData
 }
  
 //____________________________________________________________________________ 
-void AliEMCALQADataMakerRec::EndOfDetectorCycle(AliQA::TASKINDEX_t task, TObjArray ** list)
+void AliEMCALQADataMakerRec::EndOfDetectorCycle(AliQAv1::TASKINDEX_t task, TObjArray ** list)
 {
   //Detector specific actions at end of cycle
   // do the QA checking
-  AliQAChecker::Instance()->Run(AliQA::kEMCAL, task, list) ;  
+  AliQAChecker::Instance()->Run(AliQAv1::kEMCAL, task, list) ;  
 }
 
 //____________________________________________________________________________ 
 void AliEMCALQADataMakerRec::InitESDs()
 {
   //Create histograms to controll ESD
-  TH1F * h1 = new TH1F("hESDCaloClusterE",  "ESDs CaloCluster energy in EMCAL",    200, 0., 20.) ; 
+  const Bool_t expert   = kTRUE ; 
+  const Bool_t image    = kTRUE ; 
+  
+  TH1F * h1 = new TH1F("hESDCaloClusterE",  "ESDs CaloCluster energy in EMCAL;Energy [MeV];Counts",    200, 0., 20.) ; 
   h1->Sumw2() ;
-  Add2ESDsList(h1, kESDCaloClusE)  ;                                                     
+  Add2ESDsList(h1, kESDCaloClusE, !expert, image)  ;                                                     
 
-  TH1I * h2 = new TH1I("hESDCaloClusterM", "ESDs CaloCluster multiplicity in EMCAL", 100, 0,  100) ; 
+  TH1I * h2 = new TH1I("hESDCaloClusterM", "ESDs CaloCluster multiplicity in EMCAL;# of Clusters;Entries", 100, 0,  100) ; 
   h2->Sumw2() ;
-  Add2ESDsList(h2, kESDCaloClusM)  ;
+  Add2ESDsList(h2, kESDCaloClusM, !expert, image)  ;
 
-  TH1F * h3 = new TH1F("hESDCaloCellA",  "ESDs CaloCell amplitude in EMCAL",    500, 0., 250.) ; 
+  TH1F * h3 = new TH1F("hESDCaloCellA",  "ESDs CaloCell amplitude in EMCAL;Energy [MeV];Counts",    500, 0., 250.) ; 
   h3->Sumw2() ;
-  Add2ESDsList(h3, kESDCaloCellA)  ;  
+  Add2ESDsList(h3, kESDCaloCellA, !expert, image)  ;  
  
-  TH1I * h4 = new TH1I("hESDCaloCellM", "ESDs CaloCell multiplicity in EMCAL", 200, 0,  1000) ; 
+  TH1I * h4 = new TH1I("hESDCaloCellM", "ESDs CaloCell multiplicity in EMCAL;# of Clusters;Entries", 200, 0,  1000) ; 
   h4->Sumw2() ;
-  Add2ESDsList(h4, kESDCaloCellM) ;
+  Add2ESDsList(h4, kESDCaloCellM, !expert, image) ;
        
 }
 
 //____________________________________________________________________________ 
+void AliEMCALQADataMakerRec::InitDigits()
+{
+  // create Digits histograms in Digits subdir
+  const Bool_t expert   = kTRUE ; 
+  const Bool_t image    = kTRUE ; 
+  
+  TH1I * h0 = new TH1I("hEmcalDigits",    "Digits amplitude distribution in EMCAL;Amplitude [ADC counts];Counts",    500, 0, 500) ; 
+  h0->Sumw2() ;
+  Add2DigitsList(h0, 0, !expert, image) ;
+  TH1I * h1 = new TH1I("hEmcalDigitsMul", "Digits multiplicity distribution in EMCAL;# of Digits;Entries", 200, 0, 2000) ; 
+  h1->Sumw2() ;
+  Add2DigitsList(h1, 1, !expert, image) ;
+}
+
+//____________________________________________________________________________ 
 void AliEMCALQADataMakerRec::InitRecPoints()
 {
   // create Reconstructed Points histograms in RecPoints subdir
-  TH1F* h0 = new TH1F("hEMCALRpE","EMCAL RecPoint energies",200, 0.,20.); //GeV
+  const Bool_t expert   = kTRUE ; 
+  const Bool_t image    = kTRUE ; 
+  
+  TH1F* h0 = new TH1F("hEMCALRpE","EMCAL RecPoint energies;Energy [MeV];Counts",200, 0.,20.); //GeV
   h0->Sumw2();
-  Add2RecPointsList(h0,kRecPE);
+  Add2RecPointsList(h0,kRecPE, !expert, image);
 
-  TH1I* h1 = new TH1I("hEMCALRpM","EMCAL RecPoint multiplicities",100,0,100);
+  TH1I* h1 = new TH1I("hEMCALRpM","EMCAL RecPoint multiplicities;# of Clusters;Entries",100,0,100);
   h1->Sumw2();
-  Add2RecPointsList(h1,kRecPM);
+  Add2RecPointsList(h1,kRecPM, !expert, image);
 
-  TH1I* h2 = new TH1I("hEMCALRpDigM","EMCAL RecPoint Digit Multiplicities",20,0,20);
+  TH1I* h2 = new TH1I("hEMCALRpDigM","EMCAL RecPoint Digit Multiplicities;# of Digits;Entries",20,0,20);
   h2->Sumw2();
-  Add2RecPointsList(h2,kRecPDigM);
+  Add2RecPointsList(h2,kRecPDigM, !expert, image);
 
 }
 
@@ -121,43 +148,60 @@ void AliEMCALQADataMakerRec::InitRecPoints()
 void AliEMCALQADataMakerRec::InitRaws()
 {
   // create Raws histograms in Raws subdir
-  //these need more thought
-  /*
-  TH1I * h0 = new TH1I("hLowEmcalSupermodules",    "Low Gain digits in EMCAL supermodules",       12, 0, 12) ;
-  h0->Sumw2() ;
-  Add2RawsList(h0, kNsmodLG) ;
-  TH1I * h1 = new TH1I("hHighEmcalSupermodules",   "High Gain Digits in EMCAL supermodules",       12, 0, 12) ;
-  h1->Sumw2() ;
-  Add2RawsList(h1, kNsmodHG) ;
+   const Bool_t expert   = kTRUE ; 
+   const Bool_t saveCorr = kTRUE ; 
+   const Bool_t image    = kTRUE ; 
 
-  TH1F * h2 = new TH1F("hLowEmcalRawtime", "Low Gain Time of raw digits in EMCAL", 500, -50., 200.) ;
-  h2->Sumw2() ;
-  Add2RawsList(h2, kLGtime) ;
-  TH1F * h3 = new TH1F("hHighEmcalRawtime", "High Gain Time of raw digits in EMCAL", 500, -50., 200.) ;
-  h3->Sumw2() ;
-  Add2RawsList(h3, kHGtime) ;
+  int nTowersPerSM = 1152; // number of towers in a SuperModule; 24x48
+  int nTot = fSuperModules * nTowersPerSM; // max number of towers in all SuperModules
 
-  TH1F * h4 = new TH1F("hLowEmcalRawEnergy", "Low Gain Energy of raw digits in EMCAL", 500, 0., 1000.) ;
-  h4->Sumw2() ;
-  Add2RawsList(h4, kSpecLG) ;
-  TH1F * h5 = new TH1F("hHighEmcalRawEnergy", "High Gain Energy of raw digits in EMCAL",500,0., 1000.) ;
-  h5->Sumw2() ;
-  Add2RawsList(h5, kSpecHG) ;
+  // counter info: number of channels per event (bins are SM index)
+  TProfile * h0 = new TProfile("hLowEmcalSupermodules", "Low Gain EMC: # of towers vs SuperMod;SM Id;# of towers",
+                              fSuperModules, -0.5, fSuperModules-0.5) ;
+  Add2RawsList(h0, kNsmodLG, !expert, image, !saveCorr) ;
+  TProfile * h1 = new TProfile("hHighEmcalSupermodules", "High Gain EMC: # of towers vs SuperMod;SM Id;# of towers",  
+                              fSuperModules, -0.5, fSuperModules-0.5) ; 
+  Add2RawsList(h1, kNsmodHG, !expert, image, !saveCorr) ;
 
-  TH1I * h6 = new TH1I("hLowNtot", "Low Gain Total Number of raw digits in EMCAL", 500, 0, 10000) ;
+  // where did max sample occur? (bins are towers)
+  TProfile * h2 = new TProfile("hLowEmcalRawtime", "Low Gain EMC: Time at Max vs towerId;Tower Id;Time [ns]", 
+                              nTot, -0.5, nTot-0.5) ;
+  Add2RawsList(h2, kTimeLG, !expert, image, !saveCorr) ;
+  TProfile * h3 = new TProfile("hHighEmcalRawtime", "High Gain EMC: Time at Max vs towerId;Tower Id;Time [ns]", 
+                              nTot, -0.5, nTot-0.5) ;
+  Add2RawsList(h3, kTimeHG, !expert, image, !saveCorr) ;
+
+  // how much above pedestal was the max sample?  (bins are towers)
+  TProfile * h4 = new TProfile("hLowEmcalRawMaxMinusMin", "Low Gain EMC: Max - Min vs towerId;Tower Id;??", 
+                              nTot, -0.5, nTot-0.5) ;
+  Add2RawsList(h4, kSigLG, !expert, image, !saveCorr) ;
+  TProfile * h5 = new TProfile("hHighEmcalRawMaxMinusMin", "High Gain EMC: Max - Min vs towerId;Tower Id;??",
+                              nTot, -0.5, nTot-0.5) ;
+  Add2RawsList(h5, kSigHG, !expert, image, !saveCorr) ;
+
+  // total counter: channels per event
+  TH1I * h6 = new TH1I("hLowNtot", "Low Gain EMC: Total Number of found towers;# of Towers;Counts", 200, 0, nTot) ;
   h6->Sumw2() ;
-  Add2RawsList(h6, kNtotLG) ;
-  TH1I * h7 = new TH1I("hHighNtot", "High Gain Total Number of raw digits in EMCAL",500,0, 10000) ;
+  Add2RawsList(h6, kNtotLG, !expert, image, !saveCorr) ;
+  TH1I * h7 = new TH1I("hHighNtot", "High Gain EMC: Total Number of found towers;# of Towers;Counts", 200,0, nTot) ;
   h7->Sumw2() ;
-  Add2RawsList(h7, kNtotHG) ;
-
-  TH1F * h8 = new TH1F("hLowEtot", "Low Gain Total Energy of raw digits in EMCAL", 500, 0., 5000.) ;
-  h8->Sumw2() ;
-  Add2RawsList(h8, kEtotLG) ;
-  TH1F * h9 = new TH1F("hHighEtot", "High Gain Total Energy of raw digits in EMCAL",500,0., 100000.) ;
-  h9->Sumw2() ;
-  Add2RawsList(h9, kEtotHG) ;
-  */
+  Add2RawsList(h7, kNtotHG, !expert, image, !saveCorr) ;
+
+  // pedestal (bins are towers)
+  TProfile * h8 = new TProfile("hLowEmcalRawPed", "Low Gain EMC: Pedestal vs towerId;Tower Id;Pedestal [ADC counts]", 
+                              nTot, -0.5, nTot-0.5) ;
+  Add2RawsList(h8, kPedLG, !expert, image, !saveCorr) ;
+  TProfile * h9 = new TProfile("hHighEmcalRawPed", "High Gain EMC: Pedestal vs towerId;Tower Id;Pedestal [ADC counts]",
+                              nTot, -0.5, nTot-0.5) ;
+  Add2RawsList(h9, kPedHG, !expert, image, !saveCorr) ;
+
+  // pedestal rms (standard dev = sqrt of variance estimator for pedestal) (bins are towers)
+  TProfile * h10 = new TProfile("hLowEmcalRawPedRMS", "Low Gain EMC: Pedestal RMS vs towerId;Tower Id;Width [ADC counts]", 
+                               nTot, -0.5, nTot-0.5) ;
+  Add2RawsList(h10, kPedRMSLG, !expert, image, !saveCorr) ;
+  TProfile * h11 = new TProfile("hHighEmcalRawPedRMS", "High Gain EMC: Pedestal RMS vs towerId;Tower Id;Width [ADC counts]",
+                               nTot, -0.5, nTot-0.5) ;
+  Add2RawsList(h11, kPedRMSHG, !expert, image, !saveCorr) ;
   
 }
 
@@ -166,6 +210,11 @@ void AliEMCALQADataMakerRec::MakeESDs(AliESDEvent * esd)
 {
   // make QA data from ESDs
 
+  
+  // Check id histograms already created for this Event Specie
+  if ( ! GetESDsData(kESDCaloClusE) )
+    InitESDs() ;
+  
   Int_t nTot = 0 ; 
   for ( Int_t index = 0; index < esd->GetNumberOfCaloClusters() ; index++ ) {
     AliESDCaloCluster * clu = esd->GetCaloCluster(index) ;
@@ -199,7 +248,182 @@ void AliEMCALQADataMakerRec::MakeRaws(AliRawReader* rawReader)
   //emcal raw data method, but this doesn't seem to be allowed in
   //AliQADataMakerRec.h
 
-       rawReader->Reset() ;
+  // For now, to avoid redoing the expensive signal fits we just
+  // look at max vs min of the signal spextra, a la online usage in
+  // AliCaloCalibPedestal
+
+  // Check id histograms already created for this Event Specie
+  if ( ! GetRawsData(kSigLG) )
+    InitRaws() ;
+
+  rawReader->Reset() ;
+  AliCaloRawStream in(rawReader,"EMCAL"); 
+
+  // setup
+  int nTowersPerSM = 1152; // number of towers in a SuperModule; 24x48
+  int nRows = 24; // number of rows per SuperModule
+  int sampleMin = 0; 
+  int sampleMax = 0x3ff; // 1023 = 10-bit range
+
+  // for the pedestal calculation
+  Bool_t selectPedestalSamples = kTRUE;
+  int firstPedestalSample = 0;
+  int lastPedestalSample = 15;
+
+  // SM counters; decl. should be safe, assuming we don't get more than 12 SuperModules..
+  int nTotalSMLG[12] = {0};
+  int nTotalSMHG[12] = {0};
+
+  // indices for the reading
+  int iSM = 0;
+  int sample = 0;
+  int gain = 0;
+  int time = 0;
+  // counters, on sample level
+  int i = 0; // the sample number in current event.
+  int max = sampleMin, min = sampleMax;//Use these for picking the pedestal
+  int maxTime = 0;
+
+  // for the pedestal calculation
+  int sampleSum = 0; // sum of samples
+  int squaredSampleSum = 0; // sum of samples squared
+  int nSum = 0; // number of samples in sum
+  // calc. quantities
+  double meanPed = 0, squaredMean = 0, rmsPed = 0;
+  
+  while (in.Next()) { // loop over input stream
+    sample = in.GetSignal(); //Get the adc signal
+    time = in.GetTime();
+    if (sample < min) { min = sample; }
+    if (sample > max) { 
+      max = sample;
+      maxTime = time;
+    }
+    i++;
+
+    // should we add it for the pedestal calculation?
+    if ( (firstPedestalSample<=time && time<=lastPedestalSample) || // sample time in range
+        !selectPedestalSamples ) { // or we don't restrict the sample range.. - then we'll take all 
+      sampleSum += sample;
+      squaredSampleSum += sample*sample;
+      nSum++;
+    }
+
+    if ( i >= in.GetTimeLength()) {
+      // calculate pedesstal estimate: mean of possibly selected samples
+      if (nSum > 0) {
+       meanPed = sampleSum / (1.0 * nSum);
+       squaredMean = squaredSampleSum / (1.0 * nSum);
+       // The variance (rms squared) is equal to the mean of the squares minus the square of the mean..
+       rmsPed = sqrt(squaredMean - meanPed*meanPed); 
+      }
+      else {
+       meanPed = 0;
+       squaredMean = 0;
+       rmsPed  = 0;
+      }
+
+      //If we're here then we're done with this tower
+      gain = -1; // init to not valid value
+      if ( in.IsLowGain() ) {
+       gain = 0;
+      }
+      else if ( in.IsHighGain() ) {
+       gain = 1;
+      }
+      
+      iSM = in.GetModule(); //The modules are numbered starting from 0
+
+      if (iSM>=0 && iSM<fSuperModules) { // valid module reading, can go on with filling
+
+       int towerId = iSM*nTowersPerSM + in.GetColumn()*nRows + in.GetRow();
+
+       if (gain == 0) { 
+         //fill the low gain histograms, and counters
+         nTotalSMLG[iSM]++; // one more channel found
+         GetRawsData(kSigLG)->Fill(towerId, max - min);
+         GetRawsData(kTimeLG)->Fill(towerId, maxTime);
+         if (nSum>0) { // only fill pedestal info in case it could be calculated
+           GetRawsData(kPedLG)->Fill(towerId, meanPed);
+           GetRawsData(kPedRMSLG)->Fill(towerId, rmsPed);
+         }
+       } // gain==0
+       else if (gain == 1) {           
+         //fill the high gain ones
+         nTotalSMHG[iSM]++; // one more channel found
+         GetRawsData(kSigHG)->Fill(towerId, max - min);
+         GetRawsData(kTimeHG)->Fill(towerId, maxTime);
+         if (nSum>0) { // only fill pedestal info in case it could be calculated
+           GetRawsData(kPedHG)->Fill(towerId, meanPed);
+           GetRawsData(kPedRMSHG)->Fill(towerId, rmsPed);
+         }
+       }
+      } // SM index OK
+
+      // reset counters
+      max = sampleMin; min = sampleMax;
+      maxTime = 0;
+      i = 0;
+      // also pedestal calc counters
+      sampleSum = 0; // sum of samples
+      squaredSampleSum = 0; // sum of samples squared
+      nSum = 0; // number of samples in sum
+    
+    }//End if, of channel
+   
+  }//end while, of stream
+
+  // let's also fill the SM and event counter histograms
+  int nTotalHG = 0;
+  int nTotalLG = 0;
+  for (iSM=0; iSM<fSuperModules; iSM++) {  
+    nTotalLG += nTotalSMLG[iSM]; 
+    nTotalHG += nTotalSMHG[iSM]; 
+    GetRawsData(kNsmodLG)->Fill(iSM, nTotalSMLG[iSM]); 
+    GetRawsData(kNsmodHG)->Fill(iSM, nTotalSMHG[iSM]); 
+  }
+  GetRawsData(kNtotLG)->Fill(nTotalLG);
+  GetRawsData(kNtotHG)->Fill(nTotalHG);
+
+  // just in case the next rawreader consumer forgets to reset; let's do it here again..
+  rawReader->Reset() ;
+
+  return;
+}
+
+//____________________________________________________________________________
+void AliEMCALQADataMakerRec::MakeDigits(TClonesArray * digits)
+{
+  // makes data from Digits
+  // Check id histograms already created for this Event Specie
+  if ( ! GetDigitsData(0) )
+    InitDigits() ;
+
+  GetDigitsData(1)->Fill(digits->GetEntriesFast()) ; 
+  TIter next(digits) ; 
+  AliEMCALDigit * digit ; 
+  while ( (digit = dynamic_cast<AliEMCALDigit *>(next())) ) {
+    GetDigitsData(0)->Fill( digit->GetAmp()) ;
+  }  
+  
+}
+
+//____________________________________________________________________________
+void AliEMCALQADataMakerRec::MakeDigits(TTree * digitTree)
+{
+  // makes data from Digit Tree
+  TClonesArray * digits = new TClonesArray("AliEMCALDigit", 1000) ; 
+  
+  TBranch * branch = digitTree->GetBranch("EMCAL") ;
+  if ( ! branch ) {
+    AliWarning("EMCAL branch in Digit Tree not found") ; 
+  } else {
+    branch->SetAddress(&digits) ;
+    branch->GetEntry(0) ; 
+    MakeDigits(digits) ; 
+  }
+  
 }
 
 //____________________________________________________________________________
@@ -211,6 +435,11 @@ void AliEMCALQADataMakerRec::MakeRecPoints(TTree * clustersTree)
     AliError("can't get the branch with the EMCAL clusters !");
     return;
   }
+
+  // Check id histograms already created for this Event Specie
+  if ( ! GetRecPointsData(kRecPM) )
+    InitRecPoints() ;
+  
   TObjArray * emcrecpoints = new TObjArray(100) ;
   emcbranch->SetAddress(&emcrecpoints);
   emcbranch->GetEntry(0);
@@ -233,3 +462,4 @@ void AliEMCALQADataMakerRec::StartOfDetectorCycle()
   //Detector specific actions at start of cycle
   
 }
+