Added QA for digits during reconstruction (Yves)
[u/mrichter/AliRoot.git] / TPC / AliTPCQADataMakerRec.cxx
index 23b290f..c23a9ad 100644 (file)
@@ -49,6 +49,8 @@
 
 // --- ROOT system ---
 #include <TClonesArray.h>
+#include <TString.h>
+#include <TSystem.h>
 
 // --- Standard library ---
 
 #include "AliESDtrack.h"
 #include "AliLog.h"
 #include "AliTPCCalPad.h"
+#include "AliTPCCalROC.h"
 #include "AliTPCClustersRow.h"
 #include "AliTPCclusterMI.h"
+#include "AliSimDigits.h"
 
 ClassImp(AliTPCQADataMakerRec)
 
 //____________________________________________________________________________ 
 AliTPCQADataMakerRec::AliTPCQADataMakerRec() : 
-  AliQADataMakerRec(AliQA::GetDetName(AliQA::kTPC), 
+  AliQADataMakerRec(AliQAv1::GetDetName(AliQAv1::kTPC), 
                    "TPC Rec Quality Assurance Data Maker"),
-  fTPCdataQA(0),
-  fHistESDclusters(0),fHistESDratio(0), fHistESDpt(0),
-  fHistRawsOccupancy(0),
-  fHistRecPointsQmaxShort(0), fHistRecPointsQmaxMedium(0), 
-  fHistRecPointsQmaxLong(0), fHistRecPointsQShort(0), 
-  fHistRecPointsQMedium(0), fHistRecPointsQLong(0),
-  fHistRecPointsRow(0)
+  fTPCdataQA(NULL)
 {
   // ctor
+  fTPCdataQA = new AliTPCdataQA*[AliRecoParam::kNSpecies] ;
+  for (Int_t specie = 0 ; specie < AliRecoParam::kNSpecies ; specie++) 
+    fTPCdataQA[specie] = NULL ; 
+  
+  for(Int_t i = 0; i < 6; i++)
+    fMapping[i] = 0;
 }
 
 //____________________________________________________________________________ 
 AliTPCQADataMakerRec::AliTPCQADataMakerRec(const AliTPCQADataMakerRec& qadm) :
   AliQADataMakerRec(),
-  fHistESDclusters(0),  //! Clusters per ESD track
-  fHistESDratio(0),     //! Ratio of clusters to findables
-  fHistESDpt(0),        //! Pt spectrum
-  
-  fHistRawsOccupancy(0),//! Pad occupancy (1 entry per pad)
-  
-  fHistRecPointsQmaxShort(0), //! Qmax (short pads)
-  fHistRecPointsQmaxMedium(0),//! Qmax (medium pads)
-  fHistRecPointsQmaxLong(0),  //! Qmax (long pads)
-  fHistRecPointsQShort(0),    //! Q (short pads)
-  fHistRecPointsQMedium(0),   //! Q (medium pads)
-  fHistRecPointsQLong(0),     //! Q (long pads)
-  fHistRecPointsRow(0)       //! Row distribution
+  fTPCdataQA(NULL)
 {
   //copy ctor 
   // Does not copy the calibration object, instead InitRaws have to be
@@ -101,30 +93,18 @@ AliTPCQADataMakerRec::AliTPCQADataMakerRec(const AliTPCQADataMakerRec& qadm) :
   SetName((const char*)qadm.GetName()) ; 
   SetTitle((const char*)qadm.GetTitle()); 
 
+  fTPCdataQA = new AliTPCdataQA*[AliRecoParam::kNSpecies] ;
+  for (Int_t specie = 0 ; specie < AliRecoParam::kNSpecies ; specie++) 
+    fTPCdataQA[specie] = NULL ; 
+  
+  for(Int_t i = 0; i < 6; i++)
+    fMapping[i] = 0;
+
   //
   // Associate class histogram objects to the copies in the list
   // Could also be done with the indexes
   //
-  fHistESDclusters   = (TH1F*)fESDsQAList->FindObject("hESDclusters");
-  fHistESDratio             = (TH1F*)fESDsQAList->FindObject("hESDratio");
-  fHistESDpt         = (TH1F*)fESDsQAList->FindObject("hESDpt");
-
-  fHistRawsOccupancy = (TH1F*)fRawsQAList->FindObject("hRawsOccupancy");
-
-  fHistRecPointsQmaxShort  = 
-    (TH1F*)fRecPointsQAList->FindObject("hRecPointsQmaxShort");
-  fHistRecPointsQmaxMedium = 
-    (TH1F*)fRecPointsQAList->FindObject("hRecPointsQmaxMedium");
-  fHistRecPointsQmaxLong   = 
-    (TH1F*)fRecPointsQAList->FindObject("hRecPointsQmaxLong"); 
-  fHistRecPointsQShort     = 
-    (TH1F*)fRecPointsQAList->FindObject("hRecPointsQShort");
-  fHistRecPointsQMedium    = 
-    (TH1F*)fRecPointsQAList->FindObject("hRecPointsQMedium");
-  fHistRecPointsQLong      = 
-    (TH1F*)fRecPointsQAList->FindObject("hRecPointsQLong");
-  fHistRecPointsRow      = 
-    (TH1F*)fRecPointsQAList->FindObject("hRecPointsRow");
+
 }
 
 //__________________________________________________________________
@@ -135,51 +115,147 @@ AliTPCQADataMakerRec& AliTPCQADataMakerRec::operator = (const AliTPCQADataMakerR
   new(this) AliTPCQADataMakerRec(qadm);
   return *this;
 }
+
+//__________________________________________________________________
+AliTPCQADataMakerRec::~AliTPCQADataMakerRec()
+{
+  // Destructor
+  for (Int_t specie = 0 ; specie < AliRecoParam::kNSpecies ; specie++) 
+    if ( fTPCdataQA[specie] != NULL )
+      delete fTPCdataQA[specie] ; 
+  delete[] fTPCdataQA; 
+
+  for(Int_t i = 0; i < 6; i++) 
+    delete fMapping[i];
+}
  
 //____________________________________________________________________________ 
-void AliTPCQADataMakerRec::EndOfDetectorCycle(AliQA::TASKINDEX task, TObjArray * list)
+void AliTPCQADataMakerRec::EndOfDetectorCycle(AliQAv1::TASKINDEX_t task, TObjArray ** list)
 {
   //Detector specific actions at end of cycle
+  
+  for (Int_t specie = 0 ; specie < AliRecoParam::kNSpecies ; specie++) {
+    if ( !AliQAv1::Instance()->IsEventSpecieSet(specie) ) 
+      continue ; 
+    if(fTPCdataQA[specie] != NULL) { // do the final step of the QA for Raw data
 
-  if(fTPCdataQA) { // do the final step of the QA for Raw data
-
-    fTPCdataQA->Analyse(); // 31/1-08 Analyse is now protected against
+      fTPCdataQA[specie]->Analyse(); // 31/1-08 Analyse is now protected against
                            //         RAW data files with no TPC data
-    
-    //Add2RawsList(fTPCdataQA, 0);
-    // get the histograms and add them to the output
-    // 31/8-08 Histogram is only added if the Calibration class 
-    //         receives TPC data 
-    if(fTPCdataQA->GetNoThreshold()) { 
-      fHistRawsOccupancy = fTPCdataQA->GetNoThreshold()->MakeHisto1D(0, 1, -1);
-      //Add2RawsList(fHistRawsOccupancy, 1);
+
+      SetEventSpecie(AliRecoParam::ConvertIndex(specie)) ; 
+      TH1F * histRawsOccupancy                 = (TH1F*)GetRawsData(kOccupancy) ;
+      TH1F * histRawsOccupancyVsSector         = (TH1F*)GetRawsData(kOccupancyVsSector) ;
+      TH1F * histRawsNClustersPerEventVsSector = (TH1F*)GetRawsData(kNClustersPerEventVsSector) ;
+      TH1F * histRawsQVsSector                 = (TH1F*)GetRawsData(kQVsSector) ;
+      TH1F * histRawsQmaxVsSector              = (TH1F*)GetRawsData(kQmaxVsSector) ;
+      if ( !histRawsOccupancy ||
+          !histRawsOccupancyVsSector ||
+          !histRawsNClustersPerEventVsSector ||
+          !histRawsQVsSector ||
+          !histRawsQmaxVsSector) {
+        AliError("Something very wrong here, corrupted memory ?????. Please check\n") ; 
+        continue ; 
+      }
+        
+      //Add2RawsList(fTPCdataQA, 0);
+      // get the histograms and add them to the output
+      // 31/8-08 Histogram is only added if the Calibration class 
+      //         receives TPC data 
+      const Int_t eventCounter = fTPCdataQA[specie]->GetEventCounter();
+      if(eventCounter>0) { // some TPC data has been processed
+
+        // Reset histograms and refill them 
+        histRawsOccupancy->Reset();
+        histRawsOccupancyVsSector->Reset();
+        histRawsNClustersPerEventVsSector->Reset();
+        histRawsQVsSector->Reset();
+        histRawsQmaxVsSector->Reset();
+      
+        TH1F* hNorm72 = new TH1F("hNorm72", "histogram to normalize 72 sectors",
+                                 72, 0, 72);
+        hNorm72->Sumw2();
+        TH1F* hNorm108 = new TH1F("hNorm108", "histogram to normalize 108 sectors (medium and long pads are split up)",
+                                  108, 0, 108);
+        hNorm108->Sumw2();
+
+        for (Int_t iSec = 0; iSec < 72; iSec++) {
+       
+          AliTPCCalROC* occupancyROC = 
+          fTPCdataQA[specie]->GetNoThreshold()->GetCalROC(iSec); 
+          AliTPCCalROC* nclusterROC = 
+          fTPCdataQA[specie]->GetNLocalMaxima()->GetCalROC(iSec); 
+          AliTPCCalROC* qROC = 
+          fTPCdataQA[specie]->GetMeanCharge()->GetCalROC(iSec); 
+          AliTPCCalROC* qmaxROC = 
+          fTPCdataQA[specie]->GetMaxCharge()->GetCalROC(iSec); 
+
+          const Int_t nRows = occupancyROC->GetNrows(); 
+          for (Int_t iRow = 0; iRow < nRows; iRow++) {
+
+            Int_t helpSector = iSec;
+            if(iRow>=64)
+              helpSector += 36; // OROC (long pads)
+
+            const Int_t nPads = occupancyROC->GetNPads(iRow); 
+            for (Int_t iPad = 0; iPad < nPads; iPad++) {
+       
+              histRawsOccupancy->Fill(occupancyROC->GetValue(iRow, iPad));
+              hNorm72->Fill(iSec);
+              histRawsOccupancyVsSector
+              ->Fill(iSec, occupancyROC->GetValue(iRow, iPad));
+
+              const Int_t nClusters = TMath::Nint(nclusterROC->GetValue(iRow, iPad));
+           
+              if(nClusters>0) {
+             
+                histRawsNClustersPerEventVsSector->Fill(iSec, nClusters);
+                hNorm108->Fill(helpSector, nClusters);
+                histRawsQVsSector->Fill(helpSector, 
+                                       nClusters*qROC->GetValue(iRow, iPad));
+                histRawsQmaxVsSector->Fill(helpSector, 
+                                          nClusters*qmaxROC->GetValue(iRow, iPad));
+              }
+            }
+          }
+        } // end loop over sectors
+      
+        // Normalize histograms
+        histRawsOccupancyVsSector->Divide(hNorm72);
+        histRawsNClustersPerEventVsSector->Scale(1.0/Float_t(eventCounter));
+        histRawsQVsSector->Divide(hNorm108);
+        histRawsQmaxVsSector->Divide(hNorm108);
+        delete hNorm72;
+        delete hNorm108;
+      }
     }
   }
-
-  AliQAChecker::Instance()->Run(AliQA::kTPC, task, list) ;  
+  AliQAChecker::Instance()->Run(AliQAv1::kTPC, task, list) ;  
 }
 
 //____________________________________________________________________________ 
 void AliTPCQADataMakerRec::InitESDs()
 {
   //create ESDs histograms in ESDs subdir  
-  fHistESDclusters = 
+  const Bool_t expert   = kTRUE ; 
+  const Bool_t image    = kTRUE ; 
+  
+  TH1F * histESDclusters = 
     new TH1F("hESDclusters", "N TPC clusters per track; N clusters; Counts",
             160, 0, 160);
-  fHistESDclusters->Sumw2();
-  Add2ESDsList(fHistESDclusters, 0);
+  histESDclusters->Sumw2();
+  Add2ESDsList(histESDclusters, KClusters, !expert, image);
 
-  fHistESDratio = 
+  TH1F * histESDratio = 
     new TH1F("hESDratio", "Ratio: TPC clusters / findable; Ratio: cluster/findable; Counts",
             100, 0, 1);
-  fHistESDratio->Sumw2();
-  Add2ESDsList(fHistESDratio, 1);
+  histESDratio->Sumw2();
+  Add2ESDsList(histESDratio, kRatio, !expert, image);
   
-  fHistESDpt = 
+  TH1F * histESDpt = 
     new TH1F("hESDpt", "P_{T} distribution; p_{T} [GeV/c]; Counts",
             50, 0, 5);
-  fHistESDpt->Sumw2();
-  Add2ESDsList(fHistESDpt, 2);
+  histESDpt->Sumw2();
+  Add2ESDsList(histESDpt, kPt, !expert, image);
 }
 
 //____________________________________________________________________________ 
@@ -187,56 +263,114 @@ void AliTPCQADataMakerRec::InitRaws()
 {
   //
   // Adding the raw 
-  //
-  fTPCdataQA = new AliTPCdataQA();
-  fTPCdataQA->SetRangeTime(0, 999); // take all 1000 time bins 
-  Add2RawsList(fTPCdataQA, 0);
+  //  
+
+  // Modified: 7/7 - 2008
+  // Laurent Aphecetche pointed out that the mapping was read from file
+  // for each event, so now we read in the map here and set if for 
+  // the raw data qa
+  const Bool_t expert   = kTRUE ; 
+  const Bool_t saveCorr = kTRUE ; 
+  const Bool_t image    = kTRUE ; 
+  
+  for (Int_t specie = 0 ; specie < AliRecoParam::kNSpecies ; specie++) {
+    fTPCdataQA[specie] = new AliTPCdataQA(AliRecoParam::ConvertIndex(specie));
+    LoadMaps(); // Load Altro maps
+    fTPCdataQA[specie]->SetAltroMapping(fMapping); // set Altro mapping
+    fTPCdataQA[specie]->SetRangeTime(100, 920); // set time bin interval 
+//    Add2RawsList(fTPCdataQA, kTPCdataQ, !expert, image, !saveCorrA); // This is used by the AMORE monitoring <------- THIS WILL FAIL (YS)
+  }
+
+  TH1F * histRawsOccupancy = 
+    new TH1F("hRawsOccupancy", "Occupancy (all pads); Occupancy; Counts",
+            100, 0, 1);
+  histRawsOccupancy->Sumw2();
+  Add2RawsList(histRawsOccupancy, kOccupancy, !expert, image, !saveCorr);
+  
+  TH1F * histRawsOccupancyVsSector = 
+    new TH1F("hRawsOccupancyVsSector", "Occupancy vs sector; Sector; Occupancy",
+            72, 0, 72);
+  histRawsOccupancyVsSector->Sumw2();
+  Add2RawsList(histRawsOccupancyVsSector, kOccupancyVsSector, !expert, image, !saveCorr);
+
+  TH1F * histRawsNClustersPerEventVsSector = 
+    new TH1F("hRawsNClustersPerEventVsSector", "Nclusters per event vs sector; Sector; Nclusters per event",
+            72, 0, 72);
+  histRawsNClustersPerEventVsSector->Sumw2();
+  Add2RawsList(histRawsNClustersPerEventVsSector, kNClustersPerEventVsSector, !expert, image, !saveCorr);
+  
+  TH1F * histRawsQVsSector = 
+    new TH1F("hRawsQVsSector", "<Q> vs sector (OROC med: 36-71, long: 72-107); Sector; <Q>",
+            108, 0, 108);
+  histRawsQVsSector->Sumw2();
+  Add2RawsList(histRawsQVsSector, kQVsSector, !expert, image, !saveCorr);
+
+  TH1F * histRawsQmaxVsSector = 
+    new TH1F("hRawsQmaxVsSector", "<Qmax> vs sector (OROC med: 36-71, long: 72-107); Sector; <Qmax>",
+            108, 0, 108);
+  histRawsQmaxVsSector->Sumw2();
+  Add2RawsList(histRawsQmaxVsSector, kQmaxVsSector, !expert, image, !saveCorr);
+}
+
+//____________________________________________________________________________ 
+void AliTPCQADataMakerRec::InitDigits()
+{
+  const Bool_t expert   = kTRUE ; 
+  const Bool_t image    = kTRUE ; 
+  TH1F * histDigitsADC = 
+    new TH1F("hDigitsADC", "Digit ADC distribution; ADC; Counts",
+             1000, 0, 1000);
+  histDigitsADC->Sumw2();
+  Add2DigitsList(histDigitsADC, kDigitsADC, !expert, image);
 }
 
 //____________________________________________________________________________ 
 void AliTPCQADataMakerRec::InitRecPoints()
 {
-  fHistRecPointsQmaxShort = 
+  const Bool_t expert   = kTRUE ; 
+  const Bool_t image    = kTRUE ; 
+  
+  TH1F * histRecPointsQmaxShort = 
     new TH1F("hRecPointsQmaxShort", "Qmax distrbution (short pads); Qmax; Counts",
-            200, 0, 1000);
-  fHistRecPointsQmaxShort->Sumw2();
-  Add2RecPointsList(fHistRecPointsQmaxShort, 0);
+            100, 0, 300);
+  histRecPointsQmaxShort->Sumw2();
+  Add2RecPointsList(histRecPointsQmaxShort, kQmaxShort, !expert, image);
 
-  fHistRecPointsQmaxMedium = 
+  TH1F * histRecPointsQmaxMedium = 
     new TH1F("hRecPointsQmaxMedium", "Qmax distrbution (medium pads); Qmax; Counts",
-            200, 0, 1000);
-  fHistRecPointsQmaxMedium->Sumw2();
-  Add2RecPointsList(fHistRecPointsQmaxMedium, 1);
+            100, 0, 300);
+  histRecPointsQmaxMedium->Sumw2();
+  Add2RecPointsList(histRecPointsQmaxMedium, kQmaxMedium, !expert, image);
 
-  fHistRecPointsQmaxLong = 
+  TH1F * histRecPointsQmaxLong = 
     new TH1F("hRecPointsQmaxLong", "Qmax distrbution (long pads); Qmax; Counts",
-            200, 0, 1000);
-  fHistRecPointsQmaxLong->Sumw2();
-  Add2RecPointsList(fHistRecPointsQmaxLong, 2);
+            100, 0, 300);
+  histRecPointsQmaxLong->Sumw2();
+  Add2RecPointsList(histRecPointsQmaxLong, kQmaxLong, !expert, image);
 
-  fHistRecPointsQShort = 
+  TH1F * histRecPointsQShort = 
     new TH1F("hRecPointsQShort", "Q distrbution (short pads); Q; Counts",
-            200, 0, 5000);
-  fHistRecPointsQShort->Sumw2();
-  Add2RecPointsList(fHistRecPointsQShort, 3);
+            100, 0, 2000);
+  histRecPointsQShort->Sumw2();
+  Add2RecPointsList(histRecPointsQShort, kQShort, !expert, image);
 
-  fHistRecPointsQMedium = 
+  TH1F * histRecPointsQMedium = 
     new TH1F("hRecPointsQMedium", "Q distrbution (medium pads); Q; Counts",
-            200, 0, 5000);
-  fHistRecPointsQMedium->Sumw2();
-  Add2RecPointsList(fHistRecPointsQMedium, 4);
+            100, 0, 2000);
+  histRecPointsQMedium->Sumw2();
+  Add2RecPointsList(histRecPointsQMedium, kQMedium, !expert, image);
 
-  fHistRecPointsQLong = 
+  TH1F * histRecPointsQLong = 
     new TH1F("hRecPointsQLong", "Q distrbution (long pads); Q; Counts",
-            200, 0, 5000);
-  fHistRecPointsQLong->Sumw2();
-  Add2RecPointsList(fHistRecPointsQLong, 5);
+            100, 0, 2000);
+  histRecPointsQLong->Sumw2();
+  Add2RecPointsList(histRecPointsQLong, kQLong, !expert, image);
 
-  fHistRecPointsRow = 
+  TH1F * histRecPointsRow = 
     new TH1F("hRecPointsRow", "Clusters per row; Row; Counts",
             159, 0, 159);
-  fHistRecPointsRow->Sumw2();
-  Add2RecPointsList(fHistRecPointsRow, 6);
+  histRecPointsRow->Sumw2();
+  Add2RecPointsList(histRecPointsRow, kRow, !expert, image);
 }
 
 //____________________________________________________________________________
@@ -257,10 +391,10 @@ void AliTPCQADataMakerRec::MakeESDs(AliESDEvent * esd)
     
     Int_t nTPCclusters         = track->GetTPCNcls();
     Int_t nTPCclustersFindable = track->GetTPCNclsF();
-    
-    fHistESDclusters->Fill(nTPCclusters);
-    fHistESDratio->Fill(Float_t(nTPCclusters)/Float_t(nTPCclustersFindable));
-    fHistESDpt->Fill(track->Pt()); 
+    if ( nTPCclustersFindable<=0) continue;
+    GetESDsData(KClusters)->Fill(nTPCclusters);
+    GetESDsData(kRatio)->Fill(Float_t(nTPCclusters)/Float_t(nTPCclustersFindable));
+    GetESDsData(kPt)->Fill(track->Pt()); 
   }
 }
 
@@ -271,7 +405,30 @@ void AliTPCQADataMakerRec::MakeRaws(AliRawReader* rawReader)
   // To make QA for the RAW data we use the TPC Calibration framework 
   // to handle the data and then in the end extract the data
   //
-  fTPCdataQA->ProcessEvent(rawReader);  
+       rawReader->Reset() ; 
+  fTPCdataQA[AliRecoParam::AConvert(fEventSpecie)]->ProcessEvent(rawReader);  
+ }
+
+//____________________________________________________________________________
+void AliTPCQADataMakerRec::MakeDigits(TTree* digitTree)
+{
+  TBranch* branch = digitTree->GetBranch("Segment");
+  AliSimDigits* digArray = 0;
+  branch->SetAddress(&digArray);
+  
+  Int_t nEntries = Int_t(digitTree->GetEntries());
+  
+  for (Int_t n = 0; n < nEntries; n++) {
+    
+    digitTree->GetEvent(n);
+    
+    if (digArray->First())
+      do {
+        Float_t dig = digArray->CurrentDigit();
+        
+        GetDigitsData(kDigitsADC)->Fill(dig);
+      } while (digArray->Next());    
+  }
 }
 
 //____________________________________________________________________________
@@ -302,24 +459,39 @@ void AliTPCQADataMakerRec::MakeRecPoints(TTree* recTree)
 
       if(cluster->GetDetector()<36) { // IROC (short pads)
 
-       fHistRecPointsQmaxShort->Fill(Qmax);
-       fHistRecPointsQShort->Fill(Q);
+       GetRecPointsData(kQmaxShort)->Fill(Qmax);
+       GetRecPointsData(kQShort)->Fill(Q);
       } else { // OROC (medium and long pads)
        row += 63;
        if(cluster->GetRow()<64) { // medium pads
 
-         fHistRecPointsQmaxMedium->Fill(Qmax);
-         fHistRecPointsQMedium->Fill(Q);
+         GetRecPointsData(kQmaxMedium)->Fill(Qmax);
+         GetRecPointsData(kQMedium)->Fill(Q);
        } else { // long pads
 
-         fHistRecPointsQmaxLong->Fill(Qmax);
-         fHistRecPointsQLong->Fill(Q);
+         GetRecPointsData(kQmaxLong)->Fill(Qmax);
+         GetRecPointsData(kQLong)->Fill(Q);
        }
       }
       
-      fHistRecPointsRow->Fill(row);
+      GetRecPointsData(kRow)->Fill(row);
     } // end loop over clusters
   } // end loop over tree
 
   delete clrow;
 }
+
+//____________________________________________________________________________
+void AliTPCQADataMakerRec::LoadMaps()
+{
+  TString path = gSystem->Getenv("ALICE_ROOT");
+  path += "/TPC/mapping/Patch";
+
+  for(Int_t i = 0; i < 6; i++) {
+    TString path2 = path;
+    path2 += i;
+    path2 += ".data";
+    fMapping[i] = new AliTPCAltroMapping(path2.Data());
+  }
+}
+