X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=RICH%2FAliRICHClusterFinder.cxx;h=4e1a7d9a602adfbe51fbe36107372d726554d3de;hb=39853df5344e33a76cdf9d3d6f4cc7f0788a1a03;hp=d4766b721c57d970acf86a027c4d51dab44a83b2;hpb=237c933d68793192ded67c959825d4e5e881131e;p=u%2Fmrichter%2FAliRoot.git diff --git a/RICH/AliRICHClusterFinder.cxx b/RICH/AliRICHClusterFinder.cxx index d4766b721c5..4e1a7d9a602 100644 --- a/RICH/AliRICHClusterFinder.cxx +++ b/RICH/AliRICHClusterFinder.cxx @@ -13,1115 +13,323 @@ * provided "as is" without express or implied warranty. * **************************************************************************/ -/* - $Log$ - Revision 1.1 2000/04/19 13:01:48 morsch - A cluster finder and hit reconstruction class for RICH (adapted from MUON). - Cluster Finders for MUON and RICH should derive from the same class in the - future (JB, AM). - -*/ - #include "AliRICHClusterFinder.h" -#include "AliRun.h" -#include "AliRICH.h" -#include "AliRICHHit.h" -#include "AliRICHHitMapA1.h" -#include "AliRICHCerenkov.h" -#include "AliRICHPadHit.h" -#include "AliRICHDigit.h" -#include "AliRICHRawCluster.h" -#include "AliRICHRecHit.h" - -#include -#include -#include -#include -#include -#include -#include +#include "AliRICHMap.h" +#include +#include +#include +#include +#include +#include -//---------------------------------------------------------- -static AliRICHSegmentation* gSegmentation; -static AliRICHResponse* gResponse; -static Int_t gix[500]; -static Int_t giy[500]; -static Float_t gCharge[500]; -static Int_t gNbins; -static Int_t gFirst=kTRUE; -static TMinuit *gMyMinuit ; -void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag); -static Int_t gChargeTot; +void RICHMinMathieson(Int_t &npar, Double_t *gin, Double_t &chi2, Double_t *par, Int_t iflag); ClassImp(AliRICHClusterFinder) - -AliRICHClusterFinder::AliRICHClusterFinder -(AliRICHSegmentation *segmentation, AliRICHResponse *response, - TClonesArray *digits, Int_t chamber) +//__________________________________________________________________________________________________ +AliRICHClusterFinder::AliRICHClusterFinder(AliRICH *pRICH) +{//main ctor + fRICH = pRICH; + AliDebug(1,"main ctor Start."); + + fDigitMap = 0; + fRawCluster.Reset(); + fResolvedCluster.Reset(); + AliDebug(1,"main ctor Stop."); +}//main ctor +//__________________________________________________________________________________________________ +void AliRICHClusterFinder::Exec() { - -// Constructor for Cluster Finder object - - fSegmentation=segmentation; - fResponse=response; +//Main method of cluster finder. Loops on events and chambers, everything else is done in FindClusters() + AliDebug(1,"Exec Start."); - fDigits=digits; - fNdigits = fDigits->GetEntriesFast(); - fChamber=chamber; - fRawClusters=new TClonesArray("AliRICHRawCluster",10000); - fNRawClusters=0; - fCogCorr = 0; - SetNperMax(); - SetClusterSize(); - SetDeclusterFlag(); - fNPeaks=-1; -} - -AliRICHClusterFinder::AliRICHClusterFinder() -{ - -// Default constructor + R()->GetLoader() ->LoadDigits(); +// R()->GetLoader()->GetRunLoader()->LoadHeader(); + R()->GetLoader()->GetRunLoader()->LoadKinematics(); //header is already loaded - fSegmentation=0; - fResponse=0; + for(Int_t iEventN=0;iEventNGetEventsPerRun();iEventN++){//events loop + AliDebug(1,Form("Processing event %i...",iEventN)); + R()->GetLoader()->GetRunLoader()->GetEvent(iEventN); - fDigits=0; - fNdigits = 0; - fChamber=-1; - fRawClusters=new TClonesArray("AliRICHRawCluster",10000); - fNRawClusters=0; - fHitMap = 0; - fCogCorr = 0; - SetNperMax(); - SetClusterSize(); - SetDeclusterFlag(); - fNPeaks=-1; -} - -AliRICHClusterFinder::AliRICHClusterFinder(const AliRICHClusterFinder& ClusterFinder) -{ -// Copy Constructor -} - -AliRICHClusterFinder::~AliRICHClusterFinder() -{ - -// Destructor - -delete fRawClusters; -} - -void AliRICHClusterFinder::AddRawCluster(const AliRICHRawCluster c) -{ - // - // Add a raw cluster copy to the list - // - AliRICH *pRICH=(AliRICH*)gAlice->GetModule("RICH"); - pRICH->AddRawCluster(fChamber,c); - fNRawClusters++; -} - - - -void AliRICHClusterFinder::Decluster(AliRICHRawCluster *cluster) -{ - -// -// Decluster algorithm + R()->GetLoader()->MakeTree("R"); R()->MakeBranch("R"); + R()->ResetDigits(); R()->ResetClusters(); - Int_t mul = cluster->fMultiplicity; -// printf("Decluster - multiplicity %d \n",mul); - - if (mul == 1 || mul ==2) { -// -// Nothing special for 1- and 2-clusters - if (fNPeaks != 0) { - cluster->fNcluster[0]=fNPeaks; - cluster->fNcluster[1]=0; - } - AddRawCluster(*cluster); - fNPeaks++; - } else if (mul ==3) { -// -// 3-cluster, check topology -// printf("\n 3-cluster, check topology \n"); - if (fDeclusterFlag) { - if (Centered(cluster)) { - // ok, cluster is centered - } else { - // cluster is not centered, split into 2+1 - } - } else { - if (fNPeaks != 0) { - cluster->fNcluster[0]=fNPeaks; - cluster->fNcluster[1]=0; - } - AddRawCluster(*cluster); - fNPeaks++; - } - } else { -// -// 4-and more-pad clusters -// - if (mul <= fClusterSize) { - if (fDeclusterFlag) { - SplitByLocalMaxima(cluster); - } else { - if (fNPeaks != 0) { - cluster->fNcluster[0]=fNPeaks; - cluster->fNcluster[1]=0; - } - AddRawCluster(*cluster); - fNPeaks++; - } - } - } // multiplicity -} - - -Bool_t AliRICHClusterFinder::Centered(AliRICHRawCluster *cluster) + R()->GetLoader()->TreeD()->GetEntry(0); + for(Int_t iChamber=1;iChamber<=kNchambers;iChamber++){//chambers loop + FindClusters(iChamber); + }//chambers loop + R()->GetLoader()->TreeR()->Fill(); R()->GetLoader()->WriteRecPoints("OVERWRITE");//write out clusters for current event + }//events loop + + R()->ResetDigits();//reset and unload everything + R()->ResetClusters(); + R()->GetLoader() ->UnloadDigits(); + R()->GetLoader() ->UnloadRecPoints(); +// R()->GetLoader()->GetRunLoader()->UnloadHeader(); + R()->GetLoader()->GetRunLoader()->UnloadKinematics(); + + AliDebug(1,"Stop."); +}//Exec() +//__________________________________________________________________________________________________ +void AliRICHClusterFinder::FindClusters(Int_t iChamber) { - -// Is the cluster centered? - - AliRICHDigit* dig; - dig= (AliRICHDigit*)fDigits->UncheckedAt(cluster->fIndexMap[0]); - Int_t ix=dig->fPadX; - Int_t iy=dig->fPadY; - Int_t nn; - Int_t x[kMaxNeighbours], y[kMaxNeighbours], xN[kMaxNeighbours], yN[kMaxNeighbours]; - - fSegmentation->Neighbours(ix,iy,&nn,x,y); - Int_t nd=0; - for (Int_t i=0; iTestHit(x[i],y[i]) == kUsed) { - yN[nd]=x[i]; - yN[nd]=y[i]; - nd++; - - printf("Getting: %d %d %d\n",i,x[i],y[i]); - } - } - if (nd==2) { -// -// cluster is centered ! - if (fNPeaks != 0) { - cluster->fNcluster[0]=fNPeaks; - cluster->fNcluster[1]=0; - } - cluster->fCtype=0; - AddRawCluster(*cluster); - fNPeaks++; - return kTRUE; - } else if (nd ==1) { -// -// Highest signal on an edge, split cluster into 2+1 -// -// who is the neighbour ? - - printf("Calling GetIndex with x:%d y:%d\n",xN[0], yN[0]); - - Int_t nind=fHitMap->GetHitIndex(xN[0], yN[0]); - Int_t i1= (nind==cluster->fIndexMap[1]) ? 1:2; - Int_t i2= (nind==cluster->fIndexMap[1]) ? 2:1; -// -// 2-cluster - AliRICHRawCluster cnew; - if (fNPeaks == 0) { - cnew.fNcluster[0]=-1; - cnew.fNcluster[1]=fNRawClusters; - } else { - cnew.fNcluster[0]=fNPeaks; - cnew.fNcluster[1]=0; - } - cnew.fMultiplicity=2; - cnew.fIndexMap[0]=cluster->fIndexMap[0]; - cnew.fIndexMap[1]=cluster->fIndexMap[i1]; - FillCluster(&cnew); - cnew.fClusterType=cnew.PhysicsContribution(); - AddRawCluster(cnew); - fNPeaks++; -// -// 1-cluster - cluster->fMultiplicity=1; - cluster->fIndexMap[0]=cluster->fIndexMap[i2]; - cluster->fIndexMap[1]=0; - cluster->fIndexMap[2]=0; - FillCluster(cluster); - if (fNPeaks != 0) { - cluster->fNcluster[0]=fNPeaks; - cluster->fNcluster[1]=0; - } - cluster->fClusterType=cluster->PhysicsContribution(); - AddRawCluster(*cluster); - fNPeaks++; - return kFALSE; - } else { - printf("\n Completely screwed up %d !! \n",nd); +//Loops on digits for a given chamber, forms raw clusters, then tries to resolve them if requested + Int_t iNdigits=R()->Digits(iChamber)->GetEntriesFast(); + AliDebug(1,Form("Start for chamber %i with %i digits.",iChamber,iNdigits)); + + if(iNdigits==0)return;//no digits for a given chamber, nothing to do + + fDigitMap=new AliRICHMap(R()->Digits(iChamber));//create digit map for the given chamber + + for(Int_t iDigN=0;iDigNDigits(iChamber)->At(iDigN); + Int_t i=dig->X(); Int_t j=dig->Y(); + if(fDigitMap->TestHit(i,j)==kUsed) continue;//this digit is already taken, go after next digit - } - - return kFALSE; -} -void AliRICHClusterFinder::SplitByLocalMaxima(AliRICHRawCluster *c) -{ - -// -// Split the cluster according to the number of maxima inside - - - AliRICHDigit* dig[100], *digt; - Int_t ix[100], iy[100], q[100]; - Float_t x[100], y[100]; - Int_t i; // loops over digits - Int_t j; // loops over local maxima - // Float_t xPeak[2]; - // Float_t yPeak[2]; - // Int_t threshold=500; - Int_t mul=c->fMultiplicity; -// -// dump digit information into arrays -// - for (i=0; iUncheckedAt(c->fIndexMap[i]); - ix[i]= dig[i]->fPadX; - iy[i]= dig[i]->fPadY; - q[i] = dig[i]->fSignal; - fSegmentation->GetPadCxy(ix[i], iy[i], x[i], y[i]); - } -// -// Find local maxima -// - Bool_t isLocal[100]; - Int_t nLocal=0; - Int_t associatePeak[100]; - Int_t indLocal[100]; - Int_t nn; - Int_t xNei[kMaxNeighbours], yNei[kMaxNeighbours]; - for (i=0; iNeighbours(ix[i], iy[i], &nn, xNei, yNei); - isLocal[i]=kTRUE; - for (j=0; jTestHit(xNei[j], yNei[j])==kEmpty) continue; - digt=(AliRICHDigit*) fHitMap->GetHit(xNei[j], yNei[j]); - if (digt->fSignal > q[i]) { - isLocal[i]=kFALSE; - break; -// -// handle special case of neighbouring pads with equal signal - } else if (digt->fSignal == q[i]) { - if (nLocal >0) { - for (Int_t k=0; k5) { - Int_t nnew=0; - for (i=0; iSetFCN(fcn); - gMyMinuit->mninit(5,10,7); - Double_t arglist[20]; - Int_t ierflag=0; - arglist[0]=1; -// gMyMinuit->mnexcm("SET ERR",arglist,1,ierflag); -// Set starting values - static Double_t vstart[5]; - vstart[0]=x[indLocal[0]]; - vstart[1]=y[indLocal[0]]; - vstart[2]=x[indLocal[1]]; - vstart[3]=y[indLocal[1]]; - vstart[4]=Float_t(q[indLocal[0]])/ - Float_t(q[indLocal[0]]+q[indLocal[1]]); -// lower and upper limits - static Double_t lower[5], upper[5]; - Int_t isec=fSegmentation->Sector(ix[indLocal[0]], iy[indLocal[0]]); - lower[0]=vstart[0]-fSegmentation->Dpx(isec)/2; - lower[1]=vstart[1]-fSegmentation->Dpy(isec)/2; -// lower[1]=vstart[1]; - - upper[0]=lower[0]+fSegmentation->Dpx(isec); - upper[1]=lower[1]+fSegmentation->Dpy(isec); -// upper[1]=vstart[1]; - - isec=fSegmentation->Sector(ix[indLocal[1]], iy[indLocal[1]]); - lower[2]=vstart[2]-fSegmentation->Dpx(isec)/2; - lower[3]=vstart[3]-fSegmentation->Dpy(isec)/2; -// lower[3]=vstart[3]; - - upper[2]=lower[2]+fSegmentation->Dpx(isec); - upper[3]=lower[3]+fSegmentation->Dpy(isec); -// upper[3]=vstart[3]; - - lower[4]=0.; - upper[4]=1.; -// step sizes - static Double_t step[5]={0.005, 0.03, 0.005, 0.03, 0.01}; - - gMyMinuit->mnparm(0,"x1",vstart[0],step[0],lower[0],upper[0],ierflag); - gMyMinuit->mnparm(1,"y1",vstart[1],step[1],lower[1],upper[1],ierflag); - gMyMinuit->mnparm(2,"x2",vstart[2],step[2],lower[2],upper[2],ierflag); - gMyMinuit->mnparm(3,"y2",vstart[3],step[3],lower[3],upper[3],ierflag); - gMyMinuit->mnparm(4,"a0",vstart[4],step[4],lower[4],upper[4],ierflag); -// ready for minimisation - gMyMinuit->SetPrintLevel(-1); - gMyMinuit->mnexcm("SET OUT", arglist, 0, ierflag); - arglist[0]= -1; - arglist[1]= 0; - - gMyMinuit->mnexcm("SET NOGR", arglist, 0, ierflag); - gMyMinuit->mnexcm("SCAN", arglist, 0, ierflag); - gMyMinuit->mnexcm("EXIT" , arglist, 0, ierflag); -// Print results -// Double_t amin,edm,errdef; -// Int_t nvpar,nparx,icstat; -// gMyMinuit->mnstat(amin,edm,errdef,nvpar,nparx,icstat); -// gMyMinuit->mnprin(3,amin); -// Get fitted parameters - - Double_t xrec[2], yrec[2], qfrac; - TString chname; - Double_t epxz, b1, b2; - Int_t ierflg; - gMyMinuit->mnpout(0, chname, xrec[0], epxz, b1, b2, ierflg); - gMyMinuit->mnpout(1, chname, yrec[0], epxz, b1, b2, ierflg); - gMyMinuit->mnpout(2, chname, xrec[1], epxz, b1, b2, ierflg); - gMyMinuit->mnpout(3, chname, yrec[1], epxz, b1, b2, ierflg); - gMyMinuit->mnpout(4, chname, qfrac, epxz, b1, b2, ierflg); - //printf("\n %f %f %f %f %f\n", xrec[0], yrec[0], xrec[1], yrec[1],qfrac); -// delete gMyMinuit; - - - // - // One cluster for each maximum - // - for (j=0; j<2; j++) { - AliRICHRawCluster cnew; - if (fNPeaks == 0) { - cnew.fNcluster[0]=-1; - cnew.fNcluster[1]=fNRawClusters; - } else { - cnew.fNcluster[0]=fNPeaks; - cnew.fNcluster[1]=0; - } - cnew.fMultiplicity=0; - cnew.fX=Float_t(xrec[j]); - cnew.fY=Float_t(yrec[j]); - if (j==0) { - cnew.fQ=Int_t(gChargeTot*qfrac); - } else { - cnew.fQ=Int_t(gChargeTot*(1-qfrac)); - } - gSegmentation->SetHit(xrec[j],yrec[j]); - for (i=0; ifIndexMap[i]; - gSegmentation->SetPad(gix[i], giy[i]); - Float_t q1=gResponse->IntXY(gSegmentation); - cnew.fContMap[cnew.fMultiplicity]=Float_t(q[i])/(q1*cnew.fQ); - cnew.fMultiplicity++; - } - FillCluster(&cnew,0); - //printf("\n x,y %f %f ", cnew.fX, cnew.fY); - cnew.fClusterType=cnew.PhysicsContribution(); - AddRawCluster(cnew); - fNPeaks++; - } - } - - Bool_t fitted=kTRUE; - - if (nLocal !=-100 || !fitted) { - // Check if enough local clusters have been found, - // if not add global maxima to the list - // - Int_t nPerMax; - if (nLocal!=0) { - nPerMax=mul/nLocal; - } else { - printf("\n Warning, no local maximum found \n"); - nPerMax=fNperMax+1; - } - - if (nPerMax > fNperMax) { - Int_t nGlob=mul/fNperMax-nLocal+1; - if (nGlob > 0) { - Int_t nnew=0; - for (i=0; iqmax) { - dmin=d; - qmax=ql; - associatePeak[i]=j; - } - } - } - } - - - // - // One cluster for each maximum - // - for (j=0; jfIndexMap[indLocal[j]]; - cnew.fMultiplicity=1; - for (i=0; ifIndexMap[i]; - cnew.fMultiplicity++; - } - } - FillCluster(&cnew); - cnew.fClusterType=cnew.PhysicsContribution(); - AddRawCluster(cnew); - fNPeaks++; - } + if(AliRICHParam::IsResolveClusters()&&fRawCluster.Size()>1&&fRawCluster.Size()<6){ + FitCoG(); //serialization of resolved clusters will happen inside + }else{//cluster size=1 or resolving is switched off + WriteRawCluster();//simply output the formed raw cluster without deconvolution } -} - - -void AliRICHClusterFinder::FillCluster(AliRICHRawCluster* c, Int_t flag) + fRawCluster.Reset(); fResolvedCluster.Reset(); + }//digits loop for a given chamber + + delete fDigitMap; + + AliDebug(1,"Stop."); +}//FindClusters() +//__________________________________________________________________________________________________ +void AliRICHClusterFinder::FindClusterContribs(AliRICHcluster *pCluster) { -// -// Completes cluster information starting from list of digits -// - AliRICHDigit* dig; - Float_t x, y; - Int_t ix, iy; - Float_t frac=0; - - c->fPeakSignal=0; - if (flag) { - c->fX=0; - c->fY=0; - c->fQ=0; - } - //c->fQ=0; - - - for (Int_t i=0; ifMultiplicity; i++) - { - dig= (AliRICHDigit*)fDigits->UncheckedAt(c->fIndexMap[i]); - ix=dig->fPadX+c->fOffsetMap[i]; - iy=dig->fPadY; - Int_t q=dig->fSignal; - if (dig->fPhysics >= dig->fSignal) { - c->fPhysicsMap[i]=2; - } else if (dig->fPhysics == 0) { - c->fPhysicsMap[i]=0; - } else c->fPhysicsMap[i]=1; -// -// -// peak signal and track list - if (flag) { - if (q>c->fPeakSignal) { - c->fPeakSignal=q; -/* - c->fTracks[0]=dig->fTracks[0]; - c->fTracks[1]=dig->fTracks[1]; - c->fTracks[2]=dig->fTracks[2]; -*/ - //c->fTracks[0]=dig->fTrack; - c->fTracks[0]=dig->fHit; - c->fTracks[1]=dig->fTracks[0]; - c->fTracks[2]=dig->fTracks[1]; - } - } else { - if (c->fContMap[i] > frac) { - frac=c->fContMap[i]; - c->fPeakSignal=q; -/* - c->fTracks[0]=dig->fTracks[0]; - c->fTracks[1]=dig->fTracks[1]; - c->fTracks[2]=dig->fTracks[2]; -*/ - //c->fTracks[0]=dig->fTrack; - c->fTracks[0]=dig->fHit; - c->fTracks[1]=dig->fTracks[0]; - c->fTracks[2]=dig->fTracks[1]; - } - } -// - if (flag) { - fSegmentation->GetPadCxy(ix, iy, x, y); - c->fX += q*x; - c->fY += q*y; - c->fQ += q; - } - - } // loop over digits - - if (flag) { - - c->fX/=c->fQ; - c->fX=fSegmentation->GetAnod(c->fX); - c->fY/=c->fQ; -// -// apply correction to the coordinate along the anode wire -// - x=c->fX; - y=c->fY; - fSegmentation->GetPadIxy(x, y, ix, iy); - fSegmentation->GetPadCxy(ix, iy, x, y); - Int_t isec=fSegmentation->Sector(ix,iy); - TF1* cogCorr = fSegmentation->CorrFunc(isec-1); - - if (cogCorr) { - Float_t yOnPad=(c->fY-y)/fSegmentation->Dpy(isec); - c->fY=c->fY-cogCorr->Eval(yOnPad, 0, 0); - } - } -} - - -void AliRICHClusterFinder::FindCluster(Int_t i, Int_t j, AliRICHRawCluster &c){ -// -// Find clusters -// -// -// Add i,j as element of the cluster -// - - Int_t idx = fHitMap->GetHitIndex(i,j); - AliRICHDigit* dig = (AliRICHDigit*) fHitMap->GetHit(i,j); - Int_t q=dig->fSignal; - if (q > TMath::Abs(c.fPeakSignal)) { - c.fPeakSignal=q; -/* - c.fTracks[0]=dig->fTracks[0]; - c.fTracks[1]=dig->fTracks[1]; - c.fTracks[2]=dig->fTracks[2]; -*/ - //c.fTracks[0]=dig->fTrack; - c.fTracks[0]=dig->fHit; - c.fTracks[1]=dig->fTracks[0]; - c.fTracks[2]=dig->fTracks[1]; - } -// -// Make sure that list of digits is ordered -// - Int_t mu=c.fMultiplicity; - c.fIndexMap[mu]=idx; - - if (dig->fPhysics >= dig->fSignal) { - c.fPhysicsMap[mu]=2; - } else if (dig->fPhysics == 0) { - c.fPhysicsMap[mu]=0; - } else c.fPhysicsMap[mu]=1; - - if (mu > 0) { - for (Int_t ind=mu-1; ind>=0; ind--) { - Int_t ist=(c.fIndexMap)[ind]; - Int_t ql=((AliRICHDigit*)fDigits - ->UncheckedAt(ist))->fSignal; - if (q>ql) { - c.fIndexMap[ind]=idx; - c.fIndexMap[ind+1]=ist; - } else { - break; - } - } +//Finds cerenkov-feedback-mip mixture for a given cluster + AliDebug(1,"Start.");ToAliDebug(1,pCluster->Print()); + +// R()->GetLoader()->GetRunLoader()->LoadHeader(); + AliStack *pStack = R()->GetLoader()->GetRunLoader()->Stack(); + if(!pStack) + {AliInfo("No Stack found!!! No contrib to cluster found.");return;} + + TObjArray *pDigits = pCluster->Digits(); + if(!pDigits) return; //?????????? + Int_t iNmips=0,iNckovs=0,iNfeeds=0; + TArrayI contribs(3*pCluster->Size()); + Int_t *pindex = new Int_t[3*pCluster->Size()]; + for(Int_t iDigN=0;iDigNSize();iDigN++) {//loop on digits of a given cluster + contribs[3*iDigN] =((AliRICHdigit*)pDigits->At(iDigN))->GetTrack(0); + if (contribs[3*iDigN] >= 10000000) contribs[3*iDigN] = 0; + contribs[3*iDigN+1]=((AliRICHdigit*)pDigits->At(iDigN))->GetTrack(1); + if (contribs[3*iDigN+1] >= 10000000) contribs[3*iDigN+1] = 0; + contribs[3*iDigN+2]=((AliRICHdigit*)pDigits->At(iDigN))->GetTrack(2); + if (contribs[3*iDigN+2] >= 10000000) contribs[3*iDigN+2] = 0; + }//loop on digits of a given cluster + TMath::Sort(contribs.GetSize(),contribs.GetArray(),pindex); + for(Int_t iDigN=0;iDigN<3*pCluster->Size()-1;iDigN++) {//loop on digits to sort tids + AliDebug(1,Form("%4i for digit n. %4i",contribs[pindex[iDigN]],iDigN)); + if(contribs[pindex[iDigN]]!=contribs[pindex[iDigN+1]]) { + TParticle* particle = pStack->Particle(contribs[pindex[iDigN]]); + Int_t code = particle->GetPdgCode(); + Double_t charge = 0; + if(particle->GetPDG()) charge=particle->GetPDG()->Charge(); + AliDebug(1,Form(" charge of particle %f",charge)); + + if(code==50000050) iNckovs++; + if(code==50000051) iNfeeds++; + if(charge!=0) iNmips++; } + }//loop on digits to sort Tid + + if (contribs[pindex[3*pCluster->Size()-1]]!=kBad) { + + TParticle* particle = pStack->Particle(contribs[pindex[3*pCluster->Size()-1]]); + Int_t code = particle->GetPdgCode(); + Double_t charge = 0; + if(particle->GetPDG()) charge=particle->GetPDG()->Charge(); + AliDebug(1,Form(" charge of particle %f",charge)); + if(code==50000050) iNckovs++; + if(code==50000051) iNfeeds++; + if(charge!=0) iNmips++; + } - c.fMultiplicity++; - - if (c.fMultiplicity >= 50 ) { - printf("FindCluster - multiplicity >50 %d \n",c.fMultiplicity); - c.fMultiplicity=49; - } - -// Prepare center of gravity calculation - Float_t x, y; - fSegmentation->GetPadCxy(i, j, x, y); - c.fX += q*x; - c.fY += q*y; - c.fQ += q; -// Flag hit as taken - fHitMap->FlagHit(i,j); -// -// Now look recursively for all neighbours + pCluster->CFM(iNckovs,iNfeeds,iNmips); // - Int_t nn; - Int_t xList[kMaxNeighbours], yList[kMaxNeighbours]; - fSegmentation->Neighbours(i,j,&nn,xList,yList); - for (Int_t in=0; inTestHit(ix,iy)==kUnused) FindCluster(ix, iy, c); - } -} - -//_____________________________________________________________________________ - -void AliRICHClusterFinder::FindRawClusters() + delete [] pindex; + ToAliDebug(1,pCluster->Print()); + AliDebug(1,"Stop."); +}//FindClusterContribs() +//__________________________________________________________________________________________________ +void AliRICHClusterFinder::FormRawCluster(Int_t i, Int_t j) { - // - // simple RICH cluster finder from digits -- finds neighbours and - // fill the tree with raw clusters - // - if (!fNdigits) return; - - fHitMap = new AliRICHHitMapA1(fSegmentation, fDigits); - - AliRICHDigit *dig; - - //printf ("Now I'm here"); - - Int_t ndig; - Int_t nskip=0; - Int_t ncls=0; - fHitMap->FillHits(); - for (ndig=0; ndigUncheckedAt(ndig); - Int_t i=dig->fPadX; - Int_t j=dig->fPadY; - if (fHitMap->TestHit(i,j)==kUsed ||fHitMap->TestHit(i,j)==kEmpty) { - nskip++; - continue; - } - AliRICHRawCluster c; - c.fMultiplicity=0; - c.fPeakSignal=dig->fSignal; -/* - c.fTracks[0]=dig->fTracks[0]; - c.fTracks[1]=dig->fTracks[1]; - c.fTracks[2]=dig->fTracks[2]; -*/ - //c.fTracks[0]=dig->fTrack; - c.fTracks[0]=dig->fHit; - c.fTracks[1]=dig->fTracks[0]; - c.fTracks[2]=dig->fTracks[1]; - // tag the beginning of cluster list in a raw cluster - c.fNcluster[0]=-1; - FindCluster(i,j, c); - // center of gravity - c.fX /= c.fQ; - c.fX=fSegmentation->GetAnod(c.fX); - c.fY /= c.fQ; -// -// apply correction to the coordinate along the anode wire -// - Int_t ix,iy; - Float_t x=c.fX; - Float_t y=c.fY; - fSegmentation->GetPadIxy(x, y, ix, iy); - fSegmentation->GetPadCxy(ix, iy, x, y); - Int_t isec=fSegmentation->Sector(ix,iy); - TF1* cogCorr=fSegmentation->CorrFunc(isec-1); - if (cogCorr) { - Float_t yOnPad=(c.fY-y)/fSegmentation->Dpy(isec); - c.fY=c.fY-cogCorr->Eval(yOnPad,0,0); - } - -// -// Analyse cluster and decluster if necessary -// - ncls++; - c.fNcluster[1]=fNRawClusters; - c.fClusterType=c.PhysicsContribution(); - Decluster(&c); - fNPeaks=0; -// -// -// -// reset Cluster object - for (int k=0;kGetHit(i,j))->Q())); + + fRawCluster.AddDigit((AliRICHdigit*) fDigitMap->GetHit(i,j));//take this pad in cluster + fDigitMap->FlagHit(i,j);//flag this pad as taken + + Int_t listX[4], listY[4]; // Now look recursively for all neighbours + for (Int_t iNei=0;iNeiP()->PadNeighbours(i,j,listX,listY);iNei++) + if(fDigitMap->TestHit(listX[iNei],listY[iNei])==kUnused) FormRawCluster(listX[iNei],listY[iNei]); +}//FormRawCluster() +//__________________________________________________________________________________________________ +void AliRICHClusterFinder::FindLocalMaxima() { - -// Calibration - - Float_t x[5]; - Float_t y[5]; - Int_t n, i; - TF1 func; - if (fSegmentation) { - fSegmentation->GiveTestPoints(n, x, y); - for (i=0; iSetCorrFunc(i, new TF1(func)); - } +//find number of local maxima in the current raw cluster and then calculates initial center of gravity + fNlocals=0; + for(Int_t iDig1=0;iDig1At(iDig1); + TVector pad1 = pDig1->Pad(); + Int_t padQ1 = (Int_t)(pDig1->Q()+0.1); + Int_t padC1 = pDig1->ChFbMi(); + for(Int_t iDig2=0;iDig2At(iDig2); + TVector pad2 = pDig2->Pad(); + Int_t padQ2 = (Int_t)(pDig2->Q()+0.1); + if(iDig1==iDig2) continue; + Int_t diffx = TMath::Sign(Int_t(pad1[0]-pad2[0]),1); + Int_t diffy = TMath::Sign(Int_t(pad1[1]-pad2[1]),1); + if((diffx+diffy)<=1) { + if(padQ2>=padQ1) iNotMax++; + } } -} - - -void AliRICHClusterFinder:: -SinoidalFit(Float_t x, Float_t y, TF1 &func) -{ -// Sinoidal fit - - - static Int_t count=0; - char canvasname[3]; - count++; - sprintf(canvasname,"c%d",count); - - const Int_t kNs=101; - Float_t xg[kNs], yg[kNs], xrg[kNs], yrg[kNs]; - Float_t xsig[kNs], ysig[kNs]; - - AliRICHSegmentation *segmentation=fSegmentation; - - Int_t ix,iy; - segmentation->GetPadIxy(x,y,ix,iy); - segmentation->GetPadCxy(ix,iy,x,y); - Int_t isec=segmentation->Sector(ix,iy); -// Pad Limits - Float_t xmin = x-segmentation->Dpx(isec)/2; - Float_t ymin = y-segmentation->Dpy(isec)/2; -// -// Integration Limits - Float_t dxI=fResponse->SigmaIntegration()*fResponse->ChargeSpreadX(); - Float_t dyI=fResponse->SigmaIntegration()*fResponse->ChargeSpreadY(); - -// -// Scanning -// - Int_t i; - Float_t qp; -// -// y-position - Float_t yscan=ymin; - Float_t dy=segmentation->Dpy(isec)/(kNs-1); - - for (i=0; iSigGenInit(x, yscan, 0); - - for (segmentation->FirstPad(x, yscan, dxI, dyI); - segmentation->MorePads(); - segmentation->NextPad()) - { - qp=fResponse->IntXY(segmentation); - qp=TMath::Abs(qp); -// -// - if (qp > 1.e-4) { - qcheck+=qp; - Int_t ixs=segmentation->Ix(); - Int_t iys=segmentation->Iy(); - Float_t xs,ys; - segmentation->GetPadCxy(ixs,iys,xs,ys); - sum+=qp*ys; - } - } // Pad loop - Float_t ycog=sum/qcheck; - yg[i]=(yscan-y)/segmentation->Dpy(isec); - yrg[i]=(ycog-y)/segmentation->Dpy(isec); - ysig[i]=ycog-yscan; - yscan+=dy; - } // scan loop -// -// x-position - Float_t xscan=xmin; - Float_t dx=segmentation->Dpx(isec)/(kNs-1); - - for (i=0; iSigGenInit(xscan, y, 0); - - for (segmentation->FirstPad(xscan, y, dxI, dyI); - segmentation->MorePads(); - segmentation->NextPad()) - { - qp=fResponse->IntXY(segmentation); - qp=TMath::Abs(qp); -// -// - if (qp > 1.e-2) { - qcheck+=qp; - Int_t ixs=segmentation->Ix(); - Int_t iys=segmentation->Iy(); - Float_t xs,ys; - segmentation->GetPadCxy(ixs,iys,xs,ys); - sum+=qp*xs; - } - } // Pad loop - Float_t xcog=sum/qcheck; - xcog=segmentation->GetAnod(xcog); - - xg[i]=(xscan-x)/segmentation->Dpx(isec); - xrg[i]=(xcog-x)/segmentation->Dpx(isec); - xsig[i]=xcog-xscan; - xscan+=dx; + if(iNotMax==0) { + TVector2 x2=AliRICHParam::Pad2Loc(pad1); + fLocalX[fNlocals]=x2.X();fLocalY[fNlocals]=x2.Y(); + fLocalQ[fNlocals] = (Double_t)padQ1; + fLocalC[fNlocals] = padC1; + fNlocals++; } -// -// Creates a Root function based on function sinoid above -// and perform the fit -// - // TGraph *graphx = new TGraph(kNs,xg ,xsig); - // TGraph *graphxr= new TGraph(kNs,xrg,xsig); - // TGraph *graphy = new TGraph(kNs,yg ,ysig); - TGraph *graphyr= new TGraph(kNs,yrg,ysig); - - Double_t sinoid(Double_t *x, Double_t *par); - new TF1("sinoidf",sinoid,0.5,0.5,5); - graphyr->Fit("sinoidf","Q"); - func = *((TF1*)((graphyr->GetListOfFunctions())->At(0))); -/* - - TCanvas *c1=new TCanvas(canvasname,canvasname,400,10,600,700); - TPad* pad11 = new TPad("pad11"," ",0.01,0.51,0.49,0.99); - TPad* pad12 = new TPad("pad12"," ",0.51,0.51,0.99,0.99); - TPad* pad13 = new TPad("pad13"," ",0.01,0.01,0.49,0.49); - TPad* pad14 = new TPad("pad14"," ",0.51,0.01,0.99,0.49); - pad11->SetFillColor(11); - pad12->SetFillColor(11); - pad13->SetFillColor(11); - pad14->SetFillColor(11); - pad11->Draw(); - pad12->Draw(); - pad13->Draw(); - pad14->Draw(); - -// - pad11->cd(); - graphx->SetFillColor(42); - graphx->SetMarkerColor(4); - graphx->SetMarkerStyle(21); - graphx->Draw("AC"); - graphx->GetHistogram()->SetXTitle("x on pad"); - graphx->GetHistogram()->SetYTitle("xcog-x"); - - - pad12->cd(); - graphxr->SetFillColor(42); - graphxr->SetMarkerColor(4); - graphxr->SetMarkerStyle(21); - graphxr->Draw("AP"); - graphxr->GetHistogram()->SetXTitle("xcog on pad"); - graphxr->GetHistogram()->SetYTitle("xcog-x"); - - - pad13->cd(); - graphy->SetFillColor(42); - graphy->SetMarkerColor(4); - graphy->SetMarkerStyle(21); - graphy->Draw("AF"); - graphy->GetHistogram()->SetXTitle("y on pad"); - graphy->GetHistogram()->SetYTitle("ycog-y"); - - - - pad14->cd(); - graphyr->SetFillColor(42); - graphyr->SetMarkerColor(4); - graphyr->SetMarkerStyle(21); - graphyr->Draw("AF"); - graphyr->GetHistogram()->SetXTitle("ycog on pad"); - graphyr->GetHistogram()->SetYTitle("ycog-y"); - - c1->Update(); -*/ -} - -Double_t sinoid(Double_t *x, Double_t *par) + } + fRawCluster.CoG(fNlocals); //first initial approximation of the CoG...to start minimization. +}//FindLocalMaxima() +//__________________________________________________________________________________________________ +void AliRICHClusterFinder::WriteRawCluster() { - -// Sinoid function - - Double_t arg = -2*TMath::Pi()*x[0]; - Double_t fitval= par[0]*TMath::Sin(arg)+ - par[1]*TMath::Sin(2*arg)+ - par[2]*TMath::Sin(3*arg)+ - par[3]*TMath::Sin(4*arg)+ - par[4]*TMath::Sin(5*arg); - return fitval; - } - - -Double_t DoubleGauss(Double_t *x, Double_t *par) +//Add the current raw cluster to the list of clusters + AliDebug(1,"Start."); + + FindClusterContribs(&fRawCluster); + R()->AddCluster(fRawCluster); + + ToAliDebug(1,fRawCluster.Print()); AliDebug(1,"Stop."); +}//WriteRawCluster() +//__________________________________________________________________________________________________ +void AliRICHClusterFinder::WriteResolvedCluster() { - -// Doublr gaussian function - - Double_t arg1 = (x[0]-par[1])/0.18; - Double_t arg2 = (x[0]-par[3])/0.18; - Double_t fitval= par[0]*TMath::Exp(-arg1*arg1/2) - +par[2]*TMath::Exp(-arg2*arg2/2); - return fitval; - } - -Float_t DiscrCharge(Int_t i,Double_t *par) +//Add the current resolved cluster to the list of clusters + AliDebug(1,"Start."); + + FindClusterContribs(&fResolvedCluster); + R()->AddCluster(fResolvedCluster); + + ToAliDebug(1,fResolvedCluster.Print()); AliDebug(1,"Stop."); +}//WriteResolvedCluster() +//__________________________________________________________________________________________________ +void AliRICHClusterFinder::FitCoG() { -// par[0] x-position of first cluster -// par[1] y-position of first cluster -// par[2] x-position of second cluster -// par[3] y-position of second cluster -// par[4] charge fraction of first cluster -// 1-par[4] charge fraction of second cluster - - static Float_t qtot; - if (gFirst) { - qtot=0; - for (Int_t jbin=0; jbinSetPad(gix[i], giy[i]); -// First Cluster - gSegmentation->SetHit(par[0],par[1]); - Float_t q1=gResponse->IntXY(gSegmentation); - -// Second Cluster - gSegmentation->SetHit(par[2],par[3]); - Float_t q2=gResponse->IntXY(gSegmentation); - - Float_t value = qtot*(par[4]*q1+(1.-par[4])*q2); - return value; -} - -// -// Minimisation function -void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag) +//Fits cluster of size by the corresponding number of Mathieson shapes. +//This methode is only invoked in case everything is ok to start deconvolution + AliDebug(1,"Start with:"); ToAliDebug(1,fRawCluster.Print()); + + Double_t arglist; + Int_t ierflag = 0; + + TMinuit *pMinuit = new TMinuit(3*fNlocals-1); + pMinuit->mninit(5,10,7); + + arglist = -1; + pMinuit->mnexcm("SET PRI",&arglist, 1, ierflag); + pMinuit->mnexcm("SET NOW",&arglist, 0, ierflag); + + TString chname; + Int_t ierflg; + + pMinuit->SetObjectFit((TObject*)this); + pMinuit->SetFCN(RICHMinMathieson); + + Double_t vstart,lower, upper; + Double_t stepX= 0.001; + Double_t stepY= 0.001; + Double_t stepQ= 0.0001; + + for(Int_t i=0;imnparm(3*i ,Form("xCoG %i",i),vstart,stepX,lower,upper,ierflag); + vstart = fLocalY[i]; + lower = vstart - 2*AliRICHParam::PadSizeY(); + upper = vstart + 2*AliRICHParam::PadSizeY(); + pMinuit->mnparm(3*i+1,Form("yCoG %i",i),vstart,stepY,lower,upper,ierflag); + if(i==fNlocals-1) break; // last parameter is constrained + vstart = fLocalQ[i]/fRawCluster.Q(); + lower = 0; + upper = 1; + pMinuit->mnparm(3*i+2,Form("qfrac %i",i),vstart,stepQ,lower,upper,ierflag); + } + + arglist = -1; + pMinuit->mnexcm("SET NOGR",&arglist, 1, ierflag); + arglist = 1; + pMinuit->mnexcm("SET ERR", &arglist, 1,ierflg); + arglist = -1; + pMinuit->mnexcm("SIMPLEX",&arglist, 0, ierflag); + pMinuit->mnexcm("MIGRAD",&arglist, 0, ierflag); + pMinuit->mnexcm("EXIT" ,&arglist, 0, ierflag); + + Double_t xCoG[50],yCoG[50],qfracCoG[50]; + Double_t eps, b1, b2; + + Double_t qfraclast=0; + for(Int_t i=0;imnpout(3*i ,chname, xCoG[i], eps , b1, b2, ierflg); + pMinuit->mnpout(3*i+1,chname, yCoG[i], eps , b1, b2, ierflg); + if(i==fNlocals-1) break; + pMinuit->mnpout(3*i+2,chname, qfracCoG[i], eps , b1, b2, ierflg); + qfraclast+=qfracCoG[i]; + } + qfracCoG[fNlocals-1] = 1 - qfraclast; + delete pMinuit; + + for(Int_t i=0;iGetObjectFit())->GetRawCluster(); + + TVector2 centroid[50]; + Double_t q[50]; + Int_t nFunctions = (npar+1)/3; + Double_t qfract = 0; + for(Int_t i=0;iQ(); + for(Int_t i=0;iSize();i++) { + TVector pad=((AliRICHdigit *)pRawCluster->Digits()->At(i))->Pad(); + Double_t padQ = ((AliRICHdigit *)pRawCluster->Digits()->At(i))->Q(); + Double_t qfracpar=0; + for(Int_t j=0;jGetEntriesFast(); -} - -AliRICHClusterFinder& AliRICHClusterFinder::operator=(const AliRICHClusterFinder& rhs) -{ -// Assignment operator - return *this; - -} - - - - - - - - - - + chi2 += TMath::Power((qtot*qfracpar-padQ),2)/padQ; + } +}//RICHMinMathieson() +//__________________________________________________________________________________________________