//__________________________________________________________________________________________________
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)
{
//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()
//__________________________________________________________________________________________________
{
//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()
//__________________________________________________________________________________________________
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)
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;
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);
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()
//__________________________________________________________________________________________________
#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;}
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
//__________________________________________________________________________________________________
: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:
}
//__________________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();}
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;}
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
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++;
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
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;}
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
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()
{
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;
#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);
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);
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
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));
{// out the current RawCluster
Info("WriteRawCluster","Start.");
+ FindClusterContribs();
Rich()->AddCluster(fCurrentCluster);
}//WriteRawCluster()
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
#include "AliPDG.h"
#include "AliDetector.h"
#include "AliRICH.h"
-#include "AliRICHConst.h"
#include "AliRICHDisplay.h"
#include "AliRICHPoints.h"
#include "AliHeader.h"
//_____________________________________________________________________________
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;
// **************************************************************************
#include "AliRICHParam.h"
+ClassImp(AliRICHParam)
Bool_t AliRICHParam::fgIsWireSag =kTRUE;
Bool_t AliRICHParam::fgIsResolveClusters =kTRUE;
Double_t AliRICHParam::fgAngleRot =-60;
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()
{
#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;}
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;}
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
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()
//__________________________________________________________________________________________________
//__________________________________________________________________________________________________
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;
//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();
//**************************************************************************
#include "AliRICHv0.h"
-#include "AliRICHConst.h"
#include "AliRICHChamber.h"
#include <AliRun.h>
#include <AliMC.h>
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(),
{
//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();
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;
+}
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
};
}
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;
}
-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);
+
+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()
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()
-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();
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();
void SD_D()
{
Info("SD_D","Start.");
-
+ extern Int_t kBad;
r->Param()->GenSigmaThMap();
rl->LoadSDigits();
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");