[u/mrichter/AliRoot.git] / MUON / AliMUONClusterFinderv0.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$
*/

#include "AliMUONClusterFinderv0.h"
#include "AliMUONSegResV1.h"
//#include "TTree.h"
//#include "AliRun.h"
//#include <TCanvas.h>
//#include <TH1.h>
//#include <TPad.h>
//#include <TGraph.h> 

//----------------------------------------------------------
ClassImp(AliMUONClusterFinderv0)

    AliMUONClusterFinderv0::AliMUONClusterFinderv0
(AliMUONsegmentation *segmentation, AliMUONresponse *response, 
 TClonesArray *digits, Int_t chamber) : AliMUONClusterFinder(segmentation,response,digits,chamber)
{;}

    AliMUONClusterFinderv0::AliMUONClusterFinderv0():AliMUONClusterFinder()
{;}

/*
void AliMUONClusterFinder::AddRawCluster(const AliMUONRawCluster c)
{
  //
  // Add a raw cluster copy to the list
  //
    AliMUON *MUON=(AliMUON*)gAlice->GetModule("MUON");
    MUON->AddRawCluster(fChamber,c); 
    fNRawClusters++;
}
*/


void AliMUONClusterFinderv0::Decluster(AliMUONRawCluster *cluster)
{
//    AliMUONdigit *dig;
//    Int_t q;
    static int done=0;
    if (!done) {
	printf("Calling decluster\n");
	done=1;
    }
    

    Int_t mul = cluster->fMultiplicity;
//    printf("Decluster - multiplicity   %d \n",mul);

    if (mul == 1) {
//	printf("\n Nothing special for 1-clusters \n");
//
// Nothing special for 1-clusters
//
	AddRawCluster(*cluster); 
    } else if (mul ==2) {
//
// 2-cluster, compute offset
//
        SetOffset(cluster);
	FillCluster(cluster);
        AddRawCluster(*cluster); 
    } else if (mul ==3) {
//
// 3-cluster, check topology
//	printf("\n 3-cluster, check topology \n");
//
	if (Centered(cluster)) {
//
// ok, cluster is centered 
//	    printf("\n ok, cluster is centered \n");
	} else {
//
// cluster is not centered, split into 2+1
//	    printf("\n cluster is not centered, split into 2+1 \n");
	}
	    
    } else {
	if (mul >(50-5)) printf("Decluster - multiplicity %d approaching 50\n",mul);
// 
// 4-and more-pad clusters
//
	SplitByLocalMaxima(cluster);
    } // multiplicity 
}

Int_t AliMUONClusterFinderv0::PeakOffsetAndCoordinates(Int_t DigitIndex, Float_t *X, Float_t *Y)
//
// Computes for which allowed offsets the digit has the highest neighbouring charge
// Returns the value of the offset, and sets the pyisical coordinates of that pad
// Loop on physical neighbours is specific to AliMUONsegmentationV1
{
Int_t nPara, offset, returnOffset=0 ;
AliMUONdigit* dig= (AliMUONdigit*)fDigits->UncheckedAt(DigitIndex);
AliMUONsegmentationV1* seg = (AliMUONsegmentationV1*) fSegmentation;
seg->GetNParallelAndOffset(dig->fPadX,dig->fPadY,&nPara,&offset);
if (nPara>1)
  {
  Float_t qMax=0;
  for (Int_t i=0;i<nPara; i++)
    {
    // Compute the charge on the 9 neighbouring pads
    // We assume that there are no pads connected in parallel in the neighbourhood
    //
    Float_t q=0;
    for (Int_t dx=-1;dx<2;dx++)
      for (Int_t dy=-1;dy<2;dy++)
        {
        if (dx==dy && dy==0)
          continue; 
        Int_t padY=dig->fPadY+dy;
        Int_t padX=seg->Ix((Int_t) (dig->fPadX+dx+i*offset) , padY);
	if (fHitMap->TestHit(padX, padY)==empty) 
	   continue;
        AliMUONdigit* digt = (AliMUONdigit*) fHitMap->GetHit(padX,padY);
        q += digt->fSignal;
        }
    if (q>qMax)
      {
      returnOffset=i*offset;
      qMax=q;
      }
    }
  }
fSegmentation->GetPadCxy(dig->fPadX+returnOffset,dig->fPadY,*X,*Y);
return returnOffset;
}


void AliMUONClusterFinderv0::SetOffset(AliMUONRawCluster *cluster)
// compute the offsets assuming that there is only one peak !
{
//DumpCluster(cluster);
Float_t X,Y;
cluster->fOffsetMap[0]=PeakOffsetAndCoordinates(cluster->fIndexMap[0],&X,&Y);
for (Int_t i=1;i<cluster->fMultiplicity;i++) {
  AliMUONdigit* dig= (AliMUONdigit*)fDigits->UncheckedAt(cluster->fIndexMap[i]);
  fSegmentation->Distance2AndOffset(dig->fPadX,dig->fPadY,X,Y,&(cluster->fOffsetMap[i]));
  }
}

void AliMUONClusterFinderv0::DumpCluster(AliMUONRawCluster *cluster)
{    
printf ("other cluster\n");
for (Int_t i=0; i<cluster->fMultiplicity; i++)
    {
	AliMUONdigit* dig= (AliMUONdigit*)fDigits->UncheckedAt(cluster->fIndexMap[i]);
	Int_t nPara, offset;
	fSegmentation->GetNParallelAndOffset(dig->fPadX,dig->fPadY,&nPara,&offset);

	printf("X %d Y %d Q %d NPara %d \n",dig->fPadX, dig->fPadY,dig->fSignal, nPara);
    }
}

Bool_t AliMUONClusterFinderv0::Centered(AliMUONRawCluster *cluster)
{
    AliMUONdigit* dig;
    dig= (AliMUONdigit*)fDigits->UncheckedAt(cluster->fIndexMap[0]);
    Int_t ix=dig->fPadX;
    Int_t iy=dig->fPadY;
    Int_t nn;
    Int_t X[kMaxNeighbours], Y[kMaxNeighbours], XN[kMaxNeighbours], YN[kMaxNeighbours];
    fSegmentation->Neighbours(ix,iy,&nn,X,Y);

    Int_t nd=0;
    for (Int_t i=0; i<nn; i++) {
	if (fHitMap->TestHit(X[i],Y[i]) == used) {
	    XN[nd]=X[i];
	    YN[nd]=Y[i];
	    nd++;
	}
    }
    if (nd==2) {
//
// cluster is centered !
	SetOffset(cluster);
	FillCluster(cluster);
	AddRawCluster(*cluster);
	return kTRUE;
    } else if (nd ==1) {
//
// Highest signal on an edge, split cluster into 2+1
//
// who is the neighbour ?
	Int_t nind=fHitMap->GetHitIndex(XN[0], YN[0]);
	Int_t i1= (nind==cluster->fIndexMap[1]) ? 1:2;
	Int_t i2= (nind==cluster->fIndexMap[1]) ? 2:1;    
//
// 2-cluster
	AliMUONRawCluster cnew;
	cnew.fMultiplicity=2;
	cnew.fIndexMap[0]=cluster->fIndexMap[0];
	cnew.fIndexMap[1]=cluster->fIndexMap[i1];
	SetOffset(&cnew);
	FillCluster(&cnew);
	AddRawCluster(cnew);
//
// 1-cluster
	cluster->fMultiplicity=1;
	cluster->fIndexMap[0]=cluster->fIndexMap[i2];
	cluster->fIndexMap[1]=0;
	cluster->fIndexMap[2]=0;	
	FillCluster(cluster);
	AddRawCluster(*cluster);
	return kFALSE;
    } else {
	printf("\n Completely screwed up %d !! \n",nd);
	
    }
    
	return kFALSE;
}


void AliMUONClusterFinderv0::SplitByLocalMaxima(AliMUONRawCluster *c)
{
    AliMUONdigit* dig[50], *digt;
    Int_t ix[50], iy[50], q[50];
    Float_t x[50], y[50];
    Int_t i; // loops over digits
    Int_t j; // loops over local maxima
    
    Int_t mul=c->fMultiplicity;
//
//  dump digit information into arrays
//
    for (i=0; i<mul; i++)
    {
	dig[i]= (AliMUONdigit*)fDigits->UncheckedAt(c->fIndexMap[i]);
	ix[i]= dig[i]->fPadX;
	iy[i]= dig[i]->fPadY;
	q[i] = dig[i]->fSignal;
	fSegmentation->GetPadCxy(ix[i], iy[i], x[i], y[i]);
    }
//
//  Find local maxima
//
    Bool_t IsLocal[50];
    Int_t NLocal=0;
    Int_t AssocPeak[50];
    Int_t IndLocal[50];
    Int_t nn;
    Int_t X[kMaxNeighbours], Y[kMaxNeighbours];
    for (i=0; i<mul; i++) {
	fSegmentation->Neighbours(ix[i], iy[i], &nn, X, Y);
	IsLocal[i]=kTRUE;
	for (j=0; j<nn; j++) {
	    if (fHitMap->TestHit(X[j], Y[j])==empty) continue;
	    digt=(AliMUONdigit*) fHitMap->GetHit(X[j], Y[j]);
	    if (digt->fSignal > q[i]) {
		IsLocal[i]=kFALSE;
		break;
//
// handle special case of neighbouring pads with equal signal
	    } else if (digt->fSignal == q[i]) {
		if (NLocal >0) {
		    for (Int_t k=0; k<NLocal; k++) {
			if (X[j]==ix[IndLocal[k]] && Y[j]==iy[IndLocal[k]]){
			    IsLocal[i]=kFALSE;
			}
		    }
		}
	    } 
	} // loop over next neighbours
	if (IsLocal[i]) {
	    IndLocal[NLocal]=i;
            // New for LYON : we guess which is the actual position of the pad hit
	    // But this would run like that for normal chamber !
            c->fOffsetMap[i]=PeakOffsetAndCoordinates(c->fIndexMap[i], &(x[i]), &(y[i]));
	    NLocal++;
	} 
    } // loop over all digits
//    printf("Found %d local Maxima",NLocal);
// 
// Check if enough local clusters have been found,
// if not add global maxima to the list 
//
// But what the hell is that ? (Manu)
//
//    Int_t nPerMax=mul/NLocal;
//    if (nPerMax > 5) {
//	Int_t nGlob=mul/5-NLocal+1;
//	if (nGlob > 0) {
//	    Int_t nnew=0;
//	    for (i=0; i<mul; i++) {
//		if (!IsLocal[i]) {
//		    IndLocal[NLocal]=i;
//		    IsLocal[i]=kTRUE;
//		    NLocal++;
//		    nnew++;
//		}
//		if (nnew==nGlob) break;
//	    }
//	}
//   }
//
//
// Associate hits to peaks
//
    for (i=0; i<mul; i++) {
        //
        // loop on digits
        //
	Float_t dmin=1.E10;
	Float_t qmax=0;
	Int_t offset;
	if (IsLocal[i]) continue;
	for (j=0; j<NLocal; j++) {
            //
            // Loop on peaks
            //
	    Int_t il=IndLocal[j];
//	    Float_t d=TMath::Sqrt((x[i]-x[il])*(x[i]-x[il])
//				  +(y[i]-y[il])*(y[i]-y[il]));
            // Can run like that for non-Lyon chambers
	    Float_t d = fSegmentation->Distance2AndOffset(ix[i],iy[i],x[il],y[il], &offset);
	    Float_t ql=q[il];
//
// Select nearest peak
//
	    if (d<dmin) {
		dmin=d;
		qmax=ql;
		AssocPeak[i]=j;
		c->fOffsetMap[i]=offset;
	    } else if (d==dmin) {
//
// If more than one take highest peak
//
		if (ql>qmax) {
		    dmin=d;
		    qmax=ql;
		    AssocPeak[i]=j;
		    c->fOffsetMap[i]=offset;
		}
	    } // end if
	} // End loop on peaks
    } // end loop on digits
//
// One cluster for each maximum
//
    for (j=0; j<NLocal; j++) {
	AliMUONRawCluster cnew;
	cnew.fIndexMap[0]=c->fIndexMap[IndLocal[j]];
	cnew.fOffsetMap[0]=c->fOffsetMap[IndLocal[j]];
	cnew.fMultiplicity=1;
	for (i=0; i<mul; i++) {
	    if (IsLocal[i]) continue;
	    if (AssocPeak[i]==j) {
		cnew.fIndexMap[cnew.fMultiplicity]=c->fIndexMap[i];
		cnew.fOffsetMap[cnew.fMultiplicity]=c->fOffsetMap[i];
		cnew.fMultiplicity++;
	    }
	}
	FillCluster(&cnew);
	AddRawCluster(cnew);
    }
}

/*
void  AliMUONClusterFinderv0::FillCluster(AliMUONRawCluster* c)
{
//
//  Completes cluster information starting from list of digits
//
    AliMUONdigit* dig;
    Float_t x, y;
    Int_t  ix, iy;
    
    c->fPeakSignal=0;
    c->fX=0;
    c->fY=0;
    c->fQ=0;
    for (Int_t i=0; i<c->fMultiplicity; i++)
    {
	dig= (AliMUONdigit*)fDigits->UncheckedAt(c->fIndexMap[i]);
	ix=dig->fPadX + c.fOffsetMap[i]; // should be 0 for non-LYON
	iy=dig->fPadY;
	Int_t q=dig->fSignal;
//
//
// peak signal and track list
	if (q>c->fPeakSignal) {
	    c->fPeakSignal=0;
	    c->fTracks[0]=dig->fTracks[0];
	    c->fTracks[1]=dig->fTracks[1];
	    c->fTracks[2]=dig->fTracks[2];
	}
//
// centre of gravity
	fSegmentation->GetPadCxy(ix, iy, x, y);
	c->fX += q*x;
	c->fY += q*y;
	c->fQ += q;
    }
    c->fX/=c->fQ;
// Not valid for inclined tracks in X !!! (Manu)
//    c->fX=fSegmentation->GetAnod(c->fX);
    c->fY/=c->fQ; 
//
//  apply correction to the coordinate along the anode wire
//
    if (fCogCorr) {
	x=c->fX;   
	y=c->fY;
	fSegmentation->GetPadIxy(x, y, ix, iy);
	fSegmentation->GetPadCxy(ix, iy, x, y);
	Float_t YonPad=(c->fY-y)/fSegmentation->Dpy();
	c->fY=y-fCogCorr->Eval(YonPad, 0, 0);
    }

}
*/
/*
void  AliMUONClusterFinder::FindCluster(Int_t i, Int_t j, AliMUONRawCluster &c){
//
//  Find clusters
//
//
//  Add i,j as element of the cluster
//    
    
    Int_t idx = fHitMap->GetHitIndex(i,j);
    AliMUONdigit* dig = (AliMUONdigit*) fHitMap->GetHit(i,j);
    Int_t q=dig->fSignal;
    if (q > TMath::Abs(c.fPeakSignal)) {
	c.fPeakSignal=q;
	c.fTracks[0]=dig->fTracks[0];
	c.fTracks[1]=dig->fTracks[1];
	c.fTracks[2]=dig->fTracks[2];
    }
//
//  Make sure that list of digits is ordered 
// 
    Int_t mu=c.fMultiplicity;
    c.fIndexMap[mu]=idx;

    if (mu > 0) {
	for (Int_t ind=mu-1; ind>=0; ind--) {
	    Int_t ist=(c.fIndexMap)[ind];
	    Int_t ql=((AliMUONdigit*)fDigits
		      ->UncheckedAt(ist))->fSignal;
	    if (q>ql) {
		c.fIndexMap[ind]=idx;
		c.fIndexMap[ind+1]=ist;
	    } else {
		break;
	    }
	}
    }
    
    c.fMultiplicity++;
    
    if (c.fMultiplicity >= 50 ) {
	printf("FindCluster - multiplicity >50  %d \n",c.fMultiplicity);
	c.fMultiplicity=50-1;
    }

// Prepare center of gravity calculation
    Float_t x, y;
    fSegmentation->GetPadCxy(i, j, x, y);
    c.fX += q*x;
    c.fY += q*y;
    c.fQ += q;
// Flag hit as taken  
    fHitMap->FlagHit(i,j);
//
//  Now look recursively for all neighbours
//  
    Int_t nn;
    Int_t Xlist[kMaxNeighbours], Ylist[kMaxNeighbours];
    fSegmentation->Neighbours(i,j,&nn,Xlist,Ylist);
    for (Int_t in=0; in<nn; in++) {
	Int_t ix=Xlist[in];
	Int_t iy=Ylist[in];
	if (fHitMap->TestHit(ix,iy)==unused) FindCluster(ix, iy, c);
    }
}
*/

//_____________________________________________________________________________

void AliMUONClusterFinderv0::FindRawClusters()
{
  //
  // simple MUON cluster finder from digits -- finds neighbours and 
  // fill the tree with raw clusters
  //
    if (!fNdigits) return;

    fHitMap = new AliMUONHitMapA1(fSegmentation, fDigits);

    AliMUONdigit *dig;

    int ndig;
    int nskip=0;

    fHitMap->FillHits();
    for (ndig=0; ndig<fNdigits; ndig++) {
	dig = (AliMUONdigit*)fDigits->UncheckedAt(ndig);
	Int_t i=dig->fPadX;
	Int_t j=dig->fPadY;
	if (fHitMap->TestHit(i,j)==used ||fHitMap->TestHit(i,j)==empty) {
	    nskip++;
	    continue;
	}
	AliMUONRawCluster c;
	c.fMultiplicity=0;
//	c.fPeakSignal=dig->fSignal;
//	c.fTracks[0]=dig->fTracks[0];
//	c.fTracks[1]=dig->fTracks[1];
//	c.fTracks[2]=dig->fTracks[2];
	c.fPeakSignal=0;
	FindCluster(i,j, c);
	// center of gravity
	c.fX /= c.fQ;
	c.fX=fSegmentation->GetAnod(c.fX);
	c.fY /= c.fQ;
//
//  apply correction to the coordinate along the anode wire
//

	
    if (fCogCorr) {
	Int_t ix,iy;
	Float_t x=c.fX;   
	Float_t y=c.fY;
	fSegmentation->GetPadIxy(x, y, ix, iy);
	fSegmentation->GetPadCxy(ix, iy, x, y);
	Float_t YonPad=(c.fY-y)/fSegmentation->Dpy();
	c.fY=y-fCogCorr->Eval(YonPad,0,0);
    }
//
//      Analyse cluster and decluster if necessary
//	
	Decluster(&c);
//
//
//
//      reset Cluster object
	for (int k=0;k<c.fMultiplicity;k++) {
	    c.fIndexMap[k]=0;
 	    c.fOffsetMap[k]=0;
	}
	c.fMultiplicity=0;
    } // end loop ndig    
    delete fHitMap;
}

/*

void AliMUONClusterFinder::
CalibrateCOG()
{
    Float_t x[5];
    Float_t y[5];
    Int_t n, i;
    TF1 func;
    
    if (fSegmentation) {
	fSegmentation->GiveTestPoints(n, x, y);
	for (i=0; i<n; i++) {
	    Float_t xtest=x[i];
	    Float_t ytest=y[i];	    
	    SinoidalFit(xtest, ytest, func);
	}
	fCogCorr = new TF1(func);
    }
}
*/
/*

void AliMUONClusterFinder::
SinoidalFit(Float_t x, Float_t y, TF1 &func)
{
//
    static Int_t count=0;
    char canvasname[3];
    count++;
    sprintf(canvasname,"c%d",count);

// MANU : without const, error on HP
    const Int_t ns=101;
    Float_t xg[ns], yg[ns], xrg[ns], yrg[ns];
    Float_t xsig[ns], ysig[ns];
   
    AliMUONsegmentation *segmentation=fSegmentation;

    Int_t ix,iy;
    segmentation->GetPadIxy(x,y,ix,iy);   
    segmentation->GetPadCxy(ix,iy,x,y);   
    Int_t isec=segmentation->Sector(ix,iy);
// Pad Limits    
    Float_t xmin = x-segmentation->GetRealDpx(isec)/2;
    Float_t ymin = y-segmentation->Dpy()/2;
//      	
//      Integration Limits
    Float_t dxI=fResponse->Nsigma()*fResponse->ChwX();
    Float_t dyI=fResponse->Nsigma()*fResponse->ChwY();

//
//  Scanning
//
    Int_t i;
    Float_t qp;
//
//  y-position
    Float_t yscan=ymin;
    Float_t dy=segmentation->Dpy()/(ns-1);

    for (i=0; i<ns; i++) {
//
//      Pad Loop
//      
	Float_t sum=0;
	Float_t qcheck=0;
	segmentation->SigGenInit(x, yscan, 0);
	
	for (segmentation->FirstPad(x, yscan, dxI, dyI); 
	     segmentation->MorePads(); 
	     segmentation->NextPad()) 
	{
	    qp=fResponse->IntXY(segmentation);
	    qp=TMath::Abs(qp);
//
//
	    if (qp > 1.e-4) {
		qcheck+=qp;
		Int_t ixs=segmentation->Ix();
		Int_t iys=segmentation->Iy();
		Float_t xs,ys;
		segmentation->GetPadCxy(ixs,iys,xs,ys);
		sum+=qp*ys;
	    }
	} // Pad loop
	Float_t ycog=sum/qcheck;
	yg[i]=(yscan-y)/segmentation->Dpy();
	yrg[i]=(ycog-y)/segmentation->Dpy();
	ysig[i]=ycog-yscan;
	yscan+=dy;
    } // scan loop
//
//  x-position
    Float_t xscan=xmin;
    Float_t dx=segmentation->GetRealDpx(isec)/(ns-1);

    for (i=0; i<ns; i++) {
//
//      Pad Loop
//      
	Float_t sum=0;
	Float_t qcheck=0;
	segmentation->SigGenInit(xscan, y, 0);
	
	for (segmentation->FirstPad(xscan, y, dxI, dyI); 
	     segmentation->MorePads(); 
	     segmentation->NextPad()) 
	{
	    qp=fResponse->IntXY(segmentation);
	    qp=TMath::Abs(qp);
//
//
	    if (qp > 1.e-2) {
		qcheck+=qp;
		Int_t ixs=segmentation->Ix();
		Int_t iys=segmentation->Iy();
		Float_t xs,ys;
		segmentation->GetPadCxy(ixs,iys,xs,ys);
		sum+=qp*xs;
	    }
	} // Pad loop
	Float_t xcog=sum/qcheck;
	xcog=segmentation->GetAnod(xcog);
	
	xg[i]=(xscan-x)/segmentation->GetRealDpx(isec);
	xrg[i]=(xcog-x)/segmentation->GetRealDpx(isec);
	xsig[i]=xcog-xscan;
	xscan+=dx;
    }

   TCanvas *c1=new TCanvas(canvasname,canvasname,400,10,600,700);
   TPad* pad11 = new TPad("pad11"," ",0.01,0.51,0.49,0.99);
   TPad* pad12 = new TPad("pad12"," ",0.51,0.51,0.99,0.99);
   TPad* pad13 = new TPad("pad13"," ",0.01,0.01,0.49,0.49);
   TPad* pad14 = new TPad("pad14"," ",0.51,0.01,0.99,0.49);
   pad11->SetFillColor(11);
   pad12->SetFillColor(11);
   pad13->SetFillColor(11);
   pad14->SetFillColor(11);
   pad11->Draw();
   pad12->Draw();
   pad13->Draw();
   pad14->Draw();
   TGraph *graphx = new TGraph(ns,xg ,xsig);
   TGraph *graphxr= new TGraph(ns,xrg,xsig);   
   TGraph *graphy = new TGraph(ns,yg ,ysig);
   TGraph *graphyr= new TGraph(ns,yrg,ysig);
//
// Creates a Root function based on function sinoid above
// and perform the fit
//
   Double_t sinoid(Double_t *x, Double_t *par);
   TF1 *sinoidf = new TF1("sinoidf",sinoid,0.5,0.5,5);
   graphyr->Fit("sinoidf","V");
   sinoidf->Copy(func);
   func.Eval(0,0,0);
//
   pad11->cd();
   graphx->SetFillColor(42);
   graphx->SetMarkerColor(4);
   graphx->SetMarkerStyle(21);
   graphx->Draw("AC");
   graphx->GetHistogram()->SetXTitle("x on pad");
   graphx->GetHistogram()->SetYTitle("xcog-x");   


   pad12->cd();
   graphxr->SetFillColor(42);
   graphxr->SetMarkerColor(4);
   graphxr->SetMarkerStyle(21);
   graphxr->Draw("AP");
   graphxr->GetHistogram()->SetXTitle("xcog on pad");
   graphxr->GetHistogram()->SetYTitle("xcog-x");   


   pad13->cd();
   graphy->SetFillColor(42);
   graphy->SetMarkerColor(4);
   graphy->SetMarkerStyle(21);
   graphy->Draw("AF");
   graphy->GetHistogram()->SetXTitle("y on pad");
   graphy->GetHistogram()->SetYTitle("ycog-y");   


   pad14->cd();
   graphyr->SetFillColor(42);
   graphyr->SetMarkerColor(4);
   graphyr->SetMarkerStyle(21);
   graphyr->Draw("AF");
   graphyr->GetHistogram()->SetXTitle("ycog on pad");
   graphyr->GetHistogram()->SetYTitle("ycog-y");   

   c1->Update();
   
}
*/
/*
Double_t sinoid(Double_t *x, Double_t *par)
{
    Double_t arg = -2*TMath::Pi()*x[0];
    Double_t fitval= par[0]*TMath::Sin(arg)+
	par[1]*TMath::Sin(2*arg)+
	par[2]*TMath::Sin(3*arg)+
	par[3]*TMath::Sin(4*arg)+
	par[4]*TMath::Sin(5*arg);
    return fitval;
 }
*/
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
a897a37a	20	#include "AliMUONClusterFinderv0.h"
	21	#include "AliMUONSegResV1.h"
	22	//#include "TTree.h"
	23	//#include "AliRun.h"
	24	//#include <TCanvas.h>
	25	//#include <TH1.h>
	26	//#include <TPad.h>
	27	//#include <TGraph.h>
	28
	29	//----------------------------------------------------------
	30	ClassImp(AliMUONClusterFinderv0)
	31
	32	AliMUONClusterFinderv0::AliMUONClusterFinderv0
	33	(AliMUONsegmentation segmentation, AliMUONresponse response,
	34	TClonesArray *digits, Int_t chamber) : AliMUONClusterFinder(segmentation,response,digits,chamber)
	35	{;}
	36
	37	AliMUONClusterFinderv0::AliMUONClusterFinderv0():AliMUONClusterFinder()
	38	{;}
	39
	40	/*
	41	void AliMUONClusterFinder::AddRawCluster(const AliMUONRawCluster c)
	42	{
	43	//
	44	// Add a raw cluster copy to the list
	45	//
	46	AliMUON MUON=(AliMUON)gAlice->GetModule("MUON");
	47	MUON->AddRawCluster(fChamber,c);
	48	fNRawClusters++;
	49	}
	50	*/
	51
	52
	53
	54	void AliMUONClusterFinderv0::Decluster(AliMUONRawCluster *cluster)
	55	{
	56	// AliMUONdigit *dig;
	57	// Int_t q;
	58	static int done=0;
	59	if (!done) {
	60	printf("Calling decluster\n");
	61	done=1;
	62	}
	63
	64
	65
	66	Int_t mul = cluster->fMultiplicity;
	67	// printf("Decluster - multiplicity %d \n",mul);
	68
	69	if (mul == 1) {
	70	// printf("\n Nothing special for 1-clusters \n");
	71	//
	72	// Nothing special for 1-clusters
	73	//
	74	AddRawCluster(*cluster);
	75	} else if (mul ==2) {
	76	//
	77	// 2-cluster, compute offset
	78	//
	79	SetOffset(cluster);
	80	FillCluster(cluster);
	81	AddRawCluster(*cluster);
	82	} else if (mul ==3) {
	83	//
84	// 3-cluster, check topology
85	// printf("\n 3-cluster, check topology \n");
86	//
87	if (Centered(cluster)) {
88	//
89	// ok, cluster is centered
90	// printf("\n ok, cluster is centered \n");
91	} else {
92	//
93	// cluster is not centered, split into 2+1
94	// printf("\n cluster is not centered, split into 2+1 \n");
95	}
96
97	} else {
98	if (mul >(50-5)) printf("Decluster - multiplicity %d approaching 50\n",mul);
99	//
100	// 4-and more-pad clusters
101	//
102	SplitByLocalMaxima(cluster);
103	} // multiplicity
104	}
105
106	Int_t AliMUONClusterFinderv0::PeakOffsetAndCoordinates(Int_t DigitIndex, Float_t X, Float_t Y)
107	//
108	// Computes for which allowed offsets the digit has the highest neighbouring charge
109	// Returns the value of the offset, and sets the pyisical coordinates of that pad
110	// Loop on physical neighbours is specific to AliMUONsegmentationV1
111	{
112	Int_t nPara, offset, returnOffset=0 ;
113	AliMUONdigit* dig= (AliMUONdigit*)fDigits->UncheckedAt(DigitIndex);
114	AliMUONsegmentationV1* seg = (AliMUONsegmentationV1*) fSegmentation;
115	seg->GetNParallelAndOffset(dig->fPadX,dig->fPadY,&nPara,&offset);
116	if (nPara>1)
117	{
118	Float_t qMax=0;
119	for (Int_t i=0;i<nPara; i++)
120	{
121	// Compute the charge on the 9 neighbouring pads
122	// We assume that there are no pads connected in parallel in the neighbourhood
123	//
124	Float_t q=0;
125	for (Int_t dx=-1;dx<2;dx++)
126	for (Int_t dy=-1;dy<2;dy++)
127	{
128	if (dx==dy && dy==0)
129	continue;
130	Int_t padY=dig->fPadY+dy;
131	Int_t padX=seg->Ix((Int_t) (dig->fPadX+dx+i*offset) , padY);
132	if (fHitMap->TestHit(padX, padY)==empty)
133	continue;
134	AliMUONdigit* digt = (AliMUONdigit*) fHitMap->GetHit(padX,padY);
135	q += digt->fSignal;
136	}
137	if (q>qMax)
138	{
139	returnOffset=i*offset;
140	qMax=q;
141	}
142	}
143	}
144	fSegmentation->GetPadCxy(dig->fPadX+returnOffset,dig->fPadY,X,Y);
145	return returnOffset;
146	}
147
148
149	void AliMUONClusterFinderv0::SetOffset(AliMUONRawCluster *cluster)
150	// compute the offsets assuming that there is only one peak !
151	{
152	//DumpCluster(cluster);
153	Float_t X,Y;
154	cluster->fOffsetMap[0]=PeakOffsetAndCoordinates(cluster->fIndexMap[0],&X,&Y);
155	for (Int_t i=1;i<cluster->fMultiplicity;i++) {
156	AliMUONdigit* dig= (AliMUONdigit*)fDigits->UncheckedAt(cluster->fIndexMap[i]);
157	fSegmentation->Distance2AndOffset(dig->fPadX,dig->fPadY,X,Y,&(cluster->fOffsetMap[i]));
158	}
159	}
160
161	void AliMUONClusterFinderv0::DumpCluster(AliMUONRawCluster *cluster)
162	{
163	printf ("other cluster\n");
164	for (Int_t i=0; i<cluster->fMultiplicity; i++)
165	{
166	AliMUONdigit* dig= (AliMUONdigit*)fDigits->UncheckedAt(cluster->fIndexMap[i]);
167	Int_t nPara, offset;
168	fSegmentation->GetNParallelAndOffset(dig->fPadX,dig->fPadY,&nPara,&offset);
169
170	printf("X %d Y %d Q %d NPara %d \n",dig->fPadX, dig->fPadY,dig->fSignal, nPara);
171	}
172	}
173
174	Bool_t AliMUONClusterFinderv0::Centered(AliMUONRawCluster *cluster)
175	{
176	AliMUONdigit* dig;
177	dig= (AliMUONdigit*)fDigits->UncheckedAt(cluster->fIndexMap[0]);
178	Int_t ix=dig->fPadX;
179	Int_t iy=dig->fPadY;
180	Int_t nn;
181	Int_t X[kMaxNeighbours], Y[kMaxNeighbours], XN[kMaxNeighbours], YN[kMaxNeighbours];
182	fSegmentation->Neighbours(ix,iy,&nn,X,Y);
183
184	Int_t nd=0;
185	for (Int_t i=0; i<nn; i++) {
186	if (fHitMap->TestHit(X[i],Y[i]) == used) {
187	XN[nd]=X[i];
188	YN[nd]=Y[i];
189	nd++;
190	}
191	}
192	if (nd==2) {
193	//
194	// cluster is centered !
195	SetOffset(cluster);
196	FillCluster(cluster);
197	AddRawCluster(*cluster);
198	return kTRUE;
199	} else if (nd ==1) {
200	//
201	// Highest signal on an edge, split cluster into 2+1
202	//
203	// who is the neighbour ?
204	Int_t nind=fHitMap->GetHitIndex(XN[0], YN[0]);
205	Int_t i1= (nind==cluster->fIndexMap[1]) ? 1:2;
206	Int_t i2= (nind==cluster->fIndexMap[1]) ? 2:1;
207	//
208	// 2-cluster
209	AliMUONRawCluster cnew;
210	cnew.fMultiplicity=2;
211	cnew.fIndexMap[0]=cluster->fIndexMap[0];
212	cnew.fIndexMap[1]=cluster->fIndexMap[i1];
213	SetOffset(&cnew);
214	FillCluster(&cnew);
215	AddRawCluster(cnew);
216	//
217	// 1-cluster
218	cluster->fMultiplicity=1;
219	cluster->fIndexMap[0]=cluster->fIndexMap[i2];
220	cluster->fIndexMap[1]=0;
221	cluster->fIndexMap[2]=0;
222	FillCluster(cluster);
223	AddRawCluster(*cluster);
224	return kFALSE;
225	} else {
226	printf("\n Completely screwed up %d !! \n",nd);
227
228	}
229
230	return kFALSE;
231	}
232
233
234	void AliMUONClusterFinderv0::SplitByLocalMaxima(AliMUONRawCluster *c)
235	{
236	AliMUONdigit* dig[50], *digt;
237	Int_t ix[50], iy[50], q[50];
238	Float_t x[50], y[50];
239	Int_t i; // loops over digits
240	Int_t j; // loops over local maxima
241
242	Int_t mul=c->fMultiplicity;
243	//
244	// dump digit information into arrays
245	//
246	for (i=0; i<mul; i++)
247	{
248	dig[i]= (AliMUONdigit*)fDigits->UncheckedAt(c->fIndexMap[i]);
249	ix[i]= dig[i]->fPadX;
250	iy[i]= dig[i]->fPadY;
251	q[i] = dig[i]->fSignal;
252	fSegmentation->GetPadCxy(ix[i], iy[i], x[i], y[i]);
253	}
254	//
255	// Find local maxima
256	//
257	Bool_t IsLocal[50];
258	Int_t NLocal=0;
259	Int_t AssocPeak[50];
260	Int_t IndLocal[50];
261	Int_t nn;
262	Int_t X[kMaxNeighbours], Y[kMaxNeighbours];
263	for (i=0; i<mul; i++) {
264	fSegmentation->Neighbours(ix[i], iy[i], &nn, X, Y);
265	IsLocal[i]=kTRUE;
266	for (j=0; j<nn; j++) {
267	if (fHitMap->TestHit(X[j], Y[j])==empty) continue;
268	digt=(AliMUONdigit*) fHitMap->GetHit(X[j], Y[j]);
269	if (digt->fSignal > q[i]) {
270	IsLocal[i]=kFALSE;
271	break;
272	//
273	// handle special case of neighbouring pads with equal signal
274	} else if (digt->fSignal == q[i]) {
275	if (NLocal >0) {
276	for (Int_t k=0; k<NLocal; k++) {
277	if (X[j]==ix[IndLocal[k]] && Y[j]==iy[IndLocal[k]]){
278	IsLocal[i]=kFALSE;
279	}
280	}
281	}
282	}
283	} // loop over next neighbours
284	if (IsLocal[i]) {
285	IndLocal[NLocal]=i;
286	// New for LYON : we guess which is the actual position of the pad hit
287	// But this would run like that for normal chamber !
288	c->fOffsetMap[i]=PeakOffsetAndCoordinates(c->fIndexMap[i], &(x[i]), &(y[i]));
289	NLocal++;
290	}
291	} // loop over all digits
292	// printf("Found %d local Maxima",NLocal);
293	//
294	// Check if enough local clusters have been found,
295	// if not add global maxima to the list
296	//
297	// But what the hell is that ? (Manu)
298	//
299	// Int_t nPerMax=mul/NLocal;
300	// if (nPerMax > 5) {
301	// Int_t nGlob=mul/5-NLocal+1;
302	// if (nGlob > 0) {
303	// Int_t nnew=0;
304	// for (i=0; i<mul; i++) {
305	// if (!IsLocal[i]) {
306	// IndLocal[NLocal]=i;
307	// IsLocal[i]=kTRUE;
308	// NLocal++;
309	// nnew++;
310	// }
311	// if (nnew==nGlob) break;
312	// }
313	// }
314	// }
315	//
316	//
317	// Associate hits to peaks
318	//
319	for (i=0; i<mul; i++) {
320	//
321	// loop on digits
322	//
323	Float_t dmin=1.E10;
324	Float_t qmax=0;
325	Int_t offset;
326	if (IsLocal[i]) continue;
327	for (j=0; j<NLocal; j++) {
328	//
329	// Loop on peaks
330	//
331	Int_t il=IndLocal[j];
332	// Float_t d=TMath::Sqrt((x[i]-x[il])*(x[i]-x[il])
333	// +(y[i]-y[il])*(y[i]-y[il]));
334	// Can run like that for non-Lyon chambers
335	Float_t d = fSegmentation->Distance2AndOffset(ix[i],iy[i],x[il],y[il], &offset);
336	Float_t ql=q[il];
337	//
338	// Select nearest peak
339	//
340	if (d<dmin) {
341	dmin=d;
342	qmax=ql;
343	AssocPeak[i]=j;
344	c->fOffsetMap[i]=offset;
345	} else if (d==dmin) {
346	//
347	// If more than one take highest peak
348	//
349	if (ql>qmax) {
350	dmin=d;
351	qmax=ql;
352	AssocPeak[i]=j;
353	c->fOffsetMap[i]=offset;
354	}
355	} // end if
356	} // End loop on peaks
357	} // end loop on digits
358	//
359	// One cluster for each maximum
360	//
361	for (j=0; j<NLocal; j++) {
362	AliMUONRawCluster cnew;
363	cnew.fIndexMap[0]=c->fIndexMap[IndLocal[j]];
364	cnew.fOffsetMap[0]=c->fOffsetMap[IndLocal[j]];
365	cnew.fMultiplicity=1;
366	for (i=0; i<mul; i++) {
367	if (IsLocal[i]) continue;
368	if (AssocPeak[i]==j) {
369	cnew.fIndexMap[cnew.fMultiplicity]=c->fIndexMap[i];
370	cnew.fOffsetMap[cnew.fMultiplicity]=c->fOffsetMap[i];
371	cnew.fMultiplicity++;
372	}
373	}
374	FillCluster(&cnew);
375	AddRawCluster(cnew);
376	}
377	}
378
379	/*
380	void AliMUONClusterFinderv0::FillCluster(AliMUONRawCluster* c)
381	{
382	//
383	// Completes cluster information starting from list of digits
384	//
385	AliMUONdigit* dig;
386	Float_t x, y;
387	Int_t ix, iy;
388
389	c->fPeakSignal=0;
390	c->fX=0;
391	c->fY=0;
392	c->fQ=0;
393	for (Int_t i=0; i<c->fMultiplicity; i++)
394	{
395	dig= (AliMUONdigit*)fDigits->UncheckedAt(c->fIndexMap[i]);
396	ix=dig->fPadX + c.fOffsetMap[i]; // should be 0 for non-LYON
397	iy=dig->fPadY;
398	Int_t q=dig->fSignal;
399	//
400	//
401	// peak signal and track list
402	if (q>c->fPeakSignal) {
403	c->fPeakSignal=0;
404	c->fTracks[0]=dig->fTracks[0];
405	c->fTracks[1]=dig->fTracks[1];
406	c->fTracks[2]=dig->fTracks[2];
407	}
408	//
409	// centre of gravity
410	fSegmentation->GetPadCxy(ix, iy, x, y);
411	c->fX += q*x;
412	c->fY += q*y;
413	c->fQ += q;
414	}
415	c->fX/=c->fQ;
416	// Not valid for inclined tracks in X !!! (Manu)
417	// c->fX=fSegmentation->GetAnod(c->fX);
418	c->fY/=c->fQ;
419	//
420	// apply correction to the coordinate along the anode wire
421	//
422	if (fCogCorr) {
423	x=c->fX;
424	y=c->fY;
425	fSegmentation->GetPadIxy(x, y, ix, iy);
426	fSegmentation->GetPadCxy(ix, iy, x, y);
427	Float_t YonPad=(c->fY-y)/fSegmentation->Dpy();
428	c->fY=y-fCogCorr->Eval(YonPad, 0, 0);
429	}
430
431	}
432	*/
433	/*
434	void AliMUONClusterFinder::FindCluster(Int_t i, Int_t j, AliMUONRawCluster &c){
435	//
436	// Find clusters
437	//
438	//
439	// Add i,j as element of the cluster
440	//
441
442	Int_t idx = fHitMap->GetHitIndex(i,j);
443	AliMUONdigit* dig = (AliMUONdigit*) fHitMap->GetHit(i,j);
444	Int_t q=dig->fSignal;
445	if (q > TMath::Abs(c.fPeakSignal)) {
446	c.fPeakSignal=q;
447	c.fTracks[0]=dig->fTracks[0];
448	c.fTracks[1]=dig->fTracks[1];
449	c.fTracks[2]=dig->fTracks[2];
450	}
451	//
452	// Make sure that list of digits is ordered
453	//
454	Int_t mu=c.fMultiplicity;
455	c.fIndexMap[mu]=idx;
456
457	if (mu > 0) {
458	for (Int_t ind=mu-1; ind>=0; ind--) {
459	Int_t ist=(c.fIndexMap)[ind];
460	Int_t ql=((AliMUONdigit*)fDigits
461	->UncheckedAt(ist))->fSignal;
462	if (q>ql) {
463	c.fIndexMap[ind]=idx;
464	c.fIndexMap[ind+1]=ist;
465	} else {
466	break;
467	}
468	}
469	}
470
471	c.fMultiplicity++;
472
473	if (c.fMultiplicity >= 50 ) {
474	printf("FindCluster - multiplicity >50 %d \n",c.fMultiplicity);
475	c.fMultiplicity=50-1;
476	}
477
478	// Prepare center of gravity calculation
479	Float_t x, y;
480	fSegmentation->GetPadCxy(i, j, x, y);
481	c.fX += q*x;
482	c.fY += q*y;
483	c.fQ += q;
484	// Flag hit as taken
485	fHitMap->FlagHit(i,j);
486	//
487	// Now look recursively for all neighbours
488	//
489	Int_t nn;
490	Int_t Xlist[kMaxNeighbours], Ylist[kMaxNeighbours];
491	fSegmentation->Neighbours(i,j,&nn,Xlist,Ylist);
492	for (Int_t in=0; in<nn; in++) {
493	Int_t ix=Xlist[in];
494	Int_t iy=Ylist[in];
495	if (fHitMap->TestHit(ix,iy)==unused) FindCluster(ix, iy, c);
496	}
497	}
498	*/
499
500	//_____________________________________________________________________________
501
502	void AliMUONClusterFinderv0::FindRawClusters()
503	{
504	//
505	// simple MUON cluster finder from digits -- finds neighbours and
506	// fill the tree with raw clusters
507	//
508	if (!fNdigits) return;
509
510	fHitMap = new AliMUONHitMapA1(fSegmentation, fDigits);
511
512	AliMUONdigit *dig;
513
514	int ndig;
515	int nskip=0;
516
517	fHitMap->FillHits();
518	for (ndig=0; ndig<fNdigits; ndig++) {
519	dig = (AliMUONdigit*)fDigits->UncheckedAt(ndig);
520	Int_t i=dig->fPadX;
521	Int_t j=dig->fPadY;
522	if (fHitMap->TestHit(i,j)==used \|\|fHitMap->TestHit(i,j)==empty) {
523	nskip++;
524	continue;
525	}
526	AliMUONRawCluster c;
527	c.fMultiplicity=0;
528	// c.fPeakSignal=dig->fSignal;
529	// c.fTracks[0]=dig->fTracks[0];
530	// c.fTracks[1]=dig->fTracks[1];
531	// c.fTracks[2]=dig->fTracks[2];
532	c.fPeakSignal=0;
533	FindCluster(i,j, c);
534	// center of gravity
535	c.fX /= c.fQ;
536	c.fX=fSegmentation->GetAnod(c.fX);
537	c.fY /= c.fQ;
538	//
539	// apply correction to the coordinate along the anode wire
540	//
541
542
543	if (fCogCorr) {
544	Int_t ix,iy;
545	Float_t x=c.fX;
546	Float_t y=c.fY;
547	fSegmentation->GetPadIxy(x, y, ix, iy);
548	fSegmentation->GetPadCxy(ix, iy, x, y);
549	Float_t YonPad=(c.fY-y)/fSegmentation->Dpy();
550	c.fY=y-fCogCorr->Eval(YonPad,0,0);
551	}
552	//
553	// Analyse cluster and decluster if necessary
554	//
555	Decluster(&c);
556	//
557	//
558	//
559	// reset Cluster object
560	for (int k=0;k<c.fMultiplicity;k++) {
561	c.fIndexMap[k]=0;
562	c.fOffsetMap[k]=0;
563	}
564	c.fMultiplicity=0;
565	} // end loop ndig
566	delete fHitMap;
567	}
568
569	/*
570
571	void AliMUONClusterFinder::
572	CalibrateCOG()
573	{
574	Float_t x[5];
575	Float_t y[5];
576	Int_t n, i;
577	TF1 func;
578
579	if (fSegmentation) {
580	fSegmentation->GiveTestPoints(n, x, y);
581	for (i=0; i<n; i++) {
582	Float_t xtest=x[i];
583	Float_t ytest=y[i];
584	SinoidalFit(xtest, ytest, func);
585	}
586	fCogCorr = new TF1(func);
587	}
588	}
589	*/
590	/*
591
592	void AliMUONClusterFinder::
593	SinoidalFit(Float_t x, Float_t y, TF1 &func)
594	{
595	//
596	static Int_t count=0;
597	char canvasname[3];
598	count++;
599	sprintf(canvasname,"c%d",count);
600
601	// MANU : without const, error on HP
602	const Int_t ns=101;
603	Float_t xg[ns], yg[ns], xrg[ns], yrg[ns];
604	Float_t xsig[ns], ysig[ns];
605
606	AliMUONsegmentation *segmentation=fSegmentation;
607
608	Int_t ix,iy;
609	segmentation->GetPadIxy(x,y,ix,iy);
610	segmentation->GetPadCxy(ix,iy,x,y);
611	Int_t isec=segmentation->Sector(ix,iy);
612	// Pad Limits
613	Float_t xmin = x-segmentation->GetRealDpx(isec)/2;
614	Float_t ymin = y-segmentation->Dpy()/2;
615	//
616	// Integration Limits
617	Float_t dxI=fResponse->Nsigma()*fResponse->ChwX();
618	Float_t dyI=fResponse->Nsigma()*fResponse->ChwY();
619
620	//
621	// Scanning
622	//
623	Int_t i;
624	Float_t qp;
625	//
626	// y-position
627	Float_t yscan=ymin;
628	Float_t dy=segmentation->Dpy()/(ns-1);
629
630	for (i=0; i<ns; i++) {
631	//
632	// Pad Loop
633	//
634	Float_t sum=0;
635	Float_t qcheck=0;
636	segmentation->SigGenInit(x, yscan, 0);
637
638	for (segmentation->FirstPad(x, yscan, dxI, dyI);
639	segmentation->MorePads();
640	segmentation->NextPad())
641	{
642	qp=fResponse->IntXY(segmentation);
643	qp=TMath::Abs(qp);
644	//
645	//
646	if (qp > 1.e-4) {
647	qcheck+=qp;
648	Int_t ixs=segmentation->Ix();
649	Int_t iys=segmentation->Iy();
650	Float_t xs,ys;
651	segmentation->GetPadCxy(ixs,iys,xs,ys);
652	sum+=qp*ys;
653	}
654	} // Pad loop
655	Float_t ycog=sum/qcheck;
656	yg[i]=(yscan-y)/segmentation->Dpy();
657	yrg[i]=(ycog-y)/segmentation->Dpy();
658	ysig[i]=ycog-yscan;
659	yscan+=dy;
660	} // scan loop
661	//
662	// x-position
663	Float_t xscan=xmin;
664	Float_t dx=segmentation->GetRealDpx(isec)/(ns-1);
665
666	for (i=0; i<ns; i++) {
667	//
668	// Pad Loop
669	//
670	Float_t sum=0;
671	Float_t qcheck=0;
672	segmentation->SigGenInit(xscan, y, 0);
673
674	for (segmentation->FirstPad(xscan, y, dxI, dyI);
675	segmentation->MorePads();
676	segmentation->NextPad())
677	{
678	qp=fResponse->IntXY(segmentation);
679	qp=TMath::Abs(qp);
680	//
681	//
682	if (qp > 1.e-2) {
683	qcheck+=qp;
684	Int_t ixs=segmentation->Ix();
685	Int_t iys=segmentation->Iy();
686	Float_t xs,ys;
687	segmentation->GetPadCxy(ixs,iys,xs,ys);
688	sum+=qp*xs;
689	}
690	} // Pad loop
691	Float_t xcog=sum/qcheck;
692	xcog=segmentation->GetAnod(xcog);
693
694	xg[i]=(xscan-x)/segmentation->GetRealDpx(isec);
695	xrg[i]=(xcog-x)/segmentation->GetRealDpx(isec);
696	xsig[i]=xcog-xscan;
697	xscan+=dx;
698	}
699
700	TCanvas *c1=new TCanvas(canvasname,canvasname,400,10,600,700);
701	TPad* pad11 = new TPad("pad11"," ",0.01,0.51,0.49,0.99);
702	TPad* pad12 = new TPad("pad12"," ",0.51,0.51,0.99,0.99);
703	TPad* pad13 = new TPad("pad13"," ",0.01,0.01,0.49,0.49);
704	TPad* pad14 = new TPad("pad14"," ",0.51,0.01,0.99,0.49);
705	pad11->SetFillColor(11);
706	pad12->SetFillColor(11);
707	pad13->SetFillColor(11);
708	pad14->SetFillColor(11);
709	pad11->Draw();
710	pad12->Draw();
711	pad13->Draw();
712	pad14->Draw();
713	TGraph *graphx = new TGraph(ns,xg ,xsig);
714	TGraph *graphxr= new TGraph(ns,xrg,xsig);
715	TGraph *graphy = new TGraph(ns,yg ,ysig);
716	TGraph *graphyr= new TGraph(ns,yrg,ysig);
717	//
718	// Creates a Root function based on function sinoid above
719	// and perform the fit
720	//
721	Double_t sinoid(Double_t x, Double_t par);
722	TF1 *sinoidf = new TF1("sinoidf",sinoid,0.5,0.5,5);
723	graphyr->Fit("sinoidf","V");
724	sinoidf->Copy(func);
725	func.Eval(0,0,0);
726	//
727	pad11->cd();
728	graphx->SetFillColor(42);
729	graphx->SetMarkerColor(4);
730	graphx->SetMarkerStyle(21);
731	graphx->Draw("AC");
732	graphx->GetHistogram()->SetXTitle("x on pad");
733	graphx->GetHistogram()->SetYTitle("xcog-x");
734
735
736	pad12->cd();
737	graphxr->SetFillColor(42);
738	graphxr->SetMarkerColor(4);
739	graphxr->SetMarkerStyle(21);
740	graphxr->Draw("AP");
741	graphxr->GetHistogram()->SetXTitle("xcog on pad");
742	graphxr->GetHistogram()->SetYTitle("xcog-x");
743
744
745	pad13->cd();
746	graphy->SetFillColor(42);
747	graphy->SetMarkerColor(4);
748	graphy->SetMarkerStyle(21);
749	graphy->Draw("AF");
750	graphy->GetHistogram()->SetXTitle("y on pad");
751	graphy->GetHistogram()->SetYTitle("ycog-y");
752
753
754
755	pad14->cd();
756	graphyr->SetFillColor(42);
757	graphyr->SetMarkerColor(4);
758	graphyr->SetMarkerStyle(21);
759	graphyr->Draw("AF");
760	graphyr->GetHistogram()->SetXTitle("ycog on pad");
761	graphyr->GetHistogram()->SetYTitle("ycog-y");
762
763	c1->Update();
764
765	}
766	*/
767	/*
768	Double_t sinoid(Double_t x, Double_t par)
769	{
770	Double_t arg = -2TMath::Pi()x[0];
771	Double_t fitval= par[0]*TMath::Sin(arg)+
772	par[1]TMath::Sin(2arg)+
773	par[2]TMath::Sin(3arg)+
774	par[3]TMath::Sin(4arg)+
775	par[4]TMath::Sin(5arg);
776	return fitval;
777	}
778	*/
779
780
781
782
783
784
785