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

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