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

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

/*
$Log$
*/

///////////////////////////////////////////////////////////////////////////////
//                                                                           //
//  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();
}

//_____________________________________________________________________________
void AliMUONpoints::DumpDigit()
{
  //
  //   Dump digit corresponding to this point
  //
  AliMUONdigit *digit = GetDigit();
  if (digit) digit->Dump();
}

//_____________________________________________________________________________
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
  //

  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);
  Int_t nent=(Int_t)gAlice->TreeD()->GetEntries();
  gAlice->TreeD()->GetEvent(nent-2+cathode-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
  //


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

  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++;
  }
  
  if (q > TMath::Abs(c.summit->fSignal)) c.summit=b.dig;

  Float_t zpos=iChamber->ZPosition();

  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-1;
  Int_t jdown=(j-npy > 0) ? j-1 : j+1;

  Float_t x, y;
  segmentation->GetPadCxy(i-npx, j-npy,x, y);
  Int_t isec0=segmentation->Sector(i-npx,j-npy);

  Float_t dpy  = segmentation->Dpy(isec0);
  Float_t dpx  = segmentation->Dpx(isec0)/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);

  //  Int_t jtest=TMath::Abs(iyy)-npy-1;

  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);

  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);
     idown=(x >0) ? ixx+npx : -ixx+npx;
  }
 

  // 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);
  Int_t nent=(Int_t)gAlice->TreeD()->GetEntries();
  gAlice->TreeD()->GetEvent(nent-2+cathode-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);
    
   if (c.npeaks>1) {
      printf("GetCenterOfGravity -- more than one peak");
   }
   c.fX /= c.fQ;
   c.fY /= c.fQ;

   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;
   TVector *xp=new TVector(c.npoly);
   TVector *yp=new TVector(c.npoly);
   TVector *zp=new TVector(c.npoly);
   for (int i=0;i<np;i++) {
     (*xp)(i)=c.xpoly[i];
     (*yp)(i)=c.ypoly[i];
     (*zp)(i)=c.zpoly[i];
     pline->SetPoint(i,(*xp)(i),(*yp)(i),(*zp)(i));
     //printf("np, i, xp, yp, zp %d %d %f %f %f \n",np,i,(*xp)(i),(*yp)(i),(*zp)(i));
   }
   */
   pline=new TPolyLine3D(c.npoly,c.xpoly,c.ypoly,c.zpoly);
   pline->SetLineColor(kWhite);
   pline->Draw();
   /*
   delete xp;
   delete yp;
   delete zp;
   */  
   for (int k=0;k<c.npoly;k++) {
     c.xpoly[k]=c.ypoly[k]=c.zpoly[k]=0;
   }
   c.npoly=0;


}
Commit	Line	Data
4c039060	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
fe4da5cc	20	///////////////////////////////////////////////////////////////////////////////
	21	// //
	22	// This class contains the points for the ALICE event display //
	23	// //
	24	//Begin_Html
	25	/*
a897a37a	26	<img src="gif/AliMUONpointsClass.gif">
fe4da5cc	27	*/
	28	//End_Html
	29	// //
	30	// //
	31	///////////////////////////////////////////////////////////////////////////////
a897a37a	32
fe4da5cc	33	#include "AliMUONpoints.h"
a897a37a	34	#include "AliMUONdisplay.h"
fe4da5cc	35	#include "AliRun.h"
fe4da5cc	36	#include "TPad.h"
a897a37a	37	#include "TVirtualPad.h"
	38	#include "TPolyLine3D.h"
	39	#include "TPaveText.h"
fe4da5cc	40	#include "TView.h"
	41	#include "TMath.h"
	42
a897a37a	43	//const Int_t MAX_Nipx=1026, MAX_Nipy=1026;
fe4da5cc	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;
a897a37a	57	fMarker[0] = fMarker[1] = fMarker[2]=0;
a897a37a	58	fMatrix = 0;
fe4da5cc	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;
a897a37a	71	fMarker[0] = fMarker[1] = fMarker[2]=0;
a897a37a	72	fMatrix = 0;
fe4da5cc	73	}
	74
	75	//_____________________________________________________________________________
	76	AliMUONpoints::~AliMUONpoints()
	77	{
	78	//
	79	// Default destructor
	80	//
	81	fHitIndex = 0;
	82	fTrackIndex = 0;
	83	fDigitIndex = 0;
a897a37a	84	for (Int_t i=0;i<3;i++){
8eb67615	85	if (fMarker[i]) delete fMarker[i];
a897a37a	86	}
a897a37a	87	fMatrix = 0;
fe4da5cc	88	}
fe4da5cc	89
a897a37a	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	//}
fe4da5cc	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	//
a897a37a	130
	131	if (fHitIndex < 0 ) return;
	132	TVirtualPad *padsav = gPad;
fe4da5cc	133	AliMUONhit *hit = GetHit();
fe4da5cc	134	if (hit) hit->Inspect();
a897a37a	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.1dy, 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
fe4da5cc	154	}
	155
	156	//_____________________________________________________________________________
	157	void AliMUONpoints::InspectDigit()
	158	{
	159	//
	160	// Inspect digit corresponding to this point
	161	//
a897a37a	162	if (fDigitIndex < 0) return;
a897a37a	163	TVirtualPad *padsav = gPad;
fe4da5cc	164	AliMUONdigit *digit = GetDigit();
fe4da5cc	165	if (digit) digit->Inspect();
a897a37a	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.1dy, 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
fe4da5cc	190	}
	191
	192	//_____________________________________________________________________________
	193	Int_t AliMUONpoints::GetTrackIndex()
	194	{
	195	//
	196	// Dump digit corresponding to this point
	197	//
8eb67615	198
fe4da5cc	199	this->Inspect();
a897a37a	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	*/
fe4da5cc	208	return fTrackIndex;
	209	}
	210
	211	//_____________________________________________________________________________
	212	AliMUONhit *AliMUONpoints::GetHit() const
	213	{
	214	//
	215	// Returns pointer to hit index in AliRun::fParticles
	216	//
0f439796	217	AliMUON MUON = (AliMUON)gAlice->GetModule("MUON");
fe4da5cc	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
0f439796	236	AliMUON MUON = (AliMUON)gAlice->GetModule("MUON");
fe4da5cc	237	TClonesArray *MUONdigits = MUON->DigitsAddress(chamber-1);
a897a37a	238	Int_t nent=(Int_t)gAlice->TreeD()->GetEntries();
	239	gAlice->TreeD()->GetEvent(nent-2+cathode-1);
	240	//gAlice->TreeD()->GetEvent(cathode);
fe4da5cc	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;
a897a37a	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
fe4da5cc	267	};
	268	//_____________________________________________________________________________
	269
a897a37a	270	static void FindCluster(AliMUONchamber iChamber, AliMUONsegmentation segmentation, int i, int j, Bin bins[][MAX_Nipy], PreCluster &c)
fe4da5cc	271
	272	{
	273
	274	//
	275	// Find clusters
	276	//
	277
fe4da5cc	278
	279	Bin& b=bins[i][j];
	280	Int_t q=b.dig->fSignal;
	281
fe4da5cc	282	if (q<0) {
	283	q=-q;
	284	c.cut=1;
	285	}
a897a37a	286	// if (b.idx >= 0 && b.idx != c.idx) {
a897a37a	287	if (b.idx >= 0 && b.idx > c.idx) {
fe4da5cc	288	c.idx=b.idx;
	289	c.npeaks++;
	290	}
	291
	292	if (q > TMath::Abs(c.summit->fSignal)) c.summit=b.dig;
	293
8eb67615	294	Float_t zpos=iChamber->ZPosition();
a897a37a	295
fe4da5cc	296	Int_t npx = segmentation->Npx();
fe4da5cc	297	Int_t npy = segmentation->Npy();
a897a37a	298
fe4da5cc	299
a897a37a	300	// get pad coordinates and prepare the up and down steps
8eb67615	301	Int_t jup =(j-npy > 0) ? j+1 : j-1;
	302	Int_t jdown=(j-npy > 0) ? j-1 : j+1;
	303
a897a37a	304	Float_t x, y;
	305	segmentation->GetPadCxy(i-npx, j-npy,x, y);
	306	Int_t isec0=segmentation->Sector(i-npx,j-npy);
a897a37a	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);
8eb67615	314
	315	// Int_t jtest=TMath::Abs(iyy)-npy-1;
	316
a897a37a	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);
8eb67615	320
a897a37a	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);
a897a37a	324	idown=(x >0) ? ixx+npx : -ixx+npx;
	325	}
	326
a897a37a	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;
fe4da5cc	337
	338	c.fX += q*x;
	339	c.fY += q*y;
	340	c.fQ += q;
	341
	342	b.dig = 0; b.idx = c.idx;
a897a37a	343
a897a37a	344	// left and right
fe4da5cc	345	if (bins[i-1][j].dig) FindCluster(iChamber,segmentation,i-1,j,bins,c);
fe4da5cc	346	if (bins[i+1][j].dig) FindCluster(iChamber,segmentation,i+1,j,bins,c);
a897a37a	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);
fe4da5cc	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	//
a897a37a	361
	362	// const Int_t MAX_Nipx=1026, MAX_Nipy=1026;
	363
	364	Bin bins[MAX_Nipx][MAX_Nipy];
	365
fe4da5cc	366	AliMUONdisplay display=(AliMUONdisplay)gAlice->Display();
	367	Int_t chamber=display->GetChamber();
	368	Int_t cathode=display->GetCathode();
	369
0f439796	370	AliMUON MUON = (AliMUON)gAlice->GetModule("MUON");
fe4da5cc	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();
a897a37a	378
fe4da5cc	379	TClonesArray *MUONdigits = MUON->DigitsAddress(chamber-1);
a897a37a	380	Int_t nent=(Int_t)gAlice->TreeD()->GetEntries();
	381	gAlice->TreeD()->GetEvent(nent-2+cathode-1);
	382	//gAlice->TreeD()->GetEvent(cathode);
fe4da5cc	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;
fe4da5cc	390	bins[ipx+npx][ipy+npy].dig=dig;
	391
	392	int ndig;
	393	int ncls=0;
a897a37a	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;
fe4da5cc	398	}
a897a37a	399	PreCluster c; c.summit=bins[ipx+npx][ipy+npy].dig; c.idx=ncls;
a897a37a	400	FindCluster(iChamber,segmentation,ipx+npx, ipy+npy, bins, c);
a897a37a	401
a897a37a	402	if (c.npeaks>1) {
fe4da5cc	403	printf("GetCenterOfGravity -- more than one peak");
a897a37a	404	}
	405	c.fX /= c.fQ;
	406	c.fY /= c.fQ;
8eb67615	407
a897a37a	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();
fe4da5cc	420
a897a37a	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;
fe4da5cc	448
fe4da5cc	449
a897a37a	450	}
fe4da5cc	451