[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 "AliHMPIDCluster.h"  //class header
#include <TVirtualFitter.h>  //Solve()
#include <TMinuit.h>         //Solve()
#include <TClonesArray.h>    //Solve()
#include <TMarker.h>         //Draw()

#include "AliLog.h"          //FindCusterSize()

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

  fParam = AliHMPIDParam::Instance();

  //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;
  if(fParam) delete fParam; fParam=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
  fCh = -1;                                                 //init chamber
  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) continue;                                     //protection
    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
    fCh=pDig->Ch();                                         //initialize chamber number
  }//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
  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(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)->MathiesonX(par[3*j]-pClu->Dig(i)->LorsX()-0.5*AliHMPIDParam::SizePadX())-
                                           pClu->Dig(i)->MathiesonX(par[3*j]-pClu->Dig(i)->LorsX()+0.5*AliHMPIDParam::SizePadX()))*
                                           pClu->Dig(i)->IntPartMathiY(par[3*j+1]);
        derivPart[3*j+1][i] += par[3*j+2]*(pClu->Dig(i)->MathiesonY(par[3*j+1]-pClu->Dig(i)->LorsY()-0.5*AliHMPIDParam::SizePadY())-
                                           pClu->Dig(i)->MathiesonY(par[3*j+1]-pClu->Dig(i)->LorsY()+0.5*AliHMPIDParam::SizePadY()))*
                                           pClu->Dig(i)->IntPartMathiX(par[3*j]);
        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

// fit ended. Final calculations
  
  
}//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;
    case 	kBig  : status="Big Clu(>100) "   ;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,Int_t *pSigmaCut, 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
  
  Int_t rawSize = Size();                                                                //get current raw cluster size
  
  if(rawSize>100) {
    fSt = kBig;
  } else if(isTryUnfold==kFALSE) {
    fSt = kNot;
  } else if(rawSize==1) {
    fSt = kSi1;
  }
  
  if(rawSize>100 || isTryUnfold==kFALSE || rawSize==1) {                                 //No deconv if: 1 - big cluster (also avoid no zero suppression!)
                                                                                         //              2 - flag is set to FALSE
    SetClusterParams(fXX,fYY,fCh);                                                       //              3 - size = 1
    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

  delete fitter;                                                                        //temporary solution to avoid the inteference with previous instances
  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<rawSize;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<rawSize;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...
   fNlocMax = 1;
   fSt=kNoLoc;
   SetClusterParams(fXX,fYY,fCh);                                                          //need to fill the AliCluster3D part
   new ((*pCluLst)[iCluCnt++]) AliHMPIDCluster(*this);	                                   //add new unfolded cluster
   
   return fNlocMax;
 }

// case 2 -> loc max found. Check # of loc maxima 
 if ( fNlocMax >= kMaxLocMax)  { 
 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(fNlocMax>1)FindClusterSize(i,pSigmaCut);                                         //find clustersize for deconvoluted clusters
                                                                                         //after this call, fSi temporarly is the calculated size. Later is set again 
                                                                                         //to its original value
     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(fXX,fYY,fCh);                                                      //need to fill the AliCluster3D part
     new ((*pCluLst)[iCluCnt++]) AliHMPIDCluster(*this);	                         //add new unfolded cluster
     if(fNlocMax>1)SetSize(rawSize);                                                     //Original raw size is set again to its proper value
   }
 }

 return fNlocMax;
 
}//Solve()
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
void AliHMPIDCluster::FindClusterSize(Int_t i,Int_t *pSigmaCut)
{

//Estimate of the clustersize for a deconvoluted cluster
  Int_t size = 0;
  for(Int_t iDig=0;iDig<Size();iDig++) {                                               //digits loop
    AliHMPIDDigit *pDig = Dig(iDig);                                                   //take digit
    Int_t iCh = pDig->Ch();
    Double_t qPad = Q()*pDig->IntMathieson(X(),Y());                                   //pad charge
    AliDebug(1,Form("Chamber %i X %i Y %i SigmaCut %i pad %i qpadMath %8.2f qPadRaw %8.2f Qtotal %8.2f cluster n.%i",iCh,pDig->PadChX(),pDig->PadChY(),
        pSigmaCut[iCh],iDig,qPad,pDig->Q(),QRaw(),i));
    if(qPad>pSigmaCut[iCh]) size++;
   }
  AliDebug(1,Form(" Calculated size %i",size));
  if(size>0) SetSize(size);                                                            //in case of size == 0, original raw clustersize used 
}
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
9ecfba5d	16	#include "AliHMPIDCluster.h" //class header
76fd1a96	17	#include <TVirtualFitter.h> //Solve()
d3da6dc4	18	#include <TMinuit.h> //Solve()
d3da6dc4	19	#include <TClonesArray.h> //Solve()
d1bf51e1	20	#include <TMarker.h> //Draw()
27311693	21
9ecfba5d	22	#include "AliLog.h" //FindCusterSize()
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
c1afc092	42	fParam = AliHMPIDParam::Instance();
55a829a5	43
55a829a5	44	//transformation from the pad cs to local
b38ac33a	45	xL -= 0.5*fParam->SizeAllX(); //size of all pads with dead zones included
b38ac33a	46	yL -= 0.5*fParam->SizeAllY();
55a829a5	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	}
59df3af2	73	//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	74	AliHMPIDCluster::~AliHMPIDCluster()
	75	{
	76	if(fDigs) delete fDigs; fDigs=0;
	77	if(fParam) delete fParam; fParam=0;
	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
e4a3eae8	86	if(fDigs==0) return; //no digits in this cluster
aa8db9d6	87	fXX=fYY=fQRaw=0; //init summable parameters
aa8db9d6	88	fCh = -1; //init chamber
e4a3eae8	89	Int_t maxQpad=-1,maxQ=-1; //to calculate the pad with the highest charge
aa85549f	90	AliHMPIDDigit *pDig=0x0;
e4a3eae8	91	for(Int_t iDig=0;iDig<fDigs->GetEntriesFast();iDig++){ //digits loop
e4a3eae8	92	pDig=(AliHMPIDDigit*)fDigs->At(iDig); //get pointer to next digit
aa8db9d6	93	if(!pDig) continue; //protection
e4a3eae8	94	if(pDig->PadPcX() > maxPadX) maxPadX = pDig->PadPcX(); // find the minimum box that contain the cluster MaxX
	95	if(pDig->PadPcY() > maxPadY) maxPadY = pDig->PadPcY(); // MaxY
	96	if(pDig->PadPcX() < minPadX) minPadX = pDig->PadPcX(); // MinX
	97	if(pDig->PadPcY() < minPadY) minPadY = pDig->PadPcY(); // MinY
	98
	99	Float_t q=pDig->Q(); //get QDC
894758d4	100	fXX += pDig->LorsX()q;fYY +=pDig->LorsY()q; //add digit center weighted by QDC
e4a3eae8	101	fQRaw+=q; //increment total charge
e4a3eae8	102	if(q>maxQ) {maxQpad = pDig->Pad();maxQ=(Int_t)q;} // to find pad with highest charge
aa8db9d6	103	fCh=pDig->Ch(); //initialize chamber number
d3da6dc4	104	}//digits loop
e4a3eae8	105
	106	fBox=(maxPadX-minPadX+1)100+maxPadY-minPadY+1; // dimension of the box: format Xdim100+Ydim
	107
3a94c321	108	if ( fQRaw != 0 ) {fXX/=fQRaw;fYY/=fQRaw;} //final center of gravity
27311693	109
69ed32de	110	if(fDigs->GetEntriesFast()>1&&fgDoCorrSin)CorrSin(); //correct it by sinoid
e4a3eae8	111
c5c19d6a	112	fQ = fQRaw; // Before starting fit procedure, Q and QRaw must be equal
c5c19d6a	113	fMaxQpad = maxQpad; fMaxQ=maxQ; //store max charge pad to the field
	114	fChi2=0; // no Chi2 to find
	115	fNlocMax=0; // proper status from this method
d3da6dc4	116	fSt=kCoG;
55a829a5	117
59df3af2	118	SetClusterParams(fXX,fYY,fCh); //need to fill the AliCluster3D part
55a829a5	119
d3da6dc4	120	}//CoG()
	121	//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	122	void AliHMPIDCluster::CorrSin()
	123	{
	124	// Correction of cluster x position due to sinoid, see HMPID TDR page 30
	125	// Arguments: none
	126	// Returns: none
1d4857c5	127	Int_t pc,px,py;
b38ac33a	128	fParam->Lors2Pad(fXX,fYY,pc,px,py); //tmp digit to get it center
b38ac33a	129	Float_t x=fXX-fParam->LorsX(pc,px); //diff between cluster x and center of the pad contaning this cluster
894758d4	130	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	131	}
d3da6dc4	132	//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
d1bf51e1	133	void AliHMPIDCluster::Draw(Option_t*)
d1bf51e1	134	{
a1d55ff3	135	TMarker *pMark=new TMarker(X(),Y(),5); pMark->SetUniqueID(fSt);pMark->SetMarkerColor(kBlue); pMark->Draw();
d1bf51e1	136	}
d1bf51e1	137	//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
557ca324	138	void AliHMPIDCluster::FitFunc(Int_t &iNpars, Double_t* deriv, Double_t &chi2, Double_t *par, Int_t iflag)
d3da6dc4	139	{
	140	// Cluster fit function
	141	// par[0]=x par[1]=y par[2]=q for the first Mathieson shape
	142	// par[3]=x par[4]=y par[5]=q for the second Mathieson shape and so on up to iNpars/3 Mathieson shapes
76fd1a96	143	// 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	144	// Then the chi2 is calculated as the sum of this value squared for all pad in the cluster.
	145	// Arguments: iNpars - number of parameters which is number of local maxima of cluster * 3
	146	// chi2 - function result to be minimised
	147	// par - parameters array of size iNpars
	148	// Returns: none
76fd1a96	149
	150	AliHMPIDCluster pClu=(AliHMPIDCluster)TVirtualFitter::GetFitter()->GetObjectFit();
	151
	152	Int_t nPads = pClu->Size();
76fd1a96	153
d3da6dc4	154	chi2 = 0;
76fd1a96	155
	156	Int_t iNshape = iNpars/3;
	157
76fd1a96	158	for(Int_t i=0;i<nPads;i++){ //loop on all pads of the cluster
	159	Double_t dQpadMath = 0;
	160	for(Int_t j=0;j<iNshape;j++){ //Mathiesons loop as all of them may contribute to this pad
	161	Double_t fracMathi = pClu->Dig(i)->IntMathieson(par[3j],par[3j+1]);
	162	dQpadMath+=par[3j+2]fracMathi; // par[3j+2] is charge par[3j] is x par[3*j+1] is y of current Mathieson
76fd1a96	163	}
	164	if(dQpadMath>0 && pClu->Dig(i)->Q()>0) {
	165	chi2 +=TMath::Power((pClu->Dig(i)->Q()-dQpadMath),2)/pClu->Dig(i)->Q(); //chi2 function to be minimized
	166	}
	167	}
557ca324	168	//---calculate gradients...
	169	if(iflag==2) {
	170	Double_t **derivPart;
	171
	172	derivPart = new Double_t*[iNpars];
	173
	174	for(Int_t j=0;j<iNpars;j++){
	175	deriv[j] = 0;
	176	derivPart[j] = new Double_t[nPads];
	177	for(Int_t i=0;i<nPads;i++){
	178	derivPart[j][i] = 0;
	179	}
	180	}
	181
	182	for(Int_t i=0;i<nPads;i++){ //loop on all pads of the cluster
	183	for(Int_t j=0;j<iNshape;j++){ //Mathiesons loop as all of them may contribute to this pad
	184	Double_t fracMathi = pClu->Dig(i)->IntMathieson(par[3j],par[3j+1]);
c770ceb9	185	derivPart[3j ][i] += par[3j+2](pClu->Dig(i)->MathiesonX(par[3j]-pClu->Dig(i)->LorsX()-0.5*AliHMPIDParam::SizePadX())-
	186	pClu->Dig(i)->MathiesonX(par[3j]-pClu->Dig(i)->LorsX()+0.5AliHMPIDParam::SizePadX()))*
	187	pClu->Dig(i)->IntPartMathiY(par[3*j+1]);
	188	derivPart[3j+1][i] += par[3j+2](pClu->Dig(i)->MathiesonY(par[3j+1]-pClu->Dig(i)->LorsY()-0.5*AliHMPIDParam::SizePadY())-
	189	pClu->Dig(i)->MathiesonY(par[3j+1]-pClu->Dig(i)->LorsY()+0.5AliHMPIDParam::SizePadY()))*
	190	pClu->Dig(i)->IntPartMathiX(par[3*j]);
557ca324	191	derivPart[3*j+2][i] += fracMathi;
	192	}
	193	}
	194	//loop on all pads of the cluster
	195	for(Int_t i=0;i<nPads;i++){ //loop on all pads of the cluster
	196	Double_t dQpadMath = 0; //pad charge collector
	197	for(Int_t j=0;j<iNshape;j++){ //Mathiesons loop as all of them may contribute to this pad
	198	Double_t fracMathi = pClu->Dig(i)->IntMathieson(par[3j],par[3j+1]);
	199	dQpadMath+=par[3j+2]fracMathi;
	200	if(dQpadMath>0 && pClu->Dig(i)->Q()>0) {
	201	deriv[3j] += 2/pClu->Dig(i)->Q()(pClu->Dig(i)->Q()-dQpadMath)derivPart[3j ][i];
	202	deriv[3j+1] += 2/pClu->Dig(i)->Q()(pClu->Dig(i)->Q()-dQpadMath)derivPart[3j+1][i];
	203	deriv[3j+2] += 2/pClu->Dig(i)->Q()(pClu->Dig(i)->Q()-dQpadMath)derivPart[3j+2][i];
	204	}
76fd1a96	205	}
d3da6dc4	206	}
557ca324	207	//delete array...
557ca324	208	for(Int_t i=0;i<iNpars;i++) delete [] derivPart[i]; delete [] derivPart;
76fd1a96	209	}
557ca324	210	//---gradient calculations ended
54104a7c	211
	212	// fit ended. Final calculations
	213
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;
03768ab2	233	case kBig : status="Big Clu(>100) " ;break;
d1bf51e1	234
d1bf51e1	235	default: status="??????" ;break;
d3da6dc4	236	}
e4a3eae8	237	Double_t ratio=0;
e4a3eae8	238	if(Q()>0&&QRaw()>0) ratio = Q()/QRaw()*100;
86aed250	239	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	240	opt,Ch(),X(),Y(),Q(),QRaw(),ratio,Size(),fBox,fNlocMax,fChi2,status);
d1bf51e1	241	if(fDigs) fDigs->Print();
d3da6dc4	242	}//Print()
d3da6dc4	243	//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
54104a7c	244	Int_t AliHMPIDCluster::Solve(TClonesArray pCluLst,Int_t pSigmaCut, Bool_t isTryUnfold)
d3da6dc4	245	{
	246	//This methode is invoked when the cluster is formed to solve it. Solve the cluster means to try to unfold the cluster
	247	//into the local maxima number of clusters. This methode is invoked by AliHMPIDRconstructor::Dig2Clu() on cluster by cluster basis.
	248	//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,
	249	//status is preset to kFormed in AddDigit(), chamber-sector info is preseted to actual values in AddDigit()
	250	//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
	251	//Arguments: pCluLst - cluster list pointer where to add new cluster(s)
	252	// isTryUnfold - flag to switch on/off unfolding
	253	// Returns: number of local maxima of original cluster
76fd1a96	254	const Int_t kMaxLocMax=6; //max allowed number of loc max for fitting
	255	//
	256	CoG(); //First calculate CoG for the given cluster
03768ab2	257
e4a3eae8	258	Int_t iCluCnt=pCluLst->GetEntriesFast(); //get current number of clusters already stored in the list by previous operations
03768ab2	259
	260	Int_t rawSize = Size(); //get current raw cluster size
	261
416b7266	262	if(rawSize>100) {
	263	fSt = kBig;
	264	} else if(isTryUnfold==kFALSE) {
	265	fSt = kNot;
	266	} else if(rawSize==1) {
	267	fSt = kSi1;
	268	}
	269
	270	if(rawSize>100 \|\| isTryUnfold==kFALSE \|\| rawSize==1) { //No deconv if: 1 - big cluster (also avoid no zero suppression!)
	271	// 2 - flag is set to FALSE
59df3af2	272	SetClusterParams(fXX,fYY,fCh); // 3 - size = 1
d1bf51e1	273	new ((pCluLst)[iCluCnt++]) AliHMPIDCluster(this); //add this raw cluster
d1bf51e1	274	return 1;
416b7266	275
d1bf51e1	276	}
76fd1a96	277
	278	//Phase 0. Initialise Fitter
	279	Double_t arglist[10];
	280	Int_t ierflg = 0;
	281	TVirtualFitter fitter = TVirtualFitter::Fitter(this,36); //initialize Fitter
09818f05	282
	283	delete fitter; //temporary solution to avoid the inteference with previous instances
	284	fitter = TVirtualFitter::Fitter(this,3*6); //initialize Fitter
76fd1a96	285
	286	arglist[0] = -1;
	287	ierflg = fitter->ExecuteCommand("SET PRI", arglist, 1); // no printout
	288	ierflg = fitter->ExecuteCommand("SET NOW", arglist, 0); //no warning messages
	289	arglist[0] = 1;
	290	ierflg = fitter->ExecuteCommand("SET GRA", arglist, 1); //force Fitter to use my gradient
	291
	292	fitter->SetFCN(AliHMPIDCluster::FitFunc);
	293
	294	// arglist[0] = 1;
	295	// ierflg = fitter->ExecuteCommand("SET ERR", arglist ,1);
	296
	297	// Set starting values and step sizes for parameters
	298
d1bf51e1	299	//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	300	fNlocMax=0;
e4a3eae8	301
54104a7c	302	for(Int_t iDig1=0;iDig1<rawSize;iDig1++) { //first digits loop
76fd1a96	303
e4a3eae8	304	AliHMPIDDigit *pDig1 = Dig(iDig1); //take next digit
c5c19d6a	305	Int_t iCnt = 0; //counts how many neighbouring pads has QDC more then current one
76fd1a96	306
54104a7c	307	for(Int_t iDig2=0;iDig2<rawSize;iDig2++) { //loop on all digits again
76fd1a96	308
d1bf51e1	309	if(iDig1==iDig2) continue; //the same digit, no need to compare
d3da6dc4	310	AliHMPIDDigit *pDig2 = Dig(iDig2); //take second digit to compare with the first one
da08475b	311	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	312	if(dist==1) //means dig2 is a neighbour of dig1
c5c19d6a	313	if(pDig2->Q()>=pDig1->Q()) iCnt++; //count number of pads with Q more then Q of current pad
76fd1a96	314
d3da6dc4	315	}//second digits loop
76fd1a96	316
c5c19d6a	317	if(iCnt==0&&fNlocMax<kMaxLocMax){ //this pad has Q more then any neighbour so it's local maximum
76fd1a96	318
c5c19d6a	319	Double_t xStart=pDig1->LorsX();Double_t yStart=pDig1->LorsY();
b38ac33a	320	Double_t xMin=xStart-fParam->SizePadX();
	321	Double_t xMax=xStart+fParam->SizePadX();
	322	Double_t yMin=yStart-fParam->SizePadY();
	323	Double_t yMax=yStart+fParam->SizePadY();
76fd1a96	324
	325	ierflg = fitter->SetParameter(3*fNlocMax ,Form("x%i",fNlocMax),xStart,0.1,xMin,xMax); // X,Y,Q initial values of the loc max pad
	326	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	327	ierflg = fitter->SetParameter(3*fNlocMax+2,Form("q%i",fNlocMax),pDig1->Q(),0.1,0,10000); // Q constrained to be positive
76fd1a96	328
d1bf51e1	329	fNlocMax++;
76fd1a96	330
d3da6dc4	331	}//if this pad is local maximum
d3da6dc4	332	}//first digits loop
d1bf51e1	333
d3da6dc4	334	//Phase 2. Fit loc max number of Mathiesons or add this current cluster to the list
e4a3eae8	335	// case 1 -> no loc max found
76fd1a96	336	if ( fNlocMax == 0) { // case of no local maxima found: pads with same charge...
d1bf51e1	337	fNlocMax = 1;
d1bf51e1	338	fSt=kNoLoc;
59df3af2	339	SetClusterParams(fXX,fYY,fCh); //need to fill the AliCluster3D part
75c94773	340	new ((pCluLst)[iCluCnt++]) AliHMPIDCluster(this); //add new unfolded cluster
	341
	342	return fNlocMax;
d1bf51e1	343	}
e4a3eae8	344
e4a3eae8	345	// case 2 -> loc max found. Check # of loc maxima
55a829a5	346	if ( fNlocMax >= kMaxLocMax) {
59df3af2	347	SetClusterParams(fXX,fYY,fCh); // if # of local maxima exceeds kMaxLocMax...
76fd1a96	348	fSt = kMax; new ((pCluLst)[iCluCnt++]) AliHMPIDCluster(this); //...add this raw cluster
	349	} else { //or resonable number of local maxima to fit and user requested it
	350	// Now ready for minimization step
	351	arglist[0] = 500; //number of steps and sigma on pads charges
	352	arglist[1] = 1.; //
	353
	354	ierflg = fitter->ExecuteCommand("SIMPLEX",arglist,2); //start fitting with Simplex
	355	if (!ierflg)
	356	fitter->ExecuteCommand("MIGRAD" ,arglist,2); //fitting improved by Migrad
	357	if(ierflg) {
c5c19d6a	358	Double_t strategy=2;
76fd1a96	359	ierflg = fitter->ExecuteCommand("SET STR",&strategy,1); //change level of strategy
	360	if(!ierflg) {
	361	ierflg = fitter->ExecuteCommand("SIMPLEX",arglist,2); //start fitting with Simplex
	362	if (!ierflg)
	363	fitter->ExecuteCommand("MIGRAD" ,arglist,2); //fitting improved by Migrad
c5c19d6a	364	}
c5c19d6a	365	}
76fd1a96	366	if(ierflg) fSt=kAbn; //no convergence of the fit...
	367	Double_t dummy; char sName[80]; //vars to get results from Minuit
	368	Double_t edm, errdef;
	369	Int_t nvpar, nparx;
54104a7c	370
e4a3eae8	371	for(Int_t i=0;i<fNlocMax;i++){ //store the local maxima parameters
54104a7c	372	fitter->GetParameter(3*i ,sName, fXX, fErrX , dummy, dummy); // X
54104a7c	373	fitter->GetParameter(3*i+1 ,sName, fYY, fErrY , dummy, dummy); // Y
76fd1a96	374	fitter->GetParameter(3*i+2 ,sName, fQ, fErrQ , dummy, dummy); // Q
76fd1a96	375	fitter->GetStats(fChi2, edm, errdef, nvpar, nparx); //get fit infos
54104a7c	376
	377	if(fNlocMax>1)FindClusterSize(i,pSigmaCut); //find clustersize for deconvoluted clusters
	378	//after this call, fSi temporarly is the calculated size. Later is set again
	379	//to its original value
	380	if(fSt!=kAbn) {
	381	if(fNlocMax!=1)fSt=kUnf; // if unfolded
	382	if(fNlocMax==1&&fSt!=kNoLoc) fSt=kLo1; // if only 1 loc max
	383	if ( !IsInPc()) fSt = kEdg; // if Out of Pc
	384	if(fSt==kNoLoc) fNlocMax=0; // if with no loc max (pads with same charge..)
	385	}
59df3af2	386	SetClusterParams(fXX,fYY,fCh); //need to fill the AliCluster3D part
54104a7c	387	new ((pCluLst)[iCluCnt++]) AliHMPIDCluster(this); //add new unfolded cluster
54104a7c	388	if(fNlocMax>1)SetSize(rawSize); //Original raw size is set again to its proper value
e4a3eae8	389	}
e4a3eae8	390	}
d1bf51e1	391
e4a3eae8	392	return fNlocMax;
e4a3eae8	393
d3da6dc4	394	}//Solve()
d3da6dc4	395	//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
54104a7c	396	void AliHMPIDCluster::FindClusterSize(Int_t i,Int_t *pSigmaCut)
	397	{
	398
	399	//Estimate of the clustersize for a deconvoluted cluster
	400	Int_t size = 0;
	401	for(Int_t iDig=0;iDig<Size();iDig++) { //digits loop
	402	AliHMPIDDigit *pDig = Dig(iDig); //take digit
	403	Int_t iCh = pDig->Ch();
	404	Double_t qPad = Q()*pDig->IntMathieson(X(),Y()); //pad charge
	405	AliDebug(1,Form("Chamber %i X %i Y %i SigmaCut %i pad %i qpadMath %8.2f qPadRaw %8.2f Qtotal %8.2f cluster n.%i",iCh,pDig->PadChX(),pDig->PadChY(),
	406	pSigmaCut[iCh],iDig,qPad,pDig->Q(),QRaw(),i));
	407	if(qPad>pSigmaCut[iCh]) size++;
	408	}
	409	AliDebug(1,Form(" Calculated size %i",size));
9ecfba5d	410	if(size>0) SetSize(size); //in case of size == 0, original raw clustersize used
54104a7c	411	}