Add new global configuration for prod shuttle.
[u/mrichter/AliRoot.git] / 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 //====================================================================================================================================================