Method to count the physical contribution to a cluster added.
authordibari <dibari@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 25 Nov 2003 15:16:49 +0000 (15:16 +0000)
committerdibari <dibari@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 25 Nov 2003 15:16:49 +0000 (15:16 +0000)
13 files changed:
RICH/AliRICH.cxx
RICH/AliRICH.h
RICH/AliRICHClusterFinder.cxx
RICH/AliRICHClusterFinder.h
RICH/AliRICHDisplay.cxx
RICH/AliRICHParam.cxx
RICH/AliRICHParam.h
RICH/AliRICHPoints.cxx
RICH/AliRICHv0.cxx
RICH/AliRICHv1.cxx
RICH/AliRICHv1.h
RICH/ConfigRich.C
RICH/menu.C

index ec90bf8..54e6305 100644 (file)
@@ -33,23 +33,23 @@ ClassImp(AliRICHhit)
 //__________________________________________________________________________________________________
 void AliRICHhit::Print(Option_t*)const
 {
-  Info("","chamber=%2i, PID=%9i, TID=%6i, eloss=%8.3f eV",fChamber,fPid,fTrack,fEloss*1e9);
+  ::Info("hit","chamber=%2i, PID=%9i, TID=%6i, eloss=%8.3f eV",fChamber,fPid,fTrack,fEloss*1e9);
 }
 //__________________________________________________________________________________________________
 ClassImp(AliRICHdigit)
 //__________________________________________________________________________________________________
 void AliRICHdigit::Print(Option_t*)const
 {
-  Info("","ID=%6i, chamber=%2i, PadX=%3i, PadY=%3i, Q=%8.2f, TID1=%5i, TID2=%5i, TID3=%5i",
-         Id(),fChamber,fPadX,fPadY,fQdc,fTracks[0],fTracks[1],fTracks[2]);
+  ::Info("digit","ID=%6i, PID=%9i, c=%2i, x=%3i, y=%3i, q=%8.2f, TID1=%5i, TID2=%5i, TID3=%5i",
+           Id(),   fCombiPid,fChamber,fPadX,fPadY,fQdc,fTracks[0],fTracks[1],fTracks[2]);
 }
 //__________________________________________________________________________________________________
 ClassImp(AliRICHcluster)
 //__________________________________________________________________________________________________
 void AliRICHcluster::Print(Option_t*)const
 {
-  Info("","chamber=%2i,size=%3i,dim=%5i,x=%6.3f,y=%6.3f,Q=%4i,status=%i",
-            fChamber,fSize,fDimXY,fX,fY,fQdc,fStatus);
+  ::Info("cluster","CombiPid=%10i, c=%2i, size=%4i, dim=%5i, x=%7.3f, y=%7.3f, Q=%6i, st=%i",
+           fCombiPid,fChamber,fSize,fDimXY,fX,fY,fQdc,fStatus);
 }
 //__________________________________________________________________________________________________
 ClassImp(AliRICH)    
@@ -120,41 +120,6 @@ void AliRICH::Hits2SDigits()
 {
 //Create a list of sdigits corresponding to list of hits. Every hit generates one or more sdigits.
   if(GetDebug()) Info("Hit2SDigits","Start.");
-  GetLoader()->LoadHits(); 
-  
-  for(Int_t iEventN=0;iEventN<gAlice->GetEventsPerRun();iEventN++){//events loop
-    gAlice->GetRunLoader()->GetEvent(iEventN);
-    
-    if(!GetLoader()->TreeS()) GetLoader()->MakeTree("S");  
-    MakeBranch("S");
-    
-    ResetHits();  ResetSDigits();  
-    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);        
-        
-        TVector3 globX3(pHit->X(),pHit->Y(),pHit->Z());        
-        TVector3 locX3=C(pHit->C())->Glob2Loc(globX3);
-        
-        Int_t sector;
-        Int_t iTotQdc=Param()->Loc2TotQdc(locX3,pHit->Eloss(),pHit->Pid(),sector);
-        
-        Int_t iPadXmin,iPadXmax,iPadYmin,iPadYmax;
-        AliRICHParam::Loc2Area(locX3,iPadXmin,iPadYmin,iPadXmax,iPadYmax);
-        for(Int_t iPadY=iPadYmin;iPadY<=iPadYmax;iPadY++)
-          for(Int_t iPadX=iPadXmin;iPadX<=iPadXmax;iPadX++){
-            Double_t padQdc=iTotQdc*Param()->Loc2PadFrac(locX3,iPadX,iPadY);
-            if(padQdc>0.1) AddSDigit(pHit->C(),iPadX,iPadY,padQdc,pHit->GetTrack());
-          }            
-      }//specials loop
-    }//prims loop
-    GetLoader()->TreeS()->Fill();
-    GetLoader()->WriteSDigits("OVERWRITE");
-    if(GetDebug()) Info("Hit2SDigits","Event %i processed.",iEventN);
-  }//events loop  
-  GetLoader()->UnloadHits();   GetLoader()->UnloadSDigits();  
-  ResetHits();                 ResetSDigits();
   if(GetDebug()) Info("Hit2SDigits","Stop.");
 }//void AliRICH::Hits2SDigits()
 //__________________________________________________________________________________________________
@@ -162,47 +127,6 @@ void AliRICH::SDigits2Digits()
 {
 //Generate digits from sdigits.
   if(GetDebug()) Info("SDigits2Digits","Start.");
-
-  GetLoader()->LoadSDigits();
-  
-  for(Int_t iEventN=0;iEventN<gAlice->GetEventsPerRun();iEventN++){//events loop
-    gAlice->GetRunLoader()->GetEvent(iEventN);
-    
-    if(!GetLoader()->TreeD()) GetLoader()->MakeTree("D");  MakeBranch("D"); //create TreeD with RICH branches 
-    ResetSDigits();ResetDigits();//reset lists of sdigits and digits
-    GetLoader()->TreeS()->GetEntry(0);  
-    SDigits()->Sort();
-  
-    Int_t chamber,x,y,tr[3],id;
-    chamber=x=y=tr[0]=tr[1]=tr[2]=id=kBad;
-    Double_t q=kBad;
-    Int_t iNdigitsPerPad=kBad;//how many sdigits for a given pad
-        
-    for(Int_t i=0;i<SDigits()->GetEntries();i++){//sdigits loop (sorted)
-      AliRICHdigit *pSdig=(AliRICHdigit*)SDigits()->At(i);
-      if(pSdig->Id()==id){//still the same pad
-        iNdigitsPerPad++;
-        q+=pSdig->Q();
-        if(iNdigitsPerPad<=3)
-          tr[iNdigitsPerPad-1]=pSdig->T(0);
-        else
-          Info("","More then 3 sdigits for the given pad");
-      }else{//new pad, add the pevious one
-        if(id!=kBad) AddDigit(chamber,x,y,(Int_t)q,tr[0],tr[1],tr[2]);//ch-xpad-ypad-qdc-tr1-2-3
-        chamber=pSdig->C();x=pSdig->X();y=pSdig->Y();q=pSdig->Q();tr[0]=pSdig->T(0);id=pSdig->Id();
-        iNdigitsPerPad=1;tr[1]=tr[2]=kBad;
-      }
-    }//sdigits loop (sorted)
-  
-    if(SDigits()->GetEntries())AddDigit(chamber,x,y,(Int_t)q,tr[0],tr[1],tr[2]);//add the last digit
-        
-        
-    GetLoader()->TreeD()->Fill();  
-    GetLoader()->WriteDigits("OVERWRITE");
-    if(GetDebug()) Info("SDigits2Digits","Event %i processed.",iEventN);
-  }//events loop
-  GetLoader()->UnloadSDigits();  GetLoader()->UnloadDigits();  
-  ResetSDigits();                ResetDigits();
   if(GetDebug()) Info("SDigits2Digits","Stop.");
 }//SDigits2Digits()
 //__________________________________________________________________________________________________
@@ -607,7 +531,7 @@ void AliRICH::CreateGeometry()
   Float_t par[3];
     
 //External aluminium box 
-  par[0]=68.8*kcm;par[1]=13*kcm;par[2]=70.86*kcm;  gMC->Gsvolu("RICH", "BOX ", idtmed[1009], par, 3);
+  par[0]=68.8;par[1]=13;par[2]=70.86;  gMC->Gsvolu("RICH", "BOX ", idtmed[1009], par, 3);
 //Air 
   par[0]=66.3;   par[1] = 13; par[2] = 68.35;      gMC->Gsvolu("SRIC", "BOX ", idtmed[1000], par, 3); 
 //Air 2 (cutting the lower part of the box)
@@ -808,21 +732,17 @@ void AliRICH::CreateChambers()
 void AliRICH::GenerateFeedbacks(Int_t iChamber,Float_t eloss)
 {
 // Generate FeedBack photons
-  Int_t j;
-  Float_t cthf, phif, enfp = 0, sthf;
-  Float_t e1[3], e2[3], e3[3];
-  Float_t vmod, uswop;
-  Float_t dir[3], phi;
-  Float_t pol[3], mom[4];
 //Determine number of feedback photons
-  TLorentzVector globX4;
-  gMC->TrackPosition(globX4);  
+  TLorentzVector x4;
+  gMC->TrackPosition(x4);  
   Int_t sector;
-  Int_t iTotQdc=Param()->Loc2TotQdc(C(iChamber)->Glob2Loc(globX4),eloss,gMC->TrackPid(),sector);
-  Int_t iNphotons=gMC->GetRandom()->Poisson(Param()->AlphaFeedback(sector)*iTotQdc);    
-  Info("GenerateFeedbacks","N photons=%i",iNphotons);
+  Int_t iTotQdc=Param()->Loc2TotQdc(C(iChamber)->Glob2Loc(x4),eloss,gMC->TrackPid(),sector);
+  Int_t iNphotons=gMC->GetRandom()->Poisson(P()->AlphaFeedback(sector)*iTotQdc);    
+  if(GetDebug())Info("GenerateFeedbacks","N photons=%i",iNphotons);
+  Int_t j;
+  Float_t cthf, phif, enfp = 0, sthf, e1[3], e2[3], e3[3], vmod, uswop,dir[3], phi,pol[3], mom[4];
 //Generate photons
-  for(Int_t i=0;i<iNphotons;i++){
+  for(Int_t i=0;i<iNphotons;i++){//feddbacks loop
     Double_t ranf[2];
     gMC->GetRandom()->RndmArray(2,ranf);    //Sample direction
     cthf=ranf[0]*2-1.0;
@@ -863,15 +783,8 @@ void AliRICH::GenerateFeedbacks(Int_t iChamber,Float_t eloss)
       e3[j]=uswop;
     }
     
-    vmod=0;
-    for(j=0;j<3;j++) vmod+=e1[j]*e1[j];
-    vmod=TMath::Sqrt(1/vmod);
-    for(j=0;j<3;j++) e1[j]*=vmod;
-    
-    vmod=0;
-    for(j=0;j<3;j++) vmod+=e2[j]*e2[j];
-    vmod=TMath::Sqrt(1/vmod);
-    for(j=0;j<3;j++) e2[j]*=vmod;
+    vmod=0;  for(j=0;j<3;j++) vmod+=e1[j]*e1[j];  vmod=TMath::Sqrt(1/vmod);  for(j=0;j<3;j++) e1[j]*=vmod;    
+    vmod=0;  for(j=0;j<3;j++) vmod+=e2[j]*e2[j];  vmod=TMath::Sqrt(1/vmod);  for(j=0;j<3;j++) e2[j]*=vmod;
     
     phi = gMC->GetRandom()->Rndm()* 2 * TMath::Pi();
     for(j=0;j<3;j++) pol[j]=e1[j]*TMath::Sin(phi)+e2[j]*TMath::Cos(phi);
@@ -881,10 +794,11 @@ void AliRICH::GenerateFeedbacks(Int_t iChamber,Float_t eloss)
                      gAlice->GetMCApp()->GetCurrentTrackNumber(),//parent track 
                      kFeedback,                      //PID
                     mom[0],mom[1],mom[2],mom[3],    //track momentum  
-                     globX4.X(),globX4.Y(),globX4.Z(),globX4.T(),    //track origin 
+                     x4.X(),x4.Y(),x4.Z(),x4.T(),    //track origin 
                      pol[0],pol[1],pol[2],           //polarization
-                    kPFeedBackPhoton,outputNtracksStored,1.0);
-    
-  }
-}//FeedBackPhotons()
+                    kPFeedBackPhoton,
+                     outputNtracksStored,
+                     1.0);    
+  }//feedbacks loop
+}//GenerateFeedbacks()
 //__________________________________________________________________________________________________
index f875078..d92b584 100644 (file)
 #include <AliHit.h>
 #include <AliDigit.h>
 
-#include "AliRICHConst.h"
-class AliRICHChamber;
 #include "AliRICHParam.h"
 
 #include "AliRICHSDigit.h"
 
 //__________________AliRICHhit______________________________________________________________________
-//__________________________________________________________________________________________________
-//__________________________________________________________________________________________________
 class AliRICHhit : public AliHit
 {
 public:
   inline   AliRICHhit();
   inline   AliRICHhit(Int_t fIshunt, Int_t track, Int_t *vol, Float_t *hits);
-  inline   AliRICHhit(Int_t track,Int_t iPID,Int_t iChamber,TLorentzVector x4,Float_t eloss);
+  inline   AliRICHhit(Int_t tid,TVector3 x3,Double_t eloss);
   virtual ~AliRICHhit()         {;}
     
   Int_t   C()                   const{return fChamber;}
   Int_t   Chamber()             const{return fChamber;}
   Int_t   Pid()                 const{return fPid;}    
   Int_t   Particle()            const{return fPid;}
-  Float_t Theta()               const{return fTheta;}
-  Float_t Phi()                 const{return fPhi;}
-  Float_t Tlength()             const{return fTlength;}
   Float_t Eloss()               const{return fEloss;}
-  Float_t Loss()                const{return fLoss;}
-  Float_t PHfirst()             const{return fPHfirst;}
-  Float_t PHlast()              const{return fPHlast;}
   Float_t MomX()                const{return fMomX;}
   Float_t MomY()                const{return fMomY;}
   Float_t MomZ()                const{return fMomZ;}
@@ -50,17 +40,12 @@ public:
 protected:
   Int_t     fChamber;                      //chamber number
   Int_t     fPid;                          //particle code
-  Float_t   fTheta,fPhi ;                  //incident theta phi angles in degrees
-  Float_t   fTlength;                      //track length inside the chamber
-  Float_t   fEloss;                        //ionisation energy loss in gas
-  Float_t   fPHfirst;                      //first padhit
-  Float_t   fPHlast;                       //last padhit
-  Float_t   fLoss;                         // did it hit the freon?
+  Double_t  fEloss;                        //ionisation energy loss in gas
   Float_t   fMomX,fMomY,fMomZ;             //momentum at photochatode entry point
   Float_t   fNPads;                        // Pads hit
   Float_t   fCerenkovAngle;                // Dummy cerenkov angle
   Float_t   fMomFreoX,fMomFreoY,fMomFreoZ; //momentum at freon entry point
-  ClassDef(AliRICHhit,1)                   //RICH hit class
+  ClassDef(AliRICHhit,2)                   //RICH hit class
 };//class AliRICHhit
 
   //__________________________________________________________________________________________________
@@ -68,40 +53,30 @@ AliRICHhit::AliRICHhit()
            :AliHit() 
 {//default ctor  
   fChamber=fPid=kBad;
-  fTheta=fPhi=fTlength=fEloss=fPHfirst=fPHlast=fLoss=kBad;
+  fEloss=kBad;
   fMomX=fMomY=fMomZ=fNPads=fCerenkovAngle=fMomFreoX=fMomFreoY=fMomFreoZ=kBad;
 }
 //__________________________________________________________________________________________________
-AliRICHhit::AliRICHhit(Int_t shunt, Int_t track, Int_t *vol, Float_t *hit):
-            AliHit(shunt, track)
+AliRICHhit::AliRICHhit(Int_t shunt, Int_t track, Int_t *vol, Float_t *hit)
+           :AliHit(shunt, track)
 {//ctor
   fChamber=vol[0];
   fPid=(Int_t)hit[0];
   fX=hit[1];fY=hit[2];fZ=hit[3];
-  fTheta=hit[4];fPhi=hit[5];
-  fTlength=hit[6];
   fEloss=hit[7];
-  fPHfirst=(Int_t)hit[8];
-  fPHlast=(Int_t)hit[9];
-  fLoss=hit[13];
   fMomX=hit[14];fMomY=hit[15];fMomZ=hit[16];
-  fNPads=hit[17];
   fCerenkovAngle=hit[18];
   fMomFreoX=hit[19];fMomFreoY=hit[20];fMomFreoZ=hit[21];
 }
 //__________________________________________________________________________________________________
-AliRICHhit::AliRICHhit(Int_t track,Int_t iPID,Int_t iChamber,TLorentzVector x4,Float_t eloss):
-            AliHit(0, track)
+AliRICHhit::AliRICHhit(Int_t tid,TVector3 x3,Double_t eloss)
+           :AliHit(0,tid)
 {//ctor
-  fChamber=iChamber;
-  fPid=iPID;
-  fX=x4.X();fY=x4.Y();fZ=x4.Z();
+  fX=x3.X();fY=x3.Y();fZ=x3.Z();
   fEloss=eloss;
 }
 
 //__________________AliRICHCerenkov_________________________________________________________________
-//__________________________________________________________________________________________________
-//__________________________________________________________________________________________________
 class AliRICHCerenkov: public AliHit 
 {
 public:
@@ -154,16 +129,18 @@ AliRICHCerenkov::AliRICHCerenkov(Int_t shunt, Int_t track, Int_t *vol, Float_t *
 }
 
 //__________________AliRICHdigit____________________________________________________________________
-//__________________________________________________________________________________________________
-//__________________________________________________________________________________________________
 class AliRICHdigit :public AliDigit
 {
 public:
-           AliRICHdigit() {fPadX=fPadY=fChamber=fTracks[0]=fTracks[1]=fTracks[2]=kBad;fQdc=kBad;}
-  inline   AliRICHdigit(Int_t iC,Int_t iX,Int_t iY,Double_t iQdc,Int_t iT1,Int_t iT2,Int_t iT3);
+           AliRICHdigit() {fCombiPid=fChamber=fPadX=fPadY=fTracks[0]=fTracks[1]=fTracks[2]=kBad;fQdc=kBad;}
+           AliRICHdigit(Int_t c,Int_t x,Int_t y,Double_t q,Int_t cpid,Int_t tid0,Int_t tid1,Int_t tid2)
+           {fPadX=x;fPadY=y;fQdc=q;fChamber=10*c+AliRICHParam::Pad2Sec(x,y);fCombiPid=cpid;fTracks[0]=tid0;fTracks[1]=tid1;fTracks[2]=tid2;}        
   virtual ~AliRICHdigit() {;}  
-  inline   Int_t    Compare(const TObject *pObj) const;                               //virtual
+  Int_t    Compare(const TObject *pObj) const                                         //virtual
+           {if(Id()==((AliRICHdigit*)pObj)->Id()) return 0; else if(Id()>((AliRICHdigit*)pObj)->Id()) return 1;  else return -1;} 
+           
            Bool_t   IsSortable()                                  const{return kTRUE;}//virtual
+           Int_t    CombiPid()                                    const{return fCombiPid;}
            Int_t    C()                                           const{return fChamber/10;}
            Int_t    S()                                           const{return fChamber-(fChamber/10)*10;} 
            Int_t    Chamber()                                     const{return C();}
@@ -172,41 +149,22 @@ public:
            Int_t    Y()                                           const{return fPadY;}
            Int_t    Id()                                          const{return fChamber*10000000+fPadX*1000+fPadY;}
            Double_t Q()                                           const{return fQdc;}
-           Int_t    T(Int_t i)                                    const{return fTracks[i];}
+           Int_t    Tid(Int_t i)                                  const{return fTracks[i];}
            void  Print(Option_t *option="")const;                 //virtual
 protected:
-  Int_t    fChamber;  //10*module number+ sector number 
+  Int_t    fCombiPid; //1000000*Ncerenkovs+1000*Nfeedbacks+Nmips  
+  Int_t    fChamber;  //10*chamber number+ sector number 
   Int_t    fPadX;     //pad number along X
   Int_t    fPadY;     //pad number along Y
-  Double_t fQdc;      //QDC value, fractions are permitted for summable procedure
-  ClassDef(AliRICHdigit,1) //RICH digit class       
+  Double_t fQdc;      //QDC value, fractions are permitted for summable procedure  
+  ClassDef(AliRICHdigit,2) //RICH digit class       
 };//class AliRICHdigit
-//__________________________________________________________________________________________________
-AliRICHdigit::AliRICHdigit(Int_t iC,Int_t iPadX,Int_t iPadY,Double_t q,Int_t iT0,Int_t iT1,Int_t iT2)
-{
-  fPadX=iPadX;fPadY=iPadY;fQdc=q;  
-  fChamber=10*iC+AliRICHParam::Pad2Sec(iPadX,iPadY);
-  fTracks[0]=iT0;fTracks[1]=iT1;fTracks[2]=iT2;
-}
-//__________________________________________________________________________________________________
-Int_t AliRICHdigit::Compare(const TObject *pObj)const
-{
-  if(Id()==((AliRICHdigit*)pObj)->Id()) 
-    return 0;
-  else if(Id()>((AliRICHdigit*)pObj)->Id()) 
-    return 1;
-  else
-    return -1;
-} 
-//__________________AliRICHcluster__________________________________________________________________
-//__________________________________________________________________________________________________
 
-enum ClusterStatus {kOK,kEdge,kShape,kSize,kRaw};
-
-//__________________________________________________________________________________________________
+//__________________AliRICHcluster__________________________________________________________________
 class AliRICHcluster :public TObject
 {
 public:
+  enum ClusterStatus {kOK,kEdge,kShape,kSize,kRaw};
                     AliRICHcluster()                                 {fSize=fQdc=fStatus=fChamber=fDimXY=kBad;fX=fY=kBad;fDigits=0;}
   virtual          ~AliRICHcluster()                                 {delete fDigits;}  
   AliRICHcluster&   operator=(const AliRICHcluster&)                 {return *this;}
@@ -221,12 +179,20 @@ public:
          Double_t   Y()                                         const{return fY;}                          // 
          Int_t      Status()                                    const{return fStatus;}                     //
          void       SetStatus(Int_t status)                          {fStatus=status;}                     //
+         Int_t      Nmips()                                     const{return fCombiPid-1000000*Ncerenkovs()-1000*Nfeedbacks();}  //
+         Int_t      Ncerenkovs()                                const{return fCombiPid/1000000;}           //
+         Int_t      Nfeedbacks()                                const{return (fCombiPid-1000000*Ncerenkovs())/1000;}      //
+         Bool_t     IsPureMip()                                 const{return fCombiPid<1000;}   
+         Bool_t     IsPureCerenkov()                            const{return Nmips()==0&&Nfeedbacks()==0;} //
+         Bool_t     IsPureFeedback()                            const{return Nmips()==0&&Ncerenkovs()==0;} //
+         void       SetCombiPid(Int_t ckov,Int_t feeds,Int_t mips)   {fCombiPid=1000000*ckov+1000*feeds+mips;}            //
          TObjArray* Digits()                                    const{return fDigits;}                     //  
          void       Print(Option_t *option="")const;                                                       //virtual
   inline void       AddDigit(AliRICHdigit *pDig);                                                          //
   inline void       CoG();                                                                                 // 
-         void       Reset() {fSize=fQdc=fStatus=fChamber=fDimXY=kBad;fX=fY=kBad;delete fDigits;fDigits=0;} //                                                                               //
+         void       Reset() {fSize=fQdc=fStatus=fChamber=fDimXY=kBad;fX=fY=kBad;delete fDigits;fDigits=0;} //
 protected:
+  Int_t         fCombiPid;    //1000000*Ncerenkovs+1000*Nfeedbacks+Nmips  
   Int_t         fSize;        //how many digits belong to this cluster    
   Int_t         fDimXY;       //100*xdim+ydim box containing the cluster
   Int_t         fQdc;         //QDC value
@@ -235,12 +201,12 @@ protected:
   Double_t      fY;           //local y postion  
   Int_t         fStatus;      //flag to mark the quality of the cluster   
   TObjArray    *fDigits;      //! list of digits forming this cluster
-  ClassDef(AliRICHcluster,1)  //RICH cluster class       
+  ClassDef(AliRICHcluster,2)  //RICH cluster class       
 };//class AliRICHcluster
 //__________________________________________________________________________________________________
 void AliRICHcluster::AddDigit(AliRICHdigit *pDig)
 {//    
-  if(!fDigits) {fQdc=fSize=0;fDigits = new TObjArray;}
+  if(!fDigits) {fQdc=fSize=fCombiPid=0;fDigits = new TObjArray;}
   fQdc+=(Int_t)pDig->Q(); fDigits->Add(pDig);
   fChamber=10*pDig->C()+pDig->S();
   fSize++;
@@ -262,10 +228,10 @@ void AliRICHcluster::CoG()
    fDimXY = 100*(xmax-xmin+1)+ymax-ymin+1;//find box containing cluster
    fStatus=kRaw;
 }//CoG()
+
 //__________________AliRICH_________________________________________________________________________
-//__________________________________________________________________________________________________
-//__________________________________________________________________________________________________
 class AliRICHParam;
+class AliRICHChamber;
 class AliRICHSDigit;
 
 class AliRICH : public AliDetector 
@@ -285,14 +251,16 @@ public:
   inline  void    CreateSDigits();  
   inline  void    CreateDigits();  
   inline  void    CreateClusters();  
-  inline  void    AddHit(Int_t track, Int_t *vol, Float_t *hits);                                                 //virtual
-  inline  void    AddSDigit(Int_t iC,Int_t iX,Int_t iY,Double_t q,Int_t iT0,Int_t iT1=kBad,Int_t iT2=kBad);
-  inline  void    AddDigit (Int_t iC,Int_t iX,Int_t iY,Int_t iQ,Int_t iT0,Int_t iT1=kBad,Int_t iT2=kBad);     
-  inline  void    AddCluster(AliRICHcluster &clus);
-          void    ResetHits()     {AliDetector::ResetHits();fNcerenkovs=0;if(fCerenkovs)fCerenkovs->Clear();fNspecials=0;if(fSpecials)fSpecials->Clear();}  //virtual
-          void    ResetSDigits()  {fNsdigits=0;  if(fSdigits)  fSdigits ->Clear();}                                 
-          void    ResetDigits()   {if(fDigitsNew)for(int i=0;i<kNCH;i++){fDigitsNew->At(i)->Clear();fNdigitsNew[i]=0;}}
-          void    ResetClusters() {if(fClusters) for(int i=0;i<kNCH;i++){fClusters ->At(i)->Clear();fNclusters[i]=0;}}
+  void AddHit(Int_t track, Int_t *vol, Float_t *hits)     {TClonesArray &tmp=*fHits; new(tmp[fNhits++])AliRICHhit(fIshunt,track,vol,hits);}//virtual
+  void AddHit(Int_t tid,TVector3 x3,Double_t eloss=0)     {TClonesArray &tmp=*fHits;new(tmp[fNhits++])AliRICHhit(tid,x3,eloss);} 
+  inline void AddSDigit(int c,int x,int y,int q,int pid,int tid); 
+  void AddDigit(int c,int x,int y,int q,int cpid,int *tid){TClonesArray &tmp=*((TClonesArray*)fDigitsNew->At(c-1));new(tmp[fNdigitsNew[c-1]++])AliRICHdigit(c,x,y,q,cpid,tid[0],tid[1],tid[2]);}  
+  void AddCluster(AliRICHcluster &cl)                     {TClonesArray &tmp=*((TClonesArray*)fClusters->At(cl.C()-1));new(tmp[fNclusters[cl.C()-1]++])AliRICHcluster(cl);}
+  
+  void ResetHits()     {AliDetector::ResetHits();fNcerenkovs=0;if(fCerenkovs)fCerenkovs->Clear();fNspecials=0;if(fSpecials)fSpecials->Clear();}  //virtual  
+  void ResetSDigits()  {fNsdigits=0;  if(fSdigits)  fSdigits ->Clear();}                                 
+  void ResetDigits()   {if(fDigitsNew)for(int i=0;i<kNCH;i++){fDigitsNew->At(i)->Clear();fNdigitsNew[i]=0;}}
+  void ResetClusters() {if(fClusters) for(int i=0;i<kNCH;i++){fClusters ->At(i)->Clear();fNclusters[i]=0;}}
                   //Hits provided by AliDetector
   TClonesArray*   SDigits()             const{return fSdigits;}
   TClonesArray*   Digits(Int_t iC)      const{if(fDigitsNew) return (TClonesArray *)fDigitsNew->At(iC-1);else return 0;}
@@ -300,6 +268,7 @@ public:
           
   AliRICHChamber* C(Int_t iC)           const{return (AliRICHChamber*)fChambers->At(iC-1);}
   AliRICHParam*   Param()               const{return fpParam;}
+  AliRICHParam*   P()                   const{return fpParam;}
           void    CreateChambers();         
           void    CreateMaterials(); //virtual
   virtual void    BuildGeometry();   //virtual
@@ -349,31 +318,7 @@ protected:
     
   ClassDef(AliRICH,4)                        //Main RICH class 
 };//class AliRICH  
-//__________________________________________________________________________________________________
-void AliRICH::AddHit(Int_t track, Int_t *vol, Float_t *hits)
-{//Adds the current hit to the RICH hits list
-  TClonesArray &tmp=*fHits;
-  new(tmp[fNhits++])AliRICHhit(fIshunt,track,vol,hits);
-}
-//__________________________________________________________________________________________________
-void AliRICH::AddSDigit(Int_t iC,Int_t iX,Int_t iY,Double_t q,Int_t iT0,Int_t iT1,Int_t iT2)
-{//Adds the current Sdigit to the RICH list of Sdigits   
-  TClonesArray &tmp=*fSdigits;
-  new(tmp[fNsdigits++])AliRICHdigit(iC,iX,iY,q,iT0,iT1,iT2);
-} 
-//__________________________________________________________________________________________________
-void AliRICH::AddDigit(Int_t iC,Int_t iX,Int_t iY,Int_t iAdc,Int_t iT0,Int_t iT1,Int_t iT2)
-{//Adds the current digit to the corresponding RICH list of digits (individual list per chamber)
-  TClonesArray &tmp=*((TClonesArray*)fDigitsNew->At(iC-1));
-  new(tmp[fNdigitsNew[iC-1]++]) AliRICHdigit(iC,iX,iY,iAdc,iT0,iT1,iT2);
-}
-//__________________________________________________________________________________________________
-void AliRICH::AddCluster(AliRICHcluster &cluster)
-{//Adds the current cluster to the corresponding RICH list of clusters (individual list per chamber)
-  cluster.Print();
-  TClonesArray &tmp=*((TClonesArray*)fClusters->At(cluster.Chamber()-1));
-  new(tmp[fNclusters[cluster.Chamber()-1]++]) AliRICHcluster(cluster);
-}
+
 //__________________________________________________________________________________________________
 void AliRICH::CreateHits()
 {
@@ -404,10 +349,20 @@ void AliRICH::CreateClusters()
   fClusters = new TObjArray(kNCH);  
   for(Int_t i=0;i<kNCH;i++) {fClusters->AddAt(new TClonesArray("AliRICHcluster",10000), i); fNclusters[i]=0;}
 }
+//__________________________________________________________________________________________________
+void AliRICH::AddSDigit(int c,int x,int y,int q,int pid,int tid) 
+{   
+  switch(pid){
+    case 50000048: pid=1000000;break;
+    case 50000052: pid=1000;   break;
+    default:       pid=1;      break;
+  }   
+  TClonesArray &tmp=*fSdigits;
+  new(tmp[fNsdigits++])AliRICHdigit(c,x,y,q,pid,tid,kBad,kBad);
+}//AddSDigit()   
 
 
-
-//__________________________________________________________________________________________________
+//______OLD OLD OLD OLD_____________________________________________________________________________
 void AliRICH::CreateCerenkovsOld()
 {
   if(fCerenkovs) return;
index 17657ea..e59d1f2 100644 (file)
 #include "AliRICHClusterFinder.h"
 #include "AliRICHMap.h"
 #include <TMinuit.h>
+#include <TParticle.h>
 #include <TVector3.h>
 #include <AliLoader.h>
+#include <AliStack.h>
 #include <AliRun.h>
 
 void RICHMinMathieson(Int_t &npar, Double_t *gin, Double_t &chi2, Double_t *par, Int_t iflag);
@@ -73,6 +75,9 @@ void AliRICHClusterFinder::Exec()
   
   Rich()->GetLoader()->LoadDigits(); 
   
+  Rich()->GetLoader()->GetRunLoader()->LoadHeader();
+  Rich()->GetLoader()->GetRunLoader()->LoadKinematics();
+
   for(Int_t iEventN=0;iEventN<gAlice->GetEventsPerRun();iEventN++){//events loop
     gAlice->GetRunLoader()->GetEvent(iEventN);
     
@@ -81,22 +86,25 @@ void AliRICHClusterFinder::Exec()
     
     Rich()->GetLoader()->TreeD()->GetEntry(0);
     for(Int_t iChamber=1;iChamber<=kNCH;iChamber++){//chambers loop
-     FindRawClusters(iChamber);
-        
+      FindClusters(iChamber);
     }//chambers loop
-    
     Rich()->GetLoader()->TreeR()->Fill();
     Rich()->GetLoader()->WriteRecPoints("OVERWRITE");
   }//events loop  
   Rich()->GetLoader()->UnloadDigits(); Rich()->GetLoader()->UnloadRecPoints();  
   Rich()->ResetDigits();  Rich()->ResetClusters();
+  
+  Rich()->GetLoader()->GetRunLoader()->UnloadHeader();
+  Rich()->GetLoader()->GetRunLoader()->UnloadKinematics();
+
   Info("Exec","Stop.");      
 }//Exec()
 //__________________________________________________________________________________________________
-void AliRICHClusterFinder::FindRawClusters(Int_t iChamber)
-{//finds neighbours and fill the tree with raw clusters
+void AliRICHClusterFinder::FindClusters(Int_t iChamber)
+{
+  //finds neighbours and fill the tree with raw clusters
   Int_t nDigits=Rich()->Digits(iChamber)->GetEntriesFast();
-  Info("FindRawClusters","Start for Chamber %i with %i digits.",iChamber,nDigits);  
+  Info("FindClusters","Start for Chamber %i with %i digits.",iChamber,nDigits);  
   if(nDigits==0)return;
 
   fHitMap=new AliRICHMap(Rich()->Digits(iChamber));//create digit map for the given chamber
@@ -112,17 +120,52 @@ void AliRICHClusterFinder::FindRawClusters(Int_t iChamber)
       ResolveCluster(); // ResolveCluster serialization will happen inside
     } else {
       WriteRawCluster(); // simply output of the RawCluster found without deconvolution
-    }    
+    }
     fCurrentCluster.Reset();
   }//digits loop
 
   delete fHitMap;
-  Info("FindRawClusters","Stop.");
+  Info("FindClusters","Stop.");
   
-}//FindRawClusters()
+}//FindClusters()
+//__________________________________________________________________________________________________
+void AliRICHClusterFinder::FindClusterContribs()
+{
+  //finds CombiPid for a given cluster
+  Info("FindClusterContribs","Start");
+  
+  TObjArray *pDigits = fCurrentCluster.Digits();
+  Int_t iNmips=0,iNckovs=0,iNfeeds=0;
+  TArrayI contribs(3*fCurrentCluster.Size());
+  Int_t *pindex = new Int_t[3*fCurrentCluster.Size()];
+  for(Int_t iDigN=0;iDigN<fCurrentCluster.Size();iDigN++) {//loop on digits of a given cluster
+    contribs[3*iDigN]  =((AliRICHdigit*)pDigits->At(iDigN))->Tid(0);
+    contribs[3*iDigN+1]=((AliRICHdigit*)pDigits->At(iDigN))->Tid(1);
+    contribs[3*iDigN+2]=((AliRICHdigit*)pDigits->At(iDigN))->Tid(2);
+    cout << "TID1 " << contribs[3*iDigN] << " TID2 " << contribs[3*iDigN+1] << " TID3 " << contribs[3*iDigN+2] << endl;
+  }//loop on digits of a given cluster
+  TMath::Sort(contribs.GetSize(),contribs.GetArray(),pindex);
+  for(Int_t iDigN=0;iDigN<3*fCurrentCluster.Size()-1;iDigN++) {//loop on digits to sort Tid
+    cout << " contrib " << contribs[pindex[iDigN]] << " digit " << iDigN << 
+            " contrib " << contribs[pindex[iDigN+1]] << " digit " << iDigN+1 << endl;
+    if(contribs[pindex[iDigN]]!=contribs[pindex[iDigN+1]]) {
+      Int_t code   = Rich()->GetLoader()->GetRunLoader()->Stack()->Particle(contribs[pindex[iDigN]])->GetPdgCode();
+      Double_t charge = Rich()->GetLoader()->GetRunLoader()->Stack()->Particle(contribs[pindex[iDigN]])->GetPDG()->Charge();
+
+      if(code==50000050) iNckovs++;
+      else if(code==50000051) iNfeeds++;
+      else if(charge!=0) iNmips++;
+      if (contribs[pindex[iDigN+1]]==kBad) break;
+    }
+  }//loop on digits to sort Tid
+  fCurrentCluster.SetCombiPid(iNckovs,iNfeeds,iNmips);
+  fCurrentCluster.Print();
+  delete [] pindex; 
+}// FindClusterContribs()
 //__________________________________________________________________________________________________
 void  AliRICHClusterFinder::FormRawCluster(Int_t i, Int_t j)
-{// Builder of the final Raw Cluster (before deconvolution)  
+{
+  // Builder of the final Raw Cluster (before deconvolution)  
   Info("FormRawCluster","Start with digit(%i,%i)",i,j);
   
   fCurrentCluster.AddDigit((AliRICHdigit*) fHitMap->GetHit(i,j));
@@ -157,6 +200,7 @@ void AliRICHClusterFinder::WriteRawCluster()
 {// out the current RawCluster
   Info("WriteRawCluster","Start.");
   
+  FindClusterContribs();
   Rich()->AddCluster(fCurrentCluster);
   
 }//WriteRawCluster()
index 08d83ea..d1c0206 100644 (file)
@@ -17,7 +17,8 @@ public:
   
   AliRICH *Rich() {return fRICH;}                                             //Pointer to RICH  
   void     Exec();                                                            //Loop on events and chambers  
-  void     FindRawClusters(Int_t iChamber);                                   //Find raw clusters  
+  void     FindClusters(Int_t iChamber);                                      //Find raw clusters
+  void     FindClusterContribs();                                             //Find CombiPid for the current cluster
   void     FormRawCluster(Int_t i, Int_t j);                                  //form a raw cluster
   void     FindLocalMaxima();                                                 //Find local maxima in a cluster
   void     ResolveCluster();                                                  //Try to resolve a cluster with maxima > 2
index cc1c8d2..312e064 100644 (file)
@@ -54,7 +54,6 @@
 #include "AliPDG.h"
 #include "AliDetector.h"
 #include "AliRICH.h"
-#include "AliRICHConst.h"
 #include "AliRICHDisplay.h"
 #include "AliRICHPoints.h"
 #include "AliHeader.h"
@@ -82,56 +81,6 @@ AliRICHDisplay::AliRICHDisplay()
 //_____________________________________________________________________________
 AliRICHDisplay::AliRICHDisplay(Int_t size)
 {
-// Create an event display object.
-// A canvas named "edisplay" is created with a vertical size in pixels
-//
-//    A QUICK Overview of the Event Display functions
-//    ===============================================
-//
-//  The event display can ve invoked by executing the macro "display.C"
-// A canvas like in the picture below will appear.
-//
-//  On the left side of the canvas, the following buttons appear:
-//   *Next*       to move to the next event
-//   *Previous*   to move to the previous event
-
-//   *Pick*       Select this option to be able to point on a track with the
-//                mouse. Once on the track, use the right button to select
-//                an action. For example, select SetMarkerAttributes to
-//                change the marker type/color/size for the track.
-//   *Zoom*       Select this option (default) if you want to zoom.
-//                To zoom, simply select the selected area with the left button.
-//   *UnZoom*     To revert to the previous picture size.
-//
-//   slider R     On the left side, the vertical slider can be used to
-//                set the default picture size.
-//
-//    When you are in Zoom mode, you can click on the black part of the canvas
-//  to select special options with the right mouse button.
-
-//
-//  When you are in pick mode, you can "Inspect" the object pointed by the mouse.
-//  When you are on a track, select the menu item "InspectParticle"
-//  to display the current particle attributes.
-//
-//  You can activate the Root browser by selecting the Inspect menu
-//  in the canvas tool bar menu. Then select "Start Browser"
-//  This will open a new canvas with the browser. At this point, you may want
-//  to display some histograms (from the Trees). Go to the "File" menu
-//  of the browser and click on "New canvas".
-//  In the browser, click on item "ROOT files" in the left pane.
-//  Click on galice.root.
-//  Click on TH
-//  Click on TPC for example
-//  Click on any variable (eg TPC.fX) to histogram the variable.
-//
-//   If you are lost, you can click on HELP in any Root canvas or browser.
-//Begin_Html
-/*
-  <img src="gif/AliRICHDisplay.gif">
-*/
-//End_Html
-    
     
     fPad = 0;
     
index 6afae8a..3c63a82 100644 (file)
@@ -14,6 +14,7 @@
 //  **************************************************************************
 #include "AliRICHParam.h"
 
+ClassImp(AliRICHParam)
 Bool_t   AliRICHParam::fgIsWireSag            =kTRUE;
 Bool_t   AliRICHParam::fgIsResolveClusters    =kTRUE;
 Double_t AliRICHParam::fgAngleRot             =-60;
@@ -21,7 +22,7 @@ Int_t    AliRICHParam::fgHV                   =2150;
 Int_t    AliRICHParam::fgNsigmaTh             =4;
 Float_t  AliRICHParam::fgSigmaThMean          =1.5;
 Float_t  AliRICHParam::fgSigmaThSpread        =0.5;      
-ClassImp(AliRICHParam)
+Float_t  AliRICHParam::fSigmaThMap[kNCH][kNpadsX][kNpadsY];
 
 void AliRICHParam::GenSigmaThMap()
 {
index 99eb1db..fdee982 100644 (file)
@@ -1,21 +1,32 @@
 #ifndef AliRICHParam_h
 #define AliRICHParam_h
 
-#include "AliRICHConst.h"
 #include <TObject.h>
 #include <TMath.h>
 #include <TVector3.h>
 #include <TRandom.h>
 
+
+static const int kNCH=7;           //number of RICH chambers 
+static const int kNpadsX = 144;    //number of pads along X in single chamber
+static const int kNpadsY = 160;    //number of pads along Y in single chamber
+static const int kBad=-101;        //useful static const to mark initial (uninitalised) values
+
+
+static const int kadc_satm  = 4096;  //dynamic range (10 bits)
+static const int kCerenkov=50000050;  //??? go to something more general like TPDGCode
+static const int kFeedback=50000051;  //??? go to something more general like TPDGCode
+
+
 class AliRICHParam :public TObject  
 {
 public:
            AliRICHParam()                    {;}
   virtual ~AliRICHParam()                    {;}
-  static const Int_t   NpadsX()             {return kNpadsX;}
-  static const Int_t   NpadsY()             {return kNpadsY;}   
-  static Int_t   NpadsXsec()                {return NpadsX()/3;}   
-  static Int_t   NpadsYsec()                {return NpadsY()/2;}   
+  static const Int_t   NpadsX()              {return kNpadsX;}
+  static const Int_t   NpadsY()              {return kNpadsY;}   
+  static Int_t   NpadsXsec()                 {return NpadsX()/3;}   
+  static Int_t   NpadsYsec()                 {return NpadsY()/2;}   
   static Double_t DeadZone()                 {return 2.6;}
   static Double_t PadSizeX()                 {return 0.84;}
   static Double_t PadSizeY()                 {return 0.8;}
@@ -69,6 +80,7 @@ public:
   
   inline static Int_t   Loc2Sec(Double_t &x,Double_t &y); 
   inline static Int_t   Pad2Sec(Int_t &padx,Int_t &pady); 
+  static Int_t   Sector(Int_t padx,Int_t pady) {return Pad2Sec(padx,pady);}
   inline Bool_t IsOverTh(Int_t iChamber, Int_t x, Int_t y, Double_t q);
   static Int_t NsigmaTh() {return fgNsigmaTh;}
   static Float_t SigmaThMean() {return fgSigmaThMean;}
@@ -79,7 +91,7 @@ protected:
   static Bool_t  fgIsResolveClusters;                   //performs declustering or not
   static Int_t   fgHV;                                  //HV applied to anod wires
   static Double_t fgAngleRot;                           //rotation of RICH from up postion (0,0,490)cm
-  Float_t fSigmaThMap[kNCH][kNpadsX][kNpadsY];          // sigma of the pedestal distributions for all pads
+  static Float_t fSigmaThMap[kNCH][kNpadsX][kNpadsY];   // sigma of the pedestal distributions for all pads
   static Int_t fgNsigmaTh;                              // n. of sigmas to cut for zero suppression
   static Float_t fgSigmaThMean;                         // sigma threshold value
   static Float_t fgSigmaThSpread;                       // spread of sigma
@@ -89,10 +101,11 @@ protected:
 Int_t AliRICHParam::PadNeighbours(Int_t iPadX,Int_t iPadY,Int_t listX[4],Int_t listY[4])
 {
   Int_t nPads=0;
-  if(iPadY<NpadsY()){listX[nPads]=iPadX;   listY[nPads]=iPadY+1; nPads++;}       
-  if(iPadX<NpadsX()){listX[nPads]=iPadX+1; listY[nPads]=iPadY;   nPads++;}       
-  if(iPadY>1)       {listX[nPads]=iPadX;   listY[nPads]=iPadY-1; nPads++;}      
-  if(iPadX>1)       {listX[nPads]=iPadX-1; listY[nPads]=iPadY;   nPads++;}       
+  if(iPadY!=NpadsY()&&iPadY!=NpadsYsec())                      {listX[nPads]=iPadX;   listY[nPads]=iPadY+1; nPads++;}       
+  if(iPadX!=NpadsXsec()&&iPadX!=2*NpadsXsec()&&iPadX!=NpadsX()){listX[nPads]=iPadX+1; listY[nPads]=iPadY;   nPads++;}       
+  if(iPadY!=1&&iPadY!=NpadsYsec()+1)                           {listX[nPads]=iPadX;   listY[nPads]=iPadY-1; nPads++;}      
+  if(iPadX!=1&&iPadX!=NpadsXsec()+1&&iPadX!=2*NpadsXsec()+1)   {listX[nPads]=iPadX-1; listY[nPads]=iPadY;   nPads++;}
+
   return nPads;
 }//Pad2ClosePads()
 //__________________________________________________________________________________________________
@@ -165,12 +178,15 @@ void AliRICHParam::Pad2Loc(Int_t padx,Int_t pady,Double_t &x,Double_t &y)
 //__________________________________________________________________________________________________
 Double_t AliRICHParam::GainVariation(Double_t y,Int_t sector)
 {
+//returns % of gain degradation due to wire sagita  
   if(IsWireSag()){
     if(y>0) y-=SectorSizeY()/2; else  y+=SectorSizeY()/2; 
     switch(HV(sector)){
-      case 2150:
-      default:  
-        return 9e-6*TMath::Power(y,4)+2e-7*TMath::Power(y,3)-0.0316*TMath::Power(y,2)-3e-4*y+25.367;//%
+      case 2150: return 9e-6*TMath::Power(y,4)+2e-7*TMath::Power(y,3)-0.0316*TMath::Power(y,2)-3e-4*y+25.367;//%
+      case 2100: return 8e-6*TMath::Power(y,4)+2e-7*TMath::Power(y,3)-0.0283*TMath::Power(y,2)-2e-4*y+23.015;
+      case 2050: return 7e-6*TMath::Power(y,4)+1e-7*TMath::Power(y,3)-0.0254*TMath::Power(y,2)-2e-4*y+20.888;
+      case 2000: return 6e-6*TMath::Power(y,4)+8e-8*TMath::Power(y,3)-0.0227*TMath::Power(y,2)-1e-4*y+18.961;
+      default:   return 0;
     }
   }else
     return 0;
index 53e4983..b164f22 100644 (file)
@@ -210,7 +210,7 @@ void AliRICHPoints::ShowRing(Int_t highlight) {
      //printf("Particle %d belongs to ring %d \n", fTrackIndex, mdig->fTracks[1]);
 
     if (!points) continue;
-    if (fTrackIndex == mdig->T(0)) {
+    if (fTrackIndex == mdig->Tid(0)) {
 
 
       Int_t charge=(Int_t)mdig->Q();
index 0b115cd..49d569f 100644 (file)
@@ -14,7 +14,6 @@
 //**************************************************************************
 
 #include "AliRICHv0.h"
-#include "AliRICHConst.h"
 #include "AliRICHChamber.h" 
 #include <AliRun.h>
 #include <AliMC.h>
@@ -46,14 +45,15 @@ void AliRICHv0::StepManager()
       sParticle="not known";break;
   }
 
-  Info("","event=%i hunt=%i tid=%i pid=%i(%s) m=%f q=%3.1f",
+  Info("","event=%i hunt=%i tid=%i pid=%i(%s) m=%f q=%3.1f dEdX=%9.3f",
                             gMC->CurrentEvent(),
                             fIshunt,
                             gAlice->GetMCApp()->GetCurrentTrackNumber(),
                             gMC->TrackPid(),
                             sParticle,
                             gMC->TrackMass(),
-                            gMC->TrackCharge());
+                            gMC->TrackCharge(),
+                            gMC->Edep());
   Info("","Flags:alive(%i) disap(%i) enter(%i) exit(%i) inside(%i) out(%i) stop(%i) new(%i)",
                             gMC->IsTrackAlive(),
                             gMC->IsTrackDisappeared(),
index e7ce4e4..9f38988 100644 (file)
@@ -33,182 +33,49 @@ void AliRICHv1::StepManager()
 {
 //Full Step Manager
 
-  Int_t          copy, id;
+  Int_t          copy;
   static Int_t   iCurrentChamber;
-  static Int_t   vol[2];
-  Int_t          ipart;
-  static Float_t hits[22];
-  static Float_t ckovData[19];
   TLorentzVector x4,p4;
   Float_t        pos[3],mom[4],localPos[3],localMom[4];
-  Float_t        theta,phi,localTheta,localPhi;
-  Float_t        destep, step;
-  Double_t        ranf[2];
-  Int_t          nPads=kBad;
   Float_t        coscerenkov;
-  static Float_t eloss,  tlength;
-  const  Float_t kBig=1.e10;
        
   TParticle *current = (TParticle*)(*gAlice->GetMCApp()->Particles())[gAlice->GetMCApp()->GetCurrentTrackNumber()];
 
     
  
-  id=gMC->CurrentVolID(copy); iCurrentChamber=copy;
   Float_t cherenkovLoss=0;
     
-  gMC->TrackPosition(x4);    pos[0]=x4(0);    pos[1]=x4(1);    pos[2]=x4(2);
-  ckovData[1] = pos[0];   ckovData[2] = pos[1];   ckovData[3] = pos[2]; 
-  ckovData[6] = 0;                      // dummy track length
     
-    /********************Store production parameters for Cerenkov photons************************/ 
 
   if(gMC->TrackPid()==kCerenkov){//C
     Float_t ckovEnergy = current->Energy();
     if(ckovEnergy > 5.6e-09 && ckovEnergy < 7.8e-09 ){//C+E
-      if(gMC->IsTrackEntering()){//C+E+enter
-        if(gMC->VolId("FRE1")==gMC->CurrentVolID(copy)||gMC->VolId("FRE2")==gMC->CurrentVolID(copy)){//C+E+enter+FRE
-          if(gMC->IsNewTrack()){//C+E+enter+FRE+new
-            Int_t mother=current->GetFirstMother(); 
-            ckovData[10]=mother;
-            ckovData[11]=gAlice->GetMCApp()->GetCurrentTrackNumber();
-            ckovData[12]=1;             //Media where photon was produced 1->Freon, 2->Quarz
-            fCkovNumber++;
-            fFreonProd=1;
-          }//C+E+enter+FRE+new
-        }//C+E+enter+FRE
-        if(gMC->IsNewTrack()&&gMC->VolId("QUAR")==gMC->CurrentVolID(copy))  ckovData[12]=2;            
-      }//C+E+enter           
-      if(ckovData[12]==1){//C+E+produced in Freon
         if(gMC->IsTrackEntering()){                                     //is track entering?
-                   if (gMC->VolId("META")==gMC->CurrentVolID(copy)){                //is it in gap?      
-                  gMC->TrackMomentum(p4);   mom[0]=p4(0);   mom[1]=p4(1);   mom[2]=p4(2);   mom[3]=p4(3);
-                       
-                       gMC->Gmtod(mom,localMom,2);
-                       Float_t cophi = TMath::Cos(TMath::ATan2(localMom[0], localMom[1]));
-                       Float_t t = (1. - .025 / cophi) * (1. - .05 /  cophi);
-                       /**************** Photons lost in second grid have to be calculated by hand************/ 
-                       gMC->GetRandom()->RndmArray(1,ranf);
-                       if (ranf[0] > t) {
-                         gMC->StopTrack();
-                         ckovData[13] = 5;
-                         //AddCerenkov(gAlice->GetCurrentTrackNumber(),vol,ckovData);
-                       }
-                       /**********************************************************************************/
-                     }    //gap
                    
                    if (gMC->VolId("CSI ")==gMC->CurrentVolID(copy)){             //is it in csi?      
-                     
                        gMC->TrackMomentum(p4); mom[0]=p4(0);   mom[1]=p4(1);   mom[2]=p4(2);   mom[3]=p4(3);
-
                        gMC->Gmtod(mom,localMom,2);
-                       /********* Photons lost by Fresnel reflection have to be calculated by hand********/ 
-                       /***********************Cerenkov phtons (always polarised)*************************/
-                       Double_t localTc = localMom[0]*localMom[0]+localMom[2]*localMom[2];
-                       Double_t localRt = TMath::Sqrt(localTc);
-                       localTheta   = Float_t(TMath::ATan2(localRt,Double_t(localMom[1])));
-                       Double_t cotheta = TMath::Abs(cos(localTheta));
-                       Float_t t = Fresnel(ckovEnergy*1e9,cotheta,1);
-                           gMC->GetRandom()->RndmArray(1,ranf);
-                           if (ranf[0] < t) {
-                             gMC->StopTrack();
-                             ckovData[13] = 6;
-                             //AddCerenkov(gAlice->GetCurrentTrackNumber(),vol,ckovData);
+                       Double_t localTc    = localMom[0]*localMom[0]+localMom[2]*localMom[2];
+                       Double_t localTheta = TMath::ATan2(TMath::Sqrt(localTc),localMom[1]);
+                       Double_t cotheta = TMath::Abs(TMath::Cos(localTheta));
+                        if(gMC->GetRandom()->Rndm() < Fresnel(p4.E()*1e9,cotheta,1))  gMC->StopTrack();
                                
-                             //printf("Added One (2)!\n");
-                             //printf("Lost by Fresnel\n");
-                           }
-                           /**********************************************************************************/
-                     }
+                     }//C+E+produced in Freon
                  } //track entering?
-                 /********************Evaluation of losses************************/
-                 /******************still in the old fashion**********************/
-                 TArrayI procs;
-                 Int_t i1 = gMC->StepProcesses(procs);            //number of physics mechanisms acting on the particle
-                 for (Int_t i = 0; i < i1; ++i) {
-                   //        Reflection loss 
-                   if (procs[i] == kPLightReflection) {        //was it reflected
-                     ckovData[13]=10;
-                     if (gMC->VolId("FRE1")==gMC->CurrentVolID(copy) || gMC->VolId("FRE2")==gMC->CurrentVolID(copy)) 
-                       ckovData[13]=1;
-                     if (gMC->CurrentVolID(copy) == gMC->VolId("QUAR")) 
-                       ckovData[13]=2;
-                   } //reflection question
-                    
-                   //        Absorption loss 
-                   else if (procs[i] == kPLightAbsorption) {              //was it absorbed?
-                     //printf("Got in absorption\n");
-                     ckovData[13]=20;
-                     if (gMC->VolId("FRE1")==gMC->CurrentVolID(copy) || gMC->VolId("FRE2")==gMC->CurrentVolID(copy)) 
-                       ckovData[13]=11;
-                     if (gMC->CurrentVolID(copy) == gMC->VolId("QUAR")) 
-                       ckovData[13]=12;
-                     if (gMC->CurrentVolID(copy) == gMC->VolId("META")) 
-                       ckovData[13]=13;
-                     if (gMC->CurrentVolID(copy) == gMC->VolId("GAP ")) 
-                       ckovData[13]=13;
-                     
-                     if (gMC->CurrentVolID(copy) == gMC->VolId("SRIC")) 
-                       ckovData[13]=15;
-                     
-                     //        CsI inefficiency 
-                     if (gMC->CurrentVolID(copy) == gMC->VolId("CSI ")) {
-                       ckovData[13]=16;
-                     }
-                     gMC->StopTrack();
-                     //AddCerenkov(gAlice->GetCurrentTrackNumber(),vol,ckovData);
-                   } //absorption question 
-                   
-                   
-                   //        Photon goes out of tracking scope 
-                   else if (procs[i] == kPStop) {                 //is it below energy treshold?
-                     ckovData[13]=21;
-                     gMC->StopTrack();
-                     //AddCerenkov(gAlice->GetCurrentTrackNumber(),vol,ckovData);
-                   }   // energy treshold question         
-                 }  //number of mechanisms cycle
-                 /**********************End of evaluation************************/
-      }//C+E+produced in Freon
     }//C+E
   }//C
-//*************************************End of Production Parameters Storing*********************
     
-  if(gMC->TrackPid()==kCerenkov||gMC->TrackPid()==kFeedback){//CF      
-    if(gMC->VolId("CSI ")==gMC->CurrentVolID(copy)){//CF+CSI
+  if((gMC->TrackPid()==kCerenkov||gMC->TrackPid()==kFeedback)&&gMC->CurrentVolID(copy)==gMC->VolId("CSI ")){//photon in CSI     
       if(gMC->Edep()>0.){//CF+CSI+DE
         gMC->TrackPosition(x4);   pos[0]=x4(0);   pos[1]=x4(1);   pos[2]=x4(2);
         gMC->TrackMomentum(p4);   mom[0]=p4(0);   mom[1]=p4(1);   mom[2]=p4(2);   mom[3]=p4(3);
-        
-       Double_t tc = mom[0]*mom[0]+mom[1]*mom[1];
-       Double_t rt = TMath::Sqrt(tc);
-       theta   = Float_t(TMath::ATan2(rt,Double_t(mom[2])))*kRaddeg;
-       phi     = Float_t(TMath::ATan2(Double_t(mom[1]),Double_t(mom[0])))*kRaddeg;
-               
-       gMC->CurrentVolOffID(2,copy);   vol[0]=copy;    iCurrentChamber=vol[0];
+        if(IsFresnelLoss()){ gMC->StopTrack(); return;}        
+       gMC->CurrentVolOffID(2,copy);iCurrentChamber=copy;
        
         gMC->Gmtod(pos,localPos,1);     gMC->Gmtod(mom,localMom,2);
 
-       //Param()->SigGenInit(localPos[0],localPos[2]);
-        ckovData[0]=gMC->TrackPid();   
-        ckovData[1]=pos[0];      ckovData[2]=pos[1];    ckovData[3]=pos[2];            
-        ckovData[4]=theta;       ckovData[5]=phi;       //theta-phi angles of incidence 
-        ckovData[8]=(Float_t) fNsdigits;      // first sdigit
-        ckovData[9]=-1;                       // last pad hit
-        ckovData[13]=4;                       // photon was detected
-        ckovData[14]=mom[0];     ckovData[15]=mom[1];    ckovData[16]=mom[2];
-                   
-       destep = gMC->Edep();
-       gMC->SetMaxStep(kBig);
-       cherenkovLoss  += destep;
-       ckovData[7]=cherenkovLoss;
+       cherenkovLoss  += gMC->Edep();
                    
-       GenerateFeedbacks(iCurrentChamber,cherenkovLoss);//CF+CSI+DE
-                                   
-       if (fNsdigits > (Int_t)ckovData[8]) {
-         ckovData[8]= ckovData[8]+1;
-          ckovData[9]= (Float_t) fNsdigits;
-        }
-
-        ckovData[17] = nPads;
         AliRICHhit *mipHit =  (AliRICHhit*) (fHits->UncheckedAt(0));
         if(mipHit){
           mom[0] = current->Px();   mom[1] = current->Py();   mom[2] = current->Pz();
@@ -222,74 +89,29 @@ void AliRICHv1::StepManager()
             coscerenkov = (mom[0]*mipPx + mom[1]*mipPy + mom[2]*mipPz)/(rt*mipRt);
           else
             coscerenkov = 0;
-          ckovData[18]=TMath::ACos(coscerenkov);//Cerenkov angle
         }
-        AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber(),vol,ckovData);//PHOTON HIT CF+CSI+DE
-        //AddCerenkov(gAlice->GetCurrentTrackNumber(),vol,ckovData);
+        AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber(),x4.Vect());//HIT for PHOTON in conditions CF+CSI+DE
+       GenerateFeedbacks(iCurrentChamber,cherenkovLoss);//CF+CSI+DE
       }//CF+CSI+DE
-    }//CF+CSI
-  }/*CF*/else if(gMC->TrackCharge()){//MIP
-    if(gMC->VolId("FRE1")==gMC->CurrentVolID(copy)||gMC->VolId("FRE2")==gMC->CurrentVolID(copy)){//MIP+FRE
-      gMC->TrackMomentum(p4);      mom[0]=p4(0);   mom[1]=p4(1);   mom[2]=p4(2);   mom[3]=p4(3);
-      hits[19]=mom[0];     hits [20] = mom[1];     hits [21] = mom[2];     fFreonProd=1;
-    }//MIP+FRE
-    if(gMC->VolId("GAP ")==gMC->CurrentVolID(copy)){//MIP+GAP
-      gMC->CurrentVolOffID(3,copy);    vol[0]=copy;    iCurrentChamber=vol[0];
-      gMC->TrackPosition(x4);   pos[0]=x4(0);   pos[1]=x4(1);   pos[2]=x4(2);
-      gMC->TrackMomentum(p4);   mom[0]=p4(0);   mom[1]=p4(1);   mom[2]=p4(2);    mom[3]=p4(3);
-      gMC->Gmtod(pos,localPos,1);          gMC->Gmtod(mom,localMom,2);
-      ipart =gMC->TrackPid();
-      destep = gMC->Edep();step   = gMC->TrackStep();// momentum loss and steplength in last step
-      if(gMC->IsTrackEntering()){//MIP+GAP+Enter  record hit when mip enters ...
-        Double_t tc = mom[0]*mom[0]+mom[1]*mom[1];
-        Double_t rt = TMath::Sqrt(tc);
-        theta   = Float_t(TMath::ATan2(rt,Double_t(mom[2])))*kRaddeg;
-        phi     = Float_t(TMath::ATan2(Double_t(mom[1]),Double_t(mom[0])))*kRaddeg;
-       Double_t localTc = localMom[0]*localMom[0]+localMom[2]*localMom[2];
-       Double_t localRt = TMath::Sqrt(localTc);
-       localTheta   = Float_t(TMath::ATan2(localRt,Double_t(localMom[1])))*kRaddeg;                       
-       localPhi     = Float_t(TMath::ATan2(Double_t(localMom[2]),Double_t(localMom[0])))*kRaddeg;    
-        hits[0] = Float_t(ipart);         // particle type
-        hits[1] = localPos[0];    hits[2] = localPos[1];     hits[3] = localPos[2];                 
-        hits[4] = localTheta;     hits[5] = localPhi;              // theta-phi angles of incidence 
-        hits[8] = (Float_t) fNsdigits;    // first sdigit
-        hits[9] = -1;                     // last pad hit
-        hits[13] = fFreonProd;           // did id hit the freon?
-        hits[14] = mom[0];        hits[15] = mom[1];        hits[16] = mom[2];
-        hits[18] = 0;               // dummy cerenkov angle
-        tlength = 0;        eloss   = 0;        fFreonProd = 0;
-        C(iCurrentChamber)->LocaltoGlobal(localPos,hits+1);
-          
-        //Param()->SigGenInit(localPos[0], localPos[2]);
-      }/*MIP+GAP+Enter*/else if(gMC->IsTrackExiting()||gMC->IsTrackStop()||gMC->IsTrackDisappeared()){//MIP+GAP+Exit
-        gMC->SetMaxStep(kBig);
-       eloss   += destep;
-       tlength += step;
-        if (eloss > 0) {
-          if(gMC->TrackPid() == kNeutron) printf("\n\n\n\n\n Neutron Making Pad Hit!!! \n\n\n\n");
-          GenerateFeedbacks(iCurrentChamber,eloss);//MIP+GAP+Exit
-          hits[17] = nPads;
-        }
-        hits[6]=tlength;        hits[7]=eloss;
-        if(fNsdigits > (Int_t)hits[8]) {
-          hits[8]= hits[8]+1;
-          hits[9]= (Float_t) fNsdigits;
-        }
-       AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber(),vol,hits);//MIP HIT MIP+GAP+Exit
-       eloss = 0; 
-      }/*MIP+GAP+Exit*/else if(1/*Param()->SigGenCond(localPos[0], localPos[2])*/){//MIP+GAP+Spec
-        //Param()->SigGenInit(localPos[0], localPos[2]);
-        if(eloss>0){
-          if(gMC->TrackPid() == kNeutron) printf("\n\n\n\n\n Neutron Making Pad Hit!!! \n\n\n\n");
-          GenerateFeedbacks(iCurrentChamber,eloss);//MIP+GAP+Spec
-          hits[17] = nPads;
-        }
-        
-       eloss    = destep;      tlength += step ;               
-      }/*MIP+GAP+Spec*/else{//MIP+GAP+nothing special
-        eloss   += destep;
-        tlength += step ;
-      }//MIP+GAP+nothing special
-    }//MIP+GAP
-  }//MIP
+  }//CF in CSI
+  
+//Treat charged particles  
+  static Float_t eloss;
+  if(gMC->TrackCharge() && gMC->CurrentVolID(copy)==gMC->VolId("GAP ")){//MIP in GAP
+    gMC->CurrentVolOffID(3,copy); iCurrentChamber=copy;
+    if(gMC->IsTrackEntering()||gMC->IsNewTrack())//MIP in GAP Entering
+      eloss=0;                                                           
+    else if(gMC->IsTrackExiting()||gMC->IsTrackStop()||gMC->IsTrackDisappeared()){//MIP in GAP Exiting
+      eloss+=gMC->Edep();//take into account last step dEdX
+      gMC->TrackPosition(x4);  
+      AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber(),x4.Vect(),eloss);//HIT for MIP for conditions: MIP in GAP Exiting
+      GenerateFeedbacks(iCurrentChamber,eloss);//MIP+GAP+Exit
+    }else//MIP in GAP going inside
+      eloss   += gMC->Edep();
+  }//MIP in GAP
 }//void AliRICHv1::StepManager()
+
+Bool_t AliRICHv1::IsFresnelLoss()
+{
+  return kFALSE;
+}
index 7deb76a..97cc265 100644 (file)
@@ -14,7 +14,9 @@ public:
   virtual       ~AliRICHv1()                                                         {;}
   virtual void   Init()                                                              {;}
   virtual Int_t  IsVersion()                                                    const{return 1;}
+  
   virtual void   StepManager();
+          Bool_t IsFresnelLoss(); 
 private:
   ClassDef(AliRICHv1,1)//RICH full version for simulation
 };
index 044db81..4678ca1 100644 (file)
@@ -55,7 +55,7 @@ void Config()
   }   
   switch(kGen){
        case kPP7:          Pythia7(pid=kPiPlus,p=4);                break;
-        case kGun7:         Gun7(pid=kPiPlus,p=4);                   break;
+        case kGun7:         Gun7(pid=kProton);                       break;
         case kGun1:         Gun1(n=1,pid=kPiPlus,p=4,chamber=4);     break;
         case kGun0:         Gun1(n=1,pid=kNeutron,p=4,chamber=4);    break;
         default:            Fatal("Config","No generator");          break;
@@ -94,14 +94,14 @@ void Para()
 }
 
 
-void Gun7(Int_t iPID,Double_t p)
+void Gun7(Int_t iPID)
 {
-  ::Info("kir-Gun7","7 primaries of %i PID  with p=%f GeV",iPID,p);
+  ::Info("kir-Gun7","7 primaries of %i PID  with 1.5 < p < 4.5 GeV/c",iPID);
   AliGenCocktail *pCocktail=new AliGenCocktail();
   for(int i=1;i<=7;i++){
     AliGenFixed *pFixed=new AliGenFixed(1);
-    pFixed->SetMomentum(p);
-    pFixed->SetPhiRange(gRICH->C(i)->PhiD());
+    pFixed->SetMomentum(1.5+(i-1)*0.5);
+    pFixed->SetPhiRange(gRICH->C(i)->PhiD()+3);
     pFixed->SetThetaRange(gRICH->C(i)->ThetaD()+2);
     pFixed->SetOrigin(0,0,0);                 
     pFixed->SetPart(iPID);    
index 3115fa6..db32bcd 100644 (file)
@@ -1,66 +1,84 @@
+
+Int_t countContrib[7][3];
+
 void ControlPlots()
 {  
   Int_t iChamber=1;
   
-  TH1F *pCqH1=new TH1F("clusq","Cluster Charge;q [QDC]",AliRICHParam::MaxQdc(),0,AliRICHParam::MaxQdc());
-  TH1F *pDqH1=new TH1F("digq","Digit Charge;q [QDC]",AliRICHParam::MaxQdc(),0,AliRICHParam::MaxQdc());
-  al->LoadHeader();  al->LoadKinematics();
-  
-  Bool_t isHits=!rl->LoadHits();  
-    Bool_t isSDigits=!rl->LoadSDigits();  
-      Bool_t isClusters=!rl->LoadRecPoints();
-        Bool_t isDigits=!rl->LoadDigits();
-  
+  TFile *pFile = new TFile("$(HOME)/plots.root","RECREATE");   
+  TH1F *pCqH1=new TH1F("ClusQ",   "Cluster Charge all chambers;q [QDC]",r->P()->MaxQdc(),0,r->P()->MaxQdc());
+  TH1F *pCsH1=new TH1F("ClusSize","Cluster size all chambers;size [number of pads in cluster]",100,0,100);
+  TH2F *pCmH2=new TH2F("ClusMap", "Cluster map;x [cm];y [cm]",1000,-r->P()->PcSizeX()/2,r->P()->PcSizeX()/2,
+                                                             1000,-r->P()->PcSizeY()/2,r->P()->PcSizeY()/2);
+  
+  TH1F *pCqMipH1=new TH1F("MipClusQ",   "MIP Cluster Charge all chambers;q [QDC]",r->P()->MaxQdc(),0,r->P()->MaxQdc());
+  TH1F *pCsMipH1=new TH1F("MipClusSize","MIP Cluster size all chambers;size [number of pads in cluster]",100,0,100);
+  TH2F *pCmMipH2=new TH2F("MipClusMap", "MIP Cluster map;x [cm];y [cm]",1000,-r->P()->PcSizeX()/2,r->P()->PcSizeX()/2,
+                                                             1000,-r->P()->PcSizeY()/2,r->P()->PcSizeY()/2);
+  
+  TH1F *pCqCerH1=new TH1F("CerClusQ",   "Cerenkov Cluster Charge all chambers;q [QDC]",r->P()->MaxQdc(),0,r->P()->MaxQdc());
+  TH1F *pCsCerH1=new TH1F("CerClusSize","Cernekov Cluster size all chambers;size [number of pads in cluster]",100,0,100);
+  TH2F *pCmCerH2=new TH2F("CerClusMap", "Cerenkov Cluster map;x [cm];y [cm]",1000,-r->P()->PcSizeX()/2,r->P()->PcSizeX()/2,
+                                                             1000,-r->P()->PcSizeY()/2,r->P()->PcSizeY()/2);
+  TH1F *pCqFeeH1=new TH1F("FeeClusQ",   "Feedback Cluster Charge all chambers;q [QDC]",r->P()->MaxQdc(),0,r->P()->MaxQdc());
+  TH1F *pCsFeeH1=new TH1F("FeeClusSize","Feedback Cluster size all chambers;size [number of pads in cluster]",100,0,100);
+  TH2F *pCmFeeH2=new TH2F("FeeClusMap", "Feedback Cluster map;x [cm];y [cm]",1000,-r->P()->PcSizeX()/2,r->P()->PcSizeX()/2,
+                                                             1000,-r->P()->PcSizeY()/2,r->P()->PcSizeY()/2);
+  Bool_t isClusters=!rl->LoadRecPoints();
+  r->SetTreeAddress();  
   for(Int_t iEventN=0;iEventN<a->GetEventsPerRun();iEventN++){//events loop
     al->GetEvent(iEventN);    
-    if(isHits){
-      for(Int_t iPrimN=0;iPrimN<rl->TreeH()->GetEntries();iPrimN++){//prims loop
-        rl->TreeH()->GetEntry(iPrimN);      
-        Int_t iTotalHits=0,iTotalCerenkovs=0,iTotalSpecials=0;
-        iTotalHits+=r->Hits()->GetEntries();
-        iTotalCerenkovs+=r->Cerenkovs()->GetEntries();
-        iTotalSpecials+=r->Specials()->GetEntries();
-      }//prims loop
-    }//isHits
-    if(isSDigits){
-      rl->TreeS()->GetEntry(0);
-    }//isSdigits
-    if(isDigits){
-      rl->TreeD()->GetEntry(0);
-      Int_t iTotalDigits=0;
-      for(int i=1;i<=7;i++) iTotalDigits+=r->Clusters(i)->GetEntries();    
-      for(Int_t iDigitN=0;iDigitN<r->Digits(iChamber)->GetEntries();iDigitN++){
-        pDqH1->Fill(((AliRICHdigit*)r->Digits(iChamber)->At(iDigitN))->Q());
-      }//digits loop
-    }//isDigits
     if(isClusters){
       rl->TreeR()->GetEntry(0);
       Int_t iTotalClusters=0;
-      for(int i=1;i<=7;i++) iTotalClusters+=r->Clusters(i)->GetEntries();    
-      for(Int_t iClusterN=0;iClusterN<r->Clusters(iChamber)->GetEntries();iClusterN++){
-        pCqH1->Fill(((AliRICHcluster*)r->Clusters(iChamber)->At(iClusterN))->Q());
-      }//clusters loop
+      for(int i=1;i<=7;i++){//chambers loop
+        iTotalClusters+=r->Clusters(i)->GetEntries();    
+        for(Int_t iClusterN=0;iClusterN<r->Clusters(i)->GetEntries();iClusterN++){//clusters loop
+          AliRICHcluster *pClus=(AliRICHcluster*)r->Clusters(i)->At(iClusterN);
+          
+          countContrib[i-1][0] += pClus->Nmips();
+          countContrib[i-1][1] += pClus->Ncerenkovs();
+          countContrib[i-1][2] += pClus->Nfeedbacks();
+              
+          pCqH1->Fill(pClus->Q());             
+          pCsH1->Fill(pClus->Size());           
+          pCmH2->Fill(pClus->X(),pClus->Y());  
+          
+          if(pClus->IsPureMip()){ //Pure Mips
+            pCqMipH1->Fill(pClus->Q());
+            pCsMipH1->Fill(pClus->Size()); 
+            pCmMipH2->Fill(pClus->X(),pClus->Y());
+          }
+          if(pClus->IsPureCerenkov()){ //Pure Photons
+            pCqCerH1->Fill(pClus->Q());
+            pCsCerH1->Fill(pClus->Size()); 
+            pCmCerH2->Fill(pClus->X(),pClus->Y());
+          }
+          if(pClus->IsPureFeedback()){ //Pure Feedbacks
+            pCqFeeH1->Fill(pClus->Q());
+            pCsFeeH1->Fill(pClus->Size()); 
+            pCmFeeH2->Fill(pClus->X(),pClus->Y());
+          }
+        }//clusters loop
+      }//chambers loop
     }//isClusters
-    cout<<"Event "<<iEventN<<endl;
+    Info("ControlPlots","Event %i processed.",iEventN);
   }//events loop 
-  if(isHits) rl->UnloadHits();  
-    if(isSDigits) rl->UnloadSDigits(); 
-      if(isDigits) rl->UnloadDigits(); 
-        if(isClusters) rl->UnloadRecPoints();
-          al->UnloadHeader();  al->UnloadKinematics();
-  TCanvas *pC=new TCanvas("c","Control Plots",600,500);
-  pC->Divide(2,2);
-  pC->cd(1); pCqH1->Draw();
-  pC->cd(2); pDqH1->Draw();  
+  if(isClusters) rl->UnloadRecPoints();
+  
+  pFile->Write();
+  delete pFile;
+  for(Int_t i=0;i<7;i++)
+    cout <<" chamber " << i+1 << " n. mips " << countContrib[i][0] << " n. ckovs " << countContrib[i][1] << " n. fdbks " << countContrib[i][2] << endl;
 }//void ControlPlots()
 //__________________________________________________________________________________________________
 void MainTrank()
 {
   TStopwatch sw;TDatime time;
-  OLD_S_SD(); SD_D(); D_C();
+  OLD_S_SD(); SD_D();   AliRICHClusterFinder *z=new AliRICHClusterFinder(r); z->Exec();//delete z;  
   cout<<"\nInfo in <MainTrank>: Start time: ";time.Print();
-  cout<<"Info in <MainTrank>: Stop  time: ";time.Set();  time.Print();
-  cout<<"Info in <MainTrank>: Time  used: ";sw.Print();
+    cout<<"Info in <MainTrank>: Stop  time: ";time.Set();  time.Print();
+    cout<<"Info in <MainTrank>: Time  used: ";sw.Print();
 }
 //__________________________________________________________________________________________________
 void sh()
@@ -125,11 +143,7 @@ void sc()
 Double_t r2d = TMath::RadToDeg();
 Double_t d2r = TMath::DegToRad();
 
-void DisplFast()
-{
-  AliRICHDisplFast *d = new AliRICHDisplFast();
-  d->Exec();
-}  
+void DisplFast(){ AliRICHDisplFast *d = new AliRICHDisplFast();  d->Exec();}  
 
 
 void Digits2Recos()
@@ -148,20 +162,13 @@ void Digits2Recos()
 
 
 
-void D_C()
-{
-  AliRICHClusterFinder *z=new AliRICHClusterFinder(r);
-  z->Exec();  
-}
 //__________________________________________________________________________________________________
 void OLD_S_SD()
 {
-  Info("OLD_S_SD","Start.");
-  
-  rl->LoadHits(); 
-  
+  Info("OLD_S_SD","Start.");  
+  rl->LoadHits();   
   for(Int_t iEventN=0;iEventN<a->GetEventsPerRun();iEventN++){//events loop
-    al->GetEvent(iEventN);    cout<<"Event "<<iEventN<<endl;  
+    al->GetEvent(iEventN);    Info("OLD_S_SD","Processing event %i",iEventN);  
     
     rl->MakeTree("S");  r->MakeBranch("S");
     r->ResetSDigits();  r->ResetSpecialsOld();
@@ -172,13 +179,15 @@ void OLD_S_SD()
         Int_t padx=1+ ((AliRICHSDigit*)r->Specials()->At(i))->PadX()+r->Param()->NpadsX()/2;
         Int_t pady=1+ ((AliRICHSDigit*)r->Specials()->At(i))->PadY()+r->Param()->NpadsY()/2;
         Double_t q=  ((AliRICHSDigit*)r->Specials()->At(i))->QPad();
+        
         Int_t hitN= ((AliRICHSDigit*)r->Specials()->At(i))->HitNumber()-1;//!!! important -1
         Int_t chamber=((AliRICHhit*)r->Hits()->At(hitN))->C();
-        Int_t track=((AliRICHhit*)r->Hits()->At(hitN))->GetTrack();
+        Int_t tid=((AliRICHhit*)r->Hits()->At(hitN))->GetTrack();
+        Int_t pid=((AliRICHhit*)r->Hits()->At(hitN))->Pid();
         if(padx<1 || padx>r->Param()->NpadsX() ||pady<1 || pady>r->Param()->NpadsY())
           Warning("OLD_S_SD","pad is out of valid range padx= %i pady=%i event %i",padx,pady,iEventN);
         else
-          r->AddSDigit(chamber,padx,pady,q,track);
+          r->AddSDigit(chamber,padx,pady,q,pid,tid);
       }//specials loop
     }//prims loop
     rl->TreeS()->Fill();
@@ -230,7 +239,7 @@ void H_SD()
 void SD_D()
 {
   Info("SD_D","Start.");  
-
+  extern Int_t kBad; 
   r->Param()->GenSigmaThMap();
   rl->LoadSDigits();
   
@@ -240,32 +249,33 @@ void SD_D()
     r->ResetSDigits();r->ResetDigits();//reset lists of sdigits and digits
     rl->TreeS()->GetEntry(0);  
     r->SDigits()->Sort();
-  
-    Int_t kBad=-101;
-    Int_t chamber,x,y,tr[3],id;
-    Double_t q=kBad;
-    chamber=x=y=tr[0]=tr[1]=tr[2]=id=kBad;
-    Int_t iNdigitsPerPad=kBad;//how many sdigits for a given pad
-        
+      
+    Int_t combiPid,chamber,x,y,tid[3],id; Double_t q;
+    Int_t iNdigitsPerPad;//how many sdigits for a given pad
+    const int kBad=-101;//??? to be removed in code    
     for(Int_t i=0;i<r->SDigits()->GetEntries();i++){//sdigits loop (sorted)
       AliRICHdigit *pSdig=(AliRICHdigit*)r->SDigits()->At(i);
       if(pSdig->Id()==id){//still the same pad
         iNdigitsPerPad++;
         q+=pSdig->Q();
+        combiPid+=pSdig->CombiPid();
         if(iNdigitsPerPad<=3)
-          tr[iNdigitsPerPad-1]=pSdig->T(0);
+          tid[iNdigitsPerPad-1]=pSdig->Tid(0);
 //        else
 //          Info("","More then 3 sdigits for the given pad");
       }else{//new pad, add the pevious one
         if(id!=kBad&&r->Param()->IsOverTh(chamber,x,y,q))
-           r->AddDigit(chamber,x,y,q,tr[0],tr[1],tr[2]);//ch-xpad-ypad-qdc-tr1-2-3
-        chamber=pSdig->C();x=pSdig->X();y=pSdig->Y();q=pSdig->Q();tr[0]=pSdig->T(0);id=pSdig->Id();
-        iNdigitsPerPad=1;tr[1]=tr[2]=kBad;
+           r->AddDigit(chamber,x,y,q,combiPid,tid);
+        combiPid=pSdig->CombiPid();chamber=pSdig->C();id=pSdig->Id();
+        x=pSdig->X();y=pSdig->Y();
+        q=pSdig->Q();
+        tid[0]=pSdig->Tid(0);
+        iNdigitsPerPad=1;tid[1]=tid[2]=kBad;
       }
     }//sdigits loop (sorted)
   
     if(r->SDigits()->GetEntries()&&r->Param()->IsOverTh(chamber,x,y,q))
-      r->AddDigit(chamber,x,y,q,tr[0],tr[1],tr[2]);//add the last digit
+      r->AddDigit(chamber,x,y,q,combiPid,tid);//add the last digit
         
     rl->TreeD()->Fill();  
     rl->WriteDigits("OVERWRITE");