[u/mrichter/AliRoot.git] / HMPID / AliHMPIDCluster.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.                  *
//  **************************************************************************

#include <TVirtualFitter.h>  //Solve()
#include <TMinuit.h>         //Solve()
#include <TClonesArray.h>    //Solve()
#include <TMarker.h>         //Draw()

#include "AliLog.h"          //FitFunc()
#include "AliHMPIDCluster.h"  //class header

Bool_t AliHMPIDCluster::fgDoCorrSin=kTRUE;

ClassImp(AliHMPIDCluster)
    

void AliHMPIDCluster::SetClusterParams(Double_t xL,Double_t yL,Int_t iCh  )
{
  //------------------------------------------------------------------------
  //Set the cluster properties for the AliCluster3D part
  //------------------------------------------------------------------------

  //Get the volume ID from the previously set PNEntry
  UShort_t volId=AliGeomManager::LayerToVolUID(AliGeomManager::kHMPID,iCh);

  
  //get L->T cs matrix for a given chamber
  const TGeoHMatrix *t2l= AliGeomManager::GetTracking2LocalMatrix(volId);

  if(fParam->GetInstType())               //if there is no geometry we cannot retrieve the volId (only for monitoring)
  {
    new(this) AliCluster3D(); return;
  }
  

  //transformation from the pad cs to local
  xL -= 0.5*fParam->SizeAllX();      //size of all pads with dead zones included
  yL -= 0.5*fParam->SizeAllY();

  // Get the position in the tracking cs
  Double_t posL[3]={xL, yL, 0.};            //this is the LORS of HMPID
  Double_t posT[3];
  t2l->MasterToLocal(posL,posT);

 //Get the cluster covariance matrix in the tracking cs
  Double_t covL[9] = {
    0.8*0.8/12., 0.,            0.0,                 //pad size X
    0.,          0.84*0.84/12., 0.0,                 //pad size Y
    0.,          0.,            0.1,                 //just 1 , no Z dimension ???
  };
                
  TGeoHMatrix m;
  m.SetRotation(covL);
  m.Multiply(t2l);
  m.MultiplyLeft(&t2l->Inverse());
  Double_t *covT = m.GetRotationMatrix();

  new(this) AliCluster3D(volId,            // Can be done safer
       posT[0],posT[1],posT[2],
       covT[0],covT[1],covT[2],
	       covT[4],covT[5],
		       covT[8], 
			  0x0);            // No MC labels ?
}


AliHMPIDCluster::~AliHMPIDCluster(){
  if(fDigs) delete fDigs; fDigs=0;
  }

//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
void AliHMPIDCluster::CoG()
{
// Calculates naive cluster position as a center of gravity of its digits.
// Arguments: none 
//   Returns: none
  Int_t minPadX=999,minPadY=999,maxPadX=-1,maxPadY=-1;      //for box finding  
  if(fDigs==0) return;                                      //no digits in this cluster
  fXX=fYY=fQRaw=0;                                            //init summable parameters
  Int_t maxQpad=-1,maxQ=-1;                                 //to calculate the pad with the highest charge
  AliHMPIDDigit *pDig=0x0;
  for(Int_t iDig=0;iDig<fDigs->GetEntriesFast();iDig++){    //digits loop
    pDig=(AliHMPIDDigit*)fDigs->At(iDig);                   //get pointer to next digit

    if(pDig->PadPcX() > maxPadX) maxPadX = pDig->PadPcX();  // find the minimum box that contain the cluster  MaxX                            
    if(pDig->PadPcY() > maxPadY) maxPadY = pDig->PadPcY();  //                                                MaxY
    if(pDig->PadPcX() < minPadX) minPadX = pDig->PadPcX();  //                                                MinX   
    if(pDig->PadPcY() < minPadY) minPadY = pDig->PadPcY();  //                                                MinY   
    
    Float_t q=pDig->Q();                                    //get QDC 
    fXX += pDig->LorsX()*q;fYY +=pDig->LorsY()*q;             //add digit center weighted by QDC
    fQRaw+=q;                                               //increment total charge 
    if(q>maxQ) {maxQpad = pDig->Pad();maxQ=(Int_t)q;}       // to find pad with highest charge
  }//digits loop
  
  fBox=(maxPadX-minPadX+1)*100+maxPadY-minPadY+1;           // dimension of the box: format Xdim*100+Ydim
  
  if ( fQRaw != 0 ) {fXX/=fQRaw;fYY/=fQRaw;}                //final center of gravity
   
  if(fDigs->GetEntriesFast()>1&&fgDoCorrSin)CorrSin();       //correct it by sinoid   
  
  fQ  = fQRaw;                                              // Before starting fit procedure, Q and QRaw must be equal
  fCh=pDig->Ch();                                           //initialize chamber number
  fMaxQpad = maxQpad; fMaxQ=maxQ;                           //store max charge pad to the field
  fChi2=0;                                                  // no Chi2 to find
  fNlocMax=0;                                               // proper status from this method
  fSt=kCoG;
  
  if(fParam->GetInstType()) SetClusterParams(fXX,fYY,fCh);                              //need to fill the AliCluster3D part
 
}//CoG()
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
void AliHMPIDCluster::CorrSin() 
{
// Correction of cluster x position due to sinoid, see HMPID TDR  page 30
// Arguments: none
//   Returns: none
  Int_t pc,px,py;
  fParam->Lors2Pad(fXX,fYY,pc,px,py);             //tmp digit to get it center
  Float_t x=fXX-fParam->LorsX(pc,px);                    //diff between cluster x and center of the pad contaning this cluster   
  fXX+=3.31267e-2*TMath::Sin(2*TMath::Pi()/0.8*x)-2.66575e-3*TMath::Sin(4*TMath::Pi()/0.8*x)+2.80553e-3*TMath::Sin(6*TMath::Pi()/0.8*x)+0.0070;
}
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
void AliHMPIDCluster::Draw(Option_t*)
{
  TMarker *pMark=new TMarker(X(),Y(),5); pMark->SetUniqueID(fSt);pMark->SetMarkerColor(kBlue); pMark->Draw();
}
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
void AliHMPIDCluster::FitFunc(Int_t &iNpars, Double_t* deriv, Double_t &chi2, Double_t *par, Int_t iflag)
{
// Cluster fit function 
// par[0]=x par[1]=y par[2]=q for the first Mathieson shape
// par[3]=x par[4]=y par[5]=q for the second Mathieson shape and so on up to iNpars/3 Mathieson shapes
// For each pad of the cluster calculates the difference between actual pad charge and the charge induced to this pad by all Mathieson distributions
// Then the chi2 is calculated as the sum of this value squared for all pad in the cluster.  
// Arguments: iNpars - number of parameters which is number of local maxima of cluster * 3
//            chi2   - function result to be minimised 
//            par   - parameters array of size iNpars            
//   Returns: none  
  
  AliHMPIDCluster *pClu=(AliHMPIDCluster*)TVirtualFitter::GetFitter()->GetObjectFit();

  Int_t nPads = pClu->Size();  
  
  chi2 = 0;
  
  Int_t iNshape = iNpars/3;
  
  for(Int_t i=0;i<nPads;i++){                                                          //loop on all pads of the cluster
    Double_t dQpadMath = 0;
    for(Int_t j=0;j<iNshape;j++){                                                      //Mathiesons loop as all of them may contribute to this pad
      Double_t fracMathi = pClu->Dig(i)->IntMathieson(par[3*j],par[3*j+1]);
      dQpadMath+=par[3*j+2]*fracMathi;                                                 // par[3*j+2] is charge par[3*j] is x par[3*j+1] is y of current Mathieson
    }
    if(dQpadMath>0 && pClu->Dig(i)->Q()>0) {
      chi2 +=TMath::Power((pClu->Dig(i)->Q()-dQpadMath),2)/pClu->Dig(i)->Q();          //chi2 function to be minimized
    }
  }
//---calculate gradients...  
  if(iflag==2) {
    Double_t **derivPart;

    derivPart = new Double_t*[iNpars];

    for(Int_t j=0;j<iNpars;j++){                                                      
      deriv[j] = 0;
      derivPart[j] = new Double_t[nPads];
      for(Int_t i=0;i<nPads;i++){                                                          
        derivPart[j][i] = 0;
      }
    }

    for(Int_t i=0;i<nPads;i++){                                                          //loop on all pads of the cluster
      for(Int_t j=0;j<iNshape;j++){                                                      //Mathiesons loop as all of them may contribute to this pad
        Double_t fracMathi = pClu->Dig(i)->IntMathieson(par[3*j],par[3*j+1]);
        derivPart[3*j  ][i] += par[3*j+2]*(pClu->Dig(i)->Mathieson(par[3*j]-pClu->Dig(i)->LorsX()-0.5*AliHMPIDParam::SizePadX())-
                                           pClu->Dig(i)->Mathieson(par[3*j]-pClu->Dig(i)->LorsX()+0.5*AliHMPIDParam::SizePadX()))*
                                           pClu->Dig(i)->IntPartMathi(par[3*j+1],2);
        derivPart[3*j+1][i] += par[3*j+2]*(pClu->Dig(i)->Mathieson(par[3*j+1]-pClu->Dig(i)->LorsY()-0.5*AliHMPIDParam::SizePadY())-
                                           pClu->Dig(i)->Mathieson(par[3*j+1]-pClu->Dig(i)->LorsY()+0.5*AliHMPIDParam::SizePadY()))*
                                           pClu->Dig(i)->IntPartMathi(par[3*j],1);
        derivPart[3*j+2][i] += fracMathi;
      }
    }
                                                                                         //loop on all pads of the cluster     
    for(Int_t i=0;i<nPads;i++){                                                          //loop on all pads of the cluster
      Double_t dQpadMath = 0;                                                            //pad charge collector  
      for(Int_t j=0;j<iNshape;j++){                                                      //Mathiesons loop as all of them may contribute to this pad
        Double_t fracMathi = pClu->Dig(i)->IntMathieson(par[3*j],par[3*j+1]);
        dQpadMath+=par[3*j+2]*fracMathi;                                                 
        if(dQpadMath>0 && pClu->Dig(i)->Q()>0) {
          deriv[3*j]   += 2/pClu->Dig(i)->Q()*(pClu->Dig(i)->Q()-dQpadMath)*derivPart[3*j  ][i];
          deriv[3*j+1] += 2/pClu->Dig(i)->Q()*(pClu->Dig(i)->Q()-dQpadMath)*derivPart[3*j+1][i];
          deriv[3*j+2] += 2/pClu->Dig(i)->Q()*(pClu->Dig(i)->Q()-dQpadMath)*derivPart[3*j+2][i];
        }
      }
    }
    //delete array...
    for(Int_t i=0;i<iNpars;i++) delete [] derivPart[i]; delete [] derivPart;
  }
//---gradient calculations ended
  
}//FitFunction()
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
void AliHMPIDCluster::Print(Option_t* opt)const
{
//Print current cluster  
  const char *status=0;
  switch(fSt){
    case        kFrm  : status="formed        "   ;break;
    case        kUnf  : status="unfolded (fit)"   ;break;
    case        kCoG  : status="coged         "   ;break;
    case        kLo1  : status="locmax 1 (fit)"   ;break;
    case        kMax  : status="exceeded (cog)"   ;break;
    case        kNot  : status="not done (cog)"   ;break;
    case        kEmp  : status="empty         "   ;break;
    case        kEdg  : status="edge     (fit)"   ;break;
    case 	kSi1  : status="size 1   (cog)"   ;break;
    case 	kNoLoc: status="no LocMax(fit)"   ;break;
    case 	kAbn  : status="Abnormal fit  "   ;break;
    
    default:            status="??????"          ;break;   
  }
  Double_t ratio=0;
  if(Q()>0&&QRaw()>0) ratio = Q()/QRaw()*100;
  Printf("%sCLU: ch=%i                 (%7.3f,%7.3f) Q=%8.3f Qraw=%8.3f(%3.0f%%) Size=%2i DimBox=%i LocMax=%i Chi2=%7.3f   %s",
         opt,Ch(),X(),Y(),Q(),QRaw(),ratio,Size(),fBox,fNlocMax,fChi2,status);
  if(fDigs) fDigs->Print();    
}//Print()
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Int_t AliHMPIDCluster::Solve(TClonesArray *pCluLst,Bool_t isTryUnfold)
{
//This methode is invoked when the cluster is formed to solve it. Solve the cluster means to try to unfold the cluster
//into the local maxima number of clusters. This methode is invoked by AliHMPIDRconstructor::Dig2Clu() on cluster by cluster basis.  
//At this point, cluster contains a list of digits, cluster charge and size is precalculated in AddDigit(), position is preset to (-1,-1) in ctor,
//status is preset to kFormed in AddDigit(), chamber-sector info is preseted to actual values in AddDigit()
//Method first finds number of local maxima and if it's more then one tries to unfold this cluster into local maxima number of clusters
//Arguments: pCluLst     - cluster list pointer where to add new cluster(s)
//           isTryUnfold - flag to switch on/off unfolding   
//  Returns: number of local maxima of original cluster
  const Int_t kMaxLocMax=6;                                                              //max allowed number of loc max for fitting
//  
  CoG();                                                                                 //First calculate CoG for the given cluster
  Int_t iCluCnt=pCluLst->GetEntriesFast();                                               //get current number of clusters already stored in the list by previous operations
  if(isTryUnfold==kFALSE || Size()==1) {                                                 //if cluster contains single pad there is no way to improve the knowledge 
    fSt = (isTryUnfold)? kSi1: kNot;
    if(fParam->GetInstType()) SetClusterParams(fXX,fYY,fCh);  
    new ((*pCluLst)[iCluCnt++]) AliHMPIDCluster(*this);  //add this raw cluster 
    return 1;
  } 
  
//Phase 0. Initialise Fitter  
  Double_t arglist[10];
  Int_t ierflg = 0;
  TVirtualFitter *fitter = TVirtualFitter::Fitter(this,3*6);                            //initialize Fitter

  arglist[0] = -1;
  ierflg = fitter->ExecuteCommand("SET PRI", arglist, 1);                               // no printout
  ierflg = fitter->ExecuteCommand("SET NOW", arglist, 0);                               //no warning messages
  arglist[0] =  1;
  ierflg = fitter->ExecuteCommand("SET GRA", arglist, 1);                               //force Fitter to use my gradient

  fitter->SetFCN(AliHMPIDCluster::FitFunc);

//  arglist[0] = 1;
//  ierflg = fitter->ExecuteCommand("SET ERR", arglist ,1);
  
// Set starting values and step sizes for parameters
    
//Phase 1. Find number of local maxima. Strategy is to check if the current pad has QDC more then all neigbours. Also find the box contaning the cluster   
  fNlocMax=0;

  for(Int_t iDig1=0;iDig1<Size();iDig1++) {                                               //first digits loop
    
    AliHMPIDDigit *pDig1 = Dig(iDig1);                                                   //take next digit    
    Int_t iCnt = 0;                                                                      //counts how many neighbouring pads has QDC more then current one
    
    for(Int_t iDig2=0;iDig2<Size();iDig2++) {                                            //loop on all digits again
      
      if(iDig1==iDig2) continue;                                                         //the same digit, no need to compare 
      AliHMPIDDigit *pDig2 = Dig(iDig2);                                                 //take second digit to compare with the first one
      Int_t dist = TMath::Sign(Int_t(pDig1->PadChX()-pDig2->PadChX()),1)+TMath::Sign(Int_t(pDig1->PadChY()-pDig2->PadChY()),1);//distance between pads
      if(dist==1)                                                                        //means dig2 is a neighbour of dig1
         if(pDig2->Q()>=pDig1->Q()) iCnt++;                                              //count number of pads with Q more then Q of current pad
      
    }//second digits loop
    
    if(iCnt==0&&fNlocMax<kMaxLocMax){                                                    //this pad has Q more then any neighbour so it's local maximum
      
      Double_t xStart=pDig1->LorsX();Double_t yStart=pDig1->LorsY();
      Double_t xMin=xStart-fParam->SizePadX();
      Double_t xMax=xStart+fParam->SizePadX();
      Double_t yMin=yStart-fParam->SizePadY();
      Double_t yMax=yStart+fParam->SizePadY();
      
      ierflg = fitter->SetParameter(3*fNlocMax  ,Form("x%i",fNlocMax),xStart,0.1,xMin,xMax);    // X,Y,Q initial values of the loc max pad
      ierflg = fitter->SetParameter(3*fNlocMax+1,Form("y%i",fNlocMax),yStart,0.1,yMin,yMax);    // X, Y constrained to be near the loc max
      ierflg = fitter->SetParameter(3*fNlocMax+2,Form("q%i",fNlocMax),pDig1->Q(),0.1,0,10000);  // Q constrained to be positive
      
      fNlocMax++;
      
    }//if this pad is local maximum
  }//first digits loop
  
//Phase 2. Fit loc max number of Mathiesons or add this current cluster to the list
// case 1 -> no loc max found
 if ( fNlocMax == 0) {                                                                       // case of no local maxima found: pads with same charge...
   
   ierflg = fitter->SetParameter(3*fNlocMax  ,Form("x%i",fNlocMax),fXX,0.1,0,0);              // Init values taken from CoG() -> fXX,fYY,fQRaw
   ierflg = fitter->SetParameter(3*fNlocMax+1,Form("y%i",fNlocMax),fYY,0.1,0,0);              //
   ierflg = fitter->SetParameter(3*fNlocMax+2,Form("q%i",fNlocMax),fQRaw,0.1,0,10000);       //
   
   fNlocMax = 1;
   fSt=kNoLoc;
 }

// case 2 -> loc max found. Check # of loc maxima 
 if ( fNlocMax >= kMaxLocMax)  { 
 if(fParam->GetInstType()) SetClusterParams(fXX,fYY,fCh);                                                           // if # of local maxima exceeds kMaxLocMax...
   fSt = kMax;   new ((*pCluLst)[iCluCnt++]) AliHMPIDCluster(*this);                      //...add this raw cluster  
   } else {                                                                               //or resonable number of local maxima to fit and user requested it
  // Now ready for minimization step
   arglist[0] = 500;                                                                      //number of steps and sigma on pads charges
   arglist[1] = 1.;                                                                       //

   ierflg = fitter->ExecuteCommand("SIMPLEX",arglist,2);                                  //start fitting with Simplex
   if (!ierflg)
     fitter->ExecuteCommand("MIGRAD" ,arglist,2);                                         //fitting improved by Migrad
   if(ierflg) {
     Double_t strategy=2;
     ierflg = fitter->ExecuteCommand("SET STR",&strategy,1);                              //change level of strategy 
     if(!ierflg) {
       ierflg = fitter->ExecuteCommand("SIMPLEX",arglist,2);                              //start fitting with Simplex
       if (!ierflg)
         fitter->ExecuteCommand("MIGRAD" ,arglist,2);                                     //fitting improved by Migrad
     }
   }        
   if(ierflg) fSt=kAbn;                                                                   //no convergence of the fit...
   Double_t dummy; char sName[80];                                                        //vars to get results from Minuit
   Double_t edm, errdef;
   Int_t nvpar, nparx;
  
   for(Int_t i=0;i<fNlocMax;i++){                                                        //store the local maxima parameters
     fitter->GetParameter(3*i   ,sName,  fXX, fErrX , dummy, dummy);                      // X
     fitter->GetParameter(3*i+1 ,sName,  fYY, fErrY , dummy, dummy);                      // Y
     fitter->GetParameter(3*i+2 ,sName,  fQ, fErrQ , dummy, dummy);                      // Q
     fitter->GetStats(fChi2, edm, errdef, nvpar, nparx);                                 //get fit infos
      if(fSt!=kAbn) {         
        if(fNlocMax!=1)fSt=kUnf;                                                         // if unfolded
        if(fNlocMax==1&&fSt!=kNoLoc) fSt=kLo1;                                           // if only 1 loc max
        if ( !IsInPc()) fSt = kEdg;                                                      // if Out of Pc
        if(fSt==kNoLoc) fNlocMax=0;                                                      // if with no loc max (pads with same charge..)
      }
      if(fParam->GetInstType()) SetClusterParams(fXX,fYY,fCh);                                                      //need to fill the AliCluster3D part
      new ((*pCluLst)[iCluCnt++]) AliHMPIDCluster(*this);	                         //add new unfolded cluster
      
   }
 }

 return fNlocMax;
 
}//Solve()
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Commit	Line	Data
d3da6dc4	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
76fd1a96	16	#include <TVirtualFitter.h> //Solve()
d3da6dc4	17	#include <TMinuit.h> //Solve()
d3da6dc4	18	#include <TClonesArray.h> //Solve()
d1bf51e1	19	#include <TMarker.h> //Draw()
27311693	20
76fd1a96	21	#include "AliLog.h" //FitFunc()
76fd1a96	22	#include "AliHMPIDCluster.h" //class header
76fd1a96	23
69ed32de	24	Bool_t AliHMPIDCluster::fgDoCorrSin=kTRUE;
27311693	25
d3da6dc4	26	ClassImp(AliHMPIDCluster)
76fd1a96	27
55a829a5	28
	29	void AliHMPIDCluster::SetClusterParams(Double_t xL,Double_t yL,Int_t iCh )
	30	{
	31	//------------------------------------------------------------------------
	32	//Set the cluster properties for the AliCluster3D part
	33	//------------------------------------------------------------------------
	34
	35	//Get the volume ID from the previously set PNEntry
	36	UShort_t volId=AliGeomManager::LayerToVolUID(AliGeomManager::kHMPID,iCh);
	37
	38
	39	//get L->T cs matrix for a given chamber
	40	const TGeoHMatrix *t2l= AliGeomManager::GetTracking2LocalMatrix(volId);
	41
b38ac33a	42	if(fParam->GetInstType()) //if there is no geometry we cannot retrieve the volId (only for monitoring)
b87365d5	43	{
	44	new(this) AliCluster3D(); return;
	45	}
55a829a5	46
	47
	48	//transformation from the pad cs to local
b38ac33a	49	xL -= 0.5*fParam->SizeAllX(); //size of all pads with dead zones included
b38ac33a	50	yL -= 0.5*fParam->SizeAllY();
55a829a5	51
	52	// Get the position in the tracking cs
	53	Double_t posL[3]={xL, yL, 0.}; //this is the LORS of HMPID
	54	Double_t posT[3];
	55	t2l->MasterToLocal(posL,posT);
	56
	57	//Get the cluster covariance matrix in the tracking cs
	58	Double_t covL[9] = {
	59	0.8*0.8/12., 0., 0.0, //pad size X
	60	0., 0.84*0.84/12., 0.0, //pad size Y
	61	0., 0., 0.1, //just 1 , no Z dimension ???
	62	};
	63
	64	TGeoHMatrix m;
	65	m.SetRotation(covL);
	66	m.Multiply(t2l);
	67	m.MultiplyLeft(&t2l->Inverse());
	68	Double_t *covT = m.GetRotationMatrix();
	69
	70	new(this) AliCluster3D(volId, // Can be done safer
	71	posT[0],posT[1],posT[2],
	72	covT[0],covT[1],covT[2],
	73	covT[4],covT[5],
	74	covT[8],
	75	0x0); // No MC labels ?
	76	}
	77
	78
	79	AliHMPIDCluster::~AliHMPIDCluster(){
	80	if(fDigs) delete fDigs; fDigs=0;
	81	}
	82
d3da6dc4	83	//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	84	void AliHMPIDCluster::CoG()
	85	{
	86	// Calculates naive cluster position as a center of gravity of its digits.
	87	// Arguments: none
d1bf51e1	88	// Returns: none
e4a3eae8	89	Int_t minPadX=999,minPadY=999,maxPadX=-1,maxPadY=-1; //for box finding
e4a3eae8	90	if(fDigs==0) return; //no digits in this cluster
894758d4	91	fXX=fYY=fQRaw=0; //init summable parameters
e4a3eae8	92	Int_t maxQpad=-1,maxQ=-1; //to calculate the pad with the highest charge
aa85549f	93	AliHMPIDDigit *pDig=0x0;
e4a3eae8	94	for(Int_t iDig=0;iDig<fDigs->GetEntriesFast();iDig++){ //digits loop
	95	pDig=(AliHMPIDDigit*)fDigs->At(iDig); //get pointer to next digit
	96
	97	if(pDig->PadPcX() > maxPadX) maxPadX = pDig->PadPcX(); // find the minimum box that contain the cluster MaxX
	98	if(pDig->PadPcY() > maxPadY) maxPadY = pDig->PadPcY(); // MaxY
	99	if(pDig->PadPcX() < minPadX) minPadX = pDig->PadPcX(); // MinX
	100	if(pDig->PadPcY() < minPadY) minPadY = pDig->PadPcY(); // MinY
	101
	102	Float_t q=pDig->Q(); //get QDC
894758d4	103	fXX += pDig->LorsX()q;fYY +=pDig->LorsY()q; //add digit center weighted by QDC
e4a3eae8	104	fQRaw+=q; //increment total charge
e4a3eae8	105	if(q>maxQ) {maxQpad = pDig->Pad();maxQ=(Int_t)q;} // to find pad with highest charge
d3da6dc4	106	}//digits loop
e4a3eae8	107
	108	fBox=(maxPadX-minPadX+1)100+maxPadY-minPadY+1; // dimension of the box: format Xdim100+Ydim
	109
3a94c321	110	if ( fQRaw != 0 ) {fXX/=fQRaw;fYY/=fQRaw;} //final center of gravity
27311693	111
69ed32de	112	if(fDigs->GetEntriesFast()>1&&fgDoCorrSin)CorrSin(); //correct it by sinoid
e4a3eae8	113
c5c19d6a	114	fQ = fQRaw; // Before starting fit procedure, Q and QRaw must be equal
	115	fCh=pDig->Ch(); //initialize chamber number
	116	fMaxQpad = maxQpad; fMaxQ=maxQ; //store max charge pad to the field
	117	fChi2=0; // no Chi2 to find
	118	fNlocMax=0; // proper status from this method
d3da6dc4	119	fSt=kCoG;
55a829a5	120
b38ac33a	121	if(fParam->GetInstType()) SetClusterParams(fXX,fYY,fCh); //need to fill the AliCluster3D part
55a829a5	122
d3da6dc4	123	}//CoG()
	124	//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	125	void AliHMPIDCluster::CorrSin()
	126	{
	127	// Correction of cluster x position due to sinoid, see HMPID TDR page 30
	128	// Arguments: none
	129	// Returns: none
1d4857c5	130	Int_t pc,px,py;
b38ac33a	131	fParam->Lors2Pad(fXX,fYY,pc,px,py); //tmp digit to get it center
b38ac33a	132	Float_t x=fXX-fParam->LorsX(pc,px); //diff between cluster x and center of the pad contaning this cluster
894758d4	133	fXX+=3.31267e-2TMath::Sin(2TMath::Pi()/0.8x)-2.66575e-3TMath::Sin(4TMath::Pi()/0.8x)+2.80553e-3TMath::Sin(6TMath::Pi()/0.8*x)+0.0070;
d3da6dc4	134	}
d3da6dc4	135	//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
d1bf51e1	136	void AliHMPIDCluster::Draw(Option_t*)
d1bf51e1	137	{
a1d55ff3	138	TMarker *pMark=new TMarker(X(),Y(),5); pMark->SetUniqueID(fSt);pMark->SetMarkerColor(kBlue); pMark->Draw();
d1bf51e1	139	}
d1bf51e1	140	//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
557ca324	141	void AliHMPIDCluster::FitFunc(Int_t &iNpars, Double_t* deriv, Double_t &chi2, Double_t *par, Int_t iflag)
d3da6dc4	142	{
	143	// Cluster fit function
	144	// par[0]=x par[1]=y par[2]=q for the first Mathieson shape
	145	// par[3]=x par[4]=y par[5]=q for the second Mathieson shape and so on up to iNpars/3 Mathieson shapes
76fd1a96	146	// For each pad of the cluster calculates the difference between actual pad charge and the charge induced to this pad by all Mathieson distributions
d3da6dc4	147	// Then the chi2 is calculated as the sum of this value squared for all pad in the cluster.
	148	// Arguments: iNpars - number of parameters which is number of local maxima of cluster * 3
	149	// chi2 - function result to be minimised
	150	// par - parameters array of size iNpars
	151	// Returns: none
76fd1a96	152
	153	AliHMPIDCluster pClu=(AliHMPIDCluster)TVirtualFitter::GetFitter()->GetObjectFit();
	154
	155	Int_t nPads = pClu->Size();
76fd1a96	156
d3da6dc4	157	chi2 = 0;
76fd1a96	158
	159	Int_t iNshape = iNpars/3;
	160
76fd1a96	161	for(Int_t i=0;i<nPads;i++){ //loop on all pads of the cluster
	162	Double_t dQpadMath = 0;
	163	for(Int_t j=0;j<iNshape;j++){ //Mathiesons loop as all of them may contribute to this pad
	164	Double_t fracMathi = pClu->Dig(i)->IntMathieson(par[3j],par[3j+1]);
	165	dQpadMath+=par[3j+2]fracMathi; // par[3j+2] is charge par[3j] is x par[3*j+1] is y of current Mathieson
76fd1a96	166	}
	167	if(dQpadMath>0 && pClu->Dig(i)->Q()>0) {
	168	chi2 +=TMath::Power((pClu->Dig(i)->Q()-dQpadMath),2)/pClu->Dig(i)->Q(); //chi2 function to be minimized
	169	}
	170	}
557ca324	171	//---calculate gradients...
	172	if(iflag==2) {
	173	Double_t **derivPart;
	174
	175	derivPart = new Double_t*[iNpars];
	176
	177	for(Int_t j=0;j<iNpars;j++){
	178	deriv[j] = 0;
	179	derivPart[j] = new Double_t[nPads];
	180	for(Int_t i=0;i<nPads;i++){
	181	derivPart[j][i] = 0;
	182	}
	183	}
	184
	185	for(Int_t i=0;i<nPads;i++){ //loop on all pads of the cluster
	186	for(Int_t j=0;j<iNshape;j++){ //Mathiesons loop as all of them may contribute to this pad
	187	Double_t fracMathi = pClu->Dig(i)->IntMathieson(par[3j],par[3j+1]);
	188	derivPart[3j ][i] += par[3j+2](pClu->Dig(i)->Mathieson(par[3j]-pClu->Dig(i)->LorsX()-0.5*AliHMPIDParam::SizePadX())-
	189	pClu->Dig(i)->Mathieson(par[3j]-pClu->Dig(i)->LorsX()+0.5AliHMPIDParam::SizePadX()))*
	190	pClu->Dig(i)->IntPartMathi(par[3*j+1],2);
	191	derivPart[3j+1][i] += par[3j+2](pClu->Dig(i)->Mathieson(par[3j+1]-pClu->Dig(i)->LorsY()-0.5*AliHMPIDParam::SizePadY())-
	192	pClu->Dig(i)->Mathieson(par[3j+1]-pClu->Dig(i)->LorsY()+0.5AliHMPIDParam::SizePadY()))*
	193	pClu->Dig(i)->IntPartMathi(par[3*j],1);
	194	derivPart[3*j+2][i] += fracMathi;
	195	}
	196	}
	197	//loop on all pads of the cluster
	198	for(Int_t i=0;i<nPads;i++){ //loop on all pads of the cluster
	199	Double_t dQpadMath = 0; //pad charge collector
	200	for(Int_t j=0;j<iNshape;j++){ //Mathiesons loop as all of them may contribute to this pad
	201	Double_t fracMathi = pClu->Dig(i)->IntMathieson(par[3j],par[3j+1]);
	202	dQpadMath+=par[3j+2]fracMathi;
	203	if(dQpadMath>0 && pClu->Dig(i)->Q()>0) {
	204	deriv[3j] += 2/pClu->Dig(i)->Q()(pClu->Dig(i)->Q()-dQpadMath)derivPart[3j ][i];
	205	deriv[3j+1] += 2/pClu->Dig(i)->Q()(pClu->Dig(i)->Q()-dQpadMath)derivPart[3j+1][i];
	206	deriv[3j+2] += 2/pClu->Dig(i)->Q()(pClu->Dig(i)->Q()-dQpadMath)derivPart[3j+2][i];
	207	}
76fd1a96	208	}
d3da6dc4	209	}
557ca324	210	//delete array...
557ca324	211	for(Int_t i=0;i<iNpars;i++) delete [] derivPart[i]; delete [] derivPart;
76fd1a96	212	}
557ca324	213	//---gradient calculations ended
76fd1a96	214
d3da6dc4	215	}//FitFunction()
	216	//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	217	void AliHMPIDCluster::Print(Option_t* opt)const
	218	{
	219	//Print current cluster
	220	const char *status=0;
	221	switch(fSt){
d1bf51e1	222	case kFrm : status="formed " ;break;
	223	case kUnf : status="unfolded (fit)" ;break;
	224	case kCoG : status="coged " ;break;
	225	case kLo1 : status="locmax 1 (fit)" ;break;
d1bf51e1	226	case kMax : status="exceeded (cog)" ;break;
	227	case kNot : status="not done (cog)" ;break;
	228	case kEmp : status="empty " ;break;
	229	case kEdg : status="edge (fit)" ;break;
	230	case kSi1 : status="size 1 (cog)" ;break;
	231	case kNoLoc: status="no LocMax(fit)" ;break;
c5c19d6a	232	case kAbn : status="Abnormal fit " ;break;
d1bf51e1	233
d1bf51e1	234	default: status="??????" ;break;
d3da6dc4	235	}
e4a3eae8	236	Double_t ratio=0;
e4a3eae8	237	if(Q()>0&&QRaw()>0) ratio = Q()/QRaw()*100;
1d4857c5	238	Printf("%sCLU: ch=%i (%7.3f,%7.3f) Q=%8.3f Qraw=%8.3f(%3.0f%%) Size=%2i DimBox=%i LocMax=%i Chi2=%7.3f %s",
1006986d	239	opt,Ch(),X(),Y(),Q(),QRaw(),ratio,Size(),fBox,fNlocMax,fChi2,status);
d1bf51e1	240	if(fDigs) fDigs->Print();
d3da6dc4	241	}//Print()
	242	//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	243	Int_t AliHMPIDCluster::Solve(TClonesArray *pCluLst,Bool_t isTryUnfold)
	244	{
	245	//This methode is invoked when the cluster is formed to solve it. Solve the cluster means to try to unfold the cluster
	246	//into the local maxima number of clusters. This methode is invoked by AliHMPIDRconstructor::Dig2Clu() on cluster by cluster basis.
	247	//At this point, cluster contains a list of digits, cluster charge and size is precalculated in AddDigit(), position is preset to (-1,-1) in ctor,
	248	//status is preset to kFormed in AddDigit(), chamber-sector info is preseted to actual values in AddDigit()
	249	//Method first finds number of local maxima and if it's more then one tries to unfold this cluster into local maxima number of clusters
	250	//Arguments: pCluLst - cluster list pointer where to add new cluster(s)
	251	// isTryUnfold - flag to switch on/off unfolding
	252	// Returns: number of local maxima of original cluster
76fd1a96	253	const Int_t kMaxLocMax=6; //max allowed number of loc max for fitting
	254	//
	255	CoG(); //First calculate CoG for the given cluster
e4a3eae8	256	Int_t iCluCnt=pCluLst->GetEntriesFast(); //get current number of clusters already stored in the list by previous operations
e4a3eae8	257	if(isTryUnfold==kFALSE \|\| Size()==1) { //if cluster contains single pad there is no way to improve the knowledge
55a829a5	258	fSt = (isTryUnfold)? kSi1: kNot;
b38ac33a	259	if(fParam->GetInstType()) SetClusterParams(fXX,fYY,fCh);
d1bf51e1	260	new ((pCluLst)[iCluCnt++]) AliHMPIDCluster(this); //add this raw cluster
	261	return 1;
	262	}
76fd1a96	263
	264	//Phase 0. Initialise Fitter
	265	Double_t arglist[10];
	266	Int_t ierflg = 0;
	267	TVirtualFitter fitter = TVirtualFitter::Fitter(this,36); //initialize Fitter
	268
	269	arglist[0] = -1;
	270	ierflg = fitter->ExecuteCommand("SET PRI", arglist, 1); // no printout
	271	ierflg = fitter->ExecuteCommand("SET NOW", arglist, 0); //no warning messages
	272	arglist[0] = 1;
	273	ierflg = fitter->ExecuteCommand("SET GRA", arglist, 1); //force Fitter to use my gradient
	274
	275	fitter->SetFCN(AliHMPIDCluster::FitFunc);
	276
	277	// arglist[0] = 1;
	278	// ierflg = fitter->ExecuteCommand("SET ERR", arglist ,1);
	279
	280	// Set starting values and step sizes for parameters
	281
d1bf51e1	282	//Phase 1. Find number of local maxima. Strategy is to check if the current pad has QDC more then all neigbours. Also find the box contaning the cluster
d1bf51e1	283	fNlocMax=0;
e4a3eae8	284
76fd1a96	285	for(Int_t iDig1=0;iDig1<Size();iDig1++) { //first digits loop
76fd1a96	286
e4a3eae8	287	AliHMPIDDigit *pDig1 = Dig(iDig1); //take next digit
c5c19d6a	288	Int_t iCnt = 0; //counts how many neighbouring pads has QDC more then current one
76fd1a96	289
d1bf51e1	290	for(Int_t iDig2=0;iDig2<Size();iDig2++) { //loop on all digits again
76fd1a96	291
d1bf51e1	292	if(iDig1==iDig2) continue; //the same digit, no need to compare
d3da6dc4	293	AliHMPIDDigit *pDig2 = Dig(iDig2); //take second digit to compare with the first one
da08475b	294	Int_t dist = TMath::Sign(Int_t(pDig1->PadChX()-pDig2->PadChX()),1)+TMath::Sign(Int_t(pDig1->PadChY()-pDig2->PadChY()),1);//distance between pads
d1bf51e1	295	if(dist==1) //means dig2 is a neighbour of dig1
c5c19d6a	296	if(pDig2->Q()>=pDig1->Q()) iCnt++; //count number of pads with Q more then Q of current pad
76fd1a96	297
d3da6dc4	298	}//second digits loop
76fd1a96	299
c5c19d6a	300	if(iCnt==0&&fNlocMax<kMaxLocMax){ //this pad has Q more then any neighbour so it's local maximum
76fd1a96	301
c5c19d6a	302	Double_t xStart=pDig1->LorsX();Double_t yStart=pDig1->LorsY();
b38ac33a	303	Double_t xMin=xStart-fParam->SizePadX();
	304	Double_t xMax=xStart+fParam->SizePadX();
	305	Double_t yMin=yStart-fParam->SizePadY();
	306	Double_t yMax=yStart+fParam->SizePadY();
76fd1a96	307
	308	ierflg = fitter->SetParameter(3*fNlocMax ,Form("x%i",fNlocMax),xStart,0.1,xMin,xMax); // X,Y,Q initial values of the loc max pad
	309	ierflg = fitter->SetParameter(3*fNlocMax+1,Form("y%i",fNlocMax),yStart,0.1,yMin,yMax); // X, Y constrained to be near the loc max
343f4211	310	ierflg = fitter->SetParameter(3*fNlocMax+2,Form("q%i",fNlocMax),pDig1->Q(),0.1,0,10000); // Q constrained to be positive
76fd1a96	311
d1bf51e1	312	fNlocMax++;
76fd1a96	313
d3da6dc4	314	}//if this pad is local maximum
d3da6dc4	315	}//first digits loop
d1bf51e1	316
d3da6dc4	317	//Phase 2. Fit loc max number of Mathiesons or add this current cluster to the list
e4a3eae8	318	// case 1 -> no loc max found
76fd1a96	319	if ( fNlocMax == 0) { // case of no local maxima found: pads with same charge...
76fd1a96	320
894758d4	321	ierflg = fitter->SetParameter(3*fNlocMax ,Form("x%i",fNlocMax),fXX,0.1,0,0); // Init values taken from CoG() -> fXX,fYY,fQRaw
894758d4	322	ierflg = fitter->SetParameter(3*fNlocMax+1,Form("y%i",fNlocMax),fYY,0.1,0,0); //
343f4211	323	ierflg = fitter->SetParameter(3*fNlocMax+2,Form("q%i",fNlocMax),fQRaw,0.1,0,10000); //
76fd1a96	324
d1bf51e1	325	fNlocMax = 1;
	326	fSt=kNoLoc;
	327	}
e4a3eae8	328
e4a3eae8	329	// case 2 -> loc max found. Check # of loc maxima
55a829a5	330	if ( fNlocMax >= kMaxLocMax) {
b38ac33a	331	if(fParam->GetInstType()) SetClusterParams(fXX,fYY,fCh); // if # of local maxima exceeds kMaxLocMax...
76fd1a96	332	fSt = kMax; new ((pCluLst)[iCluCnt++]) AliHMPIDCluster(this); //...add this raw cluster
	333	} else { //or resonable number of local maxima to fit and user requested it
	334	// Now ready for minimization step
	335	arglist[0] = 500; //number of steps and sigma on pads charges
	336	arglist[1] = 1.; //
	337
	338	ierflg = fitter->ExecuteCommand("SIMPLEX",arglist,2); //start fitting with Simplex
	339	if (!ierflg)
	340	fitter->ExecuteCommand("MIGRAD" ,arglist,2); //fitting improved by Migrad
	341	if(ierflg) {
c5c19d6a	342	Double_t strategy=2;
76fd1a96	343	ierflg = fitter->ExecuteCommand("SET STR",&strategy,1); //change level of strategy
	344	if(!ierflg) {
	345	ierflg = fitter->ExecuteCommand("SIMPLEX",arglist,2); //start fitting with Simplex
	346	if (!ierflg)
	347	fitter->ExecuteCommand("MIGRAD" ,arglist,2); //fitting improved by Migrad
c5c19d6a	348	}
c5c19d6a	349	}
76fd1a96	350	if(ierflg) fSt=kAbn; //no convergence of the fit...
	351	Double_t dummy; char sName[80]; //vars to get results from Minuit
	352	Double_t edm, errdef;
	353	Int_t nvpar, nparx;
	354
e4a3eae8	355	for(Int_t i=0;i<fNlocMax;i++){ //store the local maxima parameters
894758d4	356	fitter->GetParameter(3*i ,sName, fXX, fErrX , dummy, dummy); // X
894758d4	357	fitter->GetParameter(3*i+1 ,sName, fYY, fErrY , dummy, dummy); // Y
76fd1a96	358	fitter->GetParameter(3*i+2 ,sName, fQ, fErrQ , dummy, dummy); // Q
76fd1a96	359	fitter->GetStats(fChi2, edm, errdef, nvpar, nparx); //get fit infos
c5c19d6a	360	if(fSt!=kAbn) {
76fd1a96	361	if(fNlocMax!=1)fSt=kUnf; // if unfolded
	362	if(fNlocMax==1&&fSt!=kNoLoc) fSt=kLo1; // if only 1 loc max
	363	if ( !IsInPc()) fSt = kEdg; // if Out of Pc
	364	if(fSt==kNoLoc) fNlocMax=0; // if with no loc max (pads with same charge..)
c5c19d6a	365	}
b38ac33a	366	if(fParam->GetInstType()) SetClusterParams(fXX,fYY,fCh); //need to fill the AliCluster3D part
e4a3eae8	367	new ((pCluLst)[iCluCnt++]) AliHMPIDCluster(this); //add new unfolded cluster
55a829a5	368
e4a3eae8	369	}
e4a3eae8	370	}
d1bf51e1	371
e4a3eae8	372	return fNlocMax;
e4a3eae8	373
d3da6dc4	374	}//Solve()
d3da6dc4	375	//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++