for(Int_t i=0;i<iNpars;i++) delete [] derivPart[i]; delete [] derivPart;
}
//---gradient calculations ended
+
+// fit ended. Final calculations
+
}//FitFunction()
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
if(fDigs) fDigs->Print();
}//Print()
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-Int_t AliHMPIDCluster::Solve(TClonesArray *pCluLst,Bool_t isTryUnfold)
+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.
//Phase 1. Find number of local maxima. Strategy is to check if the current pad has QDC more then all neigbours. Also find the box contaning the cluster
fNlocMax=0;
- for(Int_t iDig1=0;iDig1<Size();iDig1++) { //first digits loop
+ Int_t rawSize = Size();
+ 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<Size();iDig2++) { //loop on all digits again
+ 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
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 ,sName, fXX, fErrX , dummy, dummy); // X
+ fitter->GetParameter(3*i+1 ,sName, fYY, fErrY , dummy, dummy); // Y
fitter->GetParameter(3*i+2 ,sName, fQ, fErrQ , dummy, dummy); // Q
fitter->GetStats(fChi2, edm, errdef, nvpar, nparx); //get fit infos
- if(fSt!=kAbn) {
- if(fNlocMax!=1)fSt=kUnf; // if unfolded
- if(fNlocMax==1&&fSt!=kNoLoc) fSt=kLo1; // if only 1 loc max
- if ( !IsInPc()) fSt = kEdg; // if Out of Pc
- if(fSt==kNoLoc) fNlocMax=0; // if with no loc max (pads with same charge..)
- }
- if(fParam->GetInstType()) SetClusterParams(fXX,fYY,fCh); //need to fill the AliCluster3D part
- new ((*pCluLst)[iCluCnt++]) AliHMPIDCluster(*this); //add new unfolded cluster
-
+
+ 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
}
}
}//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));
+ SetSize(size);
+}
void Print (Option_t *opt="" )const; //overloaded TObject::Print() to print cluster info
static void FitFunc(Int_t &iNpars, Double_t* deriv, Double_t &chi2, Double_t *par, Int_t iflag);//fit function to be used by MINUIT
//private part
- Int_t Box ( )const{return fBox; } //Dimension of the cluster
+ Int_t Box ( )const{return fBox; } //Dimension of the cluster
void CoG ( ); //calculates center of gravity
void CorrSin ( ); //sinoidal correction
Int_t Ch ( )const{return fCh; } //chamber number
inline void Reset ( ); //cleans the cluster
void SetClusterParams(Double_t xL,Double_t yL,Int_t iCh ); //Set AliCluster3D part
Int_t Size ( )const{return fSi; } //returns number of pads in formed cluster
- Int_t Solve (TClonesArray *pCluLst,Bool_t isUnfold ); //solve cluster: MINUIT fit or CoG
+ Int_t Solve (TClonesArray *pCluLst,Int_t *pSigmaCut, Bool_t isUnfold); //solve cluster: MINUIT fit or CoG
Int_t Status ( ) const{return fSt;} //Status of cluster
Double_t QRaw ( )const{return fQRaw; } //raw cluster charge in QDC channels
Double_t Q ( )const{return fQ; } //given cluster charge in QDC channels
Double_t Qe ( )const{return fErrQ; } //Error in cluster charge in QDC channels
- Double_t X ( )const{return fXX; } //cluster x position in LRS
+ Double_t X ( )const{return fXX; } //cluster x position in LRS
Double_t Xe ( )const{return fErrX; } //cluster charge in QDC channels
- Double_t Y ( )const{return fYY; } //cluster y position in LRS
+ Double_t Y ( )const{return fYY; } //cluster y position in LRS
Double_t Ye ( )const{return fErrY; } //cluster charge in QDC channels
Double_t Chi2 ( )const{return fChi2; } //chi2 of the fit
- void DoCorrSin(Bool_t doCorrSin ){fgDoCorrSin=doCorrSin;} // Set sinoidal correction
- void SetX (Double_t x ){fXX=x;} // Setter
- void SetY (Double_t y ){fYY=y;} // Setter
+ void DoCorrSin(Bool_t doCorrSin ){fgDoCorrSin=doCorrSin;} // Set sinoidal correction
+ void SetX (Double_t x ){fXX=x;} // Setter
+ void SetY (Double_t y ){fYY=y;} // Setter
+ void SetSize (Int_t size ){fSi=size;} // Setter
+ void FindClusterSize(Int_t i,Int_t *pSigmaCut); //Find the clusterSize of deconvoluted clusters
protected:
Int_t fCh; //chamber number
}
}//AliHMPIDReconstructor
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-void AliHMPIDReconstructor::Dig2Clu(TObjArray *pDigAll,TObjArray *pCluAll,Bool_t isTryUnfold)
+void AliHMPIDReconstructor::Dig2Clu(TObjArray *pDigAll,TObjArray *pCluAll,Int_t *pUserCut,Bool_t isTryUnfold)
{
// Finds all clusters for a given digits list provided not empty. Currently digits list is a list of all digits for a single chamber.
// Puts all found clusters in separate lists, one per clusters.
AliHMPIDDigit *pDig=(AliHMPIDDigit*)pDigCur->At(iDig); //take current digit
if(!(pDig=UseDig(pDig->PadChX(),pDig->PadChY(),pDigCur,&padMap))) continue; //this digit is already taken in FormClu(), go after next digit
FormClu(&clu,pDig,pDigCur,&padMap); //form cluster starting from this digit by recursion
-
- clu.Solve(pCluCur,isTryUnfold); //solve this cluster and add all unfolded clusters to provided list
+ clu.Solve(pCluCur,pUserCut,isTryUnfold); //solve this cluster and add all unfolded clusters to provided list
clu.Reset(); //empty current cluster
}//digits loop for current chamber
}//chambers loop
pDigTree->SetBranchAddress(Form("HMPID%d",iCh),&((*fDig)[iCh]));
}
pDigTree->GetEntry(0);
- Dig2Clu(fDig,fClu); //cluster finder
+ Dig2Clu(fDig,fClu,fUserCut); //cluster finder
pCluTree->Fill(); //fill tree for current event
for(Int_t iCh=AliHMPIDParam::kMinCh;iCh<=AliHMPIDParam::kMaxCh;iCh++){
using AliReconstructor::Reconstruct; //
//private part
- static void Dig2Clu (TObjArray *pDigLst,TObjArray *pCluLst,Bool_t isUnfold=kFALSE );//digits->clusters
+ static void Dig2Clu (TObjArray *pDigLst,TObjArray *pCluLst,Int_t *pUserCut,Bool_t isUnfold=kFALSE );//digits->clusters
static void FormClu (AliHMPIDCluster *pClu,AliHMPIDDigit *pDig,TClonesArray *pDigLst,TMatrixF *pPadMap);//cluster formation recursive algorithm
static inline AliHMPIDDigit* UseDig (Int_t padX,Int_t padY, TClonesArray *pDigLst,TMatrixF *pDigMap);//use this pad's digit to form a cluster
inline Bool_t IsDigSurvive(AliHMPIDDigit *pDig )const;//check for sigma cut
-
void SetRecoParam(AliHMPIDRecoParam *recopar){ fgkRecoParam = recopar;}
static const AliHMPIDRecoParam* GetRecoParam(){ return fgkRecoParam;}
fprintf(fp," gSystem->Exec(\"touch ZZZ______finished_______SSS\");\n");
fprintf(fp," gSystem->Exec(\"aliroot rec.C &\");\n}\n");
fclose(fp);
+// rec section
char *sBatchName="rec";
FILE *fp=fopen(Form("%s.C",sBatchName),"w"); if(!fp){Info("CreateRec","Cannot open output file: %s.C",sBatchName);return;}
fprintf(fp,"void %s()\n{\n",sBatchName);
- fprintf(fp," gSystem->Exec(\"rm -rf RRR* \"); //remove garbage\n");
+ fprintf(fp," gSystem->Exec(\"rm -rf *RRR \"); //remove garbage\n");
if(fRecB->GetState()){
fprintf(fp," AliReconstruction *pRec=new AliReconstruction;\n");