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

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