Some problems with the HP compiler fixed.
[u/mrichter/AliRoot.git] / MUON / AliMUONpoints.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 /*
17 $Log$
18 */
19
20 ///////////////////////////////////////////////////////////////////////////////
21 //                                                                           //
22 //  This class contains the points for the ALICE event display               //
23 //                                                                           //
24 //Begin_Html
25 /*
26 <img src="gif/AliMUONpointsClass.gif">
27 */
28 //End_Html
29 //                                                                           //
30 //                                                                           //
31 ///////////////////////////////////////////////////////////////////////////////
32
33 #include "AliMUONpoints.h"
34 #include "AliMUONdisplay.h"
35 #include "AliRun.h"
36 #include "TPad.h"
37 #include "TVirtualPad.h"
38 #include "TPolyLine3D.h"
39 #include "TPaveText.h"
40 #include "TView.h"
41 #include "TMath.h"
42
43 //const Int_t MAX_Nipx=1026, MAX_Nipy=1026;
44 const Int_t MAX_Nipx=400, MAX_Nipy=800;
45  
46 ClassImp(AliMUONpoints)
47
48 //_____________________________________________________________________________
49 AliMUONpoints::AliMUONpoints()
50 {
51   //
52   // Default constructor
53   //
54   fHitIndex = 0;
55   fTrackIndex = 0;
56   fDigitIndex = 0;
57   fMarker[0] = fMarker[1] = fMarker[2]=0;
58   fMatrix = 0;
59 }
60
61 //_____________________________________________________________________________
62 AliMUONpoints::AliMUONpoints(Int_t npoints)
63   :AliPoints(npoints)
64 {
65   //
66   // Standard constructor
67   //
68   fHitIndex = 0;
69   fTrackIndex = 0;
70   fDigitIndex = 0;
71   fMarker[0] = fMarker[1] = fMarker[2]=0;
72   fMatrix = 0;
73 }
74          
75 //_____________________________________________________________________________
76 AliMUONpoints::~AliMUONpoints()
77 {
78   //
79   // Default destructor
80   //
81   fHitIndex = 0;
82   fTrackIndex = 0;
83   fDigitIndex = 0;
84   for (Int_t i=0;i<3;i++){
85       if (fMarker[i]) delete fMarker[i];
86   }
87   fMatrix = 0;
88 }
89
90 //_____________________________________________________________________________
91 //void AliMUONpoints::ExecuteEvent(Int_t event, Int_t px, Int_t py)
92 //{
93   //
94   //*-*-*-*-*-*-*-*-*-*Execute action corresponding to one event*-*-*-*-*-*-*-*
95   //*-*                =========================================
96   //*-*
97   //*-*  This member function must be implemented to realize the action
98   //*-*  corresponding to the mouse click on the object in the window
99   //*-*
100   //*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
101 //  gPad->SetCursor(kCross);
102   
103 //}
104 //_____________________________________________________________________________
105 void AliMUONpoints::DumpHit()
106 {
107   //
108   //   Dump hit corresponding to this point
109   //
110   AliMUONhit *hit = GetHit();
111   if (hit) hit->Dump();
112 }
113
114 //_____________________________________________________________________________
115 void AliMUONpoints::DumpDigit()
116 {
117   //
118   //   Dump digit corresponding to this point
119   //
120   AliMUONdigit *digit = GetDigit();
121   if (digit) digit->Dump();
122 }
123
124 //_____________________________________________________________________________
125 void AliMUONpoints::InspectHit()
126 {
127   //
128   //   Inspect hit corresponding to this point
129   //
130
131   if (fHitIndex < 0 ) return;
132   TVirtualPad *padsav = gPad;
133   AliMUONhit *hit = GetHit();
134   if (hit) hit->Inspect();
135   TVirtualPad *padinspect = (TVirtualPad*)(gROOT->GetListOfCanvases())->FindObject("inspect");
136    padinspect->cd();
137    Float_t xmin = gPad->GetX1();
138    Float_t xmax = gPad->GetX2();
139    Float_t ymin = gPad->GetY1();
140    Float_t ymax = gPad->GetY2();
141    Float_t dy   = ymax-ymin;
142
143       TPaveText *pad = new TPaveText(xmin, ymin+0.1*dy, xmax, ymin+0.15*dy);
144       pad->SetBit(kCanDelete);
145       pad->SetFillColor(42);
146       pad->Draw();
147       char ptitle[100];
148       sprintf(ptitle," %s , fTrack: %d  fTrackIndex: %d ",GetName(),fIndex,fTrackIndex);
149       pad->AddText(ptitle);
150       padinspect->cd();
151       padinspect->Update();
152   if (padsav) padsav->cd();
153
154 }
155
156 //_____________________________________________________________________________
157 void AliMUONpoints::InspectDigit()
158 {
159   //
160   //   Inspect digit corresponding to this point
161   //
162   if (fDigitIndex < 0) return;
163   TVirtualPad *padsav = gPad;
164   AliMUONdigit *digit = GetDigit();
165   if (digit) digit->Inspect();
166   TVirtualPad *padinspect = (TVirtualPad*)(gROOT->GetListOfCanvases())->FindObject("inspect");
167    padinspect->cd();
168    Float_t xmin = gPad->GetX1();
169    Float_t xmax = gPad->GetX2();
170    Float_t ymin = gPad->GetY1();
171    Float_t ymax = gPad->GetY2();
172    Float_t dy   = ymax-ymin;
173
174       TPaveText *pad = new TPaveText(xmin, ymin+0.1*dy, xmax, ymin+0.25*dy);
175       pad->SetBit(kCanDelete);
176       pad->SetFillColor(42);
177       pad->Draw();
178       char ptitle[11][100];
179       //      sprintf(ptitle[11],"Tracks making this digit");
180       //      pad->AddText(ptitle[11]);
181   for (int i=0;i<10;i++) {
182       if (digit->fTracks[i] == 0) continue;  
183       sprintf(ptitle[i],"fTrackIndex: %d  Charge: %d",digit->fTracks[i],digit->fTcharges[i]);
184       pad->AddText(ptitle[i]);
185   }
186       padinspect->cd();
187       padinspect->Update();
188   if (padsav) padsav->cd();
189       
190 }
191
192 //_____________________________________________________________________________
193 Int_t AliMUONpoints::GetTrackIndex()
194 {
195   //
196   //   Dump digit corresponding to this point
197   //
198
199   this->Inspect();
200   /*
201   if (fDigitIndex != 0) {
202     Int_t ncol=this->fMatrix->GetNcols();
203     for (int i=0;i<ncol;i++) {
204         printf(" track charge %f %f \n",(*(this->fMatrix))(0,i),(*(this->fMatrix))(1,i));
205     }
206   }
207   */
208   return fTrackIndex;
209 }
210
211 //_____________________________________________________________________________
212 AliMUONhit *AliMUONpoints::GetHit() const
213 {
214   //
215   //   Returns pointer to hit index in AliRun::fParticles
216   //
217   AliMUON *MUON  = (AliMUON*)gAlice->GetModule("MUON");
218   gAlice->TreeH()->GetEvent(fTrackIndex);
219   TClonesArray *MUONhits  = MUON->Hits();
220   Int_t nhits = MUONhits->GetEntriesFast();
221   if (fHitIndex < 0 || fHitIndex >= nhits) return 0;
222   return (AliMUONhit*)MUONhits->UncheckedAt(fHitIndex);
223 }
224
225 //_____________________________________________________________________________
226 AliMUONdigit *AliMUONpoints::GetDigit() const
227 {
228   //
229   //   Returns pointer to digit index in AliRun::fParticles
230   //
231
232   AliMUONdisplay *display=(AliMUONdisplay*)gAlice->Display();
233   Int_t chamber=display->GetChamber();
234   Int_t cathode=display->GetCathode();
235    
236   AliMUON *MUON  = (AliMUON*)gAlice->GetModule("MUON");
237   TClonesArray *MUONdigits  = MUON->DigitsAddress(chamber-1);
238   Int_t nent=(Int_t)gAlice->TreeD()->GetEntries();
239   gAlice->TreeD()->GetEvent(nent-2+cathode-1);
240   //gAlice->TreeD()->GetEvent(cathode);
241   Int_t ndigits = MUONdigits->GetEntriesFast();
242   if (fDigitIndex < 0 || fDigitIndex >= ndigits) return 0;
243   return (AliMUONdigit*)MUONdigits->UncheckedAt(fDigitIndex);
244 }
245 //_____________________________________________________________________________
246 struct Bin {
247    const AliMUONdigit *dig;
248    int idx;
249    Bin() {dig=0; idx=-1;}
250 };
251
252 struct PreCluster : public AliMUONreccluster {
253    const AliMUONdigit* summit;
254    int idx;
255    int cut;
256    int npeaks;
257    int npoly;
258    float xpoly[100];
259    float ypoly[100];
260    float zpoly[100];
261    PreCluster() : AliMUONreccluster() {cut=npeaks=npoly=0; 
262                                        for (int k=0;k<100;k++) {
263                                          xpoly[k]=ypoly[k]=zpoly[k]=0;
264                                        }
265    }
266                               
267 };
268 //_____________________________________________________________________________
269
270 static void FindCluster(AliMUONchamber *iChamber, AliMUONsegmentation *segmentation, int i, int j, Bin bins[][MAX_Nipy], PreCluster &c) 
271
272 {
273
274   //
275   // Find clusters
276   //
277
278
279   Bin& b=bins[i][j];
280   Int_t q=b.dig->fSignal;
281
282   if (q<0) { 
283     q=-q;
284     c.cut=1;
285   } 
286   //  if (b.idx >= 0 && b.idx != c.idx) {
287   if (b.idx >= 0 && b.idx > c.idx) {
288     c.idx=b.idx;
289     c.npeaks++;
290   }
291   
292   if (q > TMath::Abs(c.summit->fSignal)) c.summit=b.dig;
293
294   Float_t zpos=iChamber->ZPosition();
295
296   Int_t npx  = segmentation->Npx();
297   Int_t npy  = segmentation->Npy();
298
299
300   // get pad coordinates and prepare the up and down steps   
301   Int_t jup  =(j-npy > 0) ? j+1 : j-1;
302   Int_t jdown=(j-npy > 0) ? j-1 : j+1;
303
304   Float_t x, y;
305   segmentation->GetPadCxy(i-npx, j-npy,x, y);
306   Int_t isec0=segmentation->Sector(i-npx,j-npy);
307
308   Float_t dpy  = segmentation->Dpy(isec0);
309   Float_t dpx  = segmentation->Dpx(isec0)/16;
310   Int_t ixx, iyy;
311   Float_t absx=TMath::Abs(x);
312   // iup
313   (y >0) ? segmentation->GetPadIxy(absx+dpx,y+dpy,ixx,iyy) : segmentation->GetPadIxy(absx+dpx,y-dpy,ixx,iyy);
314
315   //  Int_t jtest=TMath::Abs(iyy)-npy-1;
316
317   Int_t iup=(x >0) ? ixx+npx : -ixx+npx;
318   // idown
319   (y >0) ? segmentation->GetPadIxy(absx+dpx,y-dpy,ixx,iyy) : segmentation->GetPadIxy(absx+dpx,y+dpy,ixx,iyy);
320
321   Int_t idown=(x >0) ? ixx+npx : -ixx+npx;
322   if (bins[idown][jdown].dig == 0) {
323      (y >0) ? segmentation->GetPadIxy(absx-dpx,y-dpy,ixx,iyy) : segmentation->GetPadIxy(absx-dpx,y+dpy,ixx,iyy);
324      idown=(x >0) ? ixx+npx : -ixx+npx;
325   }
326  
327
328   // calculate center of gravity
329   c.npoly++;
330   if (c.npoly > 100 ) {
331     printf("FindCluster - npoly >100,  npoly %d \n",c.npoly);
332     c.npoly=100;
333   }
334   c.xpoly[c.npoly-1]=x;
335   c.ypoly[c.npoly-1]=y;
336   c.zpoly[c.npoly-1]=zpos;
337
338   c.fX += q*x;
339   c.fY += q*y;
340   c.fQ += q;
341   
342   b.dig = 0;  b.idx = c.idx;
343
344   // left and right  
345   if (bins[i-1][j].dig) FindCluster(iChamber,segmentation,i-1,j,bins,c);
346   if (bins[i+1][j].dig) FindCluster(iChamber,segmentation,i+1,j,bins,c);
347   // up and down
348   if (bins[iup][jup].dig) FindCluster(iChamber,segmentation,iup,jup,bins,c);
349   if (bins[idown][jdown].dig) FindCluster(iChamber,segmentation,idown,jdown,bins,c);
350
351 }
352
353 //_____________________________________________________________________________
354
355 void AliMUONpoints::GetCenterOfGravity()
356 {
357   //
358   // simple MUON cluster finder from digits -- finds neighbours and 
359   // calculates center of gravity for the cluster
360   //
361
362   //  const Int_t MAX_Nipx=1026, MAX_Nipy=1026;
363
364   Bin bins[MAX_Nipx][MAX_Nipy]; 
365
366   AliMUONdisplay *display=(AliMUONdisplay*)gAlice->Display();
367   Int_t chamber=display->GetChamber();
368   Int_t cathode=display->GetCathode();
369    
370   AliMUON *MUON  = (AliMUON*)gAlice->GetModule("MUON");
371   AliMUONchamber *iChamber;
372   AliMUONsegmentation *segmentation;
373   iChamber =&(MUON->Chamber(chamber-1));
374   segmentation=iChamber->GetSegmentationModel(cathode);
375   Int_t npx  = segmentation->Npx();
376   Int_t npy  = segmentation->Npy();
377   Float_t zpos=iChamber->ZPosition();
378
379   TClonesArray *MUONdigits  = MUON->DigitsAddress(chamber-1);
380   Int_t nent=(Int_t)gAlice->TreeD()->GetEntries();
381   gAlice->TreeD()->GetEvent(nent-2+cathode-1);
382   //gAlice->TreeD()->GetEvent(cathode);
383   Int_t ndigits = MUONdigits->GetEntriesFast();
384   if (fDigitIndex < 0 || fDigitIndex >= ndigits) return;
385
386   AliMUONdigit  *dig;
387   dig=(AliMUONdigit*)MUONdigits->UncheckedAt(fDigitIndex);
388   Int_t ipx=dig->fPadX;
389   Int_t ipy=dig->fPadY;
390   bins[ipx+npx][ipy+npy].dig=dig;
391     
392   int ndig;
393   int ncls=0;
394   for (ndig=0;ndig<ndigits;ndig++) {
395     dig = (AliMUONdigit*)MUONdigits->UncheckedAt(ndig);
396     int i=dig->fPadX, j=dig->fPadY;
397     bins[i+npx][j+npy].dig=dig;
398   }
399    PreCluster c; c.summit=bins[ipx+npx][ipy+npy].dig; c.idx=ncls;
400    FindCluster(iChamber,segmentation,ipx+npx, ipy+npy, bins, c);
401     
402    if (c.npeaks>1) {
403       printf("GetCenterOfGravity -- more than one peak");
404    }
405    c.fX /= c.fQ;
406    c.fY /= c.fQ;
407
408    c.fTracks[0]=c.summit->fTracks[0];
409    c.fTracks[1]=c.summit->fTracks[1];
410    c.fTracks[2]=c.summit->fTracks[2];
411    ncls++;
412    AliMUONpoints *points = 0;
413    points = new AliMUONpoints(1);
414    points->SetMarkerColor(kYellow);
415    points->SetMarkerStyle(5);
416    points->SetMarkerSize(1.);
417    points->SetPoint(0,c.fX,c.fY,zpos);
418    points->SetParticle(-1);
419    points->Draw();
420
421    TPolyLine3D *pline=0;
422    /*
423    pline=new TPolyLine3D(c.npoly);
424    Int_t np=c.npoly;
425    TVector *xp=new TVector(c.npoly);
426    TVector *yp=new TVector(c.npoly);
427    TVector *zp=new TVector(c.npoly);
428    for (int i=0;i<np;i++) {
429      (*xp)(i)=c.xpoly[i];
430      (*yp)(i)=c.ypoly[i];
431      (*zp)(i)=c.zpoly[i];
432      pline->SetPoint(i,(*xp)(i),(*yp)(i),(*zp)(i));
433      //printf("np, i, xp, yp, zp %d %d %f %f %f \n",np,i,(*xp)(i),(*yp)(i),(*zp)(i));
434    }
435    */
436    pline=new TPolyLine3D(c.npoly,c.xpoly,c.ypoly,c.zpoly);
437    pline->SetLineColor(kWhite);
438    pline->Draw();
439    /*
440    delete xp;
441    delete yp;
442    delete zp;
443    */  
444    for (int k=0;k<c.npoly;k++) {
445      c.xpoly[k]=c.ypoly[k]=c.zpoly[k]=0;
446    }
447    c.npoly=0;
448
449
450 }
451