X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=HMPID%2FAliHMPIDCluster.cxx;h=046a6fd74d2267b1cc2cd5f2dd0ccd39edfff292;hb=93d28366943be42914694abd58d7649a9b403ead;hp=a928221ce0a1c0b3b0745fce91fc390464e258ad;hpb=1006986dcae248d3d66ae37b1f094e40f223d772;p=u%2Fmrichter%2FAliRoot.git diff --git a/HMPID/AliHMPIDCluster.cxx b/HMPID/AliHMPIDCluster.cxx index a928221ce0a..046a6fd74d2 100644 --- a/HMPID/AliHMPIDCluster.cxx +++ b/HMPID/AliHMPIDCluster.cxx @@ -17,6 +17,9 @@ #include //Solve() #include //Solve() #include //Draw() + +Bool_t AliHMPIDCluster::fgDoCorrSin=kTRUE; + ClassImp(AliHMPIDCluster) //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ void AliHMPIDCluster::CoG() @@ -46,15 +49,14 @@ void AliHMPIDCluster::CoG() fBox=(maxPadX-minPadX+1)*100+maxPadY-minPadY+1; // dimension of the box: format Xdim*100+Ydim if ( fQRaw != 0 ) fX/=fQRaw;fY/=fQRaw; //final center of gravity + + if(fDigs->GetEntriesFast()>1&&fgDoCorrSin)CorrSin(); //correct it by sinoid - - 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; // no maxima to find + 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; }//CoG() //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ @@ -64,8 +66,8 @@ void AliHMPIDCluster::CorrSin() // Arguments: none // Returns: none Int_t pc,px,py; - AliHMPIDDigit::Lors2Pad(fX,fY,pc,px,py); //tmp digit to get it center - Float_t x=fX-AliHMPIDDigit::LorsX(pc,px); //diff between cluster x and center of the pad contaning this cluster + AliHMPIDParam::Lors2Pad(fX,fY,pc,px,py); //tmp digit to get it center + Float_t x=fX-AliHMPIDParam::LorsX(pc,px); //diff between cluster x and center of the pad contaning this cluster fX+=3.31267e-2*TMath::Sin(2*TMath::Pi()/0.8*x)-2.66575e-3*TMath::Sin(4*TMath::Pi()/0.8*x)+2.80553e-3*TMath::Sin(6*TMath::Pi()/0.8*x)+0.0070; } //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ @@ -74,7 +76,7 @@ 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 *, Double_t &chi2, Double_t *par, Int_t ) +void AliHMPIDCluster::FitFunc(Int_t &iNpars, Double_t* /*deriv*/, Double_t &chi2, Double_t *par, Int_t /* */) { // Cluster fit function // par[0]=x par[1]=y par[2]=q for the first Mathieson shape @@ -92,9 +94,11 @@ void AliHMPIDCluster::FitFunc(Int_t &iNpars, Double_t *, Double_t &chi2, Double_ for(Int_t i=0;iSize();i++){ //loop on all pads of the cluster Double_t dQpadMath = 0; //pad charge collector for(Int_t j=0;jDig(i)->Mathieson(par[3*j],par[3*j+1]); // par[3*j+2] is charge par[3*j] is x par[3*j+1] is y of current Mathieson + dQpadMath+=par[3*j+2]*pClu->Dig(i)->IntMathieson(par[3*j],par[3*j+1]); // par[3*j+2] is charge par[3*j] is x par[3*j+1] is y of current Mathieson } - if(dQpadMath>0)chi2 +=TMath::Power((pClu->Dig(i)->Q()-dQpadMath),2)/dQpadMath; // +// if(dQpadMath>0)chi2 +=TMath::Power((pClu->Dig(i)->Q()-dQpadMath),2)/dQpadMath; // + if(dQpadMath>0 && pClu->Dig(i)->Q()) + chi2 +=TMath::Power((pClu->Dig(i)->Q()-dQpadMath),2)/pClu->Dig(i)->Q(); // } //loop on all pads of the cluster }//FitFunction() //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ @@ -113,6 +117,7 @@ void AliHMPIDCluster::Print(Option_t* opt)const case kEdg : status="edge (fit)" ;break; case kSi1 : status="size 1 (cog)" ;break; case kNoLoc: status="no LocMax(fit)" ;break; + case kAbn : status="Abnormal fit " ;break; default: status="??????" ;break; } @@ -142,28 +147,35 @@ Int_t AliHMPIDCluster::Solve(TClonesArray *pCluLst,Bool_t isTryUnfold) } //Phase 0. Initialise TMinuit const Int_t kMaxLocMax=6; //max allowed number of loc max for fitting - TMinuit *pMinuit = new TMinuit(3*kMaxLocMax); //init MINUIT with this number of parameters (3 params per mathieson) - pMinuit->SetObjectFit((TObject*)this); pMinuit->SetFCN(AliHMPIDCluster::FitFunc); //set fit function - Double_t aArg=-1; Int_t iErrFlg; //tmp vars for TMinuit - pMinuit->mnexcm("SET PRI",&aArg,1,iErrFlg); //suspend all printout from TMinuit - pMinuit->mnexcm("SET NOW",&aArg,0,iErrFlg); //suspend all warning printout from TMinuit + if(!gMinuit) gMinuit = new TMinuit(100); //init MINUIT with this number of parameters (3 params per mathieson) + gMinuit->mncler(); // reset Minuit list of paramters + gMinuit->SetObjectFit((TObject*)this); gMinuit->SetFCN(AliHMPIDCluster::FitFunc); //set fit function + Double_t aArg=-1; //tmp vars for TMinuit + Int_t iErrFlg; + gMinuit->mnexcm("SET PRI",&aArg,1,iErrFlg); //suspend all printout from TMinuit + gMinuit->mnexcm("SET NOW",&aArg,0,iErrFlg); //suspend all warning printout from TMinuit //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;iDig1PadChX()-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()) iHowManyMoreCnt++; //count number of pads with Q more then Q of current pad + if(pDig2->Q()>=pDig1->Q()) iCnt++; //count number of pads with Q more then Q of current pad }//second digits loop - if(iHowManyMoreCnt==0&&fNlocMaxmnparm(3*fNlocMax ,Form("x%i",fNlocMax),pDig1->LorsX(),0.1,0,0,iErrFlg); // X,Y,Q initial values of the loc max pad w - pMinuit->mnparm(3*fNlocMax+1,Form("y%i",fNlocMax),pDig1->LorsY(),0.1,0,0,iErrFlg); // - pMinuit->mnparm(3*fNlocMax+2,Form("q%i",fNlocMax),pDig1->Q(),0.1,0,100000,iErrFlg);// constrained to be positive + if(iCnt==0&&fNlocMaxLorsX();Double_t yStart=pDig1->LorsY(); + Double_t xMin=xStart-AliHMPIDParam::SizePadX(); + Double_t xMax=xStart+AliHMPIDParam::SizePadX(); + Double_t yMin=yStart-AliHMPIDParam::SizePadY(); + Double_t yMax=yStart+AliHMPIDParam::SizePadY(); + gMinuit->mnparm(3*fNlocMax ,Form("x%i",fNlocMax),xStart,0.1,xMin,xMax,iErrFlg); // X,Y,Q initial values of the loc max pad + gMinuit->mnparm(3*fNlocMax+1,Form("y%i",fNlocMax),yStart,0.1,yMin,yMax,iErrFlg); // X, Y constrained to be near the loc max + gMinuit->mnparm(3*fNlocMax+2,Form("q%i",fNlocMax),pDig1->Q(),0.1,0,100000,iErrFlg);// Q constrained to be positive fNlocMax++; }//if this pad is local maximum }//first digits loop @@ -171,9 +183,9 @@ Int_t AliHMPIDCluster::Solve(TClonesArray *pCluLst,Bool_t isTryUnfold) //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... - pMinuit->mnparm(3*fNlocMax ,Form("x%i",fNlocMax),fX,0.1,0,0,iErrFlg); // Init values taken from CoG() -> fX,fY,fQRaw - pMinuit->mnparm(3*fNlocMax+1,Form("y%i",fNlocMax),fY,0.1,0,0,iErrFlg); // - pMinuit->mnparm(3*fNlocMax+2,Form("q%i",fNlocMax),fQRaw,0.1,0,100000,iErrFlg); // + gMinuit->mnparm(3*fNlocMax ,Form("x%i",fNlocMax),fX,0.1,0,0,iErrFlg); // Init values taken from CoG() -> fX,fY,fQRaw + gMinuit->mnparm(3*fNlocMax+1,Form("y%i",fNlocMax),fY,0.1,0,0,iErrFlg); // + gMinuit->mnparm(3*fNlocMax+2,Form("q%i",fNlocMax),fQRaw,0.1,0,100000,iErrFlg); // fNlocMax = 1; fSt=kNoLoc; } @@ -185,25 +197,34 @@ Int_t AliHMPIDCluster::Solve(TClonesArray *pCluLst,Bool_t isTryUnfold) } //or... else{ //...resonable number of local maxima to fit and user requested it Double_t arglist[10];arglist[0] = 10000;arglist[1] = 1.; //number of steps and sigma on pads charges - pMinuit->mnexcm("SIMPLEX" ,arglist,2,iErrFlg); //start fitting with Simplex - pMinuit->mnexcm("MIGRAD" ,arglist,2,iErrFlg); //fitting improved by Migrad - + gMinuit->mnexcm("SIMPLEX" ,arglist,2,iErrFlg); //start fitting with Simplex + gMinuit->mnexcm("MIGRAD" ,arglist,2,iErrFlg); //fitting improved by Migrad + if(iErrFlg) { + Double_t strategy=2; + gMinuit->mnexcm("SET STR",&strategy,1,iErrFlg); //change level of strategy + if(!iErrFlg) { + gMinuit->mnexcm("SIMPLEX" ,arglist,2,iErrFlg); + gMinuit->mnexcm("MIGRAD" ,arglist,2,iErrFlg); //fitting improved by Migrad +// Printf("Try to improve fit --> err %d",iErrFlg); + } + } + if(iErrFlg) fSt=kAbn; //no convergence of the fit... Double_t dummy; TString sName; //vars to get results from Minuit for(Int_t i=0;imnpout(3*i ,sName, fX, fErrX , dummy, dummy, iErrFlg); // X - pMinuit->mnpout(3*i+1 ,sName, fY, fErrY , dummy, dummy, iErrFlg); // Y - pMinuit->mnpout(3*i+2 ,sName, fQ, fErrQ , dummy, dummy, iErrFlg); // Q - pMinuit->mnstat(fChi2,dummy,dummy,iErrFlg,iErrFlg,iErrFlg); // Chi2 of the fit - - if(fNlocMax!=1)fSt=kUnf; // - if(fNlocMax==1&&fSt!=kNoLoc) fSt=kLo1; // - if ( !IsInPc()) fSt = kEdg; // - if(fSt==kNoLoc) fNlocMax=0; // + gMinuit->mnpout(3*i ,sName, fX, fErrX , dummy, dummy, iErrFlg); // X + gMinuit->mnpout(3*i+1 ,sName, fY, fErrY , dummy, dummy, iErrFlg); // Y + gMinuit->mnpout(3*i+2 ,sName, fQ, fErrQ , dummy, dummy, iErrFlg); // Q + gMinuit->mnstat(fChi2,dummy,dummy,iErrFlg,iErrFlg,iErrFlg); // Chi2 of the fit + 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..) + } new ((*pCluLst)[iCluCnt++]) AliHMPIDCluster(*this); //add new unfolded cluster } } - delete pMinuit; return fNlocMax; }//Solve()