RICH/AliRICHCluster.cxx

   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
  16 #include "AliRICHCluster.h"
  17 #include <TMinuit.h>  //Solve()
  18
  19 ClassImp(AliRICHCluster)
  20 //__________________________________________________________________________________________________
  21 void AliRICHCluster::Print(Option_t*)const
  22 {
  23 //Print current cluster
  24   const char *status=0;
  25   switch(fStatus){
  26     case      kFormed: status="formed"     ;break;
  27     case    kUnfolded: status="unfolded"   ;break;
  28     case         kCoG: status="CoGed"      ;break;
  29     case       kEmpty: status="empty"      ;break;
  30   }
  31   Int_t iNdigs=0;  if(fDigits) iNdigs=fDigits->GetEntriesFast();
  32
  33   Printf("cfm=%10i, cs=%2i, Size=%2i Maxima=%2i, Shape=%5i, pos=(%7.3f,%7.3f) Q=%6i, %s",
  34             fCFM,fChamber,Size(),Nlocmax(),fShape,fX,fY,fQdc,status,iNdigs);
  35   for(Int_t i=0;i<iNdigs;i++) Digit(i)->Print();
  36 }//Print()
  37
  38
  39 //__________________________________________________________________________________________________
  40 TMinuit* AliRICHCluster::Solve()
  41 {
  42 //At this point, cluster contains a list of digits, cluster charge is precalculated as a sum of digits charges (in AddDigit()),
  43 //position is preset to (-1,-1) (in ctor), status is preset to kFormed in (AddDigit()), chamber-sector info is preseted to actual value (in AddDigit())
  44 //Here we decide what to do with this cluster: unfold or just calculate center of gravity
  45 //Arguments: none
  46 //  Returns: pointer to fitter or 0 if no unfolding decided
  47   TMinuit *pMinuit=0;
  48   if(Size()>=2 && AliRICHParam::IsResolveClusters())
  49     pMinuit=Unfold();
  50   else
  51     CoG(0);
  52   return pMinuit;
  53 }//Solve()
  54 //__________________________________________________________________________________________________
  55 void AliRICHCluster::FitFunc(Int_t &iNpars, Double_t *, Double_t &chi2, Double_t *aPar, Int_t )
  56 {
  57 //Cluster fit function
  58 //par[0]=x par[1]=y par[2]=q for the first Mathieson shape
  59 //par[3]=x par[4]=y par[5]=q for the second Mathieson shape and so on up to iNpars/3 Mathieson shapes
  60 //We need to calculate Qpad - Qpadmath summup over all pads of the cluster
  61 //Here Qpad is a actual charge of the pad, Qpadmath is calculated charge of the pad induced by all Mathiesons
  62 //Arguments: iNpars - number of parameters which is number of local maxima of cluster * 3
  63 //           chi2   - function result to be minimised
  64 //           aPar   - parametrs array of size iNpars
  65 //  Returns: none
  66   AliRICHCluster *pClu=(AliRICHCluster*)gMinuit->GetObjectFit();
  67   Int_t iNmathiesons = iNpars/3;
  68
  69   TVector2 curMathiesonPos;
  70   chi2 = 0;
  71   for(Int_t i=0;i<pClu->Size();i++){//digits loop
  72     TVector    pad     = pClu->Digit(i)->Pad();
  73     Double_t dQpad     = pClu->Digit(i)->Qdc();
  74     Double_t dQpadmath = 0;
  75     for(Int_t j=0;j<iNmathiesons;j++){//all Mathiesons may contribute to this pad
  76       curMathiesonPos.Set(aPar[3*j],aPar[3*j+1]);//get current position for current Mathieson
  77       dQpadmath += aPar[3*j+2]*AliRICHParam::FracQdc(curMathiesonPos,pad);//sums up contributions to the current pad from all Mathiesons
  78     }
  79     chi2 += TMath::Power((dQpadmath-dQpad),2);
  80   }//digits loop
  81 }//RichClusterFitFunction()
  82 //__________________________________________________________________________________________________
  83 TMinuit* AliRICHCluster::Unfold()
  84 {
  85 //This methode is invoked from Solve() when decided to unfold this cluster
  86 //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
  87 //Arguments: none
  88 //  Returns: pointer to fitter for retriving parameters
  89
  90   TMinuit *pMinuit = new TMinuit(15); //init MINUIT with max 15 parameters (maxim 5 mathiesons, 3 params per matheson )
  91   pMinuit->SetObjectFit((TObject*)this);
  92   pMinuit->SetFCN(AliRICHCluster::FitFunc);//set fit function
  93   Double_t aArg=-1,parStart,parStep,parLow,parHigh; Int_t iErrFlg; //tmp for MINUIT parameters definitions
  94   pMinuit->mnexcm("SET PRI" ,&aArg,1,iErrFlg); //suspend all printout from TMinuit
  95
  96   Int_t iLocMaxCnt=0;
  97 //Strategy is to check if the current pad has QDC more then all neigbours
  98   for(Int_t iDig1=0;iDig1<Size();iDig1++) {//first digits loop
  99     AliRICHDigit *pDig1 = Digit(iDig1);//take the current digit
 100     Int_t iHowManyMoreCnt = 0;//counts how many neighbouring pads has QDC more then current one
 101     for(Int_t iDig2=0;iDig2<Size();iDig2++) {//loop on all digits again
 102       AliRICHDigit *pDig2 = Digit(iDig2);
 103       if(iDig1==iDig2) continue;             //no need to compare
 104       Int_t dist = TMath::Sign(Int_t(pDig1->PadX()-pDig2->PadX()),1)+TMath::Sign(Int_t(pDig1->PadY()-pDig2->PadY()),1);//distance between pads
 105       if(dist==1)//means pads are neighbours
 106          if(pDig2->Qdc()>=pDig1->Qdc()) iHowManyMoreCnt++;//count number of pads with Q more then Q of current pad
 107     }//second digits loop
 108     if(iHowManyMoreCnt==0&&iLocMaxCnt<6){//this pad has Q more then any neighbour so it's local maximum
 109         TVector2 x2=AliRICHParam::Pad2Loc(pDig1->Pad());//take pad center position and use it as parameter for current Mathienson shape
 110         pMinuit->mnparm(3*iLocMaxCnt  ,Form("x%i",iLocMaxCnt),parStart=x2.X()      ,parStep=0.01,parLow=0,parHigh=0,iErrFlg);
 111         pMinuit->mnparm(3*iLocMaxCnt+1,Form("y%i",iLocMaxCnt),parStart=x2.Y()      ,parStep=0.01,parLow=0,parHigh=0,iErrFlg);
 112         pMinuit->mnparm(3*iLocMaxCnt+2,Form("q%i",iLocMaxCnt),parStart=pDig1->Qdc(),parStep=0.01,parLow=0,parHigh=0,iErrFlg);//
 113         iLocMaxCnt++;
 114     }//if this pad is local maximum
 115   }//first digits loop
 116
 117   fSize+=iLocMaxCnt;
 118   if(iLocMaxCnt>0&&iLocMaxCnt<6){ //resonable number of local maxima to fit
 119     Double_t aArg=0;
 120     pMinuit->mnexcm("MIGRAD",&aArg,0,iErrFlg);//start fitting
 121     fStatus=kUnfolded;
 122   }else{
 123     delete pMinuit;
 124     pMinuit=0;
 125     CoG(0);
 126   }
 127   return pMinuit;
 128 }//Unfold()
 129 //__________________________________________________________________________________________________
 130 void AliRICHCluster::CoG(Int_t nLocals)
 131 {
 132 //Calculates naive cluster position as a center of gravity of its digits.
 133 //Also determines the box fully contaning this cluster
 134 //Arguments:
 135   Float_t xmin=999,ymin=999,xmax=0,ymax=0;
 136   fX=fY=0;
 137   for(Int_t iDig=0;iDig<Size();iDig++) {
 138     AliRICHDigit *pDig=Digit(iDig);
 139     TVector pad=pDig->Pad(); Double_t q=pDig->Qdc();
 140     TVector2 x2=AliRICHParam::Pad2Loc(pad);
 141     fX += x2.X()*q;fY +=x2.Y()*q;
 142     if(pad[0]<xmin)xmin=pad[0];if(pad[0]>xmax)xmax=pad[0];if(pad[1]<ymin)ymin=pad[1];if(pad[1]>ymax)ymax=pad[1];
 143    }
 144    fX/=fQdc;fY/=fQdc;//Center of Gravity
 145
 146    TVector2 center = AliRICHParam::Pad2Loc(AliRICHParam::Loc2Pad(TVector2(fX,fY)));
 147    fX += AliRICHParam::CogCorr(fX-center.X());//correct cluster position for sinoid
 148
 149    fShape=Int_t(100*(xmax-xmin+1)+ymax-ymin+1);//find box containing cluster
 150    fSize+=nLocals;
 151    fStatus=kCoG;
 152 }//CoG()
 153 //__________________________________________________________________________________________________
 154 void AliRICHCluster::Test(const TVector2 &hitX2,Double_t dEloss)
 155 {
 156 //This is to test all cluster functionality
 157 //Method uses AddDigit() to add a predifined pad structure and then calls Solve
 158   Int_t iQtot=AliRICHParam::TotQdc(hitX2,dEloss);
 159   if(iQtot==0){
 160     Printf("Provided hit position out of sensitive area");
 161     return;
 162   }
 163   TVector area=AliRICHParam::Loc2Area(hitX2);
 164   TVector pad(2);
 165   for(pad[1]=area[1];pad[1]<=area[3];pad[1]++){//affected pads loop first y
 166     for(pad[0]=area[0];pad[0]<=area[2];pad[0]++){//then x
 167       Double_t dQpad=iQtot*AliRICHParam::FracQdc(hitX2,pad);//charge fraction from Mathieson centered at x to pad
 168       AddDigit(new AliRICHDigit(3,(Int_t)pad[0],(Int_t)pad[1],dQpad));
 169     }//affected pads loop
 170   }
 171   TMinuit *pMinuit=Solve();
 172   Print();
 173   Printf("Initial hit (%.2f,%.2f) Qtot=%i Eloss=%.2f",hitX2.X(),hitX2.Y(),iQtot,dEloss);
 174   Double_t d1,d2,d3; Int_t iErrFlg;TString sName; //tmp vars for TMinuit
 175   Double_t x,y,q;
 176   for(Int_t i=0;i<Nlocmax();i++){//retrive fitting results
 177     pMinuit->mnpout(3*i   ,sName,  x, d1 , d2, d3, iErrFlg);
 178     pMinuit->mnpout(3*i+1 ,sName,  y, d1 , d2, d3, iErrFlg);
 179     pMinuit->mnpout(3*i+2 ,sName,  q, d1 , d2, d3, iErrFlg);
 180   }
 181   Printf(" Fitted hit (%.2f,%.2f) Qfit=%.0f",x,y,q);
 182   delete pMinuit; pMinuit=0;  Reset();
 183 }//Test()