- moved AliHLTDisplay files from TPCLib to TPCLib/OnlineDisplay
[u/mrichter/AliRoot.git] / HLT / TPCLib / OnlineDisplay / AliHLTTPCDisplayResiduals.cxx
1 // $Id$
2
3 /** \class AliHLTTPCDisplayPadRow
4 <pre>
5 //_____________________________________________________________
6 // AliHLTTPCDisplayResiduals
7 //
8 // Display class for the HLT TPC-Residuals events.
9 </pre>
10 */
11 // Author: Jochen Thaeder <mailto:thaeder@kip.uni-heidelberg.de>
12 //*-- Copyright &copy ALICE HLT Group 
13
14 #include "AliHLTTPCDisplayResiduals.h"
15 #include "AliHLTTPCDisplayPadRow.h"
16
17 #include "AliHLTStdIncludes.h"
18 #include <TH2.h>
19 #include <TFile.h>
20 #include <TStyle.h>
21 #include <TGraph.h>
22 #include <TMultiGraph.h>
23 #include <TAttText.h>
24 #include <TAxis.h>
25 #include <TCanvas.h>
26
27
28 #ifdef use_aliroot
29 #include <TClonesArray.h>
30 #include <AliRun.h>
31 #include <AliSimDigits.h>
32 #include <AliTPCParam.h>
33 #endif
34
35 #include "AliHLTTPCDefinitions.h"
36 #include "AliHLTDataTypes.h"
37 #include "AliHLTTPCSpacePointData.h"
38 #include "AliHLTTPCClusterDataFormat.h"
39 #include "AliHLTTPCTrackletDataFormat.h"
40
41
42 #include "AliHLTTPCDigitReader.h"
43 #include "AliHLT_C_Component_WrapperInterface.h"
44
45 #include "AliHLTTPCDisplayMain.h"
46
47 #include "AliHLTTPCLogging.h"
48 #include "AliHLTTPCDisplay.h"
49 #include "AliHLTTPCTransform.h"
50 #include "AliHLTTPCTrack.h"
51 #include "AliHLTTPCTrackArray.h"
52 #include "AliHLTTPCMemHandler.h"
53 #include "AliHLTTPCDigitReaderPacked.h"
54
55
56 #if __GNUC__ >= 3
57 using namespace std;
58 #endif
59
60 ClassImp(AliHLTTPCDisplayResiduals)
61
62 //____________________________________________________________________________________________________
63 AliHLTTPCDisplayResiduals::AliHLTTPCDisplayResiduals(AliHLTTPCDisplayMain* display) {
64     // constructor
65     fDisplay = display;
66
67     fGraphresidualsY = NULL;
68     fGraphresidualsZ = NULL;
69     fGraphresidualsYLength = NULL;
70
71     fHistallresidualsY = new TH1F ("fHistallresidualsY","Y Residuals of all tracks in selected slices;residuals",5000,-100,100);
72     fHistallresidualsZ = new TH1F ("fHistallresidualsZ","Z Residuals of all tracks in selected slices;residuals",5000,-100,100);
73     fHistHits_S = new TH1F ("fHistHits_S","Number of cluster over track length;#Hits/s",20,0,20);
74     fHistQ_Track = new TH1F ("fHistQ_Track","Average cluster charge per track;<q>/track",500,0,500);
75     fHistQ_S   = new TH1F ("fHistQ_S","Total cluster charge over track length;Q/s",500,0,500);
76
77     fHistallresidualsY->SetTitleSize(0.03);
78     fHistallresidualsY->GetXaxis()->SetLabelSize(0.03);
79     fHistallresidualsY->GetXaxis()->SetTitleSize(0.03);
80     fHistallresidualsY->GetYaxis()->SetLabelSize(0.03);
81     fHistallresidualsY->GetYaxis()->SetTitleSize(0.03);
82
83     fHistallresidualsZ->SetTitleSize(0.03);
84     fHistallresidualsZ->GetXaxis()->SetLabelSize(0.03);
85     fHistallresidualsZ->GetXaxis()->SetTitleSize(0.03);
86     fHistallresidualsZ->GetYaxis()->SetLabelSize(0.03);
87     fHistallresidualsZ->GetYaxis()->SetTitleSize(0.03);
88
89     fHistHits_S->SetTitleSize(0.03);
90     fHistHits_S->GetXaxis()->SetLabelSize(0.03);
91     fHistHits_S->GetXaxis()->SetTitleSize(0.03);
92     fHistHits_S->GetYaxis()->SetLabelSize(0.03);
93     fHistHits_S->GetYaxis()->SetTitleSize(0.03);
94
95     fHistQ_Track->SetTitleSize(0.03);
96     fHistQ_Track->GetXaxis()->SetLabelSize(0.03);
97     fHistQ_Track->GetXaxis()->SetTitleSize(0.03);
98     fHistQ_Track->GetYaxis()->SetLabelSize(0.03);
99     fHistQ_Track->GetYaxis()->SetTitleSize(0.03);
100
101     fHistQ_S->SetTitleSize(0.03);
102     fHistQ_S->GetXaxis()->SetLabelSize(0.03);
103     fHistQ_S->GetXaxis()->SetTitleSize(0.03);
104     fHistQ_S->GetYaxis()->SetLabelSize(0.03);
105     fHistQ_S->GetYaxis()->SetTitleSize(0.03);
106 }
107
108 //____________________________________________________________________________________________________
109 AliHLTTPCDisplayResiduals::~AliHLTTPCDisplayResiduals() {
110     // destructor   
111     if ( fHistallresidualsY){
112         delete fHistallresidualsY;
113         fHistallresidualsY = NULL;
114     }
115     fHistallresidualsY = NULL;
116     if ( fHistallresidualsZ){
117         delete fHistallresidualsZ;
118         fHistallresidualsZ = NULL;
119     }
120     if ( fHistHits_S){
121         delete fHistHits_S;
122         fHistHits_S = NULL;
123     }
124     if ( fHistQ_Track){
125         delete fHistQ_Track;
126         fHistQ_Track = NULL;
127     }
128     if ( fHistQ_S){
129         delete fHistQ_S;
130         fHistQ_S = NULL;
131     }
132 }
133
134 //____________________________________________________________________________________________________
135 void AliHLTTPCDisplayResiduals::Reset(){    
136     fHistallresidualsY->Reset();
137     fHistallresidualsZ->Reset();
138     fHistHits_S->Reset();              
139     fHistQ_Track->Reset();           
140     fHistQ_S->Reset();
141
142     fDisplay->GetCanvasResiduals()->Clear();
143     fDisplay->GetCanvasHits_S()->Clear();
144     fDisplay->GetCanvasQ_Track()->Clear();
145     fDisplay->GetCanvasQ_S()->Clear();
146 }    
147
148 //____________________________________________________________________________________________________
149 void AliHLTTPCDisplayResiduals::Save(){  
150     fDisplay->GetCanvasResiduals()->SaveAs("HLT-ResidualsView.eps"); 
151     fDisplay->GetCanvasHits_S()->SaveAs("HLT-Hits_S.eps");
152     fDisplay->GetCanvasQ_Track()->SaveAs("HLT-Q_Track.eps");
153     fDisplay->GetCanvasQ_S()->SaveAs("HLT-Q_S.eps");
154 }
155
156
157 //____________________________________________________________________________________________________
158 void AliHLTTPCDisplayResiduals::Fill(){
159     // Fill Resiudals Histograms / Graphs
160    
161     AliHLTTPCTrackArray* tracks = fDisplay->GetTrackArrayPointer();
162     Int_t ntracks = tracks->GetNTracks();
163
164     Double_t maxResidualY = 0.;
165     Double_t maxResidualZ = 0.;
166
167     fGraphresidualsY = NULL;
168     fGraphresidualsYLength = NULL;
169     fGraphresidualsZ = NULL;
170
171     for(Int_t j=0; j<ntracks; j++) {    
172         
173         AliHLTTPCTrack *gtrack = tracks->GetCheckedTrack(j); 
174         if(!gtrack) continue;   
175
176         // ---------------------------------------------------------------------
177         // ++ Select if slice should be displayed or not ( selection in DISPLAY )
178
179         Int_t slice = gtrack->GetSector();
180         if (!fDisplay->GetDisplaySlice(slice)) continue;
181         
182         // -------------------
183         // Get track parameter
184
185         Int_t nHits = gtrack->GetNHits();  // Number of associated hits to track
186         Double_t radius = gtrack->GetRadius();      // radius
187         Double_t kappa = gtrack->GetKappa();        // curvature = 1/R , signed
188         Double_t lambda = atan( gtrack->GetTgl() ); // dipAngle lambda
189         
190         // -----------------------------
191         // ++ Check for straightr tracks
192
193         if (kappa == 0 && AliHLTTPCTransform::GetBFieldValue() > 0.) {
194             printf("-#- KAPPA == 0 -#-");
195             // continue;
196         }
197
198         // ------------------------------------
199         // ++ Get first/last point of the track
200
201         Double_t xyzL[3];      // lastpoint of track
202         Double_t xyzF[3];      // firstpoint of track
203         
204         xyzF[0] = gtrack->GetFirstPointX();
205         xyzF[1] = gtrack->GetFirstPointY();
206         xyzF[2] = gtrack->GetFirstPointZ();
207         
208         xyzL[0] = gtrack->GetLastPointX();
209         xyzL[1] = gtrack->GetLastPointY();
210         xyzL[2] = gtrack->GetLastPointZ();
211
212         // --------------------------
213         // ++ Calculate length of the track
214         
215         Double_t s = 0.;       // length of the track
216         if (  AliHLTTPCTransform::GetBFieldValue() == 0. || kappa == 0 ) 
217           s = sqrt ( (xyzL[0] - xyzF[0])*(xyzL[0] - xyzF[0]) + (xyzL[1] - xyzF[1])*(xyzL[1] - xyzF[1]) ); 
218         else {
219           // Calculate the length of the track. If it is to flat in in s,z plane use sxy, otherwise use sz
220           if (fabs(lambda) > 0.05){
221             // length of track calculated out of z
222             s = fabs( (xyzL[2] - xyzF[2]) / sin(lambda) ); // length of track calculated out of z
223           }
224           else {
225             Double_t d = (xyzL[0] - xyzF[0])*(xyzL[0] - xyzF[0]) + (xyzL[1] - xyzF[1])*(xyzL[1] - xyzF[1]); 
226             // length of track calculated out of xy
227             s = fabs ( acos( 0.5 * (2 - (d / (radius*radius)))) / ( kappa * cos(lambda) ) );            
228           }
229         }
230         
231         // -----------------------
232         // ++ Apply cuts on tracks
233
234         if (nHits < fDisplay->GetCutHits() ) continue;
235         if (s < fDisplay->GetCutS() ) continue;
236         if (gtrack->GetPsi() < fDisplay->GetCutPsi() ) continue;
237         if (lambda < fDisplay->GetCutLambda() ) continue; 
238         if (gtrack->GetPt() < fDisplay->GetCutPt()  &&  AliHLTTPCTransform::GetBFieldValue() != 0. ) continue;
239         if ( AliHLTTPCTransform::GetPadRow((Float_t)xyzF[0]) >   fDisplay->GetIncidentPadrow() ) continue;
240
241         // ============================================
242         // ## ROTATED Track to local coordinates BEGINN
243
244         gtrack->Rotate(slice,kTRUE);
245         
246         Int_t nRes = 0;                                            // number of resiudals
247         Double_t totalQ = 0.;                                      // total charge of track
248         UInt_t *hitnum = gtrack->GetHitNumbers();                  // hist per track
249         
250         Double_t *resY= new Double_t[nHits];                       // Y residuals of every hit
251         Double_t *resZ= new Double_t[nHits];                       // Z residuals of every hit
252         
253         Double_t *resYLength= new Double_t[2*nHits];               // length of pad in y direction
254         
255         Double_t *padrows = new Double_t[nHits];                   
256         Double_t *padrowsLength = new Double_t[2*nHits];
257
258         // ---------------------
259         // ++ Loop over all hits
260         
261         for(Int_t h=0; h<nHits; h++){
262             UInt_t id=hitnum[h];
263             Int_t patch = (id>>22) & 0x7;
264             UInt_t pos = id&0x3fffff; 
265             
266             AliHLTTPCSpacePointData *points = fDisplay->GetSpacePointDataPointer(slice,patch);
267             if (!points) continue;
268
269             Float_t xyzCtmp[3];    // cluster tmp
270             Float_t xyzTtmp[3];    // track tmp
271             
272             // -------------------------
273             // ++ Get coordinates of Hit
274             xyzCtmp[0] = points[pos].fX;
275             xyzCtmp[1] = points[pos].fY;
276             xyzCtmp[2] = points[pos].fZ;
277             totalQ += points[pos].fCharge;
278
279             // ---------------------------------
280             // ++ Get X Coordinate of the padrow
281
282             Int_t padrow = AliHLTTPCTransform::GetPadRow(points[pos].fX);
283             xyzTtmp[0] = gtrack->GetFirstPointX();
284             
285             // --------------------------------------
286             // ++ Check for CrossingPoint with Padrow
287             if(gtrack->GetCrossingPoint(padrow,xyzTtmp)) {
288                 
289                 // ----------------------
290                 // ++ Calculate Residuals
291
292                 Float_t deltaY = ( xyzCtmp[1] - xyzTtmp[1] );
293                 Float_t deltaZ = ( xyzCtmp[2] - xyzTtmp[2] );
294                 
295                 padrows[nRes] = (Double_t) padrow;
296                 resY[nRes] = (Double_t) deltaY;
297                 resZ[nRes] = (Double_t) deltaZ;
298                 
299                 resYLength[(2*nRes)] = 0.5 * AliHLTTPCTransform::GetPadLength(padrow);
300                 resYLength[(2*nRes)+1] = -0.5 * AliHLTTPCTransform::GetPadLength(padrow);
301                 padrowsLength[(2*nRes)] = (Double_t) padrow;
302                 padrowsLength[(2*nRes)+1] = (Double_t) padrow;
303                 
304                 // ---------------------------
305                 // ++ Fill residuals histogram
306
307                 fHistallresidualsY->Fill(resY[nRes]);
308                 fHistallresidualsZ->Fill(resZ[nRes]);
309
310                 if (resY[nRes] > maxResidualY ) maxResidualY = resY[nRes];
311                 if (resZ[nRes] > maxResidualZ ) maxResidualZ = resZ[nRes];
312                 nRes++;
313
314             } // END CrossingPoint
315         } // END cluster loop
316
317         gtrack->Rotate(slice,kFALSE);
318
319         // ## ROTATED Track to local coordinates END
320         // =========================================
321
322         // -------------------------------------
323         // ++ Fill Number Hits over track length
324
325         Double_t hits_S = nHits / s;
326         fHistHits_S->Fill(hits_S);
327
328         // --------------------------------
329         // ++ Fill Average charge per track
330
331         Double_t q_Track = totalQ / nHits;
332         fHistQ_Track->Fill(q_Track);
333
334         // -------------------------------------
335         // ++ Fill total charge per track length
336
337         Double_t Q_S = totalQ / s;
338         fHistQ_S->Fill(Q_S);
339
340         // --------------------------
341         // ++ Fill Graphs for 1 track
342
343         if (fDisplay->GetSelectTrackSwitch() && fDisplay->GetGlobalTrack(slice) == j){
344             // FILL Y RESIDUALS GRAPH
345             fGraphresidualsY = new TGraph(nRes-1,padrows,resY);
346             fGraphresidualsYLength = new TGraph((2*nRes)-2,padrowsLength,resYLength);
347             // FILL Z RESIDUALS GRAPH
348             fGraphresidualsZ = new TGraph(nRes-1,padrows,resZ);
349         }
350         
351         // --------------
352         // ++ Free Memory
353         
354         if ( resY ){
355             delete [] resY;
356             resY = NULL;
357         }
358         if ( resZ ){
359             delete [] resZ;
360             resZ = NULL;
361         }
362         if ( resYLength ){
363             delete [] resYLength;
364             resYLength = NULL;
365         }
366         if ( padrows ){
367             delete [] padrows;
368             padrows = NULL;
369         }
370         if ( padrowsLength ){
371             delete [] padrowsLength;
372             padrowsLength = NULL;
373         }
374         
375     } // END track loop
376
377     // ----------------------------------------
378     // ++ Set Axis Range of residual histograms
379
380     fHistallresidualsY->SetAxisRange(-maxResidualY,maxResidualY);
381     fHistallresidualsZ->SetAxisRange(-maxResidualZ,maxResidualZ);
382 }
383
384 //____________________________________________________________________________________________________
385 void AliHLTTPCDisplayResiduals::Draw(){
386
387     fDisplay->GetCanvasResiduals()->cd();
388     fDisplay->GetCanvasResiduals()->Clear();
389     
390     fDisplay->GetCanvasResiduals()->Divide(1,2);
391    
392     // ++ Y Residuals     
393     // ==============
394     fDisplay->GetCanvasResiduals()->cd(1);
395    
396     if (fDisplay->GetSelectTrackSwitch() && fGraphresidualsY){
397
398         // -------------------------
399         // ++ Y Graph for single track
400         
401         Char_t title[256];
402         sprintf(title,"Y Residuals of Track %d in Slice %d",fDisplay->GetSelectTrack(), fDisplay->GetSelectTrackSlice() );
403
404         TMultiGraph *mgY = new TMultiGraph() ;
405         mgY->SetBit(kCanDelete);
406         
407         fGraphresidualsY->GetXaxis()->SetTitle("padrow");       
408         fGraphresidualsY->GetYaxis()->SetTitle("residuals");
409         fGraphresidualsY->GetXaxis()->SetLabelSize(0.02);
410         fGraphresidualsY->GetXaxis()->SetTitleSize(0.02);
411         fGraphresidualsY->GetYaxis()->SetLabelSize(0.02);
412         fGraphresidualsY->GetYaxis()->SetTitleSize(0.02);
413
414         fGraphresidualsYLength->SetMarkerColor(2);
415         fGraphresidualsYLength->SetMarkerStyle(5);
416         fGraphresidualsY->SetMarkerColor(1);
417         fGraphresidualsY->SetMarkerStyle(3);
418
419         mgY->Add(fGraphresidualsY);
420         mgY->Add(fGraphresidualsYLength);
421         mgY->SetTitle(title);
422         mgY->Draw("AP");
423     }
424
425     else{
426
427         // -------------------------
428         // ++ Y Histogram  -> global
429         
430         fHistallresidualsY->SetStats(kFALSE);
431         fHistallresidualsY->Draw();
432     }
433     
434     // ++ Z Residuals     
435     // ==============
436     fDisplay->GetCanvasResiduals()->cd(2);
437
438     // graph for single track
439     if (fDisplay->GetSelectTrackSwitch() && fGraphresidualsZ){
440
441         // -------------------------
442         // ++ Z Graph for single track
443
444         Char_t title[256];
445         sprintf(title,"Z Residuals of Track %d in Slice %d",fDisplay->GetSelectTrack(), fDisplay->GetSelectTrackSlice() );
446         
447         TMultiGraph *mgZ = new TMultiGraph();
448         mgZ->SetBit(kCanDelete);
449
450         fGraphresidualsZ->SetTitle(title);      
451         fGraphresidualsZ->GetXaxis()->SetTitle("padrow");       
452         fGraphresidualsZ->GetYaxis()->SetTitle("residuals");
453         fGraphresidualsZ->GetXaxis()->SetLabelSize(0.02);
454         fGraphresidualsZ->GetXaxis()->SetTitleSize(0.02);
455         fGraphresidualsZ->GetYaxis()->SetLabelSize(0.02);
456         fGraphresidualsZ->GetYaxis()->SetTitleSize(0.02);
457         
458         mgZ->Add(fGraphresidualsZ);
459         mgZ->SetTitle(title);
460         mgZ->Draw("A*");
461     }
462     
463         // -------------------------
464         // ++ Z Histogram  -> global
465
466     else{
467         fHistallresidualsZ->SetStats(kFALSE);
468         fHistallresidualsZ->Draw();
469     }
470
471     fDisplay->GetCanvasResiduals()->Modified();
472     fDisplay->GetCanvasResiduals()->Update();
473
474     // -------------------------------------
475     // ++ Draw Number Hits over track length
476
477     fDisplay->GetCanvasHits_S()->cd ();
478     fDisplay->GetCanvasHits_S()->Clear();
479     
480     fHistHits_S->Draw();
481
482     fDisplay->GetCanvasHits_S()->Modified();
483     fDisplay->GetCanvasHits_S()->Update();
484
485     // --------------------------------
486     // ++ Draw Average charge per track
487
488     fDisplay->GetCanvasQ_Track()->cd();
489     fDisplay->GetCanvasQ_Track()->Clear();
490
491     fHistQ_Track->Draw();
492
493     fDisplay->GetCanvasQ_Track()->Modified();
494     fDisplay->GetCanvasQ_Track()->Update();
495
496     // -------------------------------------
497     // ++ Draw total charge per track length
498
499     fDisplay->GetCanvasQ_S()->cd();
500     fDisplay->GetCanvasQ_S()->Clear();
501     
502     fHistQ_S->Draw();
503   
504     fDisplay->GetCanvasQ_S()->Modified();
505     fDisplay->GetCanvasQ_S()->Update();
506 }
507