Introduction of the Copyright and cvs Log
[u/mrichter/AliRoot.git] / MUON / AliMUONpoints.h.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 (
86 fMarker[i]) delete fMarker[i];
87   }
88   fMatrix = 0;
89 }
90
91 //_____________________________________________________________________________
92 //void AliMUONpoints::ExecuteEvent(Int_t event, Int_t px, Int_t py)
93 //{
94   //
95   //*-*-*-*-*-*-*-*-*-*Execute action corresponding to one event*-*-*-*-*-*-*-*
96   //*-*                =========================================
97   //*-*
98   //*-*  This member function must be implemented to realize the action
99   //*-*  corresponding to the mouse click on the object in the window
100   //*-*
101   //*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
102 //  gPad->SetCursor(kCross);
103   
104 //}
105 //_____________________________________________________________________________
106 void AliMUONpoints::DumpHit()
107 {
108   //
109   //   Dump hit corresponding to this point
110   //
111   AliMUONhit *hit = GetHit();
112   if (hit) hit->Dump();
113   printf("fTrackIndex %d \n",fTrackIndex);
114 }
115
116 //_____________________________________________________________________________
117 void AliMUONpoints::DumpDigit()
118 {
119   //
120   //   Dump digit corresponding to this point
121   //
122   AliMUONdigit *digit = GetDigit();
123   if (digit) digit->Dump();
124   for (int i=0;i<10;i++) {
125         printf(" track charge %d %d \n",digit->fTracks[i],digit->fTcharges[i]);
126   }
127
128 }
129
130 //_____________________________________________________________________________
131 void AliMUONpoints::InspectHit()
132 {
133   //
134   //   Inspect hit corresponding to this point
135   //
136
137   if (fHitIndex < 0 ) return;
138   TVirtualPad *padsav = gPad;
139   AliMUONhit *hit = GetHit();
140   if (hit) hit->Inspect();
141   TVirtualPad *padinspect = (TVirtualPad*)(gROOT->GetListOfCanvases())->FindObject("inspect");
142    padinspect->cd();
143    Float_t xmin = gPad->GetX1();
144    Float_t xmax = gPad->GetX2();
145    Float_t ymin = gPad->GetY1();
146    Float_t ymax = gPad->GetY2();
147    Float_t dy   = ymax-ymin;
148
149       TPaveText *pad = new TPaveText(xmin, ymin+0.1*dy, xmax, ymin+0.15*dy);
150       pad->SetBit(kCanDelete);
151       pad->SetFillColor(42);
152       pad->Draw();
153       char ptitle[100];
154       sprintf(ptitle," %s , fTrack: %d  fTrackIndex: %d ",GetName(),fIndex,fTrackIndex);
155       pad->AddText(ptitle);
156       padinspect->cd();
157       padinspect->Update();
158   if (padsav) padsav->cd();
159
160 }
161
162 //_____________________________________________________________________________
163 void AliMUONpoints::InspectDigit()
164 {
165   //
166   //   Inspect digit corresponding to this point
167   //
168   if (fDigitIndex < 0) return;
169   TVirtualPad *padsav = gPad;
170   AliMUONdigit *digit = GetDigit();
171   if (digit) digit->Inspect();
172   TVirtualPad *padinspect = (TVirtualPad*)(gROOT->GetListOfCanvases())->FindObject("inspect");
173    padinspect->cd();
174    Float_t xmin = gPad->GetX1();
175    Float_t xmax = gPad->GetX2();
176    Float_t ymin = gPad->GetY1();
177    Float_t ymax = gPad->GetY2();
178    Float_t dy   = ymax-ymin;
179
180       TPaveText *pad = new TPaveText(xmin, ymin+0.1*dy, xmax, ymin+0.25*dy);
181       pad->SetBit(kCanDelete);
182       pad->SetFillColor(42);
183       pad->Draw();
184       char ptitle[11][100];
185       //      sprintf(ptitle[11],"Tracks making this digit");
186       //      pad->AddText(ptitle[11]);
187   for (int i=0;i<10;i++) {
188       if (digit->fTracks[i] == 0) continue;  
189       sprintf(ptitle[i],"fTrackIndex: %d  Charge: %d",digit->fTracks[i],digit->fTcharges[i]);
190       pad->AddText(ptitle[i]);
191   }
192       padinspect->cd();
193       padinspect->Update();
194   if (padsav) padsav->cd();
195       
196 }
197
198 //_____________________________________________________________________________
199 Int_t AliMUONpoints::GetTrackIndex()
200 {
201   //
202   //   Dump digit corresponding to this point
203   //
204   printf("GetTrackIndex - fTrackIndex %d \n",fTrackIndex);
205   this->Inspect();
206   /*
207   if (fDigitIndex != 0) {
208     Int_t ncol=this->fMatrix->GetNcols();
209     for (int i=0;i<ncol;i++) {
210         printf(" track charge %f %f \n",(*(this->fMatrix))(0,i),(*(this->fMatrix))(1,i));
211     }
212   }
213   */
214   return fTrackIndex;
215 }
216
217 //_____________________________________________________________________________
218 AliMUONhit *AliMUONpoints::GetHit() const
219 {
220   //
221   //   Returns pointer to hit index in AliRun::fParticles
222   //
223   AliMUON *MUON  = (AliMUON*)gAlice->GetModule("MUON");
224   gAlice->TreeH()->GetEvent(fTrackIndex);
225   TClonesArray *MUONhits  = MUON->Hits();
226   Int_t nhits = MUONhits->GetEntriesFast();
227   if (fHitIndex < 0 || fHitIndex >= nhits) return 0;
228   return (AliMUONhit*)MUONhits->UncheckedAt(fHitIndex);
229 }
230
231 //_____________________________________________________________________________
232 AliMUONdigit *AliMUONpoints::GetDigit() const
233 {
234   //
235   //   Returns pointer to digit index in AliRun::fParticles
236   //
237
238   AliMUONdisplay *display=(AliMUONdisplay*)gAlice->Display();
239   Int_t chamber=display->GetChamber();
240   Int_t cathode=display->GetCathode();
241    
242   AliMUON *MUON  = (AliMUON*)gAlice->GetModule("MUON");
243   TClonesArray *MUONdigits  = MUON->DigitsAddress(chamber-1);
244   gAlice->TreeD()->GetEvent(cathode);
245   Int_t ndigits = MUONdigits->GetEntriesFast();
246   if (fDigitIndex < 0 || fDigitIndex >= ndigits) return 0;
247   return (AliMUONdigit*)MUONdigits->UncheckedAt(fDigitIndex);
248 }
249 //_____________________________________________________________________________
250 struct Bin {
251    const AliMUONdigit *dig;
252    int idx;
253    Bin() {dig=0; idx=-1;}
254 };
255
256 struct PreCluster : public AliMUONreccluster {
257    const AliMUONdigit* summit;
258    int idx;
259    int cut;
260    int npeaks;
261    int npoly;
262    float xpoly[100];
263    float ypoly[100];
264    float zpoly[100];
265    PreCluster() : AliMUONreccluster() {cut=npeaks=npoly=0; 
266                                        for (int k=0;k<100;k++) {
267                                          xpoly[k]=ypoly[k]=zpoly[k]=0;
268                                        }
269    }
270                               
271 };
272 //_____________________________________________________________________________
273
274 static void FindCluster(AliMUONchamber *iChamber, AliMUONsegmentation *segmentation, int i, int j, Bin bins[][MAX_Nipy], PreCluster &c) 
275
276 {
277
278   //
279   // Find clusters
280   //
281
282   printf("I'm in FindCluster \n"); 
283
284   Bin& b=bins[i][j];
285   Int_t q=b.dig->fSignal;
286
287   printf("FindCluster - i j q %d %d %d\n",i,j,q);
288   
289   if (q<0) { 
290     q=-q;
291     c.cut=1;
292   } 
293   //  if (b.idx >= 0 && b.idx != c.idx) {
294   if (b.idx >= 0 && b.idx > c.idx) {
295     c.idx=b.idx;
296     c.npeaks++;
297     printf("FindCluster - one more peak \n");
298   }
299   
300   if (q > TMath::Abs(c.summit->fSignal)) c.summit=b.dig;
301
302   Float_t zpos=iChamber->ZPosition();  // check with Andreas
303
304   Int_t npx  = segmentation->Npx();
305   Int_t npy  = segmentation->Npy();
306
307   // get pad coordinates and prepare the up and down steps   
308   Int_t jup  =(j-npy > 0) ? j+1 : (j-npy-1)+npy;
309   Int_t jdown=(j-npy > 0) ? j-1 : (j-npy+1)+npy;
310   //  Int_t jup  =j+1;
311   //  Int_t jdown=j-1;
312   Float_t x, y;
313   segmentation->GetPadCxy(i-npx, j-npy,x, y);
314   Int_t isec0=segmentation->Sector(i-npx,j-npy);
315   printf("FindCluster - i-npx j-npy isec0 %d %d %d \n",i-npx,j-npy,isec0);
316   //  printf("FindCluster - x  y %f %f \n",x,y);
317
318   Float_t dpy  = segmentation->Dpy();
319   Float_t dpx  = segmentation->Dpx()/16;
320   Int_t ixx, iyy;
321   Float_t absx=TMath::Abs(x);
322   // iup
323   (y >0) ? segmentation->GetPadIxy(absx+dpx,y+dpy,ixx,iyy) : segmentation->GetPadIxy(absx+dpx,y-dpy,ixx,iyy);
324   printf(" iup: ixx iyy %d %d \n",ixx,iyy);
325   Int_t jtest=TMath::Abs(iyy)-npy-1;
326   if (j != jtest) {
327      printf(" j != jtest - something's wrong %d %d \n",j,jtest);
328   }
329   Int_t iup=(x >0) ? ixx+npx : -ixx+npx;
330   // idown
331   (y >0) ? segmentation->GetPadIxy(absx+dpx,y-dpy,ixx,iyy) : segmentation->GetPadIxy(absx+dpx,y+dpy,ixx,iyy);
332   printf(" idown: ixx iyy %d %d \n",ixx,iyy);
333   Int_t idown=(x >0) ? ixx+npx : -ixx+npx;
334   if (bins[idown][jdown].dig == 0) {
335      (y >0) ? segmentation->GetPadIxy(absx-dpx,y-dpy,ixx,iyy) : segmentation->GetPadIxy(absx-dpx,y+dpy,ixx,iyy);
336      printf(" idown: ixx iyy %d %d \n",ixx,iyy);
337      idown=(x >0) ? ixx+npx : -ixx+npx;
338   }
339  
340   printf("i, iup, idown, j, jup, jdown %d %d %d %d %d %d \n",i,iup,idown,j,jup,jdown);
341
342   // calculate center of gravity
343   c.npoly++;
344   if (c.npoly > 100 ) {
345     printf("FindCluster - npoly >100,  npoly %d \n",c.npoly);
346     c.npoly=100;
347   }
348   c.xpoly[c.npoly-1]=x;
349   c.ypoly[c.npoly-1]=y;
350   c.zpoly[c.npoly-1]=zpos;
351
352   c.fX += q*x;
353   c.fY += q*y;
354   c.fQ += q;
355   
356   b.dig = 0;  b.idx = c.idx;
357
358   // left and right  
359   if (bins[i-1][j].dig) FindCluster(iChamber,segmentation,i-1,j,bins,c);
360   if (bins[i+1][j].dig) FindCluster(iChamber,segmentation,i+1,j,bins,c);
361   // up and down
362   if (bins[iup][jup].dig) FindCluster(iChamber,segmentation,iup,jup,bins,c);
363   if (bins[idown][jdown].dig) FindCluster(iChamber,segmentation,idown,jdown,bins,c);
364
365 }
366
367 //_____________________________________________________________________________
368
369 void AliMUONpoints::GetCenterOfGravity()
370 {
371   //
372   // simple MUON cluster finder from digits -- finds neighbours and 
373   // calculates center of gravity for the cluster
374   //
375
376   //  const Int_t MAX_Nipx=1026, MAX_Nipy=1026;
377
378   Bin bins[MAX_Nipx][MAX_Nipy]; 
379
380   AliMUONdisplay *display=(AliMUONdisplay*)gAlice->Display();
381   Int_t chamber=display->GetChamber();
382   Int_t cathode=display->GetCathode();
383    
384   AliMUON *MUON  = (AliMUON*)gAlice->GetModule("MUON");
385   AliMUONchamber *iChamber;
386   AliMUONsegmentation *segmentation;
387   iChamber =&(MUON->Chamber(chamber-1));
388   segmentation=iChamber->GetSegmentationModel(cathode);
389   Int_t npx  = segmentation->Npx();
390   Int_t npy  = segmentation->Npy();
391   Float_t zpos=iChamber->ZPosition();
392
393   TClonesArray *MUONdigits  = MUON->DigitsAddress(chamber-1);
394   gAlice->TreeD()->GetEvent(cathode);
395   Int_t ndigits = MUONdigits->GetEntriesFast();
396   if (fDigitIndex < 0 || fDigitIndex >= ndigits) return;
397
398   AliMUONdigit  *dig;
399   dig=(AliMUONdigit*)MUONdigits->UncheckedAt(fDigitIndex);
400   Int_t ipx=dig->fPadX;
401   Int_t ipy=dig->fPadY;
402   bins[ipx+npx][ipy+npy].dig=dig;
403     
404   int ndig;
405   int ncls=0;
406   for (ndig=0;ndig<ndigits;ndig++) {
407     dig = (AliMUONdigit*)MUONdigits->UncheckedAt(ndig);
408     int i=dig->fPadX, j=dig->fPadY;
409     bins[i+npx][j+npy].dig=dig;
410   }
411    PreCluster c; c.summit=bins[ipx+npx][ipy+npy].dig; c.idx=ncls;
412    FindCluster(iChamber,segmentation,ipx+npx, ipy+npy, bins, c);
413    printf("GetCenterOfGravity -- npoly %d \n",c.npoly);
414     
415    if (c.npeaks>1) {
416       printf("GetCenterOfGravity -- more than one peak");
417    }
418    c.fX /= c.fQ;
419    c.fY /= c.fQ;
420    printf("GetCenterOfGravity - c.fX c.fY c.fQ c.npeaks%f %f %d %d \n",c.fX,c.fY,c.fQ,c.npeaks);
421    /*
422    c.fTracks[0]=(Int_t)(*(c.summit->fTrks))(0);
423    c.fTracks[1]=(Int_t)(*(c.summit->fTrks))(1);
424    c.fTracks[2]=(Int_t)(*(c.summit->fTrks))(2);
425    */
426    c.fTracks[0]=c.summit->fTracks[0];
427    c.fTracks[1]=c.summit->fTracks[1];
428    c.fTracks[2]=c.summit->fTracks[2];
429    ncls++;
430    AliMUONpoints *points = 0;
431    points = new AliMUONpoints(1);
432    points->SetMarkerColor(kYellow);
433    points->SetMarkerStyle(5);
434    points->SetMarkerSize(1.);
435    points->SetPoint(0,c.fX,c.fY,zpos);
436    points->SetParticle(-1);
437    points->Draw();
438
439    TPolyLine3D *pline=0;
440    /*
441    pline=new TPolyLine3D(c.npoly);
442    Int_t np=c.npoly;
443    Float_t *x=new Float_t(c.npoly);
444    Float_t *y=new Float_t(c.npoly);
445    Float_t *z=new Float_t(c.npoly);
446    for (int i=0;i<np;i++) {
447      x[i]=c.xpoly[i];
448      y[i]=c.ypoly[i];
449      x[i]=c.zpoly[i];
450      pline->SetPoint(i,x[i],y[i],z[i]);
451    }
452    */
453    pline=new TPolyLine3D(c.npoly,c.xpoly,c.ypoly,c.zpoly);
454    pline->SetLineColor(kWhite);
455    pline->Draw();
456   
457
458    printf("GetCenterOfGravity -- ncls %d \n",ncls);
459
460 }
461