[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
  Int_t maxQpad=-1,maxQ=-1;                                 //to calculate the pad with the highest charge
  AliHMPIDDigit *pDig=0x0;
  for(Int_t iDig=0;iDig<fDigs->GetEntriesFast();iDig++){    //digits loop
    pDig=(AliHMPIDDigit*)fDigs->At(iDig);                   //get pointer to next digit

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

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

    derivPart = new Double_t*[iNpars];

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

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