new class to hold event wide info relevant for TRD (Markus)
[u/mrichter/AliRoot.git] / TRD / qaRec / AliTRDcheckDetector.cxx
1 #include <TFile.h>
2 #include <TH1F.h>
3 #include <TGraph.h>
4 #include <TMath.h>
5 #include <TMap.h>
6 #include <TObjArray.h>
7 #include <TObject.h>
8 #include <TObjString.h>
9 #include <TProfile.h>
10 #include <TROOT.h>
11
12 #include "AliTRDcluster.h"
13 #include "AliESDHeader.h"
14 #include "AliESDRun.h"
15 #include "AliTRDseedV1.h"
16 #include "AliTRDtrackV1.h"
17 #include "TTreeStream.h"
18
19 #include "AliTRDtrackInfo/AliTRDtrackInfo.h"
20 #include "AliTRDtrackInfo/AliTRDeventInfo.h"
21 #include "AliTRDcheckDetector.h"
22
23 #include <cstdio>
24
25 ////////////////////////////////////////////////////////////////////////////
26 //                                                                        //
27 //  Reconstruction QA                                                     //
28 //                                                                        //
29 //  Task doing basic checks for tracking and detector performance         //
30 //                                                                        //
31 //  Authors:                                                              //
32 //    Anton Andronic <A.Andronic@gsi.de>                                  //
33 //    Markus Fasel <M.Fasel@gsi.de>                                       //
34 //                                                                        //
35 ////////////////////////////////////////////////////////////////////////////
36
37 //_______________________________________________________
38 AliTRDcheckDetector::AliTRDcheckDetector():
39   AliTRDrecoTask("DetChecker", "Basic Detector Checker")
40   ,fEventInfo(0x0)
41   ,fPHSdetector(0x0)
42   ,fPHSsector(0x0)
43   ,fQCLdetector(0x0)
44   ,fQCLsector(0x0)
45   ,fQTdetector(0x0)
46   ,fQTsector(0x0)
47   ,fTriggerNames(0x0)
48 {
49   //
50   // Default constructor
51   //
52   DefineInput(1,AliTRDeventInfo::Class());
53 }
54
55 //_______________________________________________________
56 AliTRDcheckDetector::~AliTRDcheckDetector(){
57   //
58   // Destructor
59   // 
60   if(fEventInfo) delete fEventInfo;
61   if(fPHSdetector) delete fPHSdetector;
62   if(fPHSsector) delete fPHSsector;
63   if(fQCLdetector) delete fQCLdetector;
64   if(fQCLsector) delete fQCLsector;
65   if(fQTdetector) delete fQTdetector;
66   if(fQTsector) delete fQTsector;
67   if(fTriggerNames) delete fTriggerNames;
68 }
69
70 //_______________________________________________________
71 void AliTRDcheckDetector::ConnectInputData(Option_t *opt){
72         //
73         // Connect the Input data with the task
74         //
75         AliTRDrecoTask::ConnectInputData(opt);
76         fEventInfo = dynamic_cast<AliTRDeventInfo *>(GetInputData(1));
77 }
78
79 //_______________________________________________________
80 void AliTRDcheckDetector::CreateOutputObjects(){
81   //
82   // Create Output Objects
83   //
84   fContainer = new TObjArray();
85   // Register Histograms
86   fContainer->Add(new TH1F("hNtrks", "Number of Tracks per event", 100, 0, 100));
87   fContainer->Add(new TH1F("hEventsTriggerTracks", "Trigger Class (Tracks)", 20, 0, 20));
88   fContainer->Add(new TH1F("hNcls", "Nr. of clusters per track", 30, 0, 180));
89   fContainer->Add(new TH1F("hNtls", "Nr. tracklets per track", 10, 0, 10));
90   fContainer->Add(new TH1F("hNclTls","Mean Number of clusters per tracklet", 30, 0, 30));
91   fContainer->Add(new TH1F("hChi2", "Chi2", 200, 0, 20));
92   fContainer->Add(new TH1F("hChi2n", "Norm. Chi2 (tracklets)", 50, 0, 5));
93   fContainer->Add(new TH1F("hSM", "Track Counts in Supermodule", 18, 0., 18.));
94   fContainer->Add(new TProfile("hPH", "Average PH", 31, -0.5, 30.5));
95   fContainer->Add(new TH1F("hQcl1", "Cluster charge", 200, -1200, 1200));
96   fContainer->Add(new TH1F("hQT", "Total Charge Deposit", 6000, 0, 6000));
97   fContainer->Add(new TH1F("hEventsTrigger", "Trigger Class", 20, 0, 20));
98   // Initialise the PHS, cluster charge and total charge deposit histograms for each detector
99   fPHSdetector = new TObjArray();
100   fPHSdetector->SetName("fPHSdetector");
101   fQCLdetector = new TObjArray();
102   fQCLdetector->SetName("fQCLdetector");
103   fQTdetector = new TObjArray();
104   fQTdetector->SetName("fQTdetector");
105   for(Int_t idet = 0; idet < kNDetectors; idet++){
106     fPHSdetector->Add(new TProfile(Form("hPHSdet%d", idet), Form("Average Pulse Height in Detector %d", idet), 31, -0.5, 30.5));
107     fQCLdetector->Add(new TH1F(Form("hQd%d",idet),Form("Qd%d",idet), 100,0,5000));
108     fQTdetector->Add(new TH1F(Form("hQTd%d", idet), Form("Qtotal%d", idet), 6000, 0, 6000));
109   }
110   fContainer->Add(fPHSdetector);
111   fContainer->Add(fQCLdetector);
112   fContainer->Add(fQTdetector);
113   
114   // Initialise the PHS, cluster charge and total charge deposit histograms for each sector
115   fPHSsector = new TObjArray();
116   fPHSsector->SetName("fPHSsector");
117   fQCLsector = new TObjArray();
118   fQCLsector->SetName("fQCLsector");
119   fQTsector = new TObjArray();
120   fQTsector->SetName("fQTsector");
121   for(Int_t isec = 0; isec < kNSectors; isec++){
122     fPHSsector->Add(new TProfile(Form("hPHSsec%d", isec), Form("Average Pulse Height in Sector %d", isec), 31, -0.5, 30.5));
123     fQCLsector->Add(new TH1F(Form("hQd%d",isec),Form("Qd%d",isec), 100,0,5000));
124     fQTsector->Add(new TH1F(Form("hQTd%d", isec), Form("Qtotal%d", isec), 6000, 0, 6000));
125   }
126   fContainer->Add(fPHSsector);
127   fContainer->Add(fQCLsector);
128   fContainer->Add(fQTsector);
129   fTriggerNames = new TMap();
130   fTriggerNames->SetName("fTriggerNames");
131   fContainer->Add(fTriggerNames);
132 }
133
134 //_______________________________________________________
135 void AliTRDcheckDetector::Exec(Option_t *){
136   //
137   // Execution function
138   // Filling TRD quality histos
139   //
140   if(fEventInfo->GetEventHeader()->GetEventType() != 7) return;
141   Int_t nTracks = 0;            // Count the number of tracks per event
142   Int_t triggermask = fEventInfo->GetEventHeader()->GetTriggerMask();
143   TString triggername =  fEventInfo->GetRunInfo()->GetTriggerClass(triggermask);
144   dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNEventsTrigger))->Fill(triggermask);
145   AliTRDtrackInfo *fTrackInfo = 0x0;
146   AliTRDtrackV1 *fTRDtrack = 0x0;
147   AliTRDseedV1 *fTracklet = 0x0;
148   AliTRDcluster *fTRDcluster = 0x0;
149   for(Int_t iti = 0; iti < fTracks->GetEntriesFast(); iti++){
150     fTrackInfo = dynamic_cast<AliTRDtrackInfo *>(fTracks->UncheckedAt(iti));
151     if(!fTrackInfo || !(fTRDtrack = fTrackInfo->GetTRDtrack())) continue;
152     Int_t nclusters = fTRDtrack->GetNumberOfClusters();
153     Int_t ntracklets = fTRDtrack->GetNumberOfTracklets();
154     Float_t chi2 = fTRDtrack->GetChi2();
155     // Fill Histograms
156     dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNclustersHist))->Fill(nclusters);
157     dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNtrackletsHist))->Fill(ntracklets);
158     dynamic_cast<TH1F *>(fContainer->UncheckedAt(kChi2))->Fill(chi2);
159     dynamic_cast<TH1F *>(fContainer->UncheckedAt(kChi2Normalized))->Fill(chi2/static_cast<Float_t>(ntracklets));
160     // now loop over single tracklets
161     if(ntracklets > 2){
162       Int_t sector = -1;
163       for(Int_t ilayer = 0; ilayer < kNLayers; ilayer++){
164         if(!(fTracklet = fTRDtrack->GetTracklet(ilayer)) || !fTracklet->IsOK()) continue;
165         Float_t nClustersTracklet = fTracklet->GetN();
166         dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNclusterTrackletHist))->Fill(nClustersTracklet);
167         if(nClustersTracklet){
168           Float_t Qtot = 0;
169           Int_t detector = -1;
170           for(Int_t itb = 0; itb < kNTimebins; itb++){
171             if(!(fTRDcluster = fTracklet->GetClusters(itb))) continue;
172             Int_t localtime        = fTRDcluster->GetLocalTimeBin();
173             Double_t clusterCharge = fTRDcluster->GetQ();
174             detector               = fTRDcluster->GetDetector();
175             sector                 = static_cast<Int_t>(detector/kNDetectorsSector);
176             Double_t absolute_charge = TMath::Abs(clusterCharge);
177             Qtot += absolute_charge;
178             // Fill Histograms
179             // 1st. PHS
180             dynamic_cast<TProfile *>(fContainer->UncheckedAt(kPulseHeight))->Fill(localtime, absolute_charge);
181             dynamic_cast<TProfile *>(fPHSdetector->UncheckedAt(detector))->Fill(localtime, absolute_charge);
182             dynamic_cast<TProfile *>(fPHSsector->UncheckedAt(sector))->Fill(localtime, absolute_charge);
183             // 2nd. cluster charge
184             dynamic_cast<TH1F *>(fContainer->UncheckedAt(kClusterCharge))->Fill(absolute_charge);
185             dynamic_cast<TH1F *>(fQCLdetector->UncheckedAt(detector))->Fill(absolute_charge);
186             dynamic_cast<TH1F *>(fQCLsector->UncheckedAt(sector))->Fill(absolute_charge);
187           }
188           // Fill total cluster charge histograms
189           dynamic_cast<TH1F *>(fContainer->UncheckedAt(kChargeDeposit))->Fill(Qtot);
190           dynamic_cast<TH1F *>(fQTdetector->UncheckedAt(detector))->Fill(Qtot);
191           dynamic_cast<TH1F *>(fQTsector->UncheckedAt(sector))->Fill(Qtot);
192         }
193       }
194       dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNTracksSectorHist))->Fill(sector);
195     }
196     nTracks++;
197   }
198   if(nTracks) dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNEventsTriggerTracks))->Fill(triggermask);
199   dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNTracksEventHist))->Fill(nTracks);
200   if(!fTriggerNames->FindObject(Form("%d", triggermask)))
201     fTriggerNames->Add(new TObjString(Form("%d", triggermask)), new TObjString(triggername));
202   PostData(0, fContainer);
203 }
204
205 //_______________________________________________________
206 void AliTRDcheckDetector::Terminate(Option_t *){
207   //
208   // Terminate function
209   //
210 }
211
212 //_______________________________________________________
213 Bool_t AliTRDcheckDetector::PostProcess(){
214         //
215         // Do Postprocessing (for the moment set the number of Reference histograms)
216         //
217         
218         TH1 * histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNTracksEventHist));
219         histo->GetXaxis()->SetTitle("Number of Tracks");
220         histo->GetYaxis()->SetTitle("Events");
221         // Calculate percentage of events containing tracks
222         Int_t entries = (Int_t)histo->GetEntries();
223   Float_t percTracks = entries > 0 ? (entries - histo->GetBinContent(1))/Float_t(entries) : 0.;
224         printf("Percentage of Events containing TRD Tracks: %6.3f %%\n", 1.E2*percTracks);
225         histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNclustersHist));
226         histo->GetXaxis()->SetTitle("Number of Clusters");
227         histo->GetYaxis()->SetTitle("Entries");
228         histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNtrackletsHist));
229         histo->GetXaxis()->SetTitle("Number of Tracklets");
230         histo->GetYaxis()->SetTitle("Entries");
231         histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNclusterTrackletHist));
232         histo->GetXaxis()->SetTitle("Number of Clusters");
233         histo->GetYaxis()->SetTitle("Entries");
234         histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kChi2));
235         histo->GetXaxis()->SetTitle("#chi^2");
236         histo->GetYaxis()->SetTitle("Entries");
237         histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNTracksSectorHist));
238         histo->GetXaxis()->SetTitle("Sector");
239         histo->GetYaxis()->SetTitle("Number of Tracks");
240         histo = dynamic_cast<TProfile *>(fContainer->UncheckedAt(kPulseHeight));
241         histo->GetXaxis()->SetTitle("Time / 100ns");
242         histo->GetYaxis()->SetTitle("Average Pulse Height (a. u.)");
243         histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kClusterCharge));
244         histo->GetXaxis()->SetTitle("Cluster Charge (a.u.)");
245         histo->GetYaxis()->SetTitle("Entries");
246         histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kChargeDeposit));
247         histo->GetXaxis()->SetTitle("Charge Deposit (a.u.)");
248         histo->GetYaxis()->SetTitle("Entries");
249         // Calculate the Percentage of events containing tracks as function of the trigger cluster
250         histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNEventsTrigger));
251         TH1F *histoTracks = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNEventsTriggerTracks));
252         histoTracks->Divide(histo);
253         fTriggerNames = dynamic_cast<TMap *>(fContainer->UncheckedAt(17));
254         TH1F *percentages = new TH1F("fTriggerInf", "Trigger Information", fTriggerNames->GetEntries(), 0, fTriggerNames->GetEntries());
255         TObject *triggerclass = 0x0;
256         Int_t ipt=0;
257         for(Int_t ibin = 0; ibin < histo->GetNbinsX(); ibin++){
258                 if((triggerclass = fTriggerNames->FindObject(Form("%d",ibin)))){
259                         TObjString *label = dynamic_cast<TObjString *>(dynamic_cast<TPair *>(triggerclass)->Value());
260                         //printf("Trigger Pattern for class %d: %s\n", ibin, label->String().Data());
261                         percentages->SetBinContent(percentages->FindBin(ipt), histoTracks->GetBinContent(histoTracks->FindBin(ibin)));
262                         percentages->GetXaxis()->SetBinLabel(percentages->FindBin(ipt),label->String().Data());
263                         ipt++;
264                 }       
265         }
266         percentages->GetXaxis()->SetTitle("Trigger Cluster");
267         percentages->GetYaxis()->SetTitle("%Events");
268         percentages->GetYaxis()->SetRangeUser(0,1);
269 //      percentages->SetMarkerColor(kBlue);
270 //      percentages->SetMarkerStyle(22);
271         fContainer->Add(percentages);
272         fNRefFigures = 10;
273         return kTRUE;
274 }
275
276 //_______________________________________________________
277 void AliTRDcheckDetector::GetRefFigure(Int_t ifig, Int_t &first, Int_t &last, Option_t *opt){
278         //
279         // Setting Reference Figures
280         //
281         opt = "pl";
282         switch(ifig){
283                 case 0: first = last = kNTracksEventHist;
284                                                 break;
285                 case 1: first = last = kNclustersHist;
286                                                 break;
287                 case 2: first = last = kNtrackletsHist;
288                                                 break;
289                 case 3: first = last = kNclusterTrackletHist;
290                                                 break;
291                 case 4: first = last = kChi2;
292                                                 break;
293                 case 5: first = last = kNTracksSectorHist;
294                                                 break;
295                 case 6: first = last = kPulseHeight;
296                                                 break;
297                 case 7: first = last = kClusterCharge;
298                                                 break;
299                 case 8: first = last = kChargeDeposit;
300                                                 break;
301                 case 9: first = last = kPostProcessing;
302                                                 opt="bar";
303                                                 break;
304                 default: first = last = kNTracksEventHist;
305                                                 break;
306         };
307 }
308