[u/mrichter/AliRoot.git] / MUON / AliMUONpoints.h.cxx

///////////////////////////////////////////////////////////////////////////////
//                                                                           //
//  This class contains the points for the ALICE event display               //
//                                                                           //
//Begin_Html
/*
<img src="gif/AliMUONpointsClass.gif">
*/
//End_Html
//                                                                           //
//                                                                           //
///////////////////////////////////////////////////////////////////////////////

#include "AliMUONpoints.h"
#include "AliMUONdisplay.h"
#include "AliRun.h"
#include "TPad.h"
#include "TVirtualPad.h"
#include "TPolyLine3D.h"
#include "TPaveText.h"
#include "TView.h"
#include "TMath.h"

//const Int_t MAX_Nipx=1026, MAX_Nipy=1026;
const Int_t MAX_Nipx=400, MAX_Nipy=800;
 
ClassImp(AliMUONpoints)

//_____________________________________________________________________________
AliMUONpoints::AliMUONpoints()
{
  //
  // Default constructor
  //
  fHitIndex = 0;
  fTrackIndex = 0;
  fDigitIndex = 0;
  fMarker[0] = fMarker[1] = fMarker[2]=0;
  fMatrix = 0;
}

//_____________________________________________________________________________
AliMUONpoints::AliMUONpoints(Int_t npoints)
  :AliPoints(npoints)
{
  //
  // Standard constructor
  //
  fHitIndex = 0;
  fTrackIndex = 0;
  fDigitIndex = 0;
  fMarker[0] = fMarker[1] = fMarker[2]=0;
  fMatrix = 0;
}
	 
//_____________________________________________________________________________
AliMUONpoints::~AliMUONpoints()
{
  //
  // Default destructor
  //
  fHitIndex = 0;
  fTrackIndex = 0;
  fDigitIndex = 0;
  for (Int_t i=0;i<3;i++){
      if (
fMarker[i]) delete fMarker[i];
  }
  fMatrix = 0;
}

//_____________________________________________________________________________
//void AliMUONpoints::ExecuteEvent(Int_t event, Int_t px, Int_t py)
//{
  //
  //*-*-*-*-*-*-*-*-*-*Execute action corresponding to one event*-*-*-*-*-*-*-*
  //*-*                =========================================
  //*-*
  //*-*  This member function must be implemented to realize the action
  //*-*  corresponding to the mouse click on the object in the window
  //*-*
  //*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
//  gPad->SetCursor(kCross);
  
//}
//_____________________________________________________________________________
void AliMUONpoints::DumpHit()
{
  //
  //   Dump hit corresponding to this point
  //
  AliMUONhit *hit = GetHit();
  if (hit) hit->Dump();
  printf("fTrackIndex %d \n",fTrackIndex);
}

//_____________________________________________________________________________
void AliMUONpoints::DumpDigit()
{
  //
  //   Dump digit corresponding to this point
  //
  AliMUONdigit *digit = GetDigit();
  if (digit) digit->Dump();
  for (int i=0;i<10;i++) {
        printf(" track charge %d %d \n",digit->fTracks[i],digit->fTcharges[i]);
  }

}

//_____________________________________________________________________________
void AliMUONpoints::InspectHit()
{
  //
  //   Inspect hit corresponding to this point
  //

  if (fHitIndex < 0 ) return;
  TVirtualPad *padsav = gPad;
  AliMUONhit *hit = GetHit();
  if (hit) hit->Inspect();
  TVirtualPad *padinspect = (TVirtualPad*)(gROOT->GetListOfCanvases())->FindObject("inspect");
   padinspect->cd();
   Float_t xmin = gPad->GetX1();
   Float_t xmax = gPad->GetX2();
   Float_t ymin = gPad->GetY1();
   Float_t ymax = gPad->GetY2();
   Float_t dy   = ymax-ymin;

      TPaveText *pad = new TPaveText(xmin, ymin+0.1*dy, xmax, ymin+0.15*dy);
      pad->SetBit(kCanDelete);
      pad->SetFillColor(42);
      pad->Draw();
      char ptitle[100];
      sprintf(ptitle," %s , fTrack: %d  fTrackIndex: %d ",GetName(),fIndex,fTrackIndex);
      pad->AddText(ptitle);
      padinspect->cd();
      padinspect->Update();
  if (padsav) padsav->cd();

}

//_____________________________________________________________________________
void AliMUONpoints::InspectDigit()
{
  //
  //   Inspect digit corresponding to this point
  //
  if (fDigitIndex < 0) return;
  TVirtualPad *padsav = gPad;
  AliMUONdigit *digit = GetDigit();
  if (digit) digit->Inspect();
  TVirtualPad *padinspect = (TVirtualPad*)(gROOT->GetListOfCanvases())->FindObject("inspect");
   padinspect->cd();
   Float_t xmin = gPad->GetX1();
   Float_t xmax = gPad->GetX2();
   Float_t ymin = gPad->GetY1();
   Float_t ymax = gPad->GetY2();
   Float_t dy   = ymax-ymin;

      TPaveText *pad = new TPaveText(xmin, ymin+0.1*dy, xmax, ymin+0.25*dy);
      pad->SetBit(kCanDelete);
      pad->SetFillColor(42);
      pad->Draw();
      char ptitle[11][100];
      //      sprintf(ptitle[11],"Tracks making this digit");
      //      pad->AddText(ptitle[11]);
  for (int i=0;i<10;i++) {
      if (digit->fTracks[i] == 0) continue;  
      sprintf(ptitle[i],"fTrackIndex: %d  Charge: %d",digit->fTracks[i],digit->fTcharges[i]);
      pad->AddText(ptitle[i]);
  }
      padinspect->cd();
      padinspect->Update();
  if (padsav) padsav->cd();
      
}

//_____________________________________________________________________________
Int_t AliMUONpoints::GetTrackIndex()
{
  //
  //   Dump digit corresponding to this point
  //
  printf("GetTrackIndex - fTrackIndex %d \n",fTrackIndex);
  this->Inspect();
  /*
  if (fDigitIndex != 0) {
    Int_t ncol=this->fMatrix->GetNcols();
    for (int i=0;i<ncol;i++) {
        printf(" track charge %f %f \n",(*(this->fMatrix))(0,i),(*(this->fMatrix))(1,i));
    }
  }
  */
  return fTrackIndex;
}

//_____________________________________________________________________________
AliMUONhit *AliMUONpoints::GetHit() const
{
  //
  //   Returns pointer to hit index in AliRun::fParticles
  //
  AliMUON *MUON  = (AliMUON*)gAlice->GetModule("MUON");
  gAlice->TreeH()->GetEvent(fTrackIndex);
  TClonesArray *MUONhits  = MUON->Hits();
  Int_t nhits = MUONhits->GetEntriesFast();
  if (fHitIndex < 0 || fHitIndex >= nhits) return 0;
  return (AliMUONhit*)MUONhits->UncheckedAt(fHitIndex);
}

//_____________________________________________________________________________
AliMUONdigit *AliMUONpoints::GetDigit() const
{
  //
  //   Returns pointer to digit index in AliRun::fParticles
  //

  AliMUONdisplay *display=(AliMUONdisplay*)gAlice->Display();
  Int_t chamber=display->GetChamber();
  Int_t cathode=display->GetCathode();
   
  AliMUON *MUON  = (AliMUON*)gAlice->GetModule("MUON");
  TClonesArray *MUONdigits  = MUON->DigitsAddress(chamber-1);
  gAlice->TreeD()->GetEvent(cathode);
  Int_t ndigits = MUONdigits->GetEntriesFast();
  if (fDigitIndex < 0 || fDigitIndex >= ndigits) return 0;
  return (AliMUONdigit*)MUONdigits->UncheckedAt(fDigitIndex);
}
//_____________________________________________________________________________
struct Bin {
   const AliMUONdigit *dig;
   int idx;
   Bin() {dig=0; idx=-1;}
};

struct PreCluster : public AliMUONreccluster {
   const AliMUONdigit* summit;
   int idx;
   int cut;
   int npeaks;
   int npoly;
   float xpoly[100];
   float ypoly[100];
   float zpoly[100];
   PreCluster() : AliMUONreccluster() {cut=npeaks=npoly=0; 
                                       for (int k=0;k<100;k++) {
                                         xpoly[k]=ypoly[k]=zpoly[k]=0;
				       }
   }
                              
};
//_____________________________________________________________________________

static void FindCluster(AliMUONchamber *iChamber, AliMUONsegmentation *segmentation, int i, int j, Bin bins[][MAX_Nipy], PreCluster &c) 

{

  //
  // Find clusters
  //

  printf("I'm in FindCluster \n"); 

  Bin& b=bins[i][j];
  Int_t q=b.dig->fSignal;

  printf("FindCluster - i j q %d %d %d\n",i,j,q);
  
  if (q<0) { 
    q=-q;
    c.cut=1;
  } 
  //  if (b.idx >= 0 && b.idx != c.idx) {
  if (b.idx >= 0 && b.idx > c.idx) {
    c.idx=b.idx;
    c.npeaks++;
    printf("FindCluster - one more peak \n");
  }
  
  if (q > TMath::Abs(c.summit->fSignal)) c.summit=b.dig;

  Float_t zpos=iChamber->ZPosition();  // check with Andreas

  Int_t npx  = segmentation->Npx();
  Int_t npy  = segmentation->Npy();

  // get pad coordinates and prepare the up and down steps   
  Int_t jup  =(j-npy > 0) ? j+1 : (j-npy-1)+npy;
  Int_t jdown=(j-npy > 0) ? j-1 : (j-npy+1)+npy;
  //  Int_t jup  =j+1;
  //  Int_t jdown=j-1;
  Float_t x, y;
  segmentation->GetPadCxy(i-npx, j-npy,x, y);
  Int_t isec0=segmentation->Sector(i-npx,j-npy);
  printf("FindCluster - i-npx j-npy isec0 %d %d %d \n",i-npx,j-npy,isec0);
  //  printf("FindCluster - x  y %f %f \n",x,y);

  Float_t dpy  = segmentation->Dpy();
  Float_t dpx  = segmentation->Dpx()/16;
  Int_t ixx, iyy;
  Float_t absx=TMath::Abs(x);
  // iup
  (y >0) ? segmentation->GetPadIxy(absx+dpx,y+dpy,ixx,iyy) : segmentation->GetPadIxy(absx+dpx,y-dpy,ixx,iyy);
  printf(" iup: ixx iyy %d %d \n",ixx,iyy);
  Int_t jtest=TMath::Abs(iyy)-npy-1;
  if (j != jtest) {
     printf(" j != jtest - something's wrong %d %d \n",j,jtest);
  }
  Int_t iup=(x >0) ? ixx+npx : -ixx+npx;
  // idown
  (y >0) ? segmentation->GetPadIxy(absx+dpx,y-dpy,ixx,iyy) : segmentation->GetPadIxy(absx+dpx,y+dpy,ixx,iyy);
  printf(" idown: ixx iyy %d %d \n",ixx,iyy);
  Int_t idown=(x >0) ? ixx+npx : -ixx+npx;
  if (bins[idown][jdown].dig == 0) {
     (y >0) ? segmentation->GetPadIxy(absx-dpx,y-dpy,ixx,iyy) : segmentation->GetPadIxy(absx-dpx,y+dpy,ixx,iyy);
     printf(" idown: ixx iyy %d %d \n",ixx,iyy);
     idown=(x >0) ? ixx+npx : -ixx+npx;
  }
 
  printf("i, iup, idown, j, jup, jdown %d %d %d %d %d %d \n",i,iup,idown,j,jup,jdown);

  // calculate center of gravity
  c.npoly++;
  if (c.npoly > 100 ) {
    printf("FindCluster - npoly >100,  npoly %d \n",c.npoly);
    c.npoly=100;
  }
  c.xpoly[c.npoly-1]=x;
  c.ypoly[c.npoly-1]=y;
  c.zpoly[c.npoly-1]=zpos;

  c.fX += q*x;
  c.fY += q*y;
  c.fQ += q;
  
  b.dig = 0;  b.idx = c.idx;

  // left and right  
  if (bins[i-1][j].dig) FindCluster(iChamber,segmentation,i-1,j,bins,c);
  if (bins[i+1][j].dig) FindCluster(iChamber,segmentation,i+1,j,bins,c);
  // up and down
  if (bins[iup][jup].dig) FindCluster(iChamber,segmentation,iup,jup,bins,c);
  if (bins[idown][jdown].dig) FindCluster(iChamber,segmentation,idown,jdown,bins,c);

}

//_____________________________________________________________________________

void AliMUONpoints::GetCenterOfGravity()
{
  //
  // simple MUON cluster finder from digits -- finds neighbours and 
  // calculates center of gravity for the cluster
  //

  //  const Int_t MAX_Nipx=1026, MAX_Nipy=1026;

  Bin bins[MAX_Nipx][MAX_Nipy]; 

  AliMUONdisplay *display=(AliMUONdisplay*)gAlice->Display();
  Int_t chamber=display->GetChamber();
  Int_t cathode=display->GetCathode();
   
  AliMUON *MUON  = (AliMUON*)gAlice->GetModule("MUON");
  AliMUONchamber *iChamber;
  AliMUONsegmentation *segmentation;
  iChamber =&(MUON->Chamber(chamber-1));
  segmentation=iChamber->GetSegmentationModel(cathode);
  Int_t npx  = segmentation->Npx();
  Int_t npy  = segmentation->Npy();
  Float_t zpos=iChamber->ZPosition();

  TClonesArray *MUONdigits  = MUON->DigitsAddress(chamber-1);
  gAlice->TreeD()->GetEvent(cathode);
  Int_t ndigits = MUONdigits->GetEntriesFast();
  if (fDigitIndex < 0 || fDigitIndex >= ndigits) return;

  AliMUONdigit  *dig;
  dig=(AliMUONdigit*)MUONdigits->UncheckedAt(fDigitIndex);
  Int_t ipx=dig->fPadX;
  Int_t ipy=dig->fPadY;
  bins[ipx+npx][ipy+npy].dig=dig;
    
  int ndig;
  int ncls=0;
  for (ndig=0;ndig<ndigits;ndig++) {
    dig = (AliMUONdigit*)MUONdigits->UncheckedAt(ndig);
    int i=dig->fPadX, j=dig->fPadY;
    bins[i+npx][j+npy].dig=dig;
  }
   PreCluster c; c.summit=bins[ipx+npx][ipy+npy].dig; c.idx=ncls;
   FindCluster(iChamber,segmentation,ipx+npx, ipy+npy, bins, c);
   printf("GetCenterOfGravity -- npoly %d \n",c.npoly);
    
   if (c.npeaks>1) {
      printf("GetCenterOfGravity -- more than one peak");
   }
   c.fX /= c.fQ;
   c.fY /= c.fQ;
   printf("GetCenterOfGravity - c.fX c.fY c.fQ c.npeaks%f %f %d %d \n",c.fX,c.fY,c.fQ,c.npeaks);
   /*
   c.fTracks[0]=(Int_t)(*(c.summit->fTrks))(0);
   c.fTracks[1]=(Int_t)(*(c.summit->fTrks))(1);
   c.fTracks[2]=(Int_t)(*(c.summit->fTrks))(2);
   */
   c.fTracks[0]=c.summit->fTracks[0];
   c.fTracks[1]=c.summit->fTracks[1];
   c.fTracks[2]=c.summit->fTracks[2];
   ncls++;
   AliMUONpoints *points = 0;
   points = new AliMUONpoints(1);
   points->SetMarkerColor(kYellow);
   points->SetMarkerStyle(5);
   points->SetMarkerSize(1.);
   points->SetPoint(0,c.fX,c.fY,zpos);
   points->SetParticle(-1);
   points->Draw();

   TPolyLine3D *pline=0;
   /*
   pline=new TPolyLine3D(c.npoly);
   Int_t np=c.npoly;
   Float_t *x=new Float_t(c.npoly);
   Float_t *y=new Float_t(c.npoly);
   Float_t *z=new Float_t(c.npoly);
   for (int i=0;i<np;i++) {
     x[i]=c.xpoly[i];
     y[i]=c.ypoly[i];
     x[i]=c.zpoly[i];
     pline->SetPoint(i,x[i],y[i],z[i]);
   }
   */
   pline=new TPolyLine3D(c.npoly,c.xpoly,c.ypoly,c.zpoly);
   pline->SetLineColor(kWhite);
   pline->Draw();
  

   printf("GetCenterOfGravity -- ncls %d \n",ncls);

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