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