First interation on SplitByLocalMaxima with Nico
authorkir <kir@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 21 Oct 2003 18:18:17 +0000 (18:18 +0000)
committerkir <kir@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 21 Oct 2003 18:18:17 +0000 (18:18 +0000)
RICH/AliRICHClusterFinder.cxx
RICH/AliRICHClusterFinder.h
RICH/AliRICHParam.cxx
RICH/AliRICHParam.h
RICH/AliRICHPatRec.cxx
RICH/AliRICHRawCluster.cxx
RICH/AliRICHRawCluster.h
RICH/menu.C

index 4e223eee18bd646d953e12ed785bd8683fb6ae48..7714bf5b63a961728ae68cd16d8edc8eb2109675 100644 (file)
 #include "AliRICHDigit.h"
 #include "AliRICHRawCluster.h"
 #include "AliRICHParam.h"
-
-
-
 #include <TTree.h>
 #include <TCanvas.h>
 #include <TH1.h>
 #include <TF1.h>
 #include <TPad.h>
 #include <TGraph.h> 
-#include <TPostScript.h> 
 #include <TMinuit.h>
 
-//----------------------------------------------------------
 static AliSegmentation     *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 Bool_t               gFirst=kTRUE;
 static TMinuit *gMyMinuit ;
-void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag);
+void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t);
 static Int_t                gChargeTot;
 
 ClassImp(AliRICHClusterFinder)
@@ -64,7 +59,6 @@ AliRICHClusterFinder::AliRICHClusterFinder(AliRICH *pRICH)
   fCogCorr = 0;
   SetNperMax();
   SetClusterSize();
-  SetDeclusterFlag();
   fNPeaks=-1;
 }//main ctor
 //__________________________________________________________________________________________________
@@ -76,22 +70,10 @@ void AliRICHClusterFinder::Decluster(AliRICHRawCluster *cluster)
     AddRawCluster(*cluster); 
     fNPeaks++;
   }else if(cluster->fMultiplicity==3){// 3-cluster, check topology
-    if(fDeclusterFlag)
       Centered(cluster);// ok, cluster is centered and added in Centered()
-    else{//if(fDeclusterFlag)
-      if(fNPeaks!=0){cluster->fNcluster[0]=fNPeaks;cluster->fNcluster[1]=0;} 
-      AddRawCluster(*cluster); 
-      fNPeaks++;
-    }//if(fDeclusterFlag)          
   }else{//4-and more-pad clusters
-    if(cluster->fMultiplicity<= fClusterSize){
-      if(fDeclusterFlag)
+    if(cluster->fMultiplicity<= fMaxClusterSize){
         SplitByLocalMaxima(cluster);
-      else{
-        if(fNPeaks!= 0){cluster->fNcluster[0]=fNPeaks; cluster->fNcluster[1]=0;        } 
-        AddRawCluster(*cluster);
-        fNPeaks++;
-      }//if(fDeclusterFlag)    
     }//if <= fClusterSize
   }//if multiplicity 
 }//Decluster()
@@ -168,7 +150,9 @@ Bool_t AliRICHClusterFinder::Centered(AliRICHRawCluster *cluster)
 //__________________________________________________________________________________________________
 void AliRICHClusterFinder::SplitByLocalMaxima(AliRICHRawCluster *c)
 {// Split the cluster according to the number of maxima inside
-    AliRICHDigit* dig[100], *digt;
+  Info("SplitbyLocalMaxima","Start.");
+  
+      AliRICHDigit* dig[100], *digt;
     Int_t ix[100], iy[100], q[100];
     Float_t x[100], y[100];
     Int_t i; // loops over digits
@@ -180,7 +164,7 @@ void AliRICHClusterFinder::SplitByLocalMaxima(AliRICHRawCluster *c)
     ix[i]= dig[i]->PadX();
     iy[i]= dig[i]->PadY();
     q[i] = dig[i]->Signal();
-    Rich()->Param()->Pad2Local(ix[i], iy[i], x[i], y[i]);
+    AliRICHParam::Pad2Local(ix[i], iy[i], x[i], y[i]);
   }
 //  Find local maxima
     Bool_t isLocal[100];
@@ -228,10 +212,10 @@ void AliRICHClusterFinder::SplitByLocalMaxima(AliRICHRawCluster *c)
            if (nnew==1) break;
        }
     }
-// If number of local maxima is 2 try to fit a double gaussian
-    if (nLocal==-100) {
+    if(nLocal==2) {// If number of local maxima is 2 try to fit a double gaussian
+
 //  Initialise global variables for fit
-       gFirst=1;
+       gFirst=kTRUE;
        gSegmentation=fSegmentation;
        gResponse    =fResponse;
        gNbins=mul;
@@ -241,14 +225,11 @@ void AliRICHClusterFinder::SplitByLocalMaxima(AliRICHRawCluster *c)
            giy[i]=iy[i];
            gCharge[i]=Float_t(q[i]);
        }
-       if (gFirst) {
-           gFirst=kFALSE;
-           gMyMinuit = new TMinuit(5);
-       }
+       if (gFirst)    gMyMinuit = new TMinuit(5);
+       
        gMyMinuit->SetFCN(fcn);
        gMyMinuit->mninit(5,10,7);
        Double_t arglist[20];
-       Int_t ierflag=0;
        arglist[0]=1;
 // Set starting values 
        static Double_t vstart[5];
@@ -256,55 +237,54 @@ void AliRICHClusterFinder::SplitByLocalMaxima(AliRICHRawCluster *c)
        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]]);
+       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[0]=vstart[0]-AliRICHParam::PadSizeX()/2;
+       lower[1]=vstart[1]-AliRICHParam::PadSizeY()/2;
        
-       upper[0]=lower[0]+fSegmentation->Dpx(isec);
-       upper[1]=lower[1]+fSegmentation->Dpy(isec);
+       upper[0]=vstart[0]+AliRICHParam::PadSizeX()/2;
+       upper[1]=vstart[1]+AliRICHParam::PadSizeY()/2;
        
-       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[2]=vstart[2]-AliRICHParam::PadSizeX()/2;
+       lower[3]=vstart[3]-AliRICHParam::PadSizeY()/2;
        
-       upper[2]=lower[2]+fSegmentation->Dpx(isec);
-       upper[3]=lower[3]+fSegmentation->Dpy(isec);
+       upper[2]=vstart[2]+AliRICHParam::PadSizeX()/2;
+       upper[3]=vstart[3]+AliRICHParam::PadSizeY()/2;
        
        lower[4]=0.;
        upper[4]=1.;
 // step sizes
        static Double_t step[5]={0.005, 0.03, 0.005, 0.03, 0.01};
+       Int_t iErr;
        
-       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);
+       gMyMinuit->mnparm(0,"x1",vstart[0],step[0],lower[0],upper[0],iErr);
+       gMyMinuit->mnparm(1,"y1",vstart[1],step[1],lower[1],upper[1],iErr);
+       gMyMinuit->mnparm(2,"x2",vstart[2],step[2],lower[2],upper[2],iErr);
+       gMyMinuit->mnparm(3,"y2",vstart[3],step[3],lower[3],upper[3],iErr);
+       gMyMinuit->mnparm(4,"a0",vstart[4],step[4],lower[4],upper[4],iErr);
 // ready for minimisation      
        gMyMinuit->SetPrintLevel(-1);
-       gMyMinuit->mnexcm("SET OUT", arglist, 0, ierflag);
+       gMyMinuit->mnexcm("SET OUT", arglist, 0, iErr);
        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);
+       gMyMinuit->mnexcm("SET NOGR", arglist, 0, iErr);
+       gMyMinuit->mnexcm("SIMPLEX", arglist, 0, iErr);
+       gMyMinuit->mnexcm("MIGRAD", arglist, 0, iErr);
+       gMyMinuit->mnexcm("EXIT" , arglist, 0, iErr);
 
        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);    
- // One cluster for each maximum
-       for (j=0; j<2; j++) {
+       gMyMinuit->mnpout(0, chname, xrec[0], epxz, b1, b2, iErr);      
+       gMyMinuit->mnpout(1, chname, yrec[0], epxz, b1, b2, iErr);      
+       gMyMinuit->mnpout(2, chname, xrec[1], epxz, b1, b2, iErr);      
+       gMyMinuit->mnpout(3, chname, yrec[1], epxz, b1, b2, iErr);      
+       gMyMinuit->mnpout(4, chname, qfrac,   epxz, b1, b2, iErr);      
+        
+        cout<<"xrex[0]="<<xrec[0]<<"yrec[0]="<<yrec[0]<<"xrec[1]="<<xrec[1]<<"yrec[1]="<<yrec[1]<<"qfrac="<<qfrac<<endl;
+       for (j=0; j<2; j++) { // One cluster for each maximum
            AliRICHRawCluster cnew;
            if (fNPeaks == 0) {
                cnew.fNcluster[0]=-1;
@@ -321,12 +301,9 @@ void AliRICHClusterFinder::SplitByLocalMaxima(AliRICHRawCluster *c)
            } else {
                cnew.fQ=Int_t(gChargeTot*(1-qfrac));
            }
-           gSegmentation->SetHit(xrec[j],yrec[j],0);
            for (i=0; i<mul; i++) {
                cnew.fIndexMap[cnew.fMultiplicity]=c->fIndexMap[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.fContMap[cnew.fMultiplicity]=AliRICHParam::AssignChargeToPad(xrec[j],yrec[j],gix[i], giy[i]);
                cnew.fMultiplicity++;
            }
            FillCluster(&cnew,0);
@@ -334,11 +311,10 @@ void AliRICHClusterFinder::SplitByLocalMaxima(AliRICHRawCluster *c)
            AddRawCluster(cnew);
            fNPeaks++;
        }
-    }
-
+    }//if 2 maximum in cluster
     Bool_t fitted=kTRUE;
 
-    if (nLocal !=-100 || !fitted) {
+    if (nLocal >2 || !fitted) {
        // Check if enough local clusters have been found, if not add global maxima to the list 
        Int_t nPerMax;
        if (nLocal!=0) {
@@ -415,9 +391,9 @@ void AliRICHClusterFinder::SplitByLocalMaxima(AliRICHRawCluster *c)
 void  AliRICHClusterFinder::FillCluster(AliRICHRawCluster* c, Int_t flag) 
 {//  Completes cluster information starting from list of digits
     AliRICHDigit* dig;
-    Float_t x, y, z;
+    Float_t x, y;
     Int_t  ix, iy;
-    Float_t frac=0;
+    Float_t fraction=0;
     
     c->fPeakSignal=0;
     if (flag) {
@@ -429,7 +405,7 @@ void  AliRICHClusterFinder::FillCluster(AliRICHRawCluster* c, Int_t flag)
 
     for (Int_t i=0; i<c->fMultiplicity; i++){
        dig= (AliRICHDigit*)fDigits->UncheckedAt(c->fIndexMap[i]);
-       ix=dig->PadX()+c->fOffsetMap[i];
+       ix=dig->PadX();
        iy=dig->PadY();
        Int_t q=dig->Signal();
        if (dig->Physics() >= dig->Signal()) {
@@ -446,8 +422,8 @@ void  AliRICHClusterFinder::FillCluster(AliRICHRawCluster* c, Int_t flag)
            c->fTracks[2]=dig->Track(1);
           }
        } else {
-          if (c->fContMap[i] > frac) {
-              frac=c->fContMap[i];
+          if (c->fContMap[i] > fraction) {
+              fraction=c->fContMap[i];
              c->fPeakSignal=q;
            c->fTracks[0]=dig->Hit();
            c->fTracks[1]=dig->Track(0);
@@ -455,7 +431,7 @@ void  AliRICHClusterFinder::FillCluster(AliRICHRawCluster* c, Int_t flag)
           }
        }
        if (flag) {
-           fSegmentation->GetPadC(ix, iy, x, y, z);
+           AliRICHParam::Pad2Local(ix,iy,x,y);
            c->fX += q*x;
            c->fY += q*y;
            c->fQ += q;
@@ -471,8 +447,8 @@ void  AliRICHClusterFinder::FillCluster(AliRICHRawCluster* c, Int_t flag)
 //  apply correction to the coordinate along the anode wire
      x=c->fX;   
      y=c->fY;
-     Rich()->Param()->Local2Pad(x,y,ix,iy);
-     Rich()->Param()->Pad2Local(ix,iy,x,y);
+     AliRICHParam::Local2Pad(x,y,ix,iy);
+     AliRICHParam::Pad2Local(ix,iy,x,y);
      Int_t isec=fSegmentation->Sector(ix,iy);
      TF1* cogCorr = fSegmentation->CorrFunc(isec-1);
      
@@ -481,7 +457,7 @@ void  AliRICHClusterFinder::FillCluster(AliRICHRawCluster* c, Int_t flag)
         c->fY=c->fY-cogCorr->Eval(yOnPad, 0, 0);
      }
  }
-}//FillCluster(AliRICHRawCluster* c, Int_t flag
+}//FillCluster() 
 //__________________________________________________________________________________________________
 void  AliRICHClusterFinder::AddDigit2Cluster(Int_t i, Int_t j, AliRICHRawCluster &c)
 {//Find clusters Add i,j as element of the cluster  
@@ -525,7 +501,7 @@ void  AliRICHClusterFinder::AddDigit2Cluster(Int_t i, Int_t j, AliRICHRawCluster
        c.fMultiplicity=49;
     }
   Float_t x,y;// Prepare center of gravity calculation
-  Rich()->Param()->Pad2Local(i,j,x,y);
+  AliRICHParam::Pad2Local(i,j,x,y);
   c.fX+=q*x;    c.fY+=q*y;    c.fQ += q;
   fHitMap->FlagHit(i,j);// Flag hit as taken  
 
@@ -558,7 +534,7 @@ void AliRICHClusterFinder::FindRawClusters()
     c.fX /= c.fQ;      // center of gravity
     //c.fX=fSegmentation->GetAnod(c.fX);
     c.fY /= c.fQ;
-    AddRawCluster(c);
+    //AddRawCluster(c);
     
 //    Int_t ix,iy;//  apply correction to the coordinate along the anode wire
 //    Float_t x=c.fX, y=c.fY;  
@@ -605,7 +581,6 @@ void AliRICHClusterFinder::CalibrateCOG()
 void AliRICHClusterFinder::SinoidalFit(Float_t x, Float_t y, TF1 *func)
 {//Sinoidal fit
   static Int_t count=0;
-  Float_t z;
     
     count++;
 
@@ -614,16 +589,16 @@ void AliRICHClusterFinder::SinoidalFit(Float_t x, Float_t y, TF1 *func)
     Float_t xsig[kNs], ysig[kNs];
    
     Int_t ix,iy;
-    fSegmentation->GetPadI(x,y,0,ix,iy);   
-    fSegmentation->GetPadC(ix,iy,x,y,z);   
+    AliRICHParam::Local2Pad(x,y,ix,iy);   
+    AliRICHParam::Pad2Local(ix,iy,x,y);   
     Int_t isec=fSegmentation->Sector(ix,iy);
 // Pad Limits    
-    Float_t xmin = x-fSegmentation->Dpx(isec)/2;
-    Float_t ymin = y-fSegmentation->Dpy(isec)/2;
+    Float_t xmin = x-Rich()->Param()->PadSizeX()/2;
+    Float_t ymin = y-Rich()->Param()->PadSizeY()/2;
 //             
 //      Integration Limits
-    Float_t dxI=fResponse->SigmaIntegration()*fResponse->ChargeSpreadX();
-    Float_t dyI=fResponse->SigmaIntegration()*fResponse->ChargeSpreadY();
+    Float_t dxI=Rich()->Param()->SigmaIntegration()*Rich()->Param()->ChargeSpreadX();
+    Float_t dyI=Rich()->Param()->SigmaIntegration()*Rich()->Param()->ChargeSpreadY();
 
 //
 //  Scanning
@@ -633,7 +608,7 @@ void AliRICHClusterFinder::SinoidalFit(Float_t x, Float_t y, TF1 *func)
 
 //  y-position
     Float_t yscan=ymin;
-    Float_t dy=fSegmentation->Dpy(isec)/(kNs-1);
+    Float_t dy=Rich()->Param()->PadSizeY()/(kNs-1);
 
     for (i=0; i<kNs; i++) {//      Pad Loop
        Float_t sum=0;
@@ -650,8 +625,8 @@ void AliRICHClusterFinder::SinoidalFit(Float_t x, Float_t y, TF1 *func)
                qcheck+=qp;
                Int_t ixs=fSegmentation->Ix();
                Int_t iys=fSegmentation->Iy();
-               Float_t xs,ys,zs;
-               fSegmentation->GetPadC(ixs,iys,xs,ys,zs);
+               Float_t xs,ys;
+               AliRICHParam::Pad2Local(ixs,iys,xs,ys);
                sum+=qp*ys;
            }
        } // Pad loop
@@ -680,8 +655,8 @@ void AliRICHClusterFinder::SinoidalFit(Float_t x, Float_t y, TF1 *func)
                qcheck+=qp;
                Int_t ixs=fSegmentation->Ix();
                Int_t iys=fSegmentation->Iy();
-               Float_t xs,ys,zs;
-               fSegmentation->GetPadC(ixs,iys,xs,ys,zs);
+               Float_t xs,ys;
+               AliRICHParam::Pad2Local(ixs,iys,xs,ys);
                sum+=qp*xs;
            }
        } // Pad loop
@@ -735,25 +710,20 @@ Float_t DiscrCharge(Int_t i,Double_t *par)
        for (Int_t jbin=0; jbin<gNbins; jbin++) {
            qtot+=gCharge[jbin];
        }
-       gFirst=0;
-       //printf("\n sum of charge from DiscrCharge %f\n", qtot);
+       gFirst=kFALSE;
        gChargeTot=Int_t(qtot);
        
     }
-    gSegmentation->SetPad(gix[i], giy[i]);
-//  First Cluster
-    gSegmentation->SetHit(par[0],par[1],0);
-    Float_t q1=gResponse->IntXY(gSegmentation);
+    Float_t q1=AliRICHParam::AssignChargeToPad(par[0],par[1],gix[i],giy[i]);
     
-//  Second Cluster
-    gSegmentation->SetHit(par[2],par[3],0);
-    Float_t q2=gResponse->IntXY(gSegmentation);
+    Float_t q2=AliRICHParam::AssignChargeToPad(par[2],par[3],gix[i],giy[i]);
+//    cout<<"qtot="<<gChargeTot<<" q1="<<q1<<" q2="<<q2<<" px="<<gix[i]<<" py="<<giy[i]<<endl;
     
     Float_t value = qtot*(par[4]*q1+(1.-par[4])*q2);
     return value;
 }//DiscrCharge(Int_t i,Double_t *par) 
 //__________________________________________________________________________________________________
-void fcn(Int_t &npar, Double_t */*gin*/, Double_t &f, Double_t *par, Int_t /*iflag*/)
+void fcn(Int_t &npar, Double_t */*gin*/, Double_t &f, Double_t *par, Int_t)
 {// Minimisation function
   npar=1;
     Int_t i;
@@ -770,7 +740,7 @@ void fcn(Int_t &npar, Double_t */*gin*/, Double_t &f, Double_t *par, Int_t /*ifl
        qcont+=q1;
        qtot+=q0;
     }
-    chisq=chisq+=(qtot-qcont)*(qtot-qcont)*0.5;
+//    chisq=chisq+=(qtot-qcont)*(qtot-qcont)*0.5;
     f=chisq;
 }//
 //__________________________________________________________________________________________________
index 9506660a59daa3f046dba237aa3dcaa52145eee5..2493eed1e761110a5f6b22500e5a5d44bbc1c1cc 100644 (file)
@@ -34,8 +34,7 @@ public:
   void   AddRawCluster(const AliRICHRawCluster c)  {c.Print("");Rich()->AddClusterOld(fChamber,c);fNRawClusters++;}
   void   FillCluster(AliRICHRawCluster *cluster)   {FillCluster(cluster,1);}
   void   SetNperMax(Int_t npermax=5)               {fNperMax = npermax;}     //Set max. Number of pads per local cluster
-  void   SetDeclusterFlag(Int_t flag=1)            {fDeclusterFlag =flag;}   //Decluster ?
-  void   SetClusterSize(Int_t clsize=10)           {fClusterSize = clsize;}  //Max. cluster size; bigger clusters will be rejected
+  void   SetClusterSize(Int_t clsize=100)          {fMaxClusterSize = clsize;}//Max. cluster size; bigger clusters will be rejected
   
   AliRICH * Rich() {return fRICH;}
 //protected:
@@ -49,8 +48,7 @@ public:
   Int_t                   fChamber;                      //Chamber number
   Int_t                   fNRawClusters;                 //Number of raw clusters
   Int_t                   fNperMax;                      //Number of pad hits per local maximum
-  Int_t                   fDeclusterFlag;                //Split clusters flag
-  Int_t                   fClusterSize;                  //Size of cluster 
+  Int_t                   fMaxClusterSize;               //Max size of cluster allowed
   Int_t                   fNPeaks;                       //Number of maxima in the cluster
   ClassDef(AliRICHClusterFinder,0) //Class for clustering and reconstruction of space points    
 };
index 3227aa8a9b04597da3dba07c3f37ac72ac53910c..ee28a65dcf78c3a0b61e87271b8cd73d3be1bbff 100644 (file)
@@ -1,14 +1,12 @@
 #include "AliRICHParam.h"
 #include <TMath.h>
 #include <TRandom.h>
+#include <iostream.h> 
 ClassImp(AliRICHParam)
 
 // RICH main parameters manipulator
 //__________________________________________________________________________________________________
 AliRICHParam::AliRICHParam():
-fDeadZone(0),
-fPadSizeX(0),fPadSizeY(0),
 fCurrentPadX(0),fCurrentPadY(0),fCurrentWire(0),
 fSizeZ(0),
 fAngleRot(0),fAngleYZ(0),fAngleXY(0),
@@ -31,19 +29,9 @@ fSigmaIntegration(0),
 fAlphaFeedback(0),
 fEIonisation(0),
 fMaxAdc(0),
-fSqrtKx3(0),
-fKx2(0),
-fKx4(0),
-fSqrtKy3(0),
-fKy2(0),
-fKy4(0),
-fPitch(0),
 fWireSag(0),
 fVoltage(0)
 {//defines the default parameters
-  DeadZone             (3);                 //spacer between PC planes
-  PadSize              (0.84,0.80);     
-
   Size                 (132.6*kcm,26*kcm,136.7*kcm);  //full length, not GEANT half notation
   AngleRot             (-60);                         //rotation of the whole RICH around Z, deg
   Angles               (20,19.5);                     //XY angle, YZ angle  deg  
@@ -66,18 +54,11 @@ fVoltage(0)
   MaxAdc(4096);
   AlphaFeedback(0.036);
   EIonisation(26.e-9);
-  SqrtKx3(0.77459667);
-  Kx2(0.962);
-  Kx4(0.379);
-  SqrtKy3(0.77459667);
-  Ky2(0.962);
-  Ky4(0.379);
-  Pitch(0.25);
   WireSag(1);                // 1->On, 0->Off
   Voltage(2150);             // Should only be 2000, 2050, 2100 or 2150  
 }//AliRICHParam::named ctor 
 //__________________________________________________________________________________________________
-Int_t AliRICHParam::Local2Sector(Float_t &x, Float_t &y)const
+Int_t AliRICHParam::Local2Sector(Float_t &x, Float_t &y)
 {//Determines sector for a given hit (x,y) and trasform this point to the local system of that sector.
   Int_t sector=kBad;  
   Float_t x1=-0.5*PcSizeX();      Float_t x2=-0.5*SectorSizeX()-DeadZone();  Float_t x3=-0.5*SectorSizeX();
@@ -86,29 +67,29 @@ Int_t AliRICHParam::Local2Sector(Float_t &x, Float_t &y)const
   if     (x>=x1&&x<=x2)    {sector=1;x+=0.5*PcSizeX();}
   else if(x>=x3&&x<=x4)    {sector=2;x+=0.5*SectorSizeX();}
   else if(x>=x5&&x<=x6)    {sector=3;x-=0.5*SectorSizeX()+DeadZone();}
-  else if(x< x1||x> x6)    {Error("Sector","given x position is out of PC area");return kBad;}
+  else if(x< x1||x> x6)    {return kBad;}
   else                                                        {return kBad;} //in dead zone
 
   if     (y>=-0.5*PcSizeY()   &&y<=-0.5*DeadZone())  {y+=0.5*PcSizeY();  return -sector;}
   else if(y> -0.5*DeadZone()  &&y<  0.5*DeadZone())  {return kBad;} //in dead zone
   else if(y>= 0.5*DeadZone()  &&y<= 0.5*PcSizeY())   {y-=0.5*DeadZone(); return  sector;}
-  else                                            {Error("Sector","given y position is out of PC area");return kBad;}
+  else                                               {return kBad;}
 }//Int_t AliRICHParam::Local2Sector(Float_t x, Float_t y)
 //__________________________________________________________________________________________________
-Int_t AliRICHParam::Pad2Sector(Int_t &padx, Int_t &pady)const
+Int_t AliRICHParam::Pad2Sector(Int_t &padx, Int_t &pady)
 {//Determines sector for a given pad (padx,pady) and trasform this point to the local system of that sector.
   Int_t sector=kBad;      
   if     (padx>=1            &&padx<=NpadsXsec())      {sector=1;}
   else if(padx> NpadsXsec()  &&padx<=NpadsXsec()*2)    {sector=2;padx-=NpadsXsec();}
   else if(padx> NpadsXsec()*2&&padx<=NpadsX())         {sector=3;padx-=NpadsXsec()*2;}
-  else                                   {Error("Sector","given padx position is out of PC area");return kBad;}
+  else                                                 {return kBad;}
 
   if     (pady>=1         &&pady<= NpadsYsec())     {return -sector;}
   else if(pady>NpadsYsec()&&pady<= NpadsY())        {pady-=NpadsYsec();return sector;} 
-  else                                              {Error("Sector","given y position is out of PC area");return kBad;}
+  else                                              {return kBad;}
 }//Local2Sector()
 //__________________________________________________________________________________________________
-Int_t AliRICHParam::Local2Pad(Float_t x, Float_t y, Int_t &padx, Int_t &pady)const
+Int_t AliRICHParam::Local2Pad(Float_t x, Float_t y, Int_t &padx, Int_t &pady)
 {//returns pad numbers (iPadX,iPadY) for given point in local coordinates (x,y) 
  //count starts in lower left corner from 1,1 to 144,180
   
@@ -175,3 +156,25 @@ void AliRICHParam::FirstPad(Float_t x,Float_t y)
   Int_t padx,pady;
   Local2Pad(x,y,padx,pady);
 }//void AliRICHParam::FirstPad(Float_t x,Float_t y)
+//__________________________________________________________________________________________________
+Float_t AliRICHParam::AssignChargeToPad(Float_t hitx,Float_t hity,Int_t padx,Int_t pady)
+{//
+  Float_t padXcenter=0,padYcenter=0;
+  Pad2Local(padx,pady,padXcenter,padYcenter);
+  
+  Float_t xi1=hitx-padXcenter-PadSizeX()/2;
+  Float_t xi2=hitx-padXcenter+PadSizeX()/2; 
+  Float_t yi1=hity-padYcenter-PadSizeY()/2;
+  Float_t yi2=hity-padYcenter+PadSizeY()/2;
+  xi1/=AnodeCathodeGap();
+  xi2/=AnodeCathodeGap();
+  yi1/=AnodeCathodeGap();
+  yi2/=AnodeCathodeGap();
+// The Mathieson function 
+  Double_t ux1=SqrtKx3()*TMath::TanH(Kx2()*xi1);
+  Double_t ux2=SqrtKx3()*TMath::TanH(Kx2()*xi2);    
+  Double_t uy1=SqrtKy3()*TMath::TanH(Ky2()*yi1);
+  Double_t uy2=SqrtKy3()*TMath::TanH(Ky2()*yi2);
+  return 4.*Kx4()*(TMath::ATan(ux2)-TMath::ATan(ux1))*Ky4()*(TMath::ATan(uy2)-TMath::ATan(uy1));
+}//AssignChargeToPad()
+//__________________________________________________________________________________________________
index 3baf5aaf65f1ad6f39cef543fc98bc51a10cd512..9c1f3960b358a134ce06015776acb66a4da4d87f 100644 (file)
@@ -13,32 +13,32 @@ public:
   inline  Int_t Neighbours(Int_t iPadX,Int_t iPadY,Int_t aListX[4],Int_t aListY[4])const;                      //pad->neibours
   inline  void   SigGenInit(Float_t x,Float_t y);
   inline  Bool_t SigGenCond(Float_t x,Float_t y);
-  Int_t   Local2Pad(Float_t x,Float_t y,Int_t &padx,Int_t &pady)const;               //(x,y)->(padx,pady), returns sector code 
-  Int_t   Local2PadX(Float_t x,Float_t y)    const {Int_t padx,pady;Local2Pad(x,y,padx,pady);return padx;}//(x,y)->padx
-  Int_t   Local2PadY(Float_t x,Float_t y)    const {Int_t padx,pady;Local2Pad(x,y,padx,pady);return pady;}//(x,y)->pady
-  void    Pad2Local(Int_t padx,Int_t pady,Float_t &x,Float_t &y);                                         //(padx,pady)->(x,y)
-  Int_t   LocalX2Wire(Float_t x)             const {return  Int_t((x+PcSizeX()/2)/fWirePitch)+1;}         //x->wire number
-  Float_t Wire2LocalX(Int_t iWireN)          const {return iWireN*fWirePitch-PcSizeX()/2;}                //wire number->x
+  static  Int_t   Local2Pad(Float_t x,Float_t y,Int_t &padx,Int_t &pady);               //(x,y)->(padx,pady), returns sector code 
+  static  Int_t   Local2PadX(Float_t x,Float_t y)    {Int_t padx,pady;Local2Pad(x,y,padx,pady);return padx;}//(x,y)->padx
+  static  Int_t   Local2PadY(Float_t x,Float_t y)    {Int_t padx,pady;Local2Pad(x,y,padx,pady);return pady;}//(x,y)->pady
+  static  void    Pad2Local(Int_t padx,Int_t pady,Float_t &x,Float_t &y);                                         //(padx,pady)->(x,y)
+  static  Int_t   LocalX2Wire(Float_t x)      {return  Int_t((x+PcSizeX()/2)/WirePitch())+1;}         //x->wire number
+  static  Float_t Wire2LocalX(Int_t iWireN)   {return iWireN*WirePitch()-PcSizeX()/2;}                //wire number->x
   
   Float_t Gain(Float_t y);                                 //Returns total charge induced by single photon
   Float_t TotalCharge(Int_t iPID,Float_t eloss,Float_t y); //Returns total charge induced by particle lost eloss GeV
-  Float_t PadCharge(Int_t /* iPadX */,Int_t /* iPadY */) {return 0;}   //Returns charge for a given pad
+  static  Float_t AssignChargeToPad(Float_t hx,Float_t hy, Int_t px, Int_t py); //Returns charge assigned to given pad for a given hit
   void    FirstPad(Float_t x,Float_t y);
           
+  static  Float_t AnodeCathodeGap()          {return 0.2;}
+  
   static  Int_t   NpadsX()                   {return 144;}
   static  Int_t   NpadsY()                   {return 160;}   
   static  Int_t   NpadsXsec()                {return NpadsX()/3;}   
   static  Int_t   NpadsYsec()                {return NpadsY()/2;}   
-  void    DeadZone(Float_t a)                {       fDeadZone=a;}
-  Float_t DeadZone()                    const{return fDeadZone;}
-  void    PadSize(Float_t x,Float_t y)       {       fPadSizeX=x;fPadSizeY=y;} 
-  Float_t PadSizeX()                    const{return fPadSizeX;}
-  Float_t PadSizeY()                    const{return fPadSizeY;}
-  Float_t SectorSizeX()                 const{return NpadsX()*PadSizeX()/3;}
-  Float_t SectorSizeY()                 const{return NpadsY()*PadSizeY()/2;}  
-  Float_t PcSizeX()                     const{return NpadsX()*PadSizeX()+2*DeadZone();}
-  Float_t PcSizeY()                     const{return NpadsY()*PadSizeY()+DeadZone();}
-  Float_t WirePitch()                   const{return PadSizeX()/2;}
+  static  Float_t DeadZone()                 {return 2.6;}
+  static  Float_t PadSizeX()                 {return 0.84;}
+  static  Float_t PadSizeY()                 {return 0.8;}
+  static  Float_t SectorSizeX()              {return NpadsX()*PadSizeX()/3;}
+  static  Float_t SectorSizeY()              {return NpadsY()*PadSizeY()/2;}  
+  static  Float_t PcSizeX()                  {return NpadsX()*PadSizeX()+2*DeadZone();}
+  static  Float_t PcSizeY()                  {return NpadsY()*PadSizeY()+DeadZone();}
+  static  Float_t WirePitch()                {return PadSizeX()/2;}
             
   void    Size(Float_t x,Float_t y,Float_t z){fSizeX=x;fSizeY=y;fSizeZ=z;}
   void    GeantSize(Float_t *pArr)      const{pArr[0]=fSizeX/2;pArr[1]=fSizeY/2;pArr[2]=fSizeZ/2;}  
@@ -87,28 +87,23 @@ public:
   Float_t ChargeSlope()                      {return fChargeSlope;}
   void    MaxAdc(Int_t a)                    {       fMaxAdc=a;}
   Int_t   MaxAdc()                      const{return fMaxAdc;}
-  void    Pitch(Float_t a)                   {       fPitch=a;}
-  Float_t Pitch()                       const{return fPitch;}
   void    AlphaFeedback(Float_t a)           {       fAlphaFeedback=a;}
   Float_t AlphaFeedback()               const{return fAlphaFeedback;}
   void    EIonisation(Float_t a)             {       fEIonisation=a;}
   Float_t EIonisation()                 const{return fEIonisation;}                            
-  void    SqrtKx3(Float_t a)                 {       fSqrtKx3=a;};
-  void    Kx2(Float_t a)                     {       fKx2=a;}
-  void    Kx4(Float_t a)                     {       fKx4=a;}
-  void    SqrtKy3(Float_t a)                 {       fSqrtKy3=a;}
-  void    Ky2(Float_t a)                     {       fKy2=a;}
-  void    Ky4(Float_t a)                     {       fKy4=a;}
+  static Float_t SqrtKx3()  {return 0.77459667;}
+  static Float_t Kx2()      {return 0.962;}
+  static Float_t Kx4()      {return 0.379;}
+  static Float_t SqrtKy3()  {return 0.77459667;}
+  static Float_t Ky2()      {return 0.962;}
+  static Float_t Ky4()      {return 0.379;}
+
   void    WireSag(Int_t a)                   {       fWireSag=a;}
   void    Voltage(Int_t a)                   {       fVoltage=a;}       
   Float_t Voltage()                     const{return fVoltage;}       
 protected:
-  Int_t   Local2Sector(Float_t &x,Float_t &y)const; //(x,y)->sector
-  Int_t   Pad2Sector(Int_t &padx,Int_t &pady)const; //(padx,pady)->sector
-  
-  Float_t fDeadZone;                                //space between PC sectors, cm     
-  Float_t fPadSizeX,fPadSizeY;                      //pad size, cm
-  Float_t fWirePitch;                               //distance between wires along x
+  static Int_t   Local2Sector(Float_t &x,Float_t &y); //(x,y)->sector
+  static Int_t   Pad2Sector(Int_t &padx,Int_t &pady); //(padx,pady)->sector
   
   Int_t   fCurrentPadX,fCurrentPadY;              //???
   Int_t   fCurrentWire;                           //???
@@ -133,30 +128,23 @@ protected:
   Float_t fAlphaFeedback;            //Feedback photons coefficient
   Float_t fEIonisation;              //Mean ionisation energy
   Int_t   fMaxAdc;                   //Maximum ADC channel
-  Float_t fSqrtKx3;                  //Mathieson parameters for x
-  Float_t fKx2;                      //Mathieson parameters for x
-  Float_t fKx4;                      //Mathieson parameters for x
-  Float_t fSqrtKy3;                  //Mathieson parameters for y
-  Float_t fKy2;                      //Mathieson parameters for y 
-  Float_t fKy4;                      //Mathieson parameters for y
-  Float_t fPitch;                    //Anode-cathode pitch
   Int_t   fWireSag;                  //Flag to turn on/off (0/1) wire sag
   Int_t   fVoltage;                  //Working voltage (2000, 2050, 2100, 2150)
 
-  ClassDef(AliRICHParam,1)    //RICH main parameters
+  ClassDef(AliRICHParam,2)    //RICH main parameters
 };
 //__________________________________________________________________________________________________
 void AliRICHParam::SigGenInit(Float_t x,Float_t y)
 {//Initialises pad and wire position during stepping
   Local2Pad(x,y,fCurrentPadX,fCurrentPadY);
-  fCurrentWire= (x>0) ? Int_t(x/fWirePitch)+1 : Int_t(x/fWirePitch)-1 ;
+  fCurrentWire= (x>0) ? Int_t(x/WirePitch())+1 : Int_t(x/WirePitch())-1 ;
 }
 //__________________________________________________________________________________________________
 Bool_t AliRICHParam::SigGenCond(Float_t x,Float_t y)
 {//Signal will be generated if particle crosses pad boundary or boundary between two wires.
   Int_t curPadX,curPadY;
   Local2Pad(x,y,curPadX,curPadY);
-  Int_t currentWire=(x>0) ? Int_t(x/fWirePitch)+1 : Int_t(x/fWirePitch)-1;
+  Int_t currentWire=(x>0) ? Int_t(x/WirePitch())+1 : Int_t(x/WirePitch())-1;
   if((curPadX != fCurrentPadX) || (curPadY != fCurrentPadY) || (currentWire!=fCurrentWire)) 
     return kTRUE;
   else
index 8fd4e1b8605528279e74019266ad843f6f38c115..56516d76e7651df2f8d4ca66ec6ff437b169c441 100644 (file)
@@ -95,7 +95,7 @@ void AliRICHPatRec::PatRec()
       x=pDig->PadX();
       y=pDig->PadY();
       q=pDig->Signal();
-      pRICH->Param()->Pad2Local(x,y,rx,ry);      
+      AliRICHParam::Pad2Local(x,y,rx,ry);      
 
       //printf("Pad coordinates x:%d, Real coordinates x:%f\n",x,rx);
       //printf("Pad coordinates y:%d, Real coordinates y:%f\n",y,ry);
@@ -112,7 +112,7 @@ void AliRICHPatRec::PatRec()
       Int_t xpad;
       Int_t ypad;
 
-      pRICH->Param()->Local2Pad(fXpad,fYpad,xpad,ypad);
+      AliRICHParam::Local2Pad(fXpad,fYpad,xpad,ypad);
 
       padsUsedX[goodPhotons]=xpad;
       padsUsedY[goodPhotons]=ypad;
index d0f646d35be3a92b03b754e336f340f3848bceec..40b12201328991d10f1dc390f0c40c932a31d618 100644 (file)
@@ -14,6 +14,8 @@
  **************************************************************************/
 
 #include "AliRICHRawCluster.h"
+#include <iostream.h>
+
 
 ClassImp(AliRICHRawCluster)
 //__________________________________________________________________________________________________
@@ -21,10 +23,10 @@ ClassImp(AliRICHRawCluster)
 AliRICHRawCluster :: AliRICHRawCluster() 
 {//default ctor
   fTracks[0]=fTracks[1]=fTracks[2]=-1; 
-  fQ=0; fX=fY=0; fMultiplicity=0;
+  fX=fY=0;
+  fQ=fMultiplicity=0;
   for (int k=0;k<50;k++) {
     fIndexMap[k]=-1;
-    fOffsetMap[k]=0;
     fContMap[k]=0;
     fPhysicsMap[k]=-1;
     fCtype=-1;
@@ -59,6 +61,9 @@ Int_t AliRICHRawCluster::PhysicsContribution()
 //__________________________________________________________________________________________________
 void AliRICHRawCluster::Print(Option_t*)const
 {
-  Info("","X=%7.2f, Y=%7.2f, Qdc=%4i, Peak=%4i, Multip=%2i,       T0=%5i T0=%5i T0=%5i",
-           fX,      fY,      fQ,      fPeakSignal,fMultiplicity,    fTracks[0], fTracks[1], fTracks[2]);
+  Info("","X=%7.2f, Y=%7.2f, Qdc=%4i, Peak=%4i, Multip=%2i,      T0=%5i T1=%5i T2=%5i",
+           fX,      fY,      fQ,      fPeakSignal,fMultiplicity, fTracks[0], fTracks[1], fTracks[2]);
+  
+  for(int i=0;i<fMultiplicity;i++)
+    cout<<"D"<<i<<"="<<fIndexMap[i]<<" C="<<fContMap[i]<<endl;
 }//void AliRICHRawCluster::Print(Option_t *option)const
index 5c46b7033a27b31419307dcee94e2bd72fd76b66..ce632adaac198ba2b37ebf3677795bc742ff15b5 100644 (file)
@@ -15,7 +15,7 @@ public:
   virtual ~AliRICHRawCluster(){;}
   Bool_t IsSortable() const {return kTRUE;}  //virtual
   Int_t  Compare(const TObject *obj) const;  //virtual
-  void   Print(const Option_t *option)const; //virtual  
+  void   Print(const Option_t *option="")const; //virtual  
   Int_t  PhysicsContribution();
 public:
   Int_t       fTracks[3];      //labels of overlapped tracks
@@ -24,13 +24,12 @@ public:
   Float_t     fY  ;            // Y of cluster
   Int_t       fPeakSignal;     // Charge in the peak
   Int_t       fIndexMap[50];   //indeces of digits
-  Int_t       fOffsetMap[50];  // offset map
   Float_t     fContMap[50];    //Contribution from digit
   Int_t       fPhysicsMap[50]; // physics processes
   Int_t       fMultiplicity;   //cluster multiplicity
   Int_t       fNcluster[2];    //number of clusters
   Int_t       fClusterType;    //??
   Int_t       fCtype;          //CL0, CL1, etc...    
-  ClassDef(AliRICHRawCluster,1)  //Cluster object for set:RICH
+  ClassDef(AliRICHRawCluster,2)  //Cluster object for set:RICH
 };
 #endif
index 4a1266f70d6336899918044b08ddbbe7a021a5b0..d22a9e41b359f27e1bd1d44e4abdf4b6ed149c45 100644 (file)
@@ -515,14 +515,14 @@ void TestDigitsOLD()
   d[0]=2;d[1]=2;d[2]=10;d[3]=3;d[4]=4;r->AddDigitOld(2,t,c,d); 
   d[0]=2;d[1]=3;d[2]=10;d[3]=3;d[4]=4;r->AddDigitOld(2,t,c,d);
   
-  d[0]=2;d[1]=2;d[2]=50;d[3]=3;d[4]=4;r->AddDigitOld(3,t,c,d); 
-  d[0]=2;d[1]=3;d[2]=100;d[3]=3;d[4]=4;r->AddDigitOld(3,t,c,d);
-  d[0]=3;d[1]=2;d[2]=100;d[3]=3;d[4]=4;r->AddDigitOld(3,t,c,d);
+  d[0]=2;d[1]=2;d[2]=100;d[3]=3;d[4]=4;r->AddDigitOld(3,t,c,d); 
+  d[0]=2;d[1]=3;d[2]= 50;d[3]=3;d[4]=4;r->AddDigitOld(3,t,c,d);
+  d[0]=2;d[1]=4;d[2]=200;d[3]=3;d[4]=4;r->AddDigitOld(3,t,c,d);
     
   d[0]=2;d[1]=2;d[2]=100;d[3]=3;d[4]=4;r->AddDigitOld(4,t,c,d); 
-  d[0]=2;d[1]=3;d[2]=50;d[3]=3;d[4]=4;r->AddDigitOld(4,t,c,d);
-  d[0]=3;d[1]=2;d[2]=50;d[3]=3;d[4]=4;r->AddDigitOld(4,t,c,d);
-  d[0]=3;d[1]=3;d[2]=100;d[3]=3;d[4]=4;r->AddDigitOld(4,t,c,d);
+  d[0]=2;d[1]=3;d[2]= 50;d[3]=3;d[4]=4;r->AddDigitOld(4,t,c,d);
+  d[0]=2;d[1]=4;d[2]=200;d[3]=3;d[4]=4;r->AddDigitOld(4,t,c,d);
+  d[0]=2;d[1]=5;d[2]= 50;d[3]=3;d[4]=4;r->AddDigitOld(4,t,c,d);
   
   rl->TreeD()->Fill();
   rl->WriteDigits("OVERWRITE");