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

/////////////////////////////////////////////////////
//  Segmentation and Response classes version 01   //
/////////////////////////////////////////////////////

#include <TTUBE.h>
#include <TNode.h> 
#include <TRandom.h> 

#include "AliMUONv01.h"
#include "AliRun.h"
#include "AliMC.h"
#include "iostream.h"

//___________________________________________
ClassImp(AliMUONsegmentationV01)

AliMUONsegmentationV01::AliMUONsegmentationV01() 
{
    fNsec=4;
    fRSec.Set(fNsec);    
    fNDiv.Set(fNsec);      
    fDpxD.Set(fNsec);      
}
    
void   AliMUONsegmentationV01::SetSegRadii(Float_t  r[4])
{
    for (Int_t i=0; i<4; i++) {
	fRSec[i]=r[i];
	printf("\n R %d %f",i,fRSec[i]);
	
    }
}

void AliMUONsegmentationV01::SetPadDivision(Int_t ndiv[4])
{
//
// Defines the pad size perp. to the anode wire (y) for different sectors. 
//
    for (Int_t i=0; i<4; i++) {
	fNDiv[i]=ndiv[i];
	printf("\n Ndiv %d %d",i,fNDiv[i]);
    }
    ndiv[0]=ndiv[1];
}


void AliMUONsegmentationV01::Init(AliMUONchamber*)
{
//
//  Fill the arrays fCx (x-contour) and fNpx (ix-contour) for each sector
//  These arrays help in converting from real to pad co-ordinates and
//  vice versa
//
    Int_t isec;
    printf("\n Initialise segmentation v01");
    fNpy=Int_t(fRSec[fNsec-1]/fDpy)+1;

    fDpxD[fNsec-1]=fDpx;
    if (fNsec > 1) {
	for (Int_t i=fNsec-2; i>=0; i--){
	    fDpxD[i]=fDpxD[fNsec-1]/fNDiv[i];
	    printf("\n dx %d %f",i,fDpxD[i]);
	}
    }
    fWireD=fDpxD[1]/3;
//
// fill the arrays defining the pad segmentation boundaries
    Float_t ry;
    Int_t   dnx;
//
//  loop over sections
    for(isec=0; isec<fNsec; isec++) {
//  
//  loop over pads along the aode wires
	for (Int_t iy=1; iy<=fNpy; iy++) {
//
	    Float_t x=iy*fDpy-fDpy/2;
	    if (x > fRSec[isec]) {
		fNpx[isec][iy]=0;
		 fCx[isec][iy]=0;
	    } else {
		ry=TMath::Sqrt(fRSec[isec]*fRSec[isec]-x*x);
		if (isec > 0) {
		    dnx= Int_t((ry-fCx[isec-1][iy])/fDpxD[isec]);
		    if (TMath::Odd(dnx)) dnx--;
		    fNpx[isec][iy]=fNpx[isec-1][iy]+dnx;
		    fCx[isec][iy]=fCx[isec-1][iy]+dnx*fDpxD[isec];
		} else {
		    dnx=Int_t(ry/fDpxD[isec]);
		    fNpx[isec][iy]=dnx;
		    if (TMath::Odd(dnx)) dnx--;
		    fCx[isec][iy]=dnx*fDpxD[isec];
		}
	    }
	} // y-pad loop
    } // sector loop
 
    //   
    // for debugging only 
    for (Int_t iy=0; iy<fNpy; iy++) {
	printf("\n iy %d",iy);
	for(isec=0; isec<fNsec; isec++) {
	    printf("  %d", 
		   fNpx[isec][iy]);
	}
	printf("\n iy %d",iy);
	for(isec=0; isec<fNsec; isec++) {
	    printf("  %f", 
		   fCx[isec][iy] );
	}
	printf("\n");
    }
    Float_t x,y;
    Int_t jx,jy;
    //   
    // for debugging only 
    for (Int_t ix=1; ix<100; ix++) {
	for (Int_t iy=1; iy<100; iy++) {
	    GetPadCxy(ix,iy,x,y);
	    GetPadIxy(x,y,jx,jy);
	    if ((ix != jx) && (jx!=-1)) {
		printf("\n %d %f %d %f", ix,x,iy,y);
		printf("\n %d %f %d %f", jx,x,jy,y);
		printf("\n \n **********");
	    }
	}
    }
}

Int_t AliMUONsegmentationV01::Sector(Int_t ix, Int_t iy)
{
    Int_t absix=TMath::Abs(ix);
    Int_t absiy=TMath::Abs(iy);
    Int_t isec=0;
    for (Int_t i=0; i<fNsec; i++) {
	if (absix<=fNpx[i][absiy]){
	    isec=i;
	    break;
	}
    }
    return isec;
}


Float_t AliMUONsegmentationV01::GetAnod(Float_t xhit)
{
//
// Finds anod-wire position closest to xhit
    Float_t wire= (xhit<0)? Int_t(xhit/fWireD)+0.5:Int_t(xhit/fWireD)-0.5;
//    printf("getanod: %f %f %f ",xhit, wire, fWireD*wire);
    return fWireD*wire;
}

    void AliMUONsegmentationV01::
    GetPadIxy(Float_t x, Float_t y, Int_t &ix, Int_t &iy)
{
//  returns pad coordinates (ix,iy) for given real coordinates (x,y)
//
    iy = (y>0)? Int_t(y/fDpy)+1 : Int_t(y/fDpy)-1;
    if (iy >  fNpy) iy= fNpy;
    if (iy < -fNpy) iy=-fNpy;
//
//  Find sector isec
    Int_t isec=-1;
    Float_t absx=TMath::Abs(x);
    Int_t absiy=TMath::Abs(iy);
    for (Int_t i=0; isec < fNsec; i++) {
	if (absx <= fCx[i][absiy]) {
	    isec=i;
	    break;
	}
    }
    if (isec>0) {
	ix= Int_t((absx-fCx[isec-1][absiy])/fDpxD[isec])
	    +fNpx[isec-1][absiy]+1;
    } else if (isec == 0) {
	ix= Int_t(absx/fDpxD[isec])+1;
    } else {
	ix=fNpx[fNsec-1][absiy]+1;	
    }
//    printf("\n %d %d",isec,absiy);
    
    ix = (x>0) ? ix:-ix;
}

void AliMUONsegmentationV01::
GetPadCxy(Int_t ix, Int_t iy, Float_t &x, Float_t &y)
{
//  returns real coordinates (x,y) for given pad coordinates (ix,iy)
//
    y = (iy>0) ? Float_t(iy*fDpy)-fDpy/2. : Float_t(iy*fDpy)+fDpy/2.;
//
//  Find sector isec
    Int_t isec=Sector(ix,iy);
//
    Int_t absix=TMath::Abs(ix);
    Int_t absiy=TMath::Abs(iy);
    if (isec) {
	x=fCx[isec-1][absiy]+(absix-fNpx[isec-1][absiy])*fDpxD[isec];
	x=(ix>0) ?  x-fDpxD[isec]/2 : -x+fDpxD[isec]/2;
    } else {
	x=y=0;
    }
}


void AliMUONsegmentationV01::
GetSuperPadIxy(Float_t x, Float_t y, Int_t &ix, Int_t &iy)
{
    ix = (x>0)? Int_t(x/fDpx)+1 : Int_t(x/fDpx)-1;
    iy = (y>0)? Int_t(y/fDpy)+1 : Int_t(y/fDpy)-1;
}

void AliMUONsegmentationV01::
GetSuperPadCxy(Int_t ix, Int_t iy, Float_t &x, Float_t &y) 
{
    x = (ix>0) ? Float_t(ix*fDpx)-fDpx/2. : Float_t(ix*fDpx)+fDpx/2.;
    y = (iy>0) ? Float_t(iy*fDpy)-fDpy/2. : Float_t(iy*fDpy)+fDpy/2.;
}

void AliMUONsegmentationV01::SetPADSIZ(Float_t p1, Float_t p2)
{
    fDpx=p1;
    fDpy=p2;
}

void AliMUONsegmentationV01::
FirstPad(Float_t xhit, Float_t yhit, Float_t dx, Float_t dy)
{
    //
    // Find the wire position (center of charge distribution)
    Float_t x0a=GetAnod(xhit);
    //
    // and take fNsigma*sigma around this center
    Float_t x01=x0a  - dx;
    Float_t x02=x0a  + dx;
    Float_t y01=yhit - dy;
    Float_t y02=yhit + dy;
    //
    // find the pads over which the charge distributes
    GetPadIxy(x01,y01,fixmin,fiymin);
    GetPadIxy(x02,y02,fixmax,fiymax);
    fxmin=x01;
    fxmax=x02;
    //
    // upper and lower limits should be checked
    /*
    printf("\n FirstPad called");
    printf("\n Hit Position %f %f",xhit,yhit);
    printf("\n Closest wire %f", x0a);
    printf("\n Integration limits: %i %i %i %i",fixmin,fixmax,fiymin,fiymax);
    printf("\n Integration limits: %f %f %f %f \n",x01,x02,y01,y02);
    */
    // 
    // Set current pad to lower left corner
    if (fixmax < fixmin) fixmax=fixmin;
    if (fiymax < fiymin) fiymax=fiymin;    
    fix=fixmin;
    fiy=fiymin;
    GetPadCxy(fix,fiy,fx,fy);
    
}

void AliMUONsegmentationV01::NextPad()
{
  // 
  // Step to next pad in integration region
    Float_t xc,yc;
    Int_t   iyc;
    
//  step from left to right    
    if (fx < fxmax && fx != 0) {
	fix++;
//  step up 
    } else if (fiy != fiymax) {
	fiy++;
//      get y-position of next row (yc), xc not used here 	
	GetPadCxy(fix,fiy,xc,yc);
//      get x-pad coordiante for 1 pad in row (fix)
	GetPadIxy(fxmin,yc,fix,iyc);
    } else {
	printf("\n Error: Stepping outside integration region\n ");
    }
    GetPadCxy(fix,fiy,fx,fy);
    fSector=Sector(fix,fiy);
//    printf("\n this pad %f %f %d %d",fx,fy,fix,fiy);
    
}

Int_t AliMUONsegmentationV01::MorePads()
//
// Are there more pads in the integration region
{
    if ((fx >= fxmax  && fiy >= fiymax) || fy==0) {
	return 0;
    } else {
	return 1;
    }
}

void AliMUONsegmentationV01::SigGenInit(Float_t x,Float_t y,Float_t)
{
//
//  Initialises pad and wire position during stepping
    fxt=x;
    fyt=y;
    GetPadIxy(x,y,fixt,fiyt);
    fiwt= (x>0) ? Int_t(x/fWireD)+1 : Int_t(x/fWireD)-1 ;
}

Int_t AliMUONsegmentationV01::SigGenCond(Float_t x,Float_t y,Float_t)
{
//
//  Signal will be generated if particle crosses pad boundary or
//  boundary between two wires. 
    Int_t ixt, iyt;
    GetPadIxy(x,y,ixt,iyt);

    Int_t iwt=(x>0) ? Int_t(x/fWireD)+1 : Int_t(x/fWireD)-1;
    
    if ((ixt != fixt) || (iyt !=fiyt) || (iwt != fiwt)) {
	return 1;
    } else {
	return 0;
    }
}

void AliMUONsegmentationV01::
IntegrationLimits(Float_t& x1,Float_t& x2,Float_t& y1, Float_t& y2)
{
    x1=fxt-fx-fDpxD[fSector]/2.;
    x2=x1+fDpxD[fSector];
    y1=fyt-fy-fDpy/2.;
    y2=y1+fDpy;    
}

void AliMUONsegmentationV01::
Neighbours(Int_t iX, Int_t iY, Int_t* Nlist, Int_t Xlist[10], Int_t Ylist[10])
{
    const Float_t epsilon=fDpy/1000;
    
    Float_t x,y;
    Int_t   ixx, iyy, isec1;
//
    Int_t   isec0=Sector(iX,iY);
    Int_t i=0;
//    
//  step right
    Xlist[i]=iX+1;
    Ylist[i++]=iY;
//
//  step left    
    Xlist[i]=iX-1;
    Ylist[i++]=iY;
//
//  step up 
    GetPadCxy(iX,iY,x,y);
    GetPadIxy(x+epsilon,y+fDpy,ixx,iyy);
    Xlist[i]=ixx;
    Ylist[i++]=iY+1;
    isec1=Sector(ixx,iyy);
    if (isec1==isec0) {
//
//  no sector boundary crossing
    Xlist[i]=ixx+1;
    Ylist[i++]=iY+1;
	
    Xlist[i]=ixx-1;
    Ylist[i++]=iY+1;
    } else if (isec1 < isec0) {
// finer segmentation
	Xlist[i]=ixx+1;
	Ylist[i++]=iY+1;
	
	Xlist[i]=ixx-1;
	Ylist[i++]=iY+1;
	
	Xlist[i]=ixx-2;
	Ylist[i++]=iY+1;
    } else {
// coarser segmenation	
	if (TMath::Odd(iX-fNpx[isec1-1][iY+1])) {
	    Xlist[i]=ixx-1;
	    Ylist[i++]=iY+1;
	} else {
	    Xlist[i]=ixx+1;
	    Ylist[i++]=iY+1;
	}
    }
//
// step down 
    GetPadCxy(iX,iY,x,y);
    GetPadIxy(x+epsilon,y-fDpy,ixx,iyy);
    Xlist[i]=ixx;
    Ylist[i++]=iY-1;
    isec1=Sector(ixx,iyy);
    if (isec1==isec0) {
//
//  no sector boundary crossing
    Xlist[i]=ixx+1;
    Ylist[i++]=iY-1;
	
    Xlist[i]=ixx-1;
    Ylist[i++]=iY-1;
    } else if (isec1 < isec0) {
// finer segmentation
	Xlist[i]=ixx+1;
	Ylist[i++]=iY-1;
	
	Xlist[i]=ixx-1;
	Ylist[i++]=iY-1;
	
	Xlist[i]=ixx-2;
	Ylist[i++]=iY-1;
    } else {
// coarser segmenation	
	if (TMath::Odd(iX-fNpx[isec1-1][iY-1])) {
	    Xlist[i]=ixx-1;
	    Ylist[i++]=iY-1;
	} else {
	    Xlist[i]=ixx+1;
	    Ylist[i++]=iY-1;
	}
    }
    *Nlist=i;
}

//___________________________________________
void AliMUONsegmentationV01::
FitXY(AliMUONRecCluster* Cluster,TClonesArray* MUONdigits)
    // Default : Centre of gravity method
{
    Float_t x=0;
    Float_t y=0;
    Float_t q=0;
    Float_t xToAdd;
    Float_t yToAdd;
    
    if (gAlice->TreeD()->GetReadEvent() != Cluster->GetCathod()+1)
	// next line warns if in the future cathod 1 is not event 2 !
	printf("ClusterFillXY : not reading the right cathod !\n");
    for(Int_t clusDigit=Cluster->FirstDigitIndex();
	clusDigit!=Cluster->InvalidDigitIndex();
	clusDigit=Cluster->NextDigitIndex()) {
	AliMUONdigit* pDigit=(AliMUONdigit*)MUONdigits->UncheckedAt(clusDigit);
	GetPadCxy(pDigit->fPadX,pDigit->fPadY,xToAdd,yToAdd);
	x+= xToAdd*pDigit->fSignal;
	y+= yToAdd*pDigit->fSignal;
	q+= (Float_t) pDigit->fSignal;
    }
    Cluster->fX=x/q;
    Cluster->fY=y/q;
}
Commit	Line	Data
fe4da5cc	1	/////////////////////////////////////////////////////
	2	// Segmentation and Response classes version 01 //
	3	/////////////////////////////////////////////////////
	4
	5	#include <TTUBE.h>
	6	#include <TNode.h>
	7	#include <TRandom.h>
	8
	9	#include "AliMUONv01.h"
	10	#include "AliRun.h"
	11	#include "AliMC.h"
	12	#include "iostream.h"
	13
	14	//___________________________________________
	15	ClassImp(AliMUONsegmentationV01)
	16
	17	AliMUONsegmentationV01::AliMUONsegmentationV01()
	18	{
	19	fNsec=4;
	20	fRSec.Set(fNsec);
	21	fNDiv.Set(fNsec);
	22	fDpxD.Set(fNsec);
	23	}
	24
	25	void AliMUONsegmentationV01::SetSegRadii(Float_t r[4])
	26	{
	27	for (Int_t i=0; i<4; i++) {
	28	fRSec[i]=r[i];
	29	printf("\n R %d %f",i,fRSec[i]);
	30
	31	}
	32	}
	33
	34	void AliMUONsegmentationV01::SetPadDivision(Int_t ndiv[4])
	35	{
	36	//
	37	// Defines the pad size perp. to the anode wire (y) for different sectors.
	38	//
	39	for (Int_t i=0; i<4; i++) {
	40	fNDiv[i]=ndiv[i];
	41	printf("\n Ndiv %d %d",i,fNDiv[i]);
	42	}
	43	ndiv[0]=ndiv[1];
	44	}
	45
	46
	47	void AliMUONsegmentationV01::Init(AliMUONchamber*)
	48	{
	49	//
	50	// Fill the arrays fCx (x-contour) and fNpx (ix-contour) for each sector
	51	// These arrays help in converting from real to pad co-ordinates and
	52	// vice versa
	53	//
	54	Int_t isec;
	55	printf("\n Initialise segmentation v01");
	56	fNpy=Int_t(fRSec[fNsec-1]/fDpy)+1;
	57
	58	fDpxD[fNsec-1]=fDpx;
	59	if (fNsec > 1) {
	60	for (Int_t i=fNsec-2; i>=0; i--){
	61	fDpxD[i]=fDpxD[fNsec-1]/fNDiv[i];
	62	printf("\n dx %d %f",i,fDpxD[i]);
	63	}
	64	}
65	fWireD=fDpxD[1]/3;
66	//
67	// fill the arrays defining the pad segmentation boundaries
68	Float_t ry;
69	Int_t dnx;
70	//
71	// loop over sections
72	for(isec=0; isec<fNsec; isec++) {
73	//
74	// loop over pads along the aode wires
75	for (Int_t iy=1; iy<=fNpy; iy++) {
76	//
77	Float_t x=iy*fDpy-fDpy/2;
78	if (x > fRSec[isec]) {
79	fNpx[isec][iy]=0;
80	fCx[isec][iy]=0;
81	} else {
82	ry=TMath::Sqrt(fRSec[isec]fRSec[isec]-xx);
83	if (isec > 0) {
84	dnx= Int_t((ry-fCx[isec-1][iy])/fDpxD[isec]);
85	if (TMath::Odd(dnx)) dnx--;
86	fNpx[isec][iy]=fNpx[isec-1][iy]+dnx;
87	fCx[isec][iy]=fCx[isec-1][iy]+dnx*fDpxD[isec];
88	} else {
89	dnx=Int_t(ry/fDpxD[isec]);
90	fNpx[isec][iy]=dnx;
91	if (TMath::Odd(dnx)) dnx--;
92	fCx[isec][iy]=dnx*fDpxD[isec];
93	}
94	}
95	} // y-pad loop
96	} // sector loop
97
98	//
99	// for debugging only
100	for (Int_t iy=0; iy<fNpy; iy++) {
101	printf("\n iy %d",iy);
102	for(isec=0; isec<fNsec; isec++) {
103	printf(" %d",
104	fNpx[isec][iy]);
105	}
106	printf("\n iy %d",iy);
107	for(isec=0; isec<fNsec; isec++) {
108	printf(" %f",
109	fCx[isec][iy] );
110	}
111	printf("\n");
112	}
113	Float_t x,y;
114	Int_t jx,jy;
115	//
116	// for debugging only
117	for (Int_t ix=1; ix<100; ix++) {
118	for (Int_t iy=1; iy<100; iy++) {
119	GetPadCxy(ix,iy,x,y);
120	GetPadIxy(x,y,jx,jy);
121	if ((ix != jx) && (jx!=-1)) {
122	printf("\n %d %f %d %f", ix,x,iy,y);
123	printf("\n %d %f %d %f", jx,x,jy,y);
124	printf("\n \n **********");
125	}
126	}
127	}
128	}
129
130	Int_t AliMUONsegmentationV01::Sector(Int_t ix, Int_t iy)
131	{
132	Int_t absix=TMath::Abs(ix);
133	Int_t absiy=TMath::Abs(iy);
134	Int_t isec=0;
135	for (Int_t i=0; i<fNsec; i++) {
136	if (absix<=fNpx[i][absiy]){
137	isec=i;
138	break;
139	}
140	}
141	return isec;
142	}
143
144
145	Float_t AliMUONsegmentationV01::GetAnod(Float_t xhit)
146	{
147	//
148	// Finds anod-wire position closest to xhit
149	Float_t wire= (xhit<0)? Int_t(xhit/fWireD)+0.5:Int_t(xhit/fWireD)-0.5;
150	// printf("getanod: %f %f %f ",xhit, wire, fWireD*wire);
151	return fWireD*wire;
152	}
153
154	void AliMUONsegmentationV01::
155	GetPadIxy(Float_t x, Float_t y, Int_t &ix, Int_t &iy)
156	{
157	// returns pad coordinates (ix,iy) for given real coordinates (x,y)
158	//
159	iy = (y>0)? Int_t(y/fDpy)+1 : Int_t(y/fDpy)-1;
160	if (iy > fNpy) iy= fNpy;
161	if (iy < -fNpy) iy=-fNpy;
162	//
163	// Find sector isec
164	Int_t isec=-1;
165	Float_t absx=TMath::Abs(x);
166	Int_t absiy=TMath::Abs(iy);
167	for (Int_t i=0; isec < fNsec; i++) {
168	if (absx <= fCx[i][absiy]) {
169	isec=i;
170	break;
171	}
172	}
173	if (isec>0) {
174	ix= Int_t((absx-fCx[isec-1][absiy])/fDpxD[isec])
175	+fNpx[isec-1][absiy]+1;
176	} else if (isec == 0) {
177	ix= Int_t(absx/fDpxD[isec])+1;
178	} else {
179	ix=fNpx[fNsec-1][absiy]+1;
180	}
181	// printf("\n %d %d",isec,absiy);
182
183	ix = (x>0) ? ix:-ix;
184	}
185
186	void AliMUONsegmentationV01::
187	GetPadCxy(Int_t ix, Int_t iy, Float_t &x, Float_t &y)
188	{
189	// returns real coordinates (x,y) for given pad coordinates (ix,iy)
190	//
191	y = (iy>0) ? Float_t(iyfDpy)-fDpy/2. : Float_t(iyfDpy)+fDpy/2.;
192	//
193	// Find sector isec
194	Int_t isec=Sector(ix,iy);
195	//
196	Int_t absix=TMath::Abs(ix);
197	Int_t absiy=TMath::Abs(iy);
198	if (isec) {
199	x=fCx[isec-1][absiy]+(absix-fNpx[isec-1][absiy])*fDpxD[isec];
200	x=(ix>0) ? x-fDpxD[isec]/2 : -x+fDpxD[isec]/2;
201	} else {
202	x=y=0;
203	}
204	}
205
206
207	void AliMUONsegmentationV01::
208	GetSuperPadIxy(Float_t x, Float_t y, Int_t &ix, Int_t &iy)
209	{
210	ix = (x>0)? Int_t(x/fDpx)+1 : Int_t(x/fDpx)-1;
211	iy = (y>0)? Int_t(y/fDpy)+1 : Int_t(y/fDpy)-1;
212	}
213
214	void AliMUONsegmentationV01::
215	GetSuperPadCxy(Int_t ix, Int_t iy, Float_t &x, Float_t &y)
216	{
217	x = (ix>0) ? Float_t(ixfDpx)-fDpx/2. : Float_t(ixfDpx)+fDpx/2.;
218	y = (iy>0) ? Float_t(iyfDpy)-fDpy/2. : Float_t(iyfDpy)+fDpy/2.;
219	}
220
221	void AliMUONsegmentationV01::SetPADSIZ(Float_t p1, Float_t p2)
222	{
223	fDpx=p1;
224	fDpy=p2;
225	}
226
227	void AliMUONsegmentationV01::
228	FirstPad(Float_t xhit, Float_t yhit, Float_t dx, Float_t dy)
229	{
230	//
231	// Find the wire position (center of charge distribution)
232	Float_t x0a=GetAnod(xhit);
233	//
234	// and take fNsigma*sigma around this center
235	Float_t x01=x0a - dx;
236	Float_t x02=x0a + dx;
237	Float_t y01=yhit - dy;
238	Float_t y02=yhit + dy;
239	//
240	// find the pads over which the charge distributes
241	GetPadIxy(x01,y01,fixmin,fiymin);
242	GetPadIxy(x02,y02,fixmax,fiymax);
243	fxmin=x01;
244	fxmax=x02;
245	//
246	// upper and lower limits should be checked
247	/*
248	printf("\n FirstPad called");
249	printf("\n Hit Position %f %f",xhit,yhit);
250	printf("\n Closest wire %f", x0a);
251	printf("\n Integration limits: %i %i %i %i",fixmin,fixmax,fiymin,fiymax);
252	printf("\n Integration limits: %f %f %f %f \n",x01,x02,y01,y02);
253	*/
254	//
255	// Set current pad to lower left corner
256	if (fixmax < fixmin) fixmax=fixmin;
257	if (fiymax < fiymin) fiymax=fiymin;
258	fix=fixmin;
259	fiy=fiymin;
260	GetPadCxy(fix,fiy,fx,fy);
261
262	}
263
264	void AliMUONsegmentationV01::NextPad()
265	{
266	//
267	// Step to next pad in integration region
268	Float_t xc,yc;
269	Int_t iyc;
270
271	// step from left to right
272	if (fx < fxmax && fx != 0) {
273	fix++;
274	// step up
275	} else if (fiy != fiymax) {
276	fiy++;
277	// get y-position of next row (yc), xc not used here
278	GetPadCxy(fix,fiy,xc,yc);
279	// get x-pad coordiante for 1 pad in row (fix)
280	GetPadIxy(fxmin,yc,fix,iyc);
281	} else {
282	printf("\n Error: Stepping outside integration region\n ");
283	}
284	GetPadCxy(fix,fiy,fx,fy);
285	fSector=Sector(fix,fiy);
286	// printf("\n this pad %f %f %d %d",fx,fy,fix,fiy);
287
288	}
289
290	Int_t AliMUONsegmentationV01::MorePads()
291	//
292	// Are there more pads in the integration region
293	{
294	if ((fx >= fxmax && fiy >= fiymax) \|\| fy==0) {
295	return 0;
296	} else {
297	return 1;
298	}
299	}
300
301	void AliMUONsegmentationV01::SigGenInit(Float_t x,Float_t y,Float_t)
302	{
303	//
304	// Initialises pad and wire position during stepping
305	fxt=x;
306	fyt=y;
307	GetPadIxy(x,y,fixt,fiyt);
308	fiwt= (x>0) ? Int_t(x/fWireD)+1 : Int_t(x/fWireD)-1 ;
309	}
310
311	Int_t AliMUONsegmentationV01::SigGenCond(Float_t x,Float_t y,Float_t)
312	{
313	//
314	// Signal will be generated if particle crosses pad boundary or
315	// boundary between two wires.
316	Int_t ixt, iyt;
317	GetPadIxy(x,y,ixt,iyt);
318
319	Int_t iwt=(x>0) ? Int_t(x/fWireD)+1 : Int_t(x/fWireD)-1;
320
321	if ((ixt != fixt) \|\| (iyt !=fiyt) \|\| (iwt != fiwt)) {
322	return 1;
323	} else {
324	return 0;
325	}
326	}
327
328	void AliMUONsegmentationV01::
329	IntegrationLimits(Float_t& x1,Float_t& x2,Float_t& y1, Float_t& y2)
330	{
331	x1=fxt-fx-fDpxD[fSector]/2.;
332	x2=x1+fDpxD[fSector];
333	y1=fyt-fy-fDpy/2.;
334	y2=y1+fDpy;
335	}
336
337	void AliMUONsegmentationV01::
338	Neighbours(Int_t iX, Int_t iY, Int_t* Nlist, Int_t Xlist[10], Int_t Ylist[10])
339	{
340	const Float_t epsilon=fDpy/1000;
341
342	Float_t x,y;
343	Int_t ixx, iyy, isec1;
344	//
345	Int_t isec0=Sector(iX,iY);
346	Int_t i=0;
347	//
348	// step right
349	Xlist[i]=iX+1;
350	Ylist[i++]=iY;
351	//
352	// step left
353	Xlist[i]=iX-1;
354	Ylist[i++]=iY;
355	//
356	// step up
357	GetPadCxy(iX,iY,x,y);
358	GetPadIxy(x+epsilon,y+fDpy,ixx,iyy);
359	Xlist[i]=ixx;
360	Ylist[i++]=iY+1;
361	isec1=Sector(ixx,iyy);
362	if (isec1==isec0) {
363	//
364	// no sector boundary crossing
365	Xlist[i]=ixx+1;
366	Ylist[i++]=iY+1;
367
368	Xlist[i]=ixx-1;
369	Ylist[i++]=iY+1;
370	} else if (isec1 < isec0) {
371	// finer segmentation
372	Xlist[i]=ixx+1;
373	Ylist[i++]=iY+1;
374
375	Xlist[i]=ixx-1;
376	Ylist[i++]=iY+1;
377
378	Xlist[i]=ixx-2;
379	Ylist[i++]=iY+1;
380	} else {
381	// coarser segmenation
382	if (TMath::Odd(iX-fNpx[isec1-1][iY+1])) {
383	Xlist[i]=ixx-1;
384	Ylist[i++]=iY+1;
385	} else {
386	Xlist[i]=ixx+1;
387	Ylist[i++]=iY+1;
388	}
389	}
390	//
391	// step down
392	GetPadCxy(iX,iY,x,y);
393	GetPadIxy(x+epsilon,y-fDpy,ixx,iyy);
394	Xlist[i]=ixx;
395	Ylist[i++]=iY-1;
396	isec1=Sector(ixx,iyy);
397	if (isec1==isec0) {
398	//
399	// no sector boundary crossing
400	Xlist[i]=ixx+1;
401	Ylist[i++]=iY-1;
402
403	Xlist[i]=ixx-1;
404	Ylist[i++]=iY-1;
405	} else if (isec1 < isec0) {
406	// finer segmentation
407	Xlist[i]=ixx+1;
408	Ylist[i++]=iY-1;
409
410	Xlist[i]=ixx-1;
411	Ylist[i++]=iY-1;
412
413	Xlist[i]=ixx-2;
414	Ylist[i++]=iY-1;
415	} else {
416	// coarser segmenation
417	if (TMath::Odd(iX-fNpx[isec1-1][iY-1])) {
418	Xlist[i]=ixx-1;
419	Ylist[i++]=iY-1;
420	} else {
421	Xlist[i]=ixx+1;
422	Ylist[i++]=iY-1;
423	}
424	}
425	*Nlist=i;
426	}
427
428	//___________________________________________
429	void AliMUONsegmentationV01::
430	FitXY(AliMUONRecCluster* Cluster,TClonesArray* MUONdigits)
431	// Default : Centre of gravity method
432	{
433	Float_t x=0;
434	Float_t y=0;
435	Float_t q=0;
436	Float_t xToAdd;
437	Float_t yToAdd;
438
439	if (gAlice->TreeD()->GetReadEvent() != Cluster->GetCathod()+1)
440	// next line warns if in the future cathod 1 is not event 2 !
441	printf("ClusterFillXY : not reading the right cathod !\n");
442	for(Int_t clusDigit=Cluster->FirstDigitIndex();
443	clusDigit!=Cluster->InvalidDigitIndex();
444	clusDigit=Cluster->NextDigitIndex()) {
445	AliMUONdigit* pDigit=(AliMUONdigit*)MUONdigits->UncheckedAt(clusDigit);
446	GetPadCxy(pDigit->fPadX,pDigit->fPadY,xToAdd,yToAdd);
447	x+= xToAdd*pDigit->fSignal;
448	y+= yToAdd*pDigit->fSignal;
449	q+= (Float_t) pDigit->fSignal;
450	}
451	Cluster->fX=x/q;
452	Cluster->fY=y/q;
453	}
454
455
456
457
458
459
460
461