Raw data reconstruction (K.Shileev)
authorhristov <hristov@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 8 Aug 2005 11:26:33 +0000 (11:26 +0000)
committerhristov <hristov@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 8 Aug 2005 11:26:33 +0000 (11:26 +0000)
23 files changed:
RICH/AliGenRadioactive.cxx [deleted file]
RICH/AliGenRadioactive.h [deleted file]
RICH/AliRICH.cxx
RICH/AliRICH.h
RICH/AliRICHCluster.cxx
RICH/AliRICHCluster.h
RICH/AliRICHClusterFinder.cxx
RICH/AliRICHDigit.cxx
RICH/AliRICHDigit.h
RICH/AliRICHDigitizer.cxx
RICH/AliRICHMap.cxx
RICH/AliRICHParam.h
RICH/AliRICHReconstructor.cxx
RICH/AliRICHReconstructor.h
RICH/AliRICHv1.cxx
RICH/AliRICHv1.h
RICH/RICHbaseLinkDef.h
RICH/RICHsimLinkDef.h
RICH/RichConfig.C
RICH/RichMake
RICH/api.txt
RICH/libRICHbase.pkg
RICH/libRICHsim.pkg

diff --git a/RICH/AliGenRadioactive.cxx b/RICH/AliGenRadioactive.cxx
deleted file mode 100644 (file)
index 1244b35..0000000
+++ /dev/null
@@ -1,45 +0,0 @@
-#include "AliGenRadioactive.h"
-#include <TPDGCode.h>
-#include <AliRun.h>
-#include <TDatabasePDG.h>
-#include <AliLog.h>
-
-//__________________________________________________________________________________________________
-AliGenRadioactive::AliGenRadioactive(Int_t iSrcNucleus,Int_t iNsecondaries):AliGenerator()
-{
-//Main ctor. Used to define radioactive source.
-  Double_t e[100],a[100];//arrays to store experimental points
-  Int_t nPoints; 
-  switch(iSrcNucleus){
-    case kSr90: fPartId=kElectron; nPoints=46;  //experimental part
-  a[ 0]=0.08605; a[ 1]=0.0878;  a[ 2]=0.08705; a[ 3]=0.07855; a[ 4]=0.0709; a[ 5]=0.0647;  a[ 6]=0.05015; a[ 7]=0.0372; a[ 8]=0.0268;  a[ 9]=0.0215;
-  a[10]=0.0157;  a[11]=0.01685; a[12]=0.01745; a[13]=0.01645; a[14]=0.0175; a[15]=0.01635; a[16]=0.01825; a[17]=0.0177; a[18]=0.01735; a[19]=0.0161;
-  a[20]=0.0159;  a[21]=0.0176;  a[22]=0.01605; a[23]=0.0161;  a[24]=0.01495;a[25]=0.01595; a[26]=0.01525; a[27]=0.0138; a[28]=0.0121;  a[29]=0.0101;
-  a[30]=0.01175; a[31]=0.01095; a[32]=0.0089;  a[33]=0.0091;  a[34]=0.0625; a[35]=0.0505;  a[36]=0.0475;  a[37]=0.0039; a[38]=0.0031;  a[39]=0.0028;
-  a[40]=0.0025;  a[41]=0.0017;  a[42]=4.5e-4;  a[43]=4.5e-4;  a[44]=1.5e-4; a[45]=0;     
-  break;
-    default:    AliError("Wrong source nucleus specified");  return;   
-  }
-  for(Int_t i=0;i<nPoints;i++) e[i]=0.001*(i*0.05+0.025); //kinetic energy GeV
-  fGenH1=new TH1F("Sr90","Sr90 generator hist",nPoints-1,0,e[nPoints-1]);  
-  for(Int_t i=0;i<nPoints;i++) fGenH1->Fill(e[i],a[i]);
-  fNpart=iNsecondaries;
-}
-//__________________________________________________________________________________________________
-void AliGenRadioactive::Generate()
-{
-// Generate one trigger
-  Int_t nt=0;
-  Double_t ekin=0,p=0,theta=0,phi=0,x=0,y=0,z=0,px=0,py=0,pz=0,polx=0,poly=0,polz=0;
-  Double_t m=gAlice->PDGDB()->GetParticle(fPartId)->Mass();
-  for(Int_t i=0;i<fNpart;i++){
-    x=fOrigin.At(0)+fOsigma.At(0)*(Rndm()-0.5); y=fOrigin.At(1)+fOsigma.At(1)*(Rndm()-0.5); z=fOrigin.At(2)+fOsigma.At(2)*(Rndm()-0.5);
-    ekin=fGenH1->GetRandom();  p=TMath::Sqrt(ekin*(2*m+ekin)); 
-    theta=Rndm()*fThetaMax; phi=Rndm()*fPhiMax;
-    px=p*TMath::Cos(theta)*TMath::Cos(phi);  py=p*TMath::Cos(theta)*TMath::Sin(phi);  pz=p*TMath::Sin(theta);
-//   AliDebug(1,Form("Origin=(%5.2f,%5.2f,%5.2f) Ekin=%5.3fMeV,P=%5.3fMeV (%5.2f,%5.2f,%5.2f)",x,y,z,1000*ekin,1000*p,1000*px,1000*py,1000*pz));
-    PushTrack(fTrackIt,-1,fPartId,px,py,pz,ekin+m,
-                                   x, y, z,0,
-                                    polx=0,poly=0,polz=0,kPPrimary,nt);//cm, GeV
-  }
-}
diff --git a/RICH/AliGenRadioactive.h b/RICH/AliGenRadioactive.h
deleted file mode 100644 (file)
index 0a257f4..0000000
+++ /dev/null
@@ -1,28 +0,0 @@
-#ifndef AliGenRadioactive_h
-#define AliGenRadioactive_h
-/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
- * See cxx source for full Copyright notice                               */
-
-// Simple radioactive source generator.
-// Currently Sr90 only
-
-#include <AliGenerator.h>
-#include <TH1.h>
-
-typedef enum {kSr90} RadioNuclei_t;     
-
-class AliGenRadioactive : public AliGenerator
-{
- public:
-           AliGenRadioactive() : AliGenerator(),fPartId(0),fGenH1(0) {;}//default ctor
-           AliGenRadioactive(Int_t iSource,Int_t iNparts);                  //main ctor
-  virtual ~AliGenRadioactive()                                           {if(fGenH1) delete fGenH1;}//dtor
-  virtual void Generate();                                                              //interface from AliGenerator, generates current event 
-
-
-protected:
-  Int_t    fPartId;          //ID of secondary from radioactive nucleus
-  TH1F    *fGenH1;           //Histogram containg exp spectrum to be used to sample the energy
-  ClassDef(AliGenRadioactive,1) // Radioactive source generator
-};
-#endif
index ff38f6f..a73a87a 100644 (file)
 #include <AliLog.h>
 #include <TNtupleD.h>
 #include <AliTracker.h>
-#include <AliRawDataHeader.h>
-#include <TLatex.h> //Display()
-#include <TCanvas.h> //Display()
-#include <TGraph.h> //Display()
-#include <TStyle.h> //Display()
-#include <TMarker.h> //Display()
+#include <TLatex.h>   //Display()
+#include <TCanvas.h>  //Display()
+#include <TGraph.h>   //Display()
+#include <TStyle.h>   //Display()
+#include <TMarker.h>  //Display()
  
 ClassImp(AliRICH)    
 //__________________________________________________________________________________________________
@@ -93,98 +92,6 @@ AliRICH::~AliRICH()
   AliDebug(1,"Stop.");    
 }//AliRICH::~AliRICH()
 //__________________________________________________________________________________________________
-void AliRICH::Hits2SDigits()
-{
-// Create a list of sdigits corresponding to list of hits. Every hit generates one or more sdigits.
-  AliDebug(1,"Start.");
-  for(Int_t iEventN=0;iEventN<GetLoader()->GetRunLoader()->GetAliRun()->GetEventsPerRun();iEventN++){//events loop
-    GetLoader()->GetRunLoader()->GetEvent(iEventN);//get next event
-  
-    if(!GetLoader()->TreeH()) GetLoader()->LoadHits();    GetLoader()->GetRunLoader()->LoadHeader(); 
-    if(!GetLoader()->GetRunLoader()->TreeK())             GetLoader()->GetRunLoader()->LoadKinematics();//from
-    if(!GetLoader()->TreeS()) GetLoader()->MakeTree("S"); MakeBranch("S");//to
-          
-    for(Int_t iPrimN=0;iPrimN<GetLoader()->TreeH()->GetEntries();iPrimN++){//prims loop
-      GetLoader()->TreeH()->GetEntry(iPrimN);
-      for(Int_t iHitN=0;iHitN<Hits()->GetEntries();iHitN++){//hits loop 
-        AliRICHHit *pHit=(AliRICHHit*)Hits()->At(iHitN);//get current hit                
-        TVector2 x2 = C(pHit->C())->Mrs2Anod(0.5*(pHit->InX3()+pHit->OutX3()));//hit position in the anod plane
-        Int_t iTotQdc=P()->TotQdc(x2,pHit->Eloss());//total charge produced by hit, 0 if hit in dead zone
-        if(iTotQdc==0) continue;
-        //
-        //need to quantize the anod....
-        TVector padHit=AliRICHParam::Loc2Pad(x2);
-        TVector2 padHitXY=AliRICHParam::Pad2Loc(padHit);
-        TVector2 anod;
-        if((x2.Y()-padHitXY.Y())>0) anod.Set(x2.X(),padHitXY.Y()+AliRICHParam::PitchAnod()/2);
-        else anod.Set(x2.X(),padHitXY.Y()-AliRICHParam::PitchAnod()/2);
-        //end to quantize anod
-        //
-        TVector area=P()->Loc2Area(anod);//determine affected pads, dead zones analysed inside
-        AliDebug(1,Form("hitanod(%6.2f,%6.2f)->area(%3.0f,%3.0f)-(%3.0f,%3.0f) QDC=%4i",anod.X(),anod.Y(),area[0],area[1],area[2],area[3],iTotQdc));
-        TVector pad(2);
-        for(pad[1]=area[1];pad[1]<=area[3];pad[1]++)//affected pads loop
-          for(pad[0]=area[0];pad[0]<=area[2];pad[0]++){                    
-            Double_t padQdc=iTotQdc*P()->FracQdc(anod,pad);
-            AliDebug(1,Form("current pad(%3.0f,%3.0f) with QDC  =%6.2f",pad[0],pad[1],padQdc));
-            if(padQdc>0.1) SDigitAdd(pHit->C(),pad,padQdc,GetLoader()->GetRunLoader()->Stack()->Particle(pHit->GetTrack())->GetPdgCode(),pHit->GetTrack());
-          }//affected pads loop 
-      }//hits loop
-    }//prims loop
-    GetLoader()->TreeS()->Fill();
-    GetLoader()->WriteSDigits("OVERWRITE");
-    SDigitsReset();
-  }//events loop  
-  GetLoader()->UnloadHits(); GetLoader()->GetRunLoader()->UnloadHeader(); GetLoader()->GetRunLoader()->UnloadKinematics();
-  GetLoader()->UnloadSDigits();  
-  AliDebug(1,"Stop.");
-}//Hits2SDigits()
-//__________________________________________________________________________________________________
-void AliRICH::Digits2Raw()
-{
-//Loops over all digits and creates raw data files in DDL format. GetEvent() is done outside (AliSimulation)
-//RICH has 2 DDL per chamber, even number for right part(2-4-6) odd number for left part(1-3-5) 
-//RICH has no any propriate header just uses the common one
-//Each PC is divided by 8 rows counted from 1 to 8 from top to bottom for left PCs(1-3-5) and from bottom to top for right PCc(2-4-6)     (denoted  rrrrr 5 bits 32 values)
-//Each raw is composed from 10 DILOGIC chips counted from left to right from 1 to 10                                                      (denoted   dddd 4 bits 16 values)
-//Each DILOGIC chip serves 48 channels counted from 0 to 47                                                                               (denoted aaaaaa 6 bits 64 values)
-//So RICH info word is  32 bits word with structure:   0000 0rrr rrdd ddaa aaaa qqqq qqqq qqqq (five 0 five r six a twelve q) with QDC    (denoted q...q 12 bits 4096 values)
-  AliDebug(1,"Start.");
-  GetLoader()->LoadDigits();
-  GetLoader()->TreeD()->GetEntry(0);
-  
-  Int_t kRichOffset=0x700; //currently one DDL per 3 sectors
-  
-  ofstream file135,file246;//output streams 2 DDL per chamber
-  AliRawDataHeader header;//empty DDL miniheader
-  UInt_t word32=1;        //32 bits data word 
-  
-  for(Int_t iChN=1;iChN<=kNchambers;iChN++){ //2 DDL per chamber open both in parallel   
-    file135.open(Form("RICH_%4i.ddl",kRichOffset+2*iChN-1));   //left part of chamber; sectors 1-3-5 odd DDL number
-    file246.open(Form("RICH_%4i.ddl",kRichOffset+2*iChN));     //right part of chamber; sectors 2-4-6 even DDL number
-//common DDL header defined by standart, now dummy as the total number of bytes is not yet known    
-    file135.write((char*)&header,sizeof(header)); //dummy header just place holder
-    file246.write((char*)&header,sizeof(header)); //actual will be written later
-    
-    Int_t counter135=0,counter246=0;//counts total number of records per DDL 
-    
-    for(Int_t iDigN=0;iDigN<Digits(iChN)->GetEntriesFast();iDigN++){//digits loop for a given chamber
-      AliRICHDigit *pDig=(AliRICHDigit*)Digits(iChN)->At(iDigN);
-      word32=UInt_t (pDig->Q()+0x400*pDig->X()+0x4000*pDig->Y());  //arbitrary structure
-      switch(pDig->S()){//readout by vertical sectors: 1,3,5-left DDL 2,4,6-right DDL
-        case 1: case 3: case 5: file135.write((char*)&word32,sizeof(word32)); counter135++; break;
-        case 2: case 4: case 6: file246.write((char*)&word32,sizeof(word32)); counter246++; break;
-      }//switch sector  
-    }//digits loop for a given chamber
-//now count total byte number for each DDL file and rewrite actual header    
-    header.fSize=sizeof(header)+counter135*sizeof(word32);   header.SetAttribute(0); file135.seekp(0); file135.write((char*)&header,sizeof(header));
-    header.fSize=sizeof(header)+counter246*sizeof(word32);   header.SetAttribute(0); file246.seekp(0); file246.write((char*)&header,sizeof(header));
-    file135.close(); file246.close();
-  }//chambers loop  
-  GetLoader()->UnloadDigits();
-  AliDebug(1,"Stop.");      
-}//Digits2Raw()
-//__________________________________________________________________________________________________
 void AliRICH::BuildGeometry() 
 {
 //Builds a TNode geometry for event display
@@ -1007,7 +914,7 @@ void AliRICH::DisplayEvent(Int_t iEvtNmin,Int_t iEvtNmax)const
       for(Int_t j=0;j<Digits(iChamber)->GetEntries();j++) {//digits loop
         AliRICHDigit *pDig = (AliRICHDigit*)Digits(iChamber)->At(j);
         TVector2 x2=AliRICHParam::Pad2Loc(pDig->Pad());
-        pDigitsH2[iChamber]->Fill(x2.X(),x2.Y(),pDig->Q());
+        pDigitsH2[iChamber]->Fill(x2.X(),x2.Y(),pDig->Qdc());
       }//digits loop
       if(iChamber==1) canvas->cd(7);
       if(iChamber==2) canvas->cd(8);
index 9b3f252..10cc584 100644 (file)
@@ -3,13 +3,12 @@
 /* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
  * See cxx source for full Copyright notice                               */
 
-#include <AliDetector.h>
-#include <TClonesArray.h>
+#include <AliDetector.h>  //inheritance
+#include <TClonesArray.h> 
 #include <TObjArray.h>
 #include <TVector.h>
 #include <TVector3.h>
 
-#include "AliRICHDigitizer.h"
 #include "AliRICHParam.h"
 #include "AliRICHCluster.h"
 #include "AliRICHHit.h"
@@ -30,9 +29,6 @@ public:
 //framework part  
   virtual Int_t         IsVersion()                           const =0;                                  //interface from         
   virtual void          StepManager()                               =0;                                  //interface from AliMC
-  virtual void          Hits2SDigits();                                                                  //interface from AliSimulation
-  virtual void          Digits2Raw();                                                                    //interface from AliSimulation
-  virtual AliDigitizer* CreateDigitizer(AliRunDigitizer* man) const {return new AliRICHDigitizer(man);}  //interface from AliSimulation
   virtual void          SetTreeAddress();                                                                //interface from AliLoader
   virtual void          MakeBranch(Option_t *opt=" ");                                                   //interface from AliLoader
   virtual void          CreateMaterials();                                                               //interface from AliMC
@@ -63,6 +59,7 @@ public:
     using AliDetector::Digits;  
          TClonesArray* Digits        (Int_t iC                                               )const{return fDigs ? (TClonesArray *)fDigs->At(iC-1):0;}
   inline void          DigitAdd      (Int_t c,TVector pad,int q,int cfm,int *tid             )     ;                                                 //add new digit
+  inline void          DigitAdd      (AliRICHDigit &dif                                      )     ;                                                 //add new digit
   inline void          DigitsCreate  (                                                       )     ;                                                 //create digits
          void          DigitsReset   (                                                       )     {if(fDigs)for(int i=0;i<kNchambers;i++){fDigs->At(i)->Clear();fNdigs[i]=0;}} //virtual
          void          DigitsPrint   (Int_t iEvent=0                                         )const;                                                 //prints digits
@@ -149,6 +146,13 @@ void AliRICH::DigitsCreate()
   for(Int_t i=0;i<kNchambers;i++) {fDigs->AddAt(new TClonesArray("AliRICHDigit",10000), i); fNdigs[i]=0;}
 }
 //__________________________________________________________________________________________________
+void AliRICH::DigitAdd(AliRICHDigit &dig)
+{
+//special for digit formed from raw  
+  TClonesArray &tmp=*((TClonesArray*)fDigs->At(dig.Chamber()-1));
+  new(tmp[fNdigs[dig.Chamber()-1]++])AliRICHDigit(dig);
+}    
+//__________________________________________________________________________________________________
 void AliRICH::DigitAdd(int c,TVector pad,int q,int cfm,int *tid)
 {
   TClonesArray &tmp=*((TClonesArray*)fDigs->At(c-1));
index 479c9a4..6e73674 100644 (file)
@@ -14,8 +14,7 @@
 //  **************************************************************************
 
 #include "AliRICHCluster.h"
-#include <AliLog.h>
-
+#include <TMinuit.h>  //Solve()
  
 ClassImp(AliRICHCluster)
 //__________________________________________________________________________________________________
@@ -24,16 +23,161 @@ void AliRICHCluster::Print(Option_t*)const
 //Print current cluster  
   const char *status=0;
   switch(fStatus){
-    case      kRaw: status="raw"     ;break;
-    case kResolved: status="resolved";break;
-    case    kEmpty: status="empty"   ;break;
+    case      kFormed: status="formed"     ;break;
+    case    kUnfolded: status="unfolded"   ;break;
+    case         kCoG: status="CoGed"      ;break;
+    case       kEmpty: status="empty"      ;break;
   }
-  if(fDigits)    
-    ::Info("cluster","cfm=%10i, cs=%2i, SiMa=%6i, Shape=%5i, x=%7.3f, y=%7.3f, Q=%6i, %s with %i digits",
-                             fCFM,fChamber,fSize,fShape,fX,fY,fQdc,status,fDigits->GetEntriesFast());
-  else
-    AliInfo(Form("cfm=%10i, cs=%2i, SiMa=%6i, Shape=%5i, x=%7.3f, y=%7.3f, Q=%6i, %s with %i digits",
-                             fCFM,fChamber,fSize,fShape,fX,fY,fQdc,status,0));
+  Int_t iNdigs=0;  if(fDigits) iNdigs=fDigits->GetEntriesFast();
     
+  Printf("cfm=%10i, cs=%2i, Size=%2i Maxima=%2i, Shape=%5i, pos=(%7.3f,%7.3f) Q=%6i, %s",
+            fCFM,fChamber,Size(),Nlocmax(),fShape,fX,fY,fQdc,status,iNdigs);
+  for(Int_t i=0;i<iNdigs;i++) Digit(i)->Print();    
 }//Print()
+
+
+//__________________________________________________________________________________________________
+TMinuit* AliRICHCluster::Solve()
+{
+//At this point, cluster contains a list of digits, cluster charge is precalculated as a sum of digits charges (in AddDigit()),
+//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())
+//Here we decide what to do with this cluster: unfold or just calculate center of gravity
+//Arguments: none
+//  Returns: pointer to fitter or 0 if no unfolding decided 
+  TMinuit *pMinuit=0;
+  if(Size()>=2 && AliRICHParam::IsResolveClusters())
+    pMinuit=Unfold();
+  else
+    CoG(0);
+  return pMinuit;
+}//Solve()
+//__________________________________________________________________________________________________
+void AliRICHCluster::FitFunc(Int_t &iNpars, Double_t *, Double_t &chi2, Double_t *aPar, Int_t )
+{
+//Cluster fit function 
+//par[0]=x par[1]=y par[2]=q for the first Mathieson shape
+//par[3]=x par[4]=y par[5]=q for the second Mathieson shape and so on up to iNpars/3 Mathieson shapes
+//We need to calculate Qpad - Qpadmath summup over all pads of the cluster
+//Here Qpad is a actual charge of the pad, Qpadmath is calculated charge of the pad induced by all Mathiesons
+//Arguments: iNpars - number of parameters which is number of local maxima of cluster * 3
+//           chi2   - function result to be minimised 
+//           aPar   - parametrs array of size iNpars            
+//  Returns: none  
+  AliRICHCluster *pClu=(AliRICHCluster*)gMinuit->GetObjectFit();
+  Int_t iNmathiesons = iNpars/3;
+    
+  TVector2 curMathiesonPos;
+  chi2 = 0;
+  for(Int_t i=0;i<pClu->Size();i++){//digits loop
+    TVector    pad     = pClu->Digit(i)->Pad();
+    Double_t dQpad     = pClu->Digit(i)->Qdc();
+    Double_t dQpadmath = 0;
+    for(Int_t j=0;j<iNmathiesons;j++){//all Mathiesons may contribute to this pad
+      curMathiesonPos.Set(aPar[3*j],aPar[3*j+1]);//get current position for current Mathieson
+      dQpadmath += aPar[3*j+2]*AliRICHParam::FracQdc(curMathiesonPos,pad);//sums up contributions to the current pad from all Mathiesons
+    }
+    chi2 += TMath::Power((dQpadmath-dQpad),2);
+  }//digits loop     
+}//RichClusterFitFunction()
+//__________________________________________________________________________________________________
+TMinuit* AliRICHCluster::Unfold()
+{
+//This methode is invoked from Solve() when decided to unfold this cluster
+//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
+//Arguments: none
+//  Returns: pointer to fitter for retriving parameters
+    
+  TMinuit *pMinuit = new TMinuit(15); //init MINUIT with max 15 parameters (maxim 5 mathiesons, 3 params per matheson )
+  pMinuit->SetObjectFit((TObject*)this);
+  pMinuit->SetFCN(AliRICHCluster::FitFunc);//set fit function
+  Double_t aArg=-1,parStart,parStep,parLow,parHigh; Int_t iErrFlg; //tmp for MINUIT parameters definitions
+  pMinuit->mnexcm("SET PRI" ,&aArg,1,iErrFlg); //suspend all printout from TMinuit 
+  
+  Int_t iLocMaxCnt=0;
+//Strategy is to check if the current pad has QDC more then all neigbours 
+  for(Int_t iDig1=0;iDig1<Size();iDig1++) {//first digits loop
+    AliRICHDigit *pDig1 = Digit(iDig1);//take the current digit
+    Int_t iHowManyMoreCnt = 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
+      AliRICHDigit *pDig2 = Digit(iDig2);
+      if(iDig1==iDig2) continue;             //no need to compare 
+      Int_t dist = TMath::Sign(Int_t(pDig1->PadX()-pDig2->PadX()),1)+TMath::Sign(Int_t(pDig1->PadY()-pDig2->PadY()),1);//distance between pads
+      if(dist==1)//means pads are neighbours
+         if(pDig2->Qdc()>=pDig1->Qdc()) iHowManyMoreCnt++;//count number of pads with Q more then Q of current pad
+    }//second digits loop
+    if(iHowManyMoreCnt==0&&iLocMaxCnt<6){//this pad has Q more then any neighbour so it's local maximum
+        TVector2 x2=AliRICHParam::Pad2Loc(pDig1->Pad());//take pad center position and use it as parameter for current Mathienson shape
+        pMinuit->mnparm(3*iLocMaxCnt  ,Form("x%i",iLocMaxCnt),parStart=x2.X()      ,parStep=0.01,parLow=0,parHigh=0,iErrFlg);
+        pMinuit->mnparm(3*iLocMaxCnt+1,Form("y%i",iLocMaxCnt),parStart=x2.Y()      ,parStep=0.01,parLow=0,parHigh=0,iErrFlg);
+        pMinuit->mnparm(3*iLocMaxCnt+2,Form("q%i",iLocMaxCnt),parStart=pDig1->Qdc(),parStep=0.01,parLow=0,parHigh=0,iErrFlg);//
+        iLocMaxCnt++;
+    }//if this pad is local maximum
+  }//first digits loop
+  
+  fSize+=iLocMaxCnt;
+  if(iLocMaxCnt>0&&iLocMaxCnt<6){ //resonable number of local maxima to fit
+    Double_t aArg=0;
+    pMinuit->mnexcm("MIGRAD",&aArg,0,iErrFlg);//start fitting
+    fStatus=kUnfolded;
+  }else{
+    delete pMinuit;
+    pMinuit=0;
+    CoG(0);
+  }
+  return pMinuit;  
+}//Unfold()
 //__________________________________________________________________________________________________
+void AliRICHCluster::CoG(Int_t nLocals)
+{
+//Calculates naive cluster position as a center of gravity of its digits.
+//Also determines the box fully contaning this cluster
+//Arguments:     
+  Float_t xmin=999,ymin=999,xmax=0,ymax=0;
+  fX=fY=0;
+  for(Int_t iDig=0;iDig<Size();iDig++) {
+    AliRICHDigit *pDig=Digit(iDig);
+    TVector pad=pDig->Pad(); Double_t q=pDig->Qdc();
+    TVector2 x2=AliRICHParam::Pad2Loc(pad);
+    fX += x2.X()*q;fY +=x2.Y()*q;
+    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];
+   }
+   fX/=fQdc;fY/=fQdc;//Center of Gravity
+
+   TVector2 center = AliRICHParam::Pad2Loc(AliRICHParam::Loc2Pad(TVector2(fX,fY)));
+   fX += AliRICHParam::CogCorr(fX-center.X());//correct cluster position for sinoid
+
+   fShape=Int_t(100*(xmax-xmin+1)+ymax-ymin+1);//find box containing cluster
+   fSize+=nLocals;
+   fStatus=kCoG;
+}//CoG()
+//__________________________________________________________________________________________________
+void AliRICHCluster::Test(const TVector2 &hitX2,Double_t dEloss)
+{
+//This is to test all cluster functionality
+//Method uses AddDigit() to add a predifined pad structure and then calls Solve   
+  Int_t iQtot=AliRICHParam::TotQdc(hitX2,dEloss);
+  if(iQtot==0){
+    Printf("Provided hit position out of sensitive area");
+    return;
+  }
+  TVector area=AliRICHParam::Loc2Area(hitX2);
+  TVector pad(2);
+  for(pad[1]=area[1];pad[1]<=area[3];pad[1]++){//affected pads loop first y
+    for(pad[0]=area[0];pad[0]<=area[2];pad[0]++){//then x               
+      Double_t dQpad=iQtot*AliRICHParam::FracQdc(hitX2,pad);//charge fraction from Mathieson centered at x to pad
+      AddDigit(new AliRICHDigit(3,(Int_t)pad[0],(Int_t)pad[1],dQpad));
+    }//affected pads loop 
+  }
+  TMinuit *pMinuit=Solve();
+  Print();
+  Printf("Initial hit (%.2f,%.2f) Qtot=%i Eloss=%.2f",hitX2.X(),hitX2.Y(),iQtot,dEloss);
+  Double_t d1,d2,d3; Int_t iErrFlg;TString sName; //tmp vars for TMinuit
+  Double_t x,y,q;
+  for(Int_t i=0;i<Nlocmax();i++){//retrive fitting results
+    pMinuit->mnpout(3*i   ,sName,  x, d1 , d2, d3, iErrFlg);
+    pMinuit->mnpout(3*i+1 ,sName,  y, d1 , d2, d3, iErrFlg);
+    pMinuit->mnpout(3*i+2 ,sName,  q, d1 , d2, d3, iErrFlg);
+  }
+  Printf(" Fitted hit (%.2f,%.2f) Qfit=%.0f",x,y,q);
+  delete pMinuit; pMinuit=0;  Reset();
+}//Test()
index f4b339b..74117eb 100644 (file)
@@ -3,25 +3,27 @@
 /* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
  * See cxx source for full Copyright notice                               */
 
-#include <TObject.h>
-#include <TVector.h>
-#include <TVector2.h>
+#include <TObject.h>         //base class
+#include <TVector2.h>        //DistTo
 #include "AliRICHDigit.h"
+class TMinuit;
 
 class AliRICHCluster :public TObject
 {
 public:
-  enum ClusterStatus {kEdge,kShape,kSize,kRaw,kResolved,kEmpty}; 
-                    AliRICHCluster():TObject(),fCFM(0),fSize(0),fShape(0),fQdc(0),fChamber(0),fX(0),fY(0),fStatus(kEmpty),fDigits(0) {}  //default ctor  
-  virtual          ~AliRICHCluster()                 {AliDebug(1,"Start");/*Reset();*/}                                                  //dtor
+  enum EClusterStatus {kFormed,kCoG,kUnfolded,kEmpty}; 
+  AliRICHCluster()                                               :TObject(),fCFM(-1),fSize(-1),fShape(-1),fQdc(-1),fChamber(-1),fX(-1),fY(-1),fStatus(kEmpty),fDigits(0) {}  //default ctor  
+  AliRICHCluster(Int_t cs,Double_t x,Double_t y,Int_t q,Int_t sm):TObject(),fCFM(-1),fSize(sm),fShape(-1),fQdc(q ),fChamber(cs),fX(x ),fY(y ),fStatus(kEmpty),fDigits(0) {}  //default ctor  
+  virtual          ~AliRICHCluster()                 {}                                                  //dtor
                    // AliRICHcluster(const AliRICHcluster& clus):TObject(clus)                                                         {}  //copy ctor 
-   AliRICHCluster&  operator=(const AliRICHCluster&)                 {return *this;}                                                     //copy operator
+  AliRICHCluster&  operator=(const AliRICHCluster&)                 {return *this;}                                                     //copy operator
                    
+         void      Print(Option_t *option="")const;                                                       //
   
   
-         void       Reset()                          {DeleteDigits();fCFM=fSize=fShape=fQdc=fChamber=0;fX=fY=0;fStatus=kEmpty;} //cleans the cluster
+         void       Reset()                          {DeleteDigits();fCFM=fSize=fShape=fQdc=fChamber=-1;fX=fY=-1;fStatus=kEmpty;} //cleans the cluster
          void       DeleteDigits()                   {if(fDigits) {delete fDigits;} fDigits=0;}           //deletes the list of digits  
-         Int_t      Nlocals()                   const{return fSize-10000*(fSize/10000);}                //number of local maximums
+         Int_t      Nlocmax()                   const{return fSize-10000*(fSize/10000);}                //number of local maximums
          Int_t      Size()                      const{return fSize/10000;}                              //number of digits in cluster
          Int_t      Fsize()                     const{return fSize;}                                    //
          Int_t      Shape()                     const{return fShape;}                                   //cluster shape rectangulare
@@ -47,21 +49,28 @@ public:
          Bool_t     IsFeedback()                const{return Nfeedbacks()!=0;}                             //
          Int_t      CombiPid()                  const{return fCFM;}                                        //
          void       CFM(Int_t c,Int_t f,Int_t m)     {fCFM=1000000*c+1000*f+m;}                            //cluster contributors
-         TObjArray* Digits()                    const{return fDigits;}                                     //  
-  virtual void      Print(Option_t *option="")const;                                                   //
-  inline  void      AddDigit(AliRICHDigit *pDig);                                                      //
-  inline  void      CoG(Int_t nLocals);                                                                //calculates center of gravity
+         TObjArray*    Digits()                    const{return fDigits;}                                     //  
+         
+  inline void          AddDigit(AliRICHDigit *pDig);                                                          //add new digit ot the cluster
+         AliRICHDigit* Digit  (Int_t i                      )const{return (AliRICHDigit*)fDigits->At(i); }//get pointer to i-th digit without existence check 
+         TMinuit*      Solve  (                             );                                            //calculates cluster position
+         TMinuit*      Unfold (                             );                                            //decompose cluster n. loc max clusters
+         void          Set    (Double_t x,Double_t y,Int_t q)     {fX=x;fY=y,fQdc=q;                     }//set some cluster properties
+  static void          FitFunc(Int_t &iNpars, Double_t *, Double_t &chi2, Double_t *aPar, Int_t);         //fit function to be used by MINUIT
+          
+         void      CoG(Int_t iNlocmax);                                                                   //calculates center of gravity
           void      Fill(AliRICHCluster *pRaw,Double_t x,Double_t y,Double_t q,Int_t cfm)              //form new resolved cluster from raw one
-                    {fCFM=cfm;fChamber=pRaw->Fchamber();fSize=pRaw->Fsize();fQdc=(Int_t)(q*pRaw->Q());fX=x;fY=y;fStatus=kResolved;} 
+                    {fCFM=cfm;fChamber=pRaw->Fchamber();fSize=pRaw->Fsize();fQdc=(Int_t)(q*pRaw->Q());fX=x;fY=y;fStatus=kUnfolded;} 
          Double_t   DistTo(TVector2 x)          const{return TMath::Sqrt((x.X()-fX)*(x.X()-fX)+(x.Y()-fY)*(x.Y()-fY));} //distance to given point 
          Double_t   DistX(TVector2 x)           const{return (x.X()-fX);} //distance in x to given point 
          Double_t   DistY(TVector2 x)           const{return (x.Y()-fY);} //distance to given point 
+         void       Test(const TVector2 &x,Double_t dEloss=0);            //test cluster fuctionality by provided hit with energy in eV
 protected:
   Int_t         fCFM;         //1000000*Ncerenkovs+1000*Nfeedbacks+Nmips  
-  Int_t         fSize;        //10000*(how many digits belong to this cluster) + nLocalMaxima     
+  Int_t         fSize;        //10000*(N digits) + N maxima     
   Int_t         fShape;       //100*xdim+ydim box containing the cluster
   Int_t         fQdc;         //QDC value
-  Int_t         fChamber;     //10*module number+sector number 
+  Int_t         fChamber;     //10*chamber number+sector number 
   Double_t      fX;           //local x postion 
   Double_t      fY;           //local y postion  
   Int_t         fStatus;      //flag to mark the quality of the cluster   
@@ -72,33 +81,12 @@ protected:
 void AliRICHCluster::AddDigit(AliRICHDigit *pDig)
 {
 // Adds a given digit to the list of digits belonging to this cluster    
-  if(!fDigits) {fQdc=fSize=fCFM=0;fDigits = new TObjArray;}
-  fQdc+=(Int_t)pDig->Q(); fDigits->Add(pDig);
-  fChamber=10*pDig->C()+pDig->S();
+  if(!fDigits) {fQdc=fSize=0;fDigits = new TObjArray;}
+  fDigits->Add(pDig);
+  fQdc+=(Int_t)pDig->Qdc(); 
+  fChamber=10*pDig->Chamber()+pDig->Sector();
   fSize+=10000;
-  fStatus=kRaw;
+  fStatus=kFormed;
 }
 //__________________________________________________________________________________________________
-void AliRICHCluster::CoG(Int_t nLocals)
-{
-// Calculates naive cluster position as a center of gravity of its digits.
-  Float_t xmin=999,ymin=999,xmax=0,ymax=0;   
-  fX=fY=0;
-  for(Int_t iDig=0;iDig<Size();iDig++) {
-    AliRICHDigit *pDig=(AliRICHDigit*)fDigits->At(iDig);
-    TVector pad=pDig->Pad(); Double_t q=pDig->Q();
-    TVector2 x2=AliRICHParam::Pad2Loc(pad);
-    fX += x2.X()*q;fY +=x2.Y()*q;
-    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];
-   }
-   fX/=fQdc;fY/=fQdc;//Center of Gravity
-
-   TVector2 center = AliRICHParam::Pad2Loc(AliRICHParam::Loc2Pad(TVector2(fX,fY)));
-   fX += AliRICHParam::CogCorr(fX-center.X());
-
-   fShape=Int_t(100*(xmax-xmin+1)+ymax-ymin+1);//find box containing cluster
-   fSize+=nLocals;
-   fStatus=kRaw;
-}//CoG()
-//__________________________________________________________________________________________________
 #endif
index 875cc7a..26e3004 100644 (file)
@@ -85,7 +85,7 @@ void AliRICHClusterFinder::FindClusters(Int_t iChamber)
 
   for(Int_t iDigN=0;iDigN<iNdigits;iDigN++){//digits loop for a given chamber    
     AliRICHDigit *dig=(AliRICHDigit*)R()->Digits(iChamber)->At(iDigN);
-    Int_t i=dig->X();   Int_t j=dig->Y();
+    Int_t i=dig->PadX();   Int_t j=dig->PadY();
     if(fDigitMap->TestHit(i,j)==kUsed) continue;//this digit is already taken, go after next digit
        
     FormRawCluster(i,j);//form raw cluster starting from (i,j) pad 
@@ -178,7 +178,7 @@ void AliRICHClusterFinder::FindClusterContribs(AliRICHCluster *pCluster)
 void  AliRICHClusterFinder::FormRawCluster(Int_t i, Int_t j)
 {
 //Builds the raw cluster (before deconvolution). Starts from the first pad (i,j) then calls itself recursevly  for all neighbours.
-  AliDebug(1,Form("Start with digit(%i,%i) Q=%f",i,j,((AliRICHDigit*)fDigitMap->GetHit(i,j))->Q()));
+  AliDebug(1,Form("Start with digit(%i,%i) Q=%f",i,j,((AliRICHDigit*)fDigitMap->GetHit(i,j))->Qdc()));
   
   fRawCluster.AddDigit((AliRICHDigit*) fDigitMap->GetHit(i,j));//take this pad in cluster
   fDigitMap->FlagHit(i,j);//flag this pad as taken  
@@ -198,13 +198,13 @@ void AliRICHClusterFinder::FindLocalMaxima()
     AliRICHDigit *pDig1 = (AliRICHDigit *)fRawCluster.Digits()->At(iDig1);
     if(!pDig1) {fNlocals=0;return;}
     TVector pad1 = pDig1->Pad();
-    Int_t padQ1 = (Int_t)(pDig1->Q()+0.1);
-    Int_t padC1 = pDig1->ChFbMi();
+    Int_t padQ1 = (Int_t)(pDig1->Qdc()+0.1);
+    Int_t padC1 = pDig1->Cfm();
     for(Int_t iDig2=0;iDig2<fRawCluster.Size();iDig2++) {
       AliRICHDigit *pDig2 = (AliRICHDigit *)fRawCluster.Digits()->At(iDig2);
       if(!pDig2) {fNlocals=0;return;}
       TVector pad2 = pDig2->Pad();
-      Int_t padQ2 = (Int_t)(pDig2->Q()+0.1);
+      Int_t padQ2 = (Int_t)(pDig2->Qdc()+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);
@@ -299,12 +299,9 @@ void AliRICHClusterFinder::FitCoG()
     AliDebug(1,Form("qfrac %i vstart %f lower %f upper %f ",i,vstart,lower,upper));
   }
   
-  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);
+  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);
   
@@ -355,7 +352,7 @@ void RICHMinMathieson(Int_t &npar, Double_t *, Double_t &chi2, Double_t *par, In
   Int_t qtot = pRawCluster->Q();
   for(Int_t i=0;i<pRawCluster->Size();i++) {
     TVector  pad=((AliRICHDigit *)pRawCluster->Digits()->At(i))->Pad();
-    Double_t padQ = ((AliRICHDigit *)pRawCluster->Digits()->At(i))->Q();
+    Double_t padQ = ((AliRICHDigit *)pRawCluster->Digits()->At(i))->Qdc();
     Double_t qfracpar=0;
     for(Int_t j=0;j<nFunctions;j++) {
       qfracpar += q[j]*AliRICHParam::FracQdc(centroid[j],pad);
index e3e2532..01c18cf 100644 (file)
@@ -13,8 +13,7 @@
 //  * provided "as is" without express or implied warranty.                  *
 //  **************************************************************************
 
-#include "AliRICHDigit.h"
-#include <AliLog.h>
+#include "AliRICHDigit.h"//class header
 
 ClassImp(AliRICHDigit)
 
@@ -22,7 +21,12 @@ ClassImp(AliRICHDigit)
 void AliRICHDigit::Print(Option_t*)const
 {
 //Print current digit  
-  AliInfo(Form("cfm=%9i, cs=%2i, x=%3i, y=%3i, q=%8.3f, TID1=%5i, TID2=%5i, TID3=%5i",
-                  fCFM,fChamber,fPadX,fPadY,fQdc,fTracks[0],fTracks[1],fTracks[2]));
+//Arguments: option string not used
+//  Returns: none    
+  Printf("pad=(%2i,%2i,%3i,%3i), q=%8.3f, cfm=%9i, TID=(%5i,%5i,%5i)",
+               Chamber(),Sector(),PadX(),PadY() ,  Qdc(),    Cfm()  , fTracks[0],fTracks[1],fTracks[2]);
 }
 //__________________________________________________________________________________________________
+void AliRICHDigit::Test()
+{
+}
index 20c0c1e..8063624 100644 (file)
@@ -3,32 +3,61 @@
 /* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
  * See cxx source for full Copyright notice                               */
 
-#include <AliDigit.h>
+#include <AliDigit.h>      //base class  
+#include <AliBitPacking.h> //Dig2Raw()
 #include "AliRICHParam.h"
 
+//RICH DDL ID allowed range is [0x700,0x714] or in decimal notation [1792,1812]. 20 DDL files are reserved. (kDdlOffset)
+//RICH actually uses 14 DDLs (kNddls), 2 per chamber, even number for left part(1-3-5) odd number for right part(2-4-6) 
+//So the chamber-DDL map is:
+//N  0 1L=0x700 1792          N  1 1R=0x701 1793
+//N  2 2L=0x702 1794          N  3 2R=0x703 1795
+//N  4 3L=0x704 1796          N  5 3R=0x705 1797
+//N  6 4L=0x706 1798          N  7 4R=0x707 1799
+//N  8 5L=0x708 1800          N  9 5R=0x709 1801
+//N 10 6L=0x70a 1802          N 11 6R=0x70b 1803
+//N 12 7L=0x70c 1804          N 13 7R=0x70d 1805
+//RICH has no any propriate header just uses the common one
+//RICH chamber is divide on 2 halves vertically 
+//Half chamber is divided by 24 rows counted from 1 to 24 (8 raws per sector) from top to bottom for left half chamber (sectors 1-3-5) 
+//                                                                         and from bottom to top for right half chamber (sectors 2-4-6) as seen from MARS (0,0,0)
+//Raw is composed from 10 DILOGIC chips (kNchips) counted from left to right from 1 to 10  as seen from MARS (0,0,0)
+//So each DILOGIC chip serves 48 channels for the 8x6 pads box (kChipX,kChipY). Channels counted from 0 to 47.
+//??????? Currently the exact mapping of DILOGIC addresses to pads is not known. So we invented horizontal zig-zag ???????  
+//So RICH raw word is  32 bits word with structure:   
+// 00000        rrrrr                      dddd                               aaaaaa                          qqqqqqqqqqqq 
+// 5 bits zero  5 bits raw number (1..24)  4 bits DILOGIC chip number (1..10) 6 bits DILOGIC address (0..47)  12 bits QDC value (0..4095)
+
 class AliRICHDigit :public AliDigit
 {
 public:
-  AliRICHDigit():AliDigit(),fCFM(0),fChamber(0),fPadX(0),fPadY(0),fQdc(-1){fTracks[0]=fTracks[1]=fTracks[2]=-1;}
+  enum EAbsPad {kChamber=10000000,kPadX=1000};                           //absolute pad number structure
+  enum ERawProp{kChipX=8,kChipY=6,kNchips=10,kNddls=14,kRichRawId=7,kDdlOffset=0x700};//DILOGIC is 8x6 pads
+  AliRICHDigit()                                  :AliDigit(),fCFM(-1),fChamber(-1  )  ,fPadX(-1)    ,fPadY(-1)    ,fQdc(-1) {}
+  AliRICHDigit(Int_t c,Int_t x,Int_t y,Double_t q):AliDigit(),fCFM(-1),fChamber(10*c)  ,fPadX(x )    ,fPadY(y )    ,fQdc(q ) {}
   AliRICHDigit(Int_t c,TVector pad,Double_t q,Int_t cfm,Int_t tid0,Int_t tid1,Int_t tid2):fCFM(cfm)  
        {fPadX=(Int_t)pad[0];fPadY=(Int_t)pad[1];fQdc=q;fChamber=10*c+AliRICHParam::Pad2Sec(pad);fTracks[0]=tid0;fTracks[1]=tid1;fTracks[2]=tid2;}
   virtual ~AliRICHDigit() {;}
-    
-  Int_t    Compare(const TObject *pObj) const
-                 {if(Id()==((AliRICHDigit*)pObj)->Id())return 0;else if(Id()>((AliRICHDigit*)pObj)->Id())return 1;else return -1;}  //virtual      
-  virtual Bool_t   IsSortable()                 const{return kTRUE;}                              //sort interface
-  virtual void     Print(Option_t *option="")   const;                                            //virtual
+//framework part    
+         Bool_t   IsSortable  (                   )const{return kTRUE;}                                                    //provision to use TObject::Sort()
+  inline Int_t    Compare     (const TObject *pObj)const;                                                                  //provision to use TObject::Sort()
+         void     Print       (Option_t *option="")const;                                                                  //TObject Print() overload
 //private part  
-  Int_t    ChFbMi()                     const{return fCFM;}                               //particle mixture for this digit
-  Int_t    C()                          const{return fChamber/10;}                        //chamber number
-  Int_t    S()                          const{return fChamber-(fChamber/10)*10;}          //sector number
-  Int_t    X()                          const{return fPadX;}                              //x position of the pad
-  Int_t    Y()                          const{return fPadY;}                              //y postion of the pad
-  TVector  Pad()                        const{Float_t v[2]={fPadX,fPadY}; return TVector(2,v);}
-  Int_t    Id()                         const{return fChamber*10000000+fPadX*1000+fPadY;} //absolute id of this pad
-  Double_t Q()                          const{return fQdc;}                               //charge in terms of ADC channels
-  void     AddTidOffset(Int_t offset) 
-    {for (Int_t i=0; i<3; i++) if (fTracks[i]>0) fTracks[i]+=offset;};
+         void     AddTidOffset(Int_t offset     )     {for (Int_t i=0; i<3; i++) if (fTracks[i]>0) fTracks[i]+=offset;}; //needed for merging
+         Int_t    Cfm         (                 )const{return fCFM;}                                                     //particle mixture for this digit
+         Int_t    Chamber     (                 )const{return fChamber/10;}                                              //chamber number
+         Int_t    Sector      (                 )const{return fChamber%10;}                                              //sector number
+         Int_t    PadX        (                 )const{return fPadX;}                                                    //x position of the pad
+         Int_t    PadY        (                 )const{return fPadY;}                                                    //y postion of the pad     
+         TVector  Pad         (                 )const{Float_t v[2]={fPadX,fPadY}; return TVector(2,v);}
+         Int_t    PadAbs      (                 )const{return fChamber*kChamber+fPadX*kPadX+fPadY;}                      //absolute id of this pad
+         Double_t Qdc         (                 )const{return fQdc;}                                                     //charge in terms of ADC channels
+  inline Int_t    Dig2Raw     (        UInt_t &w)const;                                                                  //returns DDL ID and fill raw 32 bits word
+  inline void     Raw2Dig     (Int_t d,UInt_t  w);                                                                       //(DDL,word32)->(ch,sec,padx,pady,QDC)
+  static Int_t    P2C         (Int_t pad        )     {return pad/kChamber;}                                             //abs pad number-> chamber number
+  static Int_t    P2X         (Int_t pad        )     {return pad%kChamber/kPadX;}                                       //abs pad number-> pad X number
+  static Int_t    P2Y         (Int_t pad        )     {return pad%kChamber%kPadX;}                                       //abs pad number-> pad Y number
+         void     Test        (                 );                                                                       //used to test all possible digit manipulations
 protected:
   Int_t    fCFM;  //1000000*Ncerenkovs+1000*Nfeedbacks+Nmips  
   Int_t    fChamber;  //10*chamber number+ sector number 
@@ -37,5 +66,58 @@ protected:
   Double_t fQdc;      //QDC value, fractions are permitted for summable procedure  
   ClassDef(AliRICHDigit,3) //RICH digit class       
 };//class AliRICHDigit
+//__________________________________________________________________________________________________
+Int_t AliRICHDigit::Compare(const TObject *pObj) const
+{
+//Used in Sort() method to compare to objects. Note that abs pad structure is first x then y, hence will be sorted on column basis.
+//This feature is used in digitizer to facilitate finding of sdigits for the same pad as they will be together after sorting.
+//Arguments: pObj - pointer to object to compare with
+//  Retunrs: -1 if AbsPad less then in pObj, 1 if more and 0 if they are the same      
+  if(PadAbs()==((AliRICHDigit*)pObj)->PadAbs())
+    return 0;
+  else if(PadAbs()>((AliRICHDigit*)pObj)->PadAbs())
+    return 1;
+  else 
+    return -1;
+}
+//__________________________________________________________________________________________________
+void AliRICHDigit::Raw2Dig(Int_t ddl,UInt_t w32)
+{
+//Reads next raw word from raw data stream and convert     
+//Arguments: w32 - 32 bits word as in raw data stream
+//           ddl - DDL file number  0 1 2 3 4 ... 13
+//  Returns: none
+      fQdc = AliBitPacking::UnpackWord(w32, 0,11);  // 0000 0rrr rrdd ddaa aaaa qqqq qqqq qqqq 
+  UInt_t a = AliBitPacking::UnpackWord(w32,12,17);  // 3322 2222 2222 1111 1111 1000 0000 0000 
+  UInt_t d = AliBitPacking::UnpackWord(w32,18,21);  // 1098 7654 3210 9876 5432 1098 7654 3210 
+  UInt_t r = AliBitPacking::UnpackWord(w32,22,26);  // r- iRawN d- iChiN a- iChiC              
+  
+  fPadY    = (r-1)*kChipY+a/kChipX+1;
+  fPadX    = (d-1)*kChipX+a%kChipX+1;      fPadX+=(ddl%2)*kChipX*kNchips;//if ddl is odd then right half of the chamber
+  TVector pad(2); pad[0]=fPadX;pad[1]=fPadY;
+  fChamber = ((ddl+2)/2)*10+AliRICHParam::Pad2Sec(pad);  // ddl 0..13 to chamber 1..7
+}  
+//__________________________________________________________________________________________________
+Int_t AliRICHDigit::Dig2Raw(UInt_t &w32)const
+{
+//Convert digit structure to raw word format
+//Arguments: 32 bits raw word to fill
+//  Returns: DDL ID where to write this digit
+  Int_t ddl=2*Chamber()-1;              //chamber 1..7 -> DDL 0..13, this idDdl is for right half (sectors 2 4 6), to be decremented if d < kNchips
+  UInt_t a =  (PadY()-1)%kChipY*kChipX+(PadX()-1)%kChipX;          //invented to be horizontal zig-zag
+  UInt_t r =1+(PadY()-1)/kChipY;      
+  UInt_t d =1+(PadX()-1)/kChipX;    
+  if(d>kNchips)
+    d-=kNchips;              //chip number more then kNchips means right half of chamber, goes to this ddl
+  else
+    ddl--;                   //chip number less then kNchips means left half of the chamber, goes to ddl-1 
+  
+  w32=0;
+  AliBitPacking::PackWord((UInt_t)fQdc,w32, 0,11);  // 0000 0rrr rrdd ddaa aaaa qqqq qqqq qqqq
+  AliBitPacking::PackWord(           a,w32,12,17);  // 3322 2222 2222 1111 1111 1000 0000 0000
+  AliBitPacking::PackWord(           d,w32,18,21);  // 1098 7654 3210 9876 5432 1098 7654 3210 
+  AliBitPacking::PackWord(           r,w32,22,26);  
+  return ddl; //ddl 0..13 where to write this digit 
+}
 
 #endif
index 160ea92..4e85828 100644 (file)
@@ -61,14 +61,14 @@ void AliRICHDigitizer::Exec(Option_t*)
   Int_t iNdigsPerPad=0;                   //how many sdigits for a given pad
   for(Int_t i=0;i<tmpCA.GetEntries();i++){//sdigits loop (sorted)
     AliRICHDigit *pSdig=(AliRICHDigit*)tmpCA.At(i);//get new sdigit
-    if(pSdig->Id()==iId){//still the same pad
-      iNdigsPerPad++;         dQdc+=pSdig->Q();      iCfm+=pSdig->ChFbMi();//sum up charge and cfm
-      if(pSdig->ChFbMi()==1) aTids[0] = pSdig->GetTrack(0); // force the first tid to be mip's tid if it exists in the current pad
+    if(pSdig->PadAbs()==iId){//still the same pad
+      iNdigsPerPad++;         dQdc+=pSdig->Qdc();      iCfm+=pSdig->Cfm();//sum up charge and cfm
+      if(pSdig->Cfm()==1) aTids[0] = pSdig->GetTrack(0); // force the first tid to be mip's tid if it exists in the current pad
       if(iNdigsPerPad<=3)        aTids[iNdigsPerPad-1]=pSdig->GetTrack(0);
       else                         AliWarning(Form("More then 3 sdigits for the given pad X:%3.0f Y:%3.0f",pSdig->Pad()(0),pSdig->Pad()(1)));
     }else{//new pad, add the pevious one
         if(iId!=-1 && AliRICHParam::IsOverTh(iChamber,pad,dQdc)) pOutRich->DigitAdd(iChamber,pad,(Int_t)dQdc,iCfm,aTids); //add newly created dig
-        iChamber=pSdig->C(); pad=pSdig->Pad(); iCfm=pSdig->ChFbMi(); dQdc=pSdig->Q();  iId=pSdig->Id();                    //init all values by current sdig
+        iChamber=pSdig->Chamber(); pad=pSdig->Pad(); iCfm=pSdig->Cfm(); dQdc=pSdig->Qdc();  iId=pSdig->PadAbs();                    //init all values by current sdig
         iNdigsPerPad=1; aTids[0]=pSdig->GetTrack(0); aTids[1]=aTids[2]=-1;
       }
   }//sdigits loop (sorted)
index b6837a9..1236fc0 100644 (file)
@@ -34,7 +34,7 @@ void  AliRICHMap::FillHits()
   if(!fNdigits) return;    
   for(Int_t iDigN=0;iDigN<fNdigits;iDigN++){
     AliRICHDigit *pDig= (AliRICHDigit*)fDigits->At(iDigN);
-    SetHit(pDig->X(),pDig->Y(),iDigN);
+    SetHit(pDig->PadX(),pDig->PadY(),iDigN);
   }
 }
 //__________________________________________________________________________________________________
index b41b68d..432bb9b 100644 (file)
@@ -16,7 +16,6 @@
 #include <AliLog.h>
 #include <TClass.h>
 
-
 static const int kNchambers=7;     //number of RICH chambers 
 static const int kNpadsX = 160;    //number of pads along X in single chamber
 static const int kNpadsY = 144;    //number of pads along Y in single chamber
@@ -62,29 +61,55 @@ public:
   AliRICHChamber* C(Int_t i)                 {return (AliRICHChamber*)fpChambers->UncheckedAt(i-1);}      //returns pointer to chamber i
   Int_t           Nchambers()                {return fpChambers->GetEntriesFast();}      //returns number of chambers 
 //Geometrical properties  
-  static Int_t    NpadsX()                   {return kNpadsX;}                           //pads along X in chamber
-  static Int_t    NpadsY()                   {return kNpadsY;}                           //pads along Y in chamber
-  static Int_t    NpadsXsec()                {return NpadsX()/2;}                        //pads along X in sector
-  static Int_t    NpadsYsec()                {return NpadsY()/3;}                        //pads along Y in sector
-  static Double_t DeadZone()                 {return 2.6;}                               //dead zone size in cm  
-  static Double_t PadSizeX()                 {return 0.8;}                               //pad size x, cm 
-  static Double_t PadSizeY()                 {return 0.84;}                              //pad size y, cm   
-  
-  static Double_t SectorSizeX()              {return NpadsX()*PadSizeX()/2;}             //sector size x, cm
-  static Double_t SectorSizeY()              {return NpadsY()*PadSizeY()/3;}             //sector size y, cm 
-  static Double_t PcSizeX()                  {return NpadsX()*PadSizeX()+DeadZone();}    //PC size x, cm
-  static Double_t PcSizeY()                  {return NpadsY()*PadSizeY()+2*DeadZone();}  //PC size y, cm
-//   
-  static Double_t Zfreon()                   {return 1.5;}                               //freon thinkness, cm
-  static Double_t Zwin()                     {return 0.5;}                               //radiator quartz window, cm   
-  static Double_t Pc2Win()                   {return 8.0;}                               //cm between CsI PC and radiator quartz window
-  static Double_t Pc2Coll()                  {return 7.0;}                               //cm between CsI PC and third wire grid (collection wires)     
-  static Double_t Pc2Anod()                  {return 0.204;}                             //cm between CsI PC and first wire grid (anod wires)     
-  static Double_t Pc2Cath()                  {return 0.445;}                             //cm between CsI PC and second wire grid (cathode wires)
-  static Double_t Freon2Pc()                 {return Zfreon()+Zwin()+Pc2Win();}          //cm between CsI PC and entrance to freon
-  static Double_t PitchAnod()                {return PadSizeY()/2;}                      //cm between anode wires
-  static Double_t PitchCath()                {return PadSizeY()/4;}                      //cm between cathode wires
-  static Double_t PitchColl()                {return 0.5;}                               //cm between collection wires
+  static        Int_t      NpadsX()                   {return kNpadsX;}                           //pads along X in chamber
+  static        Int_t      NpadsY()                   {return kNpadsY;}                           //pads along Y in chamber
+  static        Int_t      NpadsXsec()                {return NpadsX()/2;}                        //pads along X in sector
+  static        Int_t      NpadsYsec()                {return NpadsY()/3;}                        //pads along Y in sector
+  static        Double_t   DeadZone()                 {return 2.6;}                               //dead zone size in cm  
+  static        Double_t   SectorSizeX()              {return NpadsX()*PadSizeX()/2;}             //sector size x, cm
+  static        Double_t   SectorSizeY()              {return NpadsY()*PadSizeY()/3;}             //sector size y, cm 
+  static        Double_t   PcSizeX()                  {return NpadsX()*PadSizeX()+DeadZone();}    //PC size x, cm
+  static        Double_t   PcSizeY()                  {return NpadsY()*PadSizeY()+2*DeadZone();}  //PC size y, cm
+  static        Double_t   Zfreon()                   {return 1.5;}                               //freon thinkness, cm
+  static        Double_t   Zwin()                     {return 0.5;}                               //radiator quartz window, cm   
+  static        Double_t   Pc2Win()                   {return 8.0;}                               //cm between CsI PC and radiator quartz window
+  static        Double_t   Pc2Coll()                  {return 7.0;}                               //cm between CsI PC and third wire grid (collection wires)     
+  static        Double_t   Pc2Anod()                  {return 0.204;}                             //cm between CsI PC and first wire grid (anod wires)     
+  static        Double_t   Pc2Cath()                  {return 0.445;}                             //cm between CsI PC and second wire grid (cathode wires)
+  static        Double_t   Freon2Pc()                 {return Zfreon()+Zwin()+Pc2Win();}          //cm between CsI PC and entrance to freon
+  static        Double_t   PitchAnod()                {return PadSizeY()/2;}                      //cm between anode wires
+  static        Double_t   PitchCath()                {return PadSizeY()/4;}                      //cm between cathode wires
+  static        Double_t   PitchColl()                {return 0.5;}                               //cm between collection wires
+  static        Double_t   PadSizeX()                 {return 0.8;}                               //pad size x, cm 
+  static        Double_t   PadSizeY      (                               ){return 0.84;}                              //pad size y, cm   
+//trasformation methodes
+  static        Int_t      Pad2Cha       (Int_t pad                      ){return pad/100000000;                     }//abs pad -> chamber
+  static        Int_t      Pad2Sec       (Int_t pad                      ){return pad%100000000/1000000;             }//abs pad -> sector
+  static        Int_t      Pad2PadX      (Int_t pad                      ){return pad%1000000/1000;                  }//abs pad -> pad x 
+  static        Int_t      Pad2PadY      (Int_t pad                      ){return pad%1000000%100;                   }//abs pad -> pad y
+  static        Int_t      PadAbs        (Int_t c,Int_t s,Int_t x,Int_t y){return 100000000*c+1000000*s+1000*x+y;    }//(c,s,x,y) -> abs pad
+  static inline TVector2   Pad2Loc       (Int_t pad                      );                                           //abs pad ->LORS
+  static inline TVector2   Pad2Loc       (TVector pad                    );                                           //pad  -> LORS returns center of the pad
+  static        TVector2   Pad2Loc       (Int_t x,Int_t y                ){TVector pad(2);pad[0]=x;pad[1]=y;return Pad2Loc(pad);}//return center of the pad (x,y)
+  static inline TVector    Loc2Area      (const TVector2 &x2             );                                           //pads area affected by hit x2. area is LeftDown-RightUp pad numbers
+  static inline Int_t      Loc2Sec       (const TVector2 &x2             );                                           //LORS -> sector
+  static        Int_t      Loc2Sec       (Double_t x,Double_t y          ){return Loc2Sec(TVector2(x,y));}            //LORS -> sector
+  static inline TVector    Loc2Pad       (const TVector2 &x2             );                                           //LORS -> pad
+  static        TVector    Loc2Pad       (Double_t x,Double_t y          ){return Loc2Pad(TVector2(x,y));}            //LORS -> pad
+  static inline Int_t      Pad2Sec       (const TVector &pad             );                                           //pad  -> sector
+  static inline Int_t      PadNeighbours (Int_t iPadX,Int_t iPadY,Int_t aListX[4],Int_t aListY[4]);                   //pad -> list of it neighbours
+  static        Bool_t     IsAccepted    (const TVector2 &x2             ){return ( x2.X()>=0 && x2.X()<=PcSizeX() && x2.Y()>=0 && x2.Y()<=PcSizeY() ) ? kTRUE:kFALSE;}
+//optical properties methodes  
+  static        Double_t   MeanCkovEnergy(                               ){return 6.766;}                 //mean Ckov energy according to the total trasmission curve
+  static        Float_t    PhotonEnergy  (Int_t i                        ){return 0.1*i+5.5;}                         //photon energy (eV) for i-th point
+  static        Float_t    AbsCH4        (Float_t ev                     );                                           //CH4 abs len (cm) 
+  static        Float_t    AbsGel        (Float_t                        ){return 500;}                               //Aerogel abs len (cm)
+  static        Float_t    RefIdxC6F14   (Float_t eV                     ){return eV*0.0172+1.177;}                   //Freon ref idx
+  static        Float_t    RefIdxCH4     (Float_t                        ){return 1.000444;}                          //Methane ref idx 
+  static        Float_t    RefIdxSiO2    (Float_t eV                     ){Float_t e1=10.666,e2=18.125,f1=46.411,f2= 228.71; return TMath::Sqrt(1.+f1/(e1*e1-eV*eV)+f2/(e2*e2-eV*eV));}//Quartz window ref index from TDR p.35
+  static        Float_t    RefIdxGel     (Float_t                        ){return 1.05;}                              //aerogel ref index 
+  static        Float_t    DenGel        (                               ){return (RefIdxGel(0)-1)/0.21;}             //aerogel density gr/cm^3 parametrization by E.Nappi
+
   
   static Double_t IonisationPotential()      {return 26.0e-9;}                            //for CH4 in GeV taken from ????
   static TVector2 MathiesonDelta()           {return TVector2(5*0.18,5*0.18);}            //area of 5 sigmas of Mathieson distribution (cm)
@@ -92,28 +117,6 @@ public:
   
   static Int_t    HV(Int_t sector)           {if (sector>=1 && sector <=6) return fgHV[sector-1];  else return -1;} //high voltage for this sector
   static void     SetHV(Int_t sector,Int_t hv){fgHV[sector-1]=hv;}  
-//optical properties methodes  
-  static Double_t MeanCkovEnergy()         {return 6.766;}                 //mean Ckov energy according to the total trasmission curve
-  static Float_t  PhotonEnergy(Int_t i)    {return 0.1*i+5.5;}             //photon energy (eV) for i-th point
-  static Float_t  AbsCH4(Float_t ev);                                      //CH4 absorption length (cm) for photon with given energy (eV)
-  static Float_t  AbsGel(Float_t)          {return 500;}                   //Aerogel absorption length (cm) for photon with given energy (eV)
-  static Float_t  RefIdxC6F14(Float_t eV)  {return eV*0.0172+1.177;}       //Freon ref index for photon with given energy (eV)
-  static Float_t  RefIdxCH4(Float_t)       {return 1.000444;}              //Methane ref index for photon with given energy (eV)
-  static Float_t  RefIdxSiO2(Float_t eV)   {Float_t e1=10.666,e2=18.125,f1=46.411,f2= 228.71;
-                                     return TMath::Sqrt(1.+f1/(e1*e1-eV*eV)+f2/(e2*e2-eV*eV));}//Quartz window ref index from TDR p.35
-  static Float_t  RefIdxGel(Float_t)       {return 1.05;}                  //aerogel ref index 
-  static Float_t  DenGel()                 {return (RefIdxGel(0)-1)/0.21;} //aerogel density gr/cm^3 parametrization by E.Nappi
-//trasformation methodes
-  inline static TVector  Loc2Area(const TVector2 &x2);                                             //return area affected by hit x2
-  inline static Int_t    Loc2Sec(const TVector2 &x2);                                              //return sector containing given position
-         static Int_t    Loc2Sec(Double_t x,Double_t y) {return Loc2Sec(TVector2(x,y));}           //return sector containing given position
-  inline static TVector  Loc2Pad(const TVector2 &x2);                                              //return pad containing given position
-         static TVector  Loc2Pad(Double_t x,Double_t y) {return Loc2Pad(TVector2(x,y));}           //return pad containing given position
-  inline static TVector2 Pad2Loc(TVector pad);                                                     //return center of the pad
-         static TVector2 Pad2Loc(Int_t x,Int_t y) {TVector pad(2);pad[0]=x;pad[1]=y;return Pad2Loc(pad);}//return center of the pad (x,y)
-  inline static Int_t    Pad2Sec(const TVector &pad);                                              //return sector of given pad
-  inline static Int_t    PadNeighbours(Int_t iPadX,Int_t iPadY,Int_t aListX[4],Int_t aListY[4]);   //number of neighbours for this pad
-         static Bool_t   IsAccepted(const TVector2 &x2) {return ( x2.X()>=0 && x2.X()<=PcSizeX() && x2.Y()>=0 && x2.Y()<=PcSizeY() ) ? kTRUE:kFALSE;}
 //charge response methodes  
   inline static Double_t Mathieson(Double_t x1,Double_t x2,Double_t y1,Double_t y2);               //Mathienson integral over given limits
   inline static Double_t GainSag(Double_t x,Int_t sector);                                         //gain variations in %
@@ -122,7 +125,7 @@ public:
                           if(fgIsWireSag) return QdcSlope(Loc2Sec(x2))*(1+GainSag(x2.X(),Loc2Sec(x2))/100);
                           else            return QdcSlope(Loc2Sec(x2));}
   inline static Double_t FracQdc(const TVector2 &x2,const TVector &pad);                           //charge fraction to pad from hit
-  inline static Int_t    TotQdc(TVector2 x2,Double_t eloss);                                       //total charge for hit eloss=0 for photons
+  inline static Int_t    TotQdc(TVector2 x2,Double_t eloss);                                       //total charge for Eloss (GeV) 0 for photons
   inline static Bool_t   IsOverTh(Int_t c,TVector pad,Double_t q);                                 //is QDC of the pad registered by FEE  
          static Int_t    NsigmaTh()                    {return fgNsigmaTh;}                        //
          static Float_t  SigmaThMean()                 {return fgSigmaThMean;}                     //QDC electronic noise mean
@@ -140,7 +143,7 @@ public:
          static Double_t SnellAngle(Float_t n1, Float_t n2, Float_t theta1);                              // Snell law
          static void     AnglesInDRS(Double_t trackTheta,Double_t trackPhi,Double_t thetaCerenkov,Double_t phiCerenkov,Double_t &tout,Double_t &pout);//It finds photon angles in 
                                                                                                                                                       //Detector Reference System
-  
+
   static Bool_t     fgIsAerogel;                            //aerogel geometry instead of normal RICH flag
 protected:
   static Bool_t     fgIsRadioSrc;                           //radioactive source instead of radiators flag
@@ -266,6 +269,17 @@ TVector2 AliRICHParam::Pad2Loc(TVector pad)
   return TVector2(x,y);
 }
 //__________________________________________________________________________________________________
+TVector2 AliRICHParam::Pad2Loc(Int_t pad)
+{
+//Converts absolute pad number to local position in LORS
+//LORS is a chamber  reference system with origin in left-down coner looking from IP
+//Arguments: pad- absolute pad number
+//  Returns: pad center position as TVector2 in PCRS  
+  TVector2 pos;
+  pos.Set((Pad2PadX(pad)-0.5)*PadSizeX() , (Pad2PadY(pad)-0.5)*PadSizeY());//set to sector LORS
+  return pos;
+}
+//__________________________________________________________________________________________________
 Double_t AliRICHParam::GainSag(Double_t x,Int_t sector)
 {
 //Returns % of gain variation due to wire sagita.
@@ -308,19 +322,21 @@ Double_t AliRICHParam::FracQdc(const TVector2 &x2,const TVector &pad)
   return (Loc2Sec(x2)!=Pad2Sec(pad)) ? 0:Mathieson(normXmin, normYmin, normXmax, normYmax);
 }
 //__________________________________________________________________________________________________
-Double_t AliRICHParam::Mathieson(Double_t xMin,Double_t yMin,Double_t xMax,Double_t yMax)
+Double_t AliRICHParam::Mathieson(Double_t x1,Double_t y1,Double_t x2,Double_t y2)
 {
-//All arguments are parametrised according to NIM A370(1988)602-603
-//Returns a charge fraction.   
+//This is the answer to electrostatic problem of charge distrubution in MWPC described elsewhere. (NIM A370(1988)602-603)
+//Arguments: x1- diff between center of distribution and left margin of interested pad divided by anod-cathode distance
+//           x2,y1,y2- analogically  
+//  Returns: a charge fraction [0-1].
   const Double_t kSqrtKx3=0.77459667;const Double_t kX2=0.962;const Double_t kX4=0.379;
   const Double_t kSqrtKy3=0.77459667;const Double_t kY2=0.962;const Double_t kY4=0.379;
 
-  Double_t ux1=kSqrtKx3*TMath::TanH(kX2*xMin);
-  Double_t ux2=kSqrtKx3*TMath::TanH(kX2*xMax);    
-  Double_t uy1=kSqrtKy3*TMath::TanH(kY2*yMin);
-  Double_t uy2=kSqrtKy3*TMath::TanH(kY2*yMax);
+  Double_t ux1=kSqrtKx3*TMath::TanH(kX2*x1);
+  Double_t ux2=kSqrtKx3*TMath::TanH(kX2*x2);
+  Double_t uy1=kSqrtKy3*TMath::TanH(kY2*y1);
+  Double_t uy2=kSqrtKy3*TMath::TanH(kY2*y2);
   return 4*kX4*(TMath::ATan(ux2)-TMath::ATan(ux1))*kY4*(TMath::ATan(uy2)-TMath::ATan(uy1));
-}  
+}
 //__________________________________________________________________________________________________
 TVector AliRICHParam::Loc2Area(const TVector2 &x2)
 {
@@ -341,6 +357,5 @@ Bool_t AliRICHParam::IsOverTh(Int_t ,TVector ,Double_t q)
 //Checks if the current q is over threshold and FEE will save this value to data concentrator.
   return (q>NsigmaTh()*(SigmaThMean()+(1.-2*gRandom->Rndm())*SigmaThSpread()));
 }
+//__________________________________________________________________________________________________
 #endif //AliRICHParam_h
-
-
index bc5074f..833e921 100644 (file)
  * provided "as is" without express or implied warranty.                  *
  **************************************************************************/
 
-#include "AliRICHReconstructor.h"
-
+#include "AliRICHReconstructor.h" //class header
+#include <AliRawReader.h>         //Reconstruct(...)
+#include "AliRICHv1.h"         //Reconstruct(...) 
+#include <AliRun.h>               //ConvertDigits uses gAlice
+#include <TMinuit.h>              //Dig2Clu()
 ClassImp(AliRICHReconstructor)
 
+void AliRICHReconstructor::Reconstruct(AliRunLoader *pAL,AliRawReader* pRR)const
+{
+  AliLoader *pRL=pAL->GetDetectorLoader("RICH");
+  AliRICH *pRich=(AliRICH*)pAL->GetAliRun()->GetDetector("RICH");
+  
+  AliRICHDigit dig; //tmp digit, raw digit will be converted to it
+  TClonesArray *pDigList=new TClonesArray("AliRICHDigit"); Int_t iDigCnt=0; //tmp list of digits for single chamber only
+  
+  Int_t iEvtN=0;
+  while(pRR->NextEvent()){//events loop
+    pAL->GetEvent(iEvtN++);
+    pRL->MakeTree("R");  pRich->MakeBranch("R");
+    
+    for(Int_t iChN=1;iChN<=7;iChN++){//chambers loop
+      pRR->Select(AliRICHDigit::kRichRawId,2*iChN-2,2*iChN-1);//select only DDL files for the current chamber      
+      UInt_t w32=0;
+      while(pRR->ReadNextInt(w32)){//raw records loop (in selected DDL files)
+        UInt_t ddl=pRR->GetDDLID(); //returns 0,1,2 ... 13
+        dig.Raw2Dig(ddl,w32);  
+        AliDebug(1,Form("Ch=%i DDL=%i raw=0x%x digit=(%3i,%3i,%3i,%3i) Q=%5.2f",iChN,ddl,w32,dig.Chamber(),dig.Sector(),dig.PadX(),dig.PadY(),dig.Qdc()));
+        new((*pDigList)[iDigCnt++]) AliRICHDigit(dig); //add this digit to the tmp list
+      }//raw records loop
+      if(iDigCnt) Dig2Clu(pDigList,pRich->Clusters(iChN));//cluster finder for the current chamber if any digits present
+      pRR->Reset();        
+      pDigList->Delete();  iDigCnt=0;//clean up list of digits for the current chamber
+    }//chambers loop
+    pRL->TreeR()->Fill();            //fill tree for current event
+    pRL->WriteRecPoints("OVERWRITE");//write out clusters for current event
+    pRich->ClustersReset();
+  }//events loop  
+  pRL->UnloadRecPoints();  
+}//Reconstruct raw data
+//__________________________________________________________________________________________________
+void AliRICHReconstructor::Dig2Clu(TClonesArray *pDigList,TClonesArray *pCluList)const
+{
+//Finds all clusters for a given digits list provided not empty. Currently digits list provided is a list of all digits for a single chamber.
+//Puts all found clusters in the given clusters list. 
+//Arguments: pDigList - list of digits provided not empty
+//  Returns: none    
+  TMatrixF digMap(1,AliRICHParam::NpadsX(),1,AliRICHParam::NpadsY());  digMap=(Float_t)-1; //digit map for one chamber reseted to -1
+  for(Int_t iDigN=0 ; iDigN < pDigList->GetEntriesFast() ; iDigN++){ //digits loop to fill digits map
+    AliRICHDigit *pDig= (AliRICHDigit*)pDigList->At(iDigN);
+    digMap( pDig->PadX(), pDig->PadY() )=iDigN;                     //(padx,pady) cell takes digit number
+  }
+  
+  AliRICHCluster clu; Int_t iCluCnt=0; //tmp cluster and cluster counter
+  
+  for(Int_t iDigN=0;iDigN<pDigList->GetEntriesFast();iDigN++){//digits loop    
+    AliRICHDigit *pDig=(AliRICHDigit*)pDigList->At(iDigN);
+    if(!(pDig=UseDig(pDig->PadX(),pDig->PadY(),pDigList,&digMap))) continue;  //this digit is already taken in FormCluster(), go after next digit
+    FormCluster(&clu,pDig,pDigList,&digMap);                            //form cluster starting from this digit
+    TMinuit *pMinuit=clu.Solve();                              //solve this cluster
+    
+    if(pMinuit){//means cluster is solved into local maxima number of clusters, so add all of them in loop
+      Double_t x=-1,y=-1,q=-1;TString str; Double_t b1,b2,b3; Int_t iErrFlg;//tmp to withdraw resulting parameters
+      for(Int_t i=0;i<clu.Nlocmax();i++){//take resulting parameters 
+        pMinuit->mnpout(3*i  ,str, x, b1, b2, b3, iErrFlg);
+        pMinuit->mnpout(3*i+1,str, y, b1, b2, b3, iErrFlg);
+        pMinuit->mnpout(3*i+2,str, q, b1, b2, b3, iErrFlg);
+        clu.Set(x,y,(Int_t)q);
+        new((*pCluList)[iCluCnt++]) AliRICHCluster(clu);
+      }
+      delete pMinuit;
+    }else//means cluster is solved as simple center of gravity cluster, add it
+      new((*pCluList)[iCluCnt++]) AliRICHCluster(clu); 
+    clu.Reset();//make current cluster empty
+  }//digits loop
+}//Dig2Clu()
+//__________________________________________________________________________________________________
+void  AliRICHReconstructor::FormCluster(AliRICHCluster *pClu,AliRICHDigit *pDig,TClonesArray *pDigList,TMatrixF *pDigMap)const
+{
+//Forms the initial cluster as a sum of all adjascent digits. Starts from the given digit
+//then calls itself recursevly  for all neighbours.
+//Arguments: pClu - pointer to cluster being formed
+//  Returns: none
+  pClu->AddDigit(pDig);//take this digit in cluster and mark it as taken
+
+  Int_t x[4],y[4];
+  
+  Int_t iNnei=AliRICHParam::PadNeighbours(pDig->PadX(),pDig->PadY(),x,y);//returns in x,y all possible neighbours of the given one
+  for (Int_t i=0;i<iNnei;i++)
+    if((pDig=UseDig(x[i],y[i],pDigList,pDigMap))) FormCluster(pClu,pDig,pDigList,pDigMap);    
+}//FormCluster()
index f3f4830..4d28e88 100644 (file)
@@ -3,22 +3,44 @@
 /* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
  * See cxx source for full Copyright notice                               */
 
-#include <AliReconstructor.h>
-#include "AliRICHTracker.h"
-#include "AliRICHClusterFinder.h"
+#include <AliReconstructor.h>       //base class
+#include "AliRICHTracker.h"         //CreateTracker()
+#include "AliRICHClusterFinder.h"   //Reconstruct()
+
+class AliRawReader;
+class TTree;
 
 class AliRICHReconstructor: public AliReconstructor 
 {
 public:
            AliRICHReconstructor(): AliReconstructor()              {}//default ctor
   virtual ~AliRICHReconstructor()                                  {}//dtor  
-  virtual AliTracker*  CreateTracker(AliRunLoader*           )const{return new AliRICHTracker;}                         //interface from AliReconstructor
-  virtual void         Reconstruct  (AliRunLoader* pRunLoader)const{AliRICHClusterFinder clus(pRunLoader); clus.Exec();}//interface from AliReconstructor
-  using AliReconstructor::Reconstruct;                                                                                  //to get rid of virtual hidden warning 
-//  virtual void         FillESD(AliRunLoader* pAL, AliESD* pESD) const;    //virtual from AliReconstructor
-//  using AliReconstructor::FillESD;                                        //to get rid of virtual hidden warning 
+//framework part  
+  AliTracker*  CreateTracker         (AliRunLoader*                      )const{return new AliRICHTracker;}                 //interface from AliReconstructor
+  void         Reconstruct           (AliRunLoader* pAL                  )const{AliRICHClusterFinder cf(pAL);    cf.Exec();}//from AliReconstruction for digits->clusters
+  void         Reconstruct           (AliRunLoader* pAL,AliRawReader *pRR)const;                                            //from AliReconstruction for raw->clusters
+  using AliReconstructor::Reconstruct;                                                                            //to get rid of virtual hidden warning 
+//private part  
+         void          Dig2Clu   (TClonesArray*pDigList,TClonesArray *pCluList                                     )const;//form clusters out of provided digits list
+         void          FormCluster(AliRICHCluster *pClu,AliRICHDigit *pDig,TClonesArray *pDigList,TMatrixF *pDigMap)const;//form cluster recursive algorithm
+  inline AliRICHDigit *UseDig    (Int_t padX,Int_t padY,TClonesArray *pDigList,TMatrixF *pDigMap                   )const;//use this pad's digit to form a cluster
 protected:
   ClassDef(AliRICHReconstructor, 0)   //class for the RICH reconstruction
 };
+//__________________________________________________________________________________________________
+AliRICHDigit* AliRICHReconstructor::UseDig(Int_t padX,Int_t padY,TClonesArray *pDigList,TMatrixF *pDigMap)const
+{
+//Digit map contains a matrix if digit numbers.
+//Main operation in forming initial cluster is done here. Requested digit pointer is returned and this digit marked as taken.
+//Arguments: padX,padY - pad number
+//           pDigList  - list of digits for one sector
+//           pDigMap   - map of those digits
+//  Returns: pointer to digit if not yet used or 0 if used
+  Int_t iDig=(Int_t)(*pDigMap)(padX,padY);(*pDigMap)(padX,padY)=-1;//take digit number from the map and reset this map cell to -1
+  if(iDig!=-1)
+    return (AliRICHDigit*)pDigList->At(iDig);    //digit pointer
+  else      
+    return 0;
+}
 
 #endif
index ce4a41a..cfe2ed6 100644 (file)
 // **************************************************************************
 
 
-#include "AliRICHv1.h"
+#include "AliRICHv1.h"     //class header
 #include "AliRICHParam.h"
 #include "AliRICHChamber.h"
 #include <TParticle.h> 
 #include <TRandom.h> 
 #include <TVirtualMC.h>
 #include <TPDGCode.h>
-
+#include <AliStack.h>      //Hits2SDigits()
 #include <AliConst.h>
 #include <AliPDG.h>
 #include <AliRun.h>
 #include <AliMC.h>
+#include <AliRawDataHeader.h> //Digits2Raw()
 
 ClassImp(AliRICHv1)    
 //__________________________________________________________________________________________________
@@ -180,3 +181,89 @@ void AliRICHv1::GenerateFeedbacks(Int_t iChamber,Float_t eloss)
   }//feedbacks loop
   AliDebug(1,"Stop.");
 }//GenerateFeedbacks()
+//__________________________________________________________________________________________________
+void AliRICHv1::Hits2SDigits()
+{
+// Create a list of sdigits corresponding to list of hits. Every hit generates one or more sdigits.
+  AliDebug(1,"Start.");
+  for(Int_t iEventN=0;iEventN<GetLoader()->GetRunLoader()->GetAliRun()->GetEventsPerRun();iEventN++){//events loop
+    GetLoader()->GetRunLoader()->GetEvent(iEventN);//get next event
+  
+    if(!GetLoader()->TreeH()) GetLoader()->LoadHits();    GetLoader()->GetRunLoader()->LoadHeader(); 
+    if(!GetLoader()->GetRunLoader()->TreeK())             GetLoader()->GetRunLoader()->LoadKinematics();//from
+    if(!GetLoader()->TreeS()) GetLoader()->MakeTree("S"); MakeBranch("S");//to
+          
+    for(Int_t iPrimN=0;iPrimN<GetLoader()->TreeH()->GetEntries();iPrimN++){//prims loop
+      GetLoader()->TreeH()->GetEntry(iPrimN);
+      for(Int_t iHitN=0;iHitN<Hits()->GetEntries();iHitN++){//hits loop 
+        AliRICHHit *pHit=(AliRICHHit*)Hits()->At(iHitN);//get current hit                
+        TVector2 x2 = C(pHit->C())->Mrs2Anod(0.5*(pHit->InX3()+pHit->OutX3()));//hit position in the anod plane
+        Int_t iTotQdc=P()->TotQdc(x2,pHit->Eloss());//total charge produced by hit, 0 if hit in dead zone
+        if(iTotQdc==0) continue;
+        //
+        //need to quantize the anod....
+        TVector padHit=AliRICHParam::Loc2Pad(x2);
+        TVector2 padHitXY=AliRICHParam::Pad2Loc(padHit);
+        TVector2 anod;
+        if((x2.Y()-padHitXY.Y())>0) anod.Set(x2.X(),padHitXY.Y()+AliRICHParam::PitchAnod()/2);
+        else anod.Set(x2.X(),padHitXY.Y()-AliRICHParam::PitchAnod()/2);
+        //end to quantize anod
+        //
+        TVector area=P()->Loc2Area(anod);//determine affected pads, dead zones analysed inside
+        AliDebug(1,Form("hitanod(%6.2f,%6.2f)->area(%3.0f,%3.0f)-(%3.0f,%3.0f) QDC=%4i",anod.X(),anod.Y(),area[0],area[1],area[2],area[3],iTotQdc));
+        TVector pad(2);
+        for(pad[1]=area[1];pad[1]<=area[3];pad[1]++)//affected pads loop
+          for(pad[0]=area[0];pad[0]<=area[2];pad[0]++){                    
+            Double_t padQdc=iTotQdc*P()->FracQdc(anod,pad);
+            AliDebug(1,Form("current pad(%3.0f,%3.0f) with QDC  =%6.2f",pad[0],pad[1],padQdc));
+            if(padQdc>0.1) SDigitAdd(pHit->C(),pad,padQdc,GetLoader()->GetRunLoader()->Stack()->Particle(pHit->GetTrack())->GetPdgCode(),pHit->GetTrack());
+          }//affected pads loop 
+      }//hits loop
+    }//prims loop
+    GetLoader()->TreeS()->Fill();
+    GetLoader()->WriteSDigits("OVERWRITE");
+    SDigitsReset();
+  }//events loop  
+  GetLoader()->UnloadHits(); GetLoader()->GetRunLoader()->UnloadHeader(); GetLoader()->GetRunLoader()->UnloadKinematics();
+  GetLoader()->UnloadSDigits();  
+  AliDebug(1,"Stop.");
+}//Hits2SDigits()
+//__________________________________________________________________________________________________
+void AliRICHv1::Digits2Raw()
+{
+//Creates raw data files in DDL format. Invoked by AliSimulation
+//loop over events is done outside in AliSimulation
+//Arguments: none
+//  Returns: none    
+  AliDebug(1,"Start.");
+  GetLoader()->LoadDigits();
+  GetLoader()->TreeD()->GetEntry(0);
+  
+  ofstream file[AliRICHDigit::kNddls];   //output streams
+  Int_t    cnt[AliRICHDigit::kNddls];        //data words counters for DDLs
+  AliRawDataHeader header; //empty DDL header
+  UInt_t w32=0;            //32 bits data word 
+  
+  for(Int_t i=0;i<AliRICHDigit::kNddls;i++){//open all 14 DDL in parallel
+    file[i].open(Form("RICH_%4i.ddl",AliRICHDigit::kDdlOffset+i));  //first is 0x700
+    file[i].write((char*)&header,sizeof(header));                //write dummy header as place holder, actual will be written later when total size of DDL is known
+    cnt[i]=0; //reset counters
+  }
+  
+  for(Int_t iChN=1;iChN<=kNchambers;iChN++){ //digits are stored on chamber by chamber basis   
+    for(Int_t iDigN=0;iDigN<Digits(iChN)->GetEntriesFast();iDigN++){//digits loop for a given chamber
+      AliRICHDigit *pDig=(AliRICHDigit*)Digits(iChN)->At(iDigN);
+      Int_t ddl=pDig->Dig2Raw(w32);  //ddl is 0..13 
+      file[ddl].write((char*)&w32,sizeof(w32));  cnt[ddl]++;//write formated digit to the propriate file (as decided in Dig2Raw) and increment corresponding counter
+    }//digits loop for a given chamber
+  }//chambers loop    
+  for(Int_t i=0;i<AliRICHDigit::kNddls;i++){
+    header.fSize=sizeof(header)+cnt[i]*sizeof(w32); //now calculate total number of bytes for each DDL file  
+    header.SetAttribute(0); 
+    file[i].seekp(0); file[i].write((char*)&header,sizeof(header));//rewrite DDL header with fSize field properly set
+    file[i].close();                                               //close DDL file
+  }
+  GetLoader()->UnloadDigits();
+  AliDebug(1,"Stop.");      
+}//Digits2Raw()
+//__________________________________________________________________________________________________
index 1457920..8065a36 100644 (file)
@@ -4,7 +4,8 @@
 /* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
  * See cxx source for full Copyright notice                               */
 
-#include "AliRICH.h"
+#include "AliRICH.h"             //base class 
+#include "AliRICHDigitizer.h"  //CreateDigitizer()
 
 class AliRICHv1 : public AliRICH 
 {
@@ -12,11 +13,16 @@ public:
                  AliRICHv1():AliRICH()                                               {;}          //default ctor
                  AliRICHv1(const char *name, const char *title):AliRICH(name,title)  {;}          //named ctor
   virtual       ~AliRICHv1()                                                         {;}          //dtor
-  virtual void   Init()                                                              {;}
-  virtual Int_t  IsVersion()                                                    const{return 1;}  
-  virtual void   StepManager();                                                                   //full slow step manager
-          Bool_t IsLostByFresnel();                                                               //checks if the photon lost on Fresnel reflection  
-          void   GenerateFeedbacks(Int_t iChamber,Float_t eloss=0);                               //generates feedback photons; eloss=0 for photon
+//framework part  
+  void   Init()                                                              {;}
+  Int_t          IsVersion        (                          )const{return 1;}  
+  void           StepManager      (                          );                                           //called from AliSimulation or AliRun when transport particles
+  void           Hits2SDigits     (                          );                                           //called from AliSimulation for Hits->SDigits
+  AliDigitizer*  CreateDigitizer  (AliRunDigitizer *pMan     )const{return new AliRICHDigitizer(pMan);}   //called from AliSimulation for SDigits->Digits
+  void           Digits2Raw       (                          );                                           //called from AliSimulation for Digits->Raw
+//private part  
+          Bool_t IsLostByFresnel  ();                                                                     //checks if the photon lost on Fresnel reflection  
+          void   GenerateFeedbacks(Int_t iChamber,Float_t eloss=0);                                       //generates feedback photons; eloss=0 for photon
 protected:
   ClassDef(AliRICHv1,1)                                                 //RICH full version for simulation
 };
index 545324f..f4e5e42 100644 (file)
@@ -8,7 +8,6 @@
 #pragma link C++ class  AliRICHHit+;
 #pragma link C++ class  AliRICHDigit+;
 #pragma link C++ class  AliRICHCluster+;
-#pragma link C++ class  AliRICHDigitizer+;
 #pragma link C++ class  AliRICHTracker+;
 #pragma link C++ class  AliRICHRecon+;
 #pragma link C++ class  AliRICHHelix+;
index 7d4450f..0e87282 100644 (file)
@@ -4,6 +4,5 @@
 #pragma link off all functions;
 #pragma link C++ class  AliRICHv0+;
 #pragma link C++ class  AliRICHv1+;
-#pragma link C++ class  AliGenRadioactive+;
-#pragma link C++ enum   RadioNuclei_t;
+#pragma link C++ class  AliRICHDigitizer+;
 #endif
index 4ef4223..1512a42 100644 (file)
@@ -1,3 +1,5 @@
+#include <TPDGCode.h>
+
 class RichConfig : public TGMainFrame
 {
 RQ_OBJECT()
@@ -22,8 +24,9 @@ public:
   void    VerSlot(Int_t ver)             {if(ver==kTestBeam){ fMagCB->SetState(kButtonUp); fGenTypeCO->Select(kGunAlongZ);}}        
   
   Float_t Eta2Theta(Float_t arg)    const{return (180./TMath::Pi())*2.*TMath::ATan(TMath::Exp(-arg));}
-  void    CreateConfig();
-  void    CreateRichBatch();
+  void    CreateConfig();      //create Config.C
+  void    Generator(FILE *pF); //write generatro part of Config.C
+  void    CreateRichBatch();   //create RichBatch.C
   void    Exit();
      
   TGComboBox   *fRichCO;  TGButtonGroup *fRichBG;                                                                              //RICH
@@ -87,7 +90,6 @@ RichConfig::RichConfig(const char *sFileName):TGMainFrame(gClient->GetRoot(),700
   fGenTypeCO->AddEntry("HIJING"                       ,kHijing);
   fGenTypeCO->AddEntry("HIJING para"                  ,kHijingPara);
   fGenTypeCO->AddEntry("2 p+HIJING"                   ,kHijingPara2Proton);
-  fGenTypeCO->AddEntry("Sr90 source"                  ,kSr90);
   fGenTypeCO->AddEntry("RICH lib"                     ,kRichLib);
   fGenTypeCO->AddEntry("RICH lib+HIJING"              ,kSignalHijing);
   fGenTypeCO->Select(kHijing);
@@ -236,7 +238,66 @@ void RichConfig::CreateConfig()
   fprintf(fp,"  pDecayer->SetForceDecay(kAll);\n"); 
   fprintf(fp,"  pDecayer->Init();\n"); 
   fprintf(fp,"  gMC->SetExternalDecayer(pDecayer);\n\n");
-//Physics
+  Physics(fp);
+//Field
+  if(fMagCB->GetState()==kButtonDown) fprintf(fp,"  gAlice->SetField();\n\n");else fprintf(fp,"  gAlice->SetField(0);\n\n");
+  fprintf(fp,"  pAL->CdGAFile();\n\n");                                 //????       
+//BODY-ALIC 
+  fprintf(fp,"  new AliBODY(\"BODY\",\"Alice envelop\");\n\n");
+//RICH
+  if(!fRichBG->GetButton(kDeclust)->GetState()) fprintf(fp,"  AliRICHParam::SetDeclustering(kFALSE);\n");
+  if(!fRichBG->GetButton(kSagita)->GetState())  fprintf(fp,"  AliRICHParam::SetWireSag(kFALSE);\n")     ;
+  switch(fRichCO->GetSelected()){
+    case kNo:                                                                                                  break;
+    case kVer0:  fprintf(fp,"  AliRICH *pRICH=new AliRICHv0(\"RICH\",\"Rich with debug StepManager\");\n\n"); break;      
+    case kVer1:  fprintf(fp,"  AliRICH *pRICH=new AliRICHv1(\"RICH\",\"Normal\");\n\n");                 break;      
+  }//switch RICH
+//Generator
+  Generator(fp);
+//central before RICH detectors                  
+  if(fDetBG->GetButton(kPIPE )->GetState()) fprintf(fp,"\n  new AliPIPEv0(\"PIPE\",\"Beam Pipe\");\n");
+  if(fDetBG->GetButton(kSHILD)->GetState()) fprintf(fp,"\n  new AliSHILv2(\"SHIL\",\"Shielding Version 2\");\n");  
+  if(fDetBG->GetButton(kITS  )->GetState()){
+    fprintf(fp,"\n  AliITSvPPRasymmFMD *pIts =new AliITSvPPRasymmFMD(\"ITS\",\"ITS PPR detailed version\");\n");
+    fprintf(fp,"  pIts->SetMinorVersion(2); pIts->SetReadDet(kTRUE);\n");
+    fprintf(fp,"  pIts->SetThicknessDet1(200.); pIts->SetThicknessDet2(200.);\n");
+    fprintf(fp,"  pIts->SetThicknessChip1(200.); pIts->SetThicknessChip2(200.);\n");
+    fprintf(fp,"  pIts->SetRails(0); pIts->SetCoolingFluid(1);\n");
+    fprintf(fp,"  pIts->SetEUCLID(0);\n");
+  }  
+  if(fDetBG->GetButton(kTPC  )->GetState()) {fprintf(fp,"\n  AliTPC *pTpc=new AliTPCv2(\"TPC\",\"Default\"); pTpc->SetSecAU(-1);pTpc->SetSecAL(-1);\n");}  
+  if(fDetBG->GetButton(kFRAME)->GetState())  fprintf(fp,"\n  AliFRAMEv2 *pFrame=new AliFRAMEv2(\"FRAME\",\"Space Frame\"); pFrame->SetHoles(1);\n");
+  if(fDetBG->GetButton(kTRD  )->GetState()) {
+    fprintf(fp,"\n  AliTRD *pTrd=new AliTRDv1(\"TRD\",\"TRD slow simulator\");\n");
+    fprintf(fp,"  pTrd->SetGasMix(1); pTrd->SetPHOShole();pTrd->CreateTR();\n");
+  }  
+  if(fDetBG->GetButton(kTOF  )->GetState()) fprintf(fp,"\n  new AliTOFv4T0(\"TOF\", \"normal TOF\");\n");
+//central after RICH detectors  
+  if(fDetBG->GetButton(kMAG  )->GetState()) fprintf(fp,"\n  new AliMAG(\"MAG\",\"Magnet\");\n");  
+  if(fDetBG->GetButton(kHALL )->GetState()) fprintf(fp,"\n  new AliHALL(\"HALL\",\"Alice Hall\");\n");
+//forward detectors  
+  if(fDetBG->GetButton(kFMD  )->GetState()) fprintf(fp,"\n  new AliFMDv1(\"FMD\",\"normal FMD\");\n");
+  if(fDetBG->GetButton(kABSO )->GetState()) fprintf(fp,"\n  new AliABSOv0(\"ABSO\",\"Muon absorber\");\n");
+  if(fDetBG->GetButton(kDIPO )->GetState()) fprintf(fp,"\n  new AliDIPOv2(\"DIPO\",\"Dipole version 2\");\n");
+  if(fDetBG->GetButton(kMUON )->GetState()) fprintf(fp,"\n  new AliMUONv1(\"MUON\",\"default\");\n");
+  if(fDetBG->GetButton(kPMD  )->GetState()) fprintf(fp,"\n  new AliPMDv1(\"PMD\",\"normal PMD\");\n");
+  if(fDetBG->GetButton(kSTART)->GetState()) fprintf(fp,"\n  new AliSTARTv1(\"START\",\"START Detector\");\n");
+  if(fDetBG->GetButton(kVZERO)->GetState()) fprintf(fp,"\n  new AliVZEROv2(\"VZERO\",\"normal VZERO\");\n");
+  if(fDetBG->GetButton(kZDC  )->GetState()) fprintf(fp,"\n  new AliZDCv2(\"ZDC\",\"normal ZDC\");\n");
+//different phase space detectors  
+  if(fDetBG->GetButton(kPHOS )->GetState()) fprintf(fp,"\n  new AliPHOSv1(\"PHOS\",\"IHEP\");\n");
+  if(fDetBG->GetButton(kEMCAL)->GetState()) fprintf(fp,"\n  new AliEMCALv1(\"EMCAL\",\"G56_2_55_19_104_14\");\n");
+  if(fDetBG->GetButton(kCRT  )->GetState()) fprintf(fp,"\n  new AliCRTv0(\"CRT\",\"normal ACORDE\");\n");
+
+  fprintf(fp,"\n  ::Info(\"RICH private config\",\"Stop\");\n"); 
+  fprintf(fp,"}\n");
+out:  
+  fclose(fp);  
+//  CloseWindow();
+}//CreateConfig
+//__________________________________________________________________________________________________
+void RichConfig::Physics(FILE *fp)
+{
   if(fProcBG->GetButton(kDCAY)->GetState()) fprintf(fp,"  gMC->SetProcess(\"DCAY\",1);");  else fprintf(fp,"  gMC->SetProcess(\"DCAY\",0);");
   if(fProcBG->GetButton(kPAIR)->GetState()) fprintf(fp,"  gMC->SetProcess(\"PAIR\",1);\n");else fprintf(fp,"  gMC->SetProcess(\"PAIR\",0);\n");
   if(fProcBG->GetButton(kCOMP)->GetState()) fprintf(fp,"  gMC->SetProcess(\"COMP\",1);");  else fprintf(fp,"  gMC->SetProcess(\"COMP\",0);");
@@ -259,20 +320,10 @@ void RichConfig::CreateConfig()
   fprintf(fp,"  gMC->SetCut(\"BCUTM\" ,0.001); ");  fprintf(fp,"  gMC->SetCut(\"DCUTE\" ,0.001); ");
   fprintf(fp,"  gMC->SetCut(\"DCUTM\" ,0.001);\n"); fprintf(fp,"  gMC->SetCut(\"PPCUTM\",0.001); ");
   fprintf(fp,"  gMC->SetCut(\"TOFMAX\",1e10);\n\n"); 
-//Field
-  if(fMagCB->GetState()==kButtonDown) fprintf(fp,"  gAlice->SetField();\n\n");else fprintf(fp,"  gAlice->SetField(0);\n\n");
-  fprintf(fp,"  pAL->CdGAFile();\n\n");                                 //????       
-//BODY-ALIC 
-  fprintf(fp,"  new AliBODY(\"BODY\",\"Alice envelop\");\n\n");
-//RICH
-  if(!fRichBG->GetButton(kDeclust)->GetState()) fprintf(fp,"  AliRICHParam::SetDeclustering(kFALSE);\n");
-  if(!fRichBG->GetButton(kSagita)->GetState())  fprintf(fp,"  AliRICHParam::SetWireSag(kFALSE);\n")     ;
-  switch(fRichCO->GetSelected()){
-    case kNo:                                                                                                  break;
-    case kVer0:  fprintf(fp,"  AliRICH *pRICH=new AliRICHv0(\"RICH\",\"Rich with debug StepManager\");\n\n"); break;      
-    case kVer1:  fprintf(fp,"  AliRICH *pRICH=new AliRICHv1(\"RICH\",\"Normal\");\n\n");                 break;      
-  }//switch RICH
-//Generator
+}//Physics()
+//__________________________________________________________________________________________________
+void RichConfig::Generator(FILE *fp)
+{
   switch(fGenTypeCO->GetSelected()){
     case kHijingPara: 
       fprintf(fp,"  AliGenHIJINGpara *pGen=new AliGenHIJINGpara(%i);\n",(int)fGenPrimNE->GetNumber());
@@ -326,11 +377,6 @@ void RichConfig::CreateConfig()
       fprintf(fp,"  pPythia->SetProcess(kPyMb);  pPythia->SetEnergyCMS(14000);\n");      
       fprintf(fp,"  pCocktail->AddGenerator(pPythia,\"Pythia\",1);\n");  
       fprintf(fp,"  pCocktail->Init();\n");
-    case kSr90:  
-      fprintf(fp,"  AliGenRadioactive *pGen=new AliGenRadioactive(kSr90,%i);\n",(int)fGenNprimEntry->GetNumber());
-      fprintf(fp,"  pGen->SetOrigin(0,0,0);  pGen->SetSigma(0.1,0.1,0.05);\n");      
-      fprintf(fp,"  pGen->Init();\n");
-    break;  
     case kHijing:  
       fprintf(fp,"  AliGenHijing *pGen=new AliGenHijing(-1); pGen->SetEnergyCMS(5500); pGen->SetReferenceFrame(\"CMS\");\n");
       fprintf(fp,"  pGen->SetProjectile(\"A\", 208, 82); pGen->SetTarget(\"A\", 208, 82);\n");      
@@ -384,47 +430,9 @@ void RichConfig::CreateConfig()
       fprintf(fp,"  pCocktail->Init();\n");
     break;
   }
-//central before RICH detectors                  
-  if(fDetBG->GetButton(kPIPE )->GetState()) fprintf(fp,"\n  new AliPIPEv0(\"PIPE\",\"Beam Pipe\");\n");
-  if(fDetBG->GetButton(kSHILD)->GetState()) fprintf(fp,"\n  new AliSHILv2(\"SHIL\",\"Shielding Version 2\");\n");  
-  if(fDetBG->GetButton(kITS  )->GetState()){
-    fprintf(fp,"\n  AliITSvPPRasymmFMD *pIts =new AliITSvPPRasymmFMD(\"ITS\",\"ITS PPR detailed version\");\n");
-    fprintf(fp,"  pIts->SetMinorVersion(2); pIts->SetReadDet(kTRUE);\n");
-    fprintf(fp,"  pIts->SetThicknessDet1(200.); pIts->SetThicknessDet2(200.);\n");
-    fprintf(fp,"  pIts->SetThicknessChip1(200.); pIts->SetThicknessChip2(200.);\n");
-    fprintf(fp,"  pIts->SetRails(0); pIts->SetCoolingFluid(1);\n");
-    fprintf(fp,"  pIts->SetEUCLID(0);\n");
-  }  
-  if(fDetBG->GetButton(kTPC  )->GetState()) {fprintf(fp,"\n  AliTPC *pTpc=new AliTPCv2(\"TPC\",\"Default\"); pTpc->SetSecAU(-1);pTpc->SetSecAL(-1);\n");}  
-  if(fDetBG->GetButton(kFRAME)->GetState())  fprintf(fp,"\n  AliFRAMEv2 *pFrame=new AliFRAMEv2(\"FRAME\",\"Space Frame\"); pFrame->SetHoles(1);\n");
-  if(fDetBG->GetButton(kTRD  )->GetState()) {
-    fprintf(fp,"\n  AliTRD *pTrd=new AliTRDv1(\"TRD\",\"TRD slow simulator\");\n");
-    fprintf(fp,"  pTrd->SetGasMix(1); pTrd->SetPHOShole();pTrd->CreateTR();\n");
-  }  
-  if(fDetBG->GetButton(kTOF  )->GetState()) fprintf(fp,"\n  new AliTOFv4T0(\"TOF\", \"normal TOF\");\n");
-//central after RICH detectors  
-  if(fDetBG->GetButton(kMAG  )->GetState()) fprintf(fp,"\n  new AliMAG(\"MAG\",\"Magnet\");\n");  
-  if(fDetBG->GetButton(kHALL )->GetState()) fprintf(fp,"\n  new AliHALL(\"HALL\",\"Alice Hall\");\n");
-//forward detectors  
-  if(fDetBG->GetButton(kFMD  )->GetState()) fprintf(fp,"\n  new AliFMDv1(\"FMD\",\"normal FMD\");\n");
-  if(fDetBG->GetButton(kABSO )->GetState()) fprintf(fp,"\n  new AliABSOv0(\"ABSO\",\"Muon absorber\");\n");
-  if(fDetBG->GetButton(kDIPO )->GetState()) fprintf(fp,"\n  new AliDIPOv2(\"DIPO\",\"Dipole version 2\");\n");
-  if(fDetBG->GetButton(kMUON )->GetState()) fprintf(fp,"\n  new AliMUONv1(\"MUON\",\"default\");\n");
-  if(fDetBG->GetButton(kPMD  )->GetState()) fprintf(fp,"\n  new AliPMDv1(\"PMD\",\"normal PMD\");\n");
-  if(fDetBG->GetButton(kSTART)->GetState()) fprintf(fp,"\n  new AliSTARTv1(\"START\",\"START Detector\");\n");
-  if(fDetBG->GetButton(kVZERO)->GetState()) fprintf(fp,"\n  new AliVZEROv2(\"VZERO\",\"normal VZERO\");\n");
-  if(fDetBG->GetButton(kZDC  )->GetState()) fprintf(fp,"\n  new AliZDCv2(\"ZDC\",\"normal ZDC\");\n");
-//different phase space detectors  
-  if(fDetBG->GetButton(kPHOS )->GetState()) fprintf(fp,"\n  new AliPHOSv1(\"PHOS\",\"IHEP\");\n");
-  if(fDetBG->GetButton(kEMCAL)->GetState()) fprintf(fp,"\n  new AliEMCALv1(\"EMCAL\",\"G56_2_55_19_104_14\");\n");
-  if(fDetBG->GetButton(kCRT  )->GetState()) fprintf(fp,"\n  new AliCRTv0(\"CRT\",\"normal ACORDE\");\n");
+}//Generator()
+
 
-  fprintf(fp,"\n  ::Info(\"RICH private config\",\"Stop\");\n"); 
-  fprintf(fp,"}\n");
-out:  
-  fclose(fp);  
-//  CloseWindow();
-}//CreateConfig
 //__________________________________________________________________________________________________
 void RichConfig::CreateRichBatch()
 {//creates RichBatch.C file
index 6c35345..e19bc30 100644 (file)
@@ -10,7 +10,7 @@ include       lib$(Module)rec.pkg
 SrcRec         :=$(SRCS)
 
 RootTarget      :=$(shell root-config --arch)
-DirOut         :=/tmp/$(Module)
+DirOut         :=/tmp/$(Module)_$(USER)
 LibBase                :=$(LIB_MY)/lib$(Module)base.so
 LibSim         :=$(LIB_MY)/lib$(Module)sim.so
 LibRec         :=$(LIB_MY)/lib$(Module)rec.so
index 5f1a697..1254a3a 100644 (file)
@@ -63,14 +63,6 @@ What are the meanings of different VMC flags:
        gMC->IsTrackAlive()
        gMC->IsTrackStop()
        gMC->IsTrackDisappeared()
-How is sdigit produced from hit:
-       Responsible method is AliRICH::Hits2SDigits
-       One hit may affect one or more pads.
-        Hit position is taken on the anode wires plane as the most of avalanche is developed there.
-        This position is not directly available, track intersections with entrance and exit of amplification gap are only stored.
-        So the position in the middle of the gap is calculated as average out of pHit->In() and pHit->Out() positions.
-        Then, total charge collected for this hit is calculated by AliRICHParam::Hit2Qdc.    
-        Area of disintegration is a list of pads affected by current hit. This is a parameter of Mathienson    
 How to get pad number for a local position:
        use static TVector AliRICHParam::Loc2Pad(TVector2 position);
 Why list of chambers belongs to AliRICHParam:
@@ -87,3 +79,70 @@ How to loop over all possible object:
       }//hits loop
     }//TreeH loop
   }//events loop
+
+
+RICH full simulation-reconstruction sequence
+
+hits->sdigit:
+       Responsible method is AliRICH::Hits2SDigits
+       One hit may affect one or more pads.
+        Hit position is taken on the anode wires plane as the most of avalanche is developed there.
+        This position is not directly available, track intersections with entrance and exit of amplification gap are only stored.
+        So the position in the middle of the gap is calculated as average out of pHit->In() and pHit->Out() positions.
+        Then, total charge collected for this hit is calculated by AliRICHParam::Hit2Qdc.    
+        Area of disintegration is a list of pads affected by current hit. This is a parameter of Mathienson    
+sdigits->digits:
+       The necessety of sdigits is dictated by the fact that trasport engine transports track by track. It means that it may happen that the same
+        pad is affected by few tracks. But this might be known after the trasport of full event only. 
+
+Generalized structure of AliReconstruction:
+
+Run()
+{
+  if(there is galice.root)                                          |
+    AliRunLoader::Open(....)                                        | this is done in InitRunLoader()
+  else                                                              |
+    if(raw data process requested)                                  |
+      create galice.root on the base of AliRawReader::NextEvent     |   
+  
+  for(all detectors){                                               |
+   if(detector not selected to run) skip this detector              |
+   reconstructor=get detector's reconstructor                       | this is done in RunLocalReconstruction()
+   if(detector HasLocalReconstruction) skip this detector           | IMPORTANT! if HasLocalReconstruction() returns YES use RunLocalEventReconstruction instead
+   if(run upon raw data)                                            |  
+     reconstructor->Reconstruct(fRunLoader, fRawReader);            |   
+   else                                                             |
+     reconstructor->Reconstruct(fRunLoader);                        |
+  }
+
+  for(all events){                                                      
+  
+    for(all detectors){                                                 |            
+      if(detector not selected to run) skip this detector               |
+      reconstructor=get detector's reconstructor                        | 
+      loader=get detector's loader                                      | this is done in RunLocalEventReconstruction()
+      if(raw data process requested and detector HasDigitConversion){   | 
+        loader->LoadDigits("update");                                   | open file and invoke  detector->SetTreeAddress();
+        loader->CleanDigits();                                          |   
+        loader->MakeDigitsContainer();                                  | create tree
+        reconstructor->Reconstruct(fRawReader,loader->TreeD());         | expected to fill TreeD out of raw reader
+        loader->WriteDigits("overwrite");                               | 
+        loader->UnloadDigits();                                         |  
+      }                                                                 |
+      if(detector do not HasLocalReconstruction) skip this detector     | IMPORTANT! assumed that this detector is already processed in RunLocalReconstruction()
+      loader->LoadRecPoints("update");                                  |
+      loader->CleanRecPoints();                                         |
+      loader->MakeRecPointsContainer();                                 | 
+      if(fRawReader && reconstructor do not HasDigitConversion()){      | 
+        reconstructor->Reconstruct(fRawReader, loader->TreeR());        | expected to fill TreeR out of raw reader
+      }else{                                                            |
+        loader->LoadDigits("read");                                     |
+        reconstructor->Reconstruct(loader->TreeD(),loader->TreeR());    | the only operations inside are pDigTree->GetEntry(0) and pCluTree->Fill();
+        loader->UnloadDigits();                                         | 
+      }                                                                 |
+      loader->WriteRecPoints("OVERWRITE");                              |
+      loader->UnloadRecPoints();                                        |
+    }//detectors loop                                                   |
+    
+  }//events loop 
+}
index b3c6a07..3c971d2 100644 (file)
@@ -1,4 +1,4 @@
-SRCS:= AliRICHParam.cxx AliRICHChamber.cxx AliRICHHit.cxx AliRICHDigit.cxx AliRICHCluster.cxx AliRICH.cxx AliRICHDigitizer.cxx  AliRICHTracker.cxx AliRICHRecon.cxx AliRICHHelix.cxx
+SRCS:= AliRICHParam.cxx AliRICHChamber.cxx AliRICHHit.cxx AliRICHDigit.cxx AliRICHCluster.cxx AliRICH.cxx AliRICHTracker.cxx AliRICHRecon.cxx AliRICHHelix.cxx
 
 HDRS:= $(SRCS:.cxx=.h)
 DHDR:= RICHbaseLinkDef.h
index 18a8ee1..57c0d65 100644 (file)
@@ -1,4 +1,4 @@
-SRCS:= AliRICHv0.cxx AliRICHv1.cxx AliGenRadioactive.cxx 
+SRCS:= AliRICHv0.cxx AliRICHv1.cxx AliRICHDigitizer.cxx
 
 HDRS:= $(SRCS:.cxx=.h)
 DHDR:= RICHsimLinkDef.h