MFT/AliMFTClusterQA.cxx

   1 #include "TObject.h"
   2 #include "AliRunLoader.h"
   3 #include "AliRun.h"
   4 #include "AliLoader.h"
   5 #include "AliMFT.h"
   6 #include "TClonesArray.h"
   7 #include "AliMFTCluster.h"
   8 #include "AliMFTSegmentation.h"
   9 #include "TFile.h"
  10 #include "TH1D.h"
  11 #include "TH2D.h"
  12 #include "AliLog.h"
  13 #include "TString.h"
  14
  15 #include "AliMFTClusterQA.h"
  16
  17 //====================================================================================================================================================
  18 //
  19 // Class for the analysis of the MFT clusters (a.k.a. rec points). Few QA histograms are created
  20 //
  21 // Contact author: antonio.uras@cern.ch
  22 //
  23 //====================================================================================================================================================
  24
  25 ClassImp(AliMFTClusterQA)
  26
  27 //====================================================================================================================================================
  28
  29 AliMFTClusterQA::AliMFTClusterQA():
  30   TObject(),
  31   fMFTLoader(0),
  32   fRunLoader(0),
  33   fMFT(0),
  34   fNPlanes(0),
  35   fNEvents(0),
  36   fEv(0),
  37   fFileOut(0),
  38   fReadDir(0),
  39   fOutDir(0)
  40 {
  41
  42   // default constructor
  43
  44   for (Int_t iPlane=0; iPlane<fNMaxPlanes; iPlane++) {
  45     fHistNClustersPerEvent[iPlane]     = 0;
  46     fHistNPixelsPerCluster[iPlane]     = 0;
  47     fHistClusterSizeX[iPlane]          = 0;
  48     fHistClusterSizeY[iPlane]          = 0;
  49     fHistClusterRadialPosition[iPlane] = 0;
  50     fClusterScatterPlotXY[iPlane]      = 0;
  51   }
  52
  53 }
  54
  55 //====================================================================================================================================================
  56
  57 void AliMFTClusterQA::Init(Char_t *readDir, Char_t *outDir, Int_t nEventsToAnalyze) {
  58
  59   fReadDir = readDir;
  60   fOutDir  = outDir;
  61
  62   fRunLoader = AliRunLoader::Open(Form("%s/galice.root", fReadDir.Data()));
  63   gAlice = fRunLoader->GetAliRun();
  64   if (!gAlice) fRunLoader->LoadgAlice();
  65   fNEvents = fRunLoader->GetNumberOfEvents();
  66   if (nEventsToAnalyze>0) fNEvents = TMath::Min(fNEvents, nEventsToAnalyze);
  67
  68   fMFT = (AliMFT*) gAlice->GetDetector("MFT");
  69   fNPlanes = fMFT->GetSegmentation()->GetNPlanes();
  70
  71   BookHistos();
  72
  73   fMFTLoader = fRunLoader->GetDetectorLoader("MFT");
  74   fMFTLoader -> LoadRecPoints("READ");
  75
  76 }
  77
  78 //====================================================================================================================================================
  79
  80 Bool_t AliMFTClusterQA::LoadNextEvent() {
  81
  82   if (fEv>=fNEvents) return kFALSE;
  83   AliDebug(1, Form("event %5d",fEv));
  84
  85   fRunLoader->GetEvent(fEv);
  86   fEv++;
  87
  88   if (!fMFTLoader->TreeR()->GetEvent()) return kTRUE;
  89
  90   for (Int_t iPlane=0; iPlane<fNPlanes; iPlane++) {
  91     Int_t nClusters = fMFT->GetRecPointsList(iPlane)->GetEntries();
  92     fHistNClustersPerEvent[iPlane] -> Fill(nClusters);
  93     fClusterScatterPlotXY[iPlane]  -> Fill(0., 0.);    // "scaler" bin
  94     AliDebug(1,Form("nClusters = %5d", nClusters));
  95     for (Int_t iCluster=0; iCluster<nClusters; iCluster++) {
  96       AliMFTCluster *cluster = (AliMFTCluster*) (fMFT->GetRecPointsList(iPlane))->At(iCluster);
  97       fHistNPixelsPerCluster[iPlane] -> Fill(cluster->GetSize());
  98       fHistClusterSizeX[iPlane]      -> Fill(cluster->GetErrX()*1.e4);   // converted in microns
  99       fHistClusterSizeY[iPlane]      -> Fill(cluster->GetErrY()*1.e4);   // converted in microns
 100       fHistClusterRadialPosition[iPlane] -> Fill(TMath::Sqrt(cluster->GetX()*cluster->GetX()+cluster->GetY()*cluster->GetY()));
 101       fClusterScatterPlotXY[iPlane]      -> Fill(cluster->GetX(), cluster->GetY());
 102     }
 103   }
 104
 105   return kTRUE;
 106
 107 }
 108
 109 //====================================================================================================================================================
 110
 111 void AliMFTClusterQA::BookHistos() {
 112
 113   fFileOut = new TFile(Form("%s/MFT.RecPoints.QA.root",fOutDir.Data()), "recreate");
 114
 115   for (Int_t iPlane=0; iPlane<fNPlanes; iPlane++) {
 116
 117     fHistNClustersPerEvent[iPlane] = new TH1D(Form("fHistNClustersPerEvent_Pl%02d",iPlane),
 118                                               Form("Number of clusters per event in Plane%02d",iPlane),
 119                                               25000, -0.5, 24999.5);
 120
 121     fHistNPixelsPerCluster[iPlane] = new TH1D(Form("fHistNPixelsPerCluster_Pl%02d",iPlane),
 122                                               Form("Number of pixels per cluster in Plane%02d",iPlane),
 123                                               15, -0.5, 14.5);
 124
 125     fHistClusterSizeX[iPlane]      = new TH1D(Form("fHistClusterSizeX_Pl%02d",iPlane),
 126                                               Form("#Deltax for clusters in Plane%02d",iPlane),
 127                                               100, 0., 100.);
 128
 129     fHistClusterSizeY[iPlane]      = new TH1D(Form("fHistClusterSizeY_Pl%02d",iPlane),
 130                                               Form("#Deltay for clusters in Plane%02d",iPlane),
 131                                               100, 0., 100.);
 132
 133     fHistNClustersPerEvent[iPlane] -> SetXTitle("N_{clusters} per Event");
 134     fHistNPixelsPerCluster[iPlane] -> SetXTitle("N_{pixels} per Cluster");
 135     fHistClusterSizeX[iPlane]      -> SetXTitle("#Deltax  [#mum]");
 136     fHistClusterSizeY[iPlane]      -> SetXTitle("#Deltay  [#mum]");
 137
 138     fHistNClustersPerEvent[iPlane] -> Sumw2();
 139     fHistNPixelsPerCluster[iPlane] -> Sumw2();
 140     fHistClusterSizeX[iPlane]      -> Sumw2();
 141     fHistClusterSizeY[iPlane]      -> Sumw2();
 142
 143     //------------------------------------------------------------
 144
 145     Int_t rMax = Int_t(10.*(fMFT->GetSegmentation()->GetPlane(iPlane)->GetRMaxSupport()));  // expressed in mm
 146
 147     fHistClusterRadialPosition[iPlane] = new TH1D(Form("fHistClusterRadialPosition_Pl%02d",iPlane),
 148                                                   Form("Cluster radial position (Plane%02d)",iPlane),
 149                                                   rMax, 0, Double_t(rMax)/10.);
 150
 151     fClusterScatterPlotXY[iPlane] = new TH2D(Form("fClusterScatterPlotXY_Pl%02d",iPlane),
 152                                              Form("Cluster scatter plot (Plane%02d)",iPlane),
 153                                              2*rMax+1, (-rMax-0.5)/10., (rMax+0.5)/10., 2*rMax+1, (-rMax-0.5)/10., (rMax+0.5)/10.);
 154
 155     fHistClusterRadialPosition[iPlane] -> SetXTitle("R  [cm]");
 156     fClusterScatterPlotXY[iPlane]      -> SetXTitle("X  [cm]");
 157     fClusterScatterPlotXY[iPlane]      -> SetYTitle("Y  [cm]");
 158
 159     fHistClusterRadialPosition[iPlane] -> Sumw2();
 160     fClusterScatterPlotXY[iPlane]      -> Sumw2();
 161
 162   }
 163
 164 }
 165
 166 //====================================================================================================================================================
 167
 168 void AliMFTClusterQA::Terminate() {
 169
 170   AliInfo("Writing QA histos...");
 171
 172   // ---- equalize radial clusters distribution ----------------------
 173
 174   for (Int_t iPlane=0; iPlane<fNPlanes; iPlane++) {
 175     for (Int_t iBin=1; iBin<=fHistClusterRadialPosition[iPlane]->GetNbinsX(); iBin++) {
 176       Double_t rMin = fHistClusterRadialPosition[iPlane]->GetBinLowEdge(iBin);        // in cm
 177       Double_t rMax = fHistClusterRadialPosition[iPlane]->GetBinWidth(iBin) + rMin;   // in cm
 178       Double_t area = 100.*TMath::Pi()*(rMax*rMax - rMin*rMin);                       // in mm^2
 179       Double_t density = fHistClusterRadialPosition[iPlane]->GetBinContent(iBin)/area;
 180       fHistClusterRadialPosition[iPlane]->SetBinContent(iBin, density);
 181       fHistClusterRadialPosition[iPlane]->SetBinError(iBin, fHistClusterRadialPosition[iPlane]->GetBinError(iBin)/area);
 182     }
 183     fHistClusterRadialPosition[iPlane] -> SetBinContent(1, fEv);      // "scaler" bin
 184   }
 185
 186   // -----------------------------------------------------------------
 187
 188   fFileOut->cd();
 189
 190   for (Int_t iPlane=0; iPlane<fNPlanes; iPlane++) {
 191     fHistNClustersPerEvent[iPlane]     -> Write();
 192     fHistNPixelsPerCluster[iPlane]     -> Write();
 193     fHistClusterSizeX[iPlane]          -> Write();
 194     fHistClusterSizeY[iPlane]          -> Write();
 195     fHistClusterRadialPosition[iPlane] -> Write();
 196     fClusterScatterPlotXY[iPlane]      -> Write();
 197   }
 198
 199   fFileOut -> Close();
 200
 201 }
 202
 203 //====================================================================================================================================================