c60c1031b850fdbda227202e4667fff904445716
[u/mrichter/AliRoot.git] / TRD / AliTRDdisplayDigits3D.C
1 //_____________________________________________________________________________
2 Int_t AliTRDdisplayDigits3D(Int_t event = 0, Int_t thresh = 0
3                           , Bool_t sdigits = kFALSE) 
4 {
5   //  
6   //  TRD digits display
7   //
8   //  Input parameter:
9   //    <event>   : Event number 
10   //    <thresh>  : Threshold to suppress the noise
11   //    <sdigits> : If kTRUE it will display summable digits, normal digits otherwise.
12   //                The signal event is displayed in yellow.
13   //
14
15   Char_t *inputFile = "galice.root";
16
17   const Int_t kNdict = 3;
18
19   // Define the objects
20   AliTRDv1       *trd;
21   AliTRDgeometry *geo;
22
23   Int_t           track;
24   Int_t           idict;
25
26   // Connect the AliRoot file containing Geometry, Kine, Hits, and digits
27   TFile *gafl = (TFile*) gROOT->GetListOfFiles()->FindObject(inputFile);
28   if (!gafl) {
29     printf("Open the ALIROOT-file %s\n",inputFile);
30     gafl = new TFile(inputFile);
31   }
32   else {
33     printf("%s is already open",inputFile);
34   }
35
36   // Get AliRun object from file or create it if not on file
37   gAlice = (AliRun *) gafl->Get("gAlice");
38   if (!gAlice) return 1;
39
40   // Import the Trees for the event nEvent in the file
41   Int_t nparticles = gAlice->GetEvent(event);
42   if (nparticles <= 0) return 1;
43   
44   // Get the pointer to the detector object
45   trd = (AliTRDv1*) gAlice->GetDetector("TRD");
46
47   // Get the pointer to the geometry object
48   if (trd) {
49     geo = trd->GetGeometry();
50   }
51   else {
52     printf("Cannot find the geometry\n");
53     return 1;
54   }
55
56   TCanvas *c1 = new TCanvas("digits","TRD digits display",0,0,700,730);
57   TView   *v  = new TView(1);
58   v->SetRange(-430,-560,-430,430,560,1710);
59   c1->Clear();
60   c1->SetFillColor(1);
61   c1->SetTheta(90.0);
62   c1->SetPhi(0.0);
63
64   Int_t markerColorSignal = 2;
65   Int_t markerColorBgnd   = 7;
66   Int_t markerColorMerged = 5;
67   Int_t mask              = 10000000;
68
69   // Create the digits manager
70   AliTRDdigitsManager *digitsManager = new AliTRDdigitsManager();
71   digitsManager->SetDebug(1);
72   digitsManager->SetSDigits(sdigits);
73
74   // Read the digits from the file
75   digitsManager->Open(inputFile);
76   digitsManager->ReadDigits();
77
78   Int_t totalsignal = 0;
79   Int_t totalbgnd   = 0;
80   Int_t totalmerged = 0;
81
82   // Loop through all detectors
83   for (Int_t idet = 0; idet < geo->Ndet(); idet++) {
84
85     printf("<AliTRDdisplayDigits3D> Loading detector %d\n",idet);
86     AliTRDdataArrayI *digits  = digitsManager->GetDigits(idet);
87     digits->Expand();
88     AliTRDdataArrayI *tracks[kNdict];
89     for (Int_t idict = 0; idict < kNdict; idict++) {
90       tracks[idict] = digitsManager->GetDictionary(idet,idict);
91       tracks[idict]->Expand();
92     }
93
94     Int_t isec    = geo->GetSector(idet);
95     Int_t icha    = geo->GetChamber(idet);
96     Int_t ipla    = geo->GetPlane(idet);
97     Int_t  rowMax = geo->GetRowMax(ipla,icha,isec);
98     Int_t  colMax = geo->GetColMax(ipla);
99     Int_t timeMax = geo->GetTimeMax();
100
101     Int_t ndigits = digits->GetOverThreshold(thresh);
102
103     if (ndigits > 0) {
104
105       TPolyMarker3D *pmSignal = new TPolyMarker3D(ndigits);
106       TPolyMarker3D *pmBgnd   = new TPolyMarker3D(ndigits);
107       TPolyMarker3D *pmMerged = new TPolyMarker3D(ndigits);
108  
109       Int_t ibgnd   = 0;
110       Int_t isignal = 0;
111       Int_t imerged = 0;
112
113       for (Int_t time = 0; time < timeMax; time++) {
114         for (Int_t  col = 0;  col <  colMax;  col++) {
115           for (Int_t  row = 0;  row <  rowMax;  row++) {
116
117             Int_t type = 1;
118
119             Int_t amp = digits->GetDataUnchecked(row,col,time);
120             for (idict = 0; idict < kNdict; idict++) {
121               Int_t trk = tracks[idict]->GetDataUnchecked(row,col,time) - 1;
122               if ((idict == 0) && (trk >= mask)) {
123                 type = 2;
124               }
125               if ((type  == 1) && (trk >= mask)) {
126                 type = 3;
127               }              
128             }
129
130             if (amp > thresh) {
131           
132               Float_t glb[3];
133               Float_t loc[3];
134
135               loc[0] = row;
136               loc[1] = col;
137               loc[2] = time;
138               geo->Local2Global(idet,loc,glb);
139               Float_t x = glb[0];
140               Float_t y = glb[1];
141               Float_t z = glb[2];
142
143               if      (type == 1) {
144                 pmSignal->SetPoint(isignal,x,y,z);
145                 isignal++;
146                 totalsignal++;
147               }
148               else if (type == 2) {
149                 pmBgnd->SetPoint(ibgnd,x,y,z);
150                 ibgnd++;
151                 totalbgnd++;
152               }
153               else if (type == 3) {
154                 pmMerged->SetPoint(imerged,x,y,z);
155                 imerged++;
156                 totalmerged++;
157               }
158
159             }
160
161           }
162         }
163       }
164
165       digits->Compress(1,0);
166       for (idict = 0; idict < kNdict; idict++) {
167         tracks[idict]->Compress(1,0);
168       }
169
170       pmMerged->SetMarkerSize(1); 
171       pmMerged->SetMarkerColor(markerColorMerged);
172       pmMerged->SetMarkerStyle(1);
173       pmMerged->Draw();
174
175       pmBgnd->SetMarkerSize(1); 
176       pmBgnd->SetMarkerColor(markerColorBgnd);
177       pmBgnd->SetMarkerStyle(1);
178       pmBgnd->Draw();
179
180       pmSignal->SetMarkerSize(1); 
181       pmSignal->SetMarkerColor(markerColorSignal);
182       pmSignal->SetMarkerStyle(1);
183       pmSignal->Draw();
184    
185     }
186
187   }
188
189   TGeometry *geoAlice = (TGeometry *) gafl->Get("AliceGeom");
190   TNode     *main     = (TNode *) ((geoAlice->GetListOfNodes())->First());
191   TIter      next(main->GetListOfNodes());
192   TNode     *module   = 0;
193   while ((module = (TNode *) next())) {
194     Char_t ch[100];
195     sprintf(ch,"%s\n",module->GetTitle());
196     if ((ch[0] == 'T') && ((ch[1] == 'R') || (ch[1] == 'P'))) {
197       module->SetVisibility( 3);
198     }
199     else {
200       module->SetVisibility(-1);
201     }
202   }
203       
204   geoAlice->Draw("same");
205
206   c1->Modified(); 
207   c1->Update(); 
208
209   gafl->Close();
210
211   printf("<AliTRDdisplayDigits3D> Number of digits:\n");
212   printf("                        signal = %d, bgnd = %d, merged = %d\n"
213         ,totalsignal,totalbgnd,totalmerged);
214
215   return 0;
216
217 }
218