--- /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 <TH1F.h>
#include <TH2F.h>
#include <TStopwatch.h>
+#include <AliLog.h>
ClassImp(AliRICHhit)
//__________________________________________________________________________________________________
void AliRICHhit::Print(Option_t*)const
{
- ::Info("hit","Ch=%1i, TID=%6i, eloss=%9.3f eV, in-out dist=%9.4f, OUT(%7.2f,%7.2f,%7.2f)"
- ,fChamber,fTrack,fEloss*1e9,Length(),fOutX3.X(),fOutX3.Y(),fOutX3.Z());
+ AliInfo(Form("Ch=%1i, TID=%6i, eloss=%9.3f eV, in-out dist=%9.4f, OUT(%7.2f,%7.2f,%7.2f)"
+ ,fChamber,fTrack,fEloss*1e9,Length(),fOutX3.X(),fOutX3.Y(),fOutX3.Z()));
}
//__________________________________________________________________________________________________
ClassImp(AliRICHdigit)
//__________________________________________________________________________________________________
void AliRICHdigit::Print(Option_t*)const
{
- ::Info("digit","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]);
+ 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]));
}
//__________________________________________________________________________________________________
ClassImp(AliRICHcluster)
::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
- ::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,0);
+ 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));
}
//__________________________________________________________________________________________________
:AliDetector(name,title),fpParam(new AliRICHParam),fSdigits(0),fNsdigits(0),fDigitsNew(0),fClusters(0)
{
//Named ctor
- if(GetDebug())Info("named ctor","Start.");
+ AliDebug(1,"Start.");
//AliDetector ctor deals with Hits and Digits (reset them to 0, does not create them)
CreateHits(); gAlice->GetMCApp()->AddHitList(fHits);
fCounters.ResizeTo(20); fCounters.Zero();
- if(GetDebug())Info("named ctor","Stop.");
+ AliDebug(1,"Stop.");
}//AliRICH::AliRICH(const char *name, const char *title)
//__________________________________________________________________________________________________
AliRICH::~AliRICH()
{
//dtor
- if(GetDebug()) Info("dtor","Start.");
+ AliDebug(1,"Start.");
if(fpParam) delete fpParam;
if(fDigits) delete fDigits;
if(fDigitsNew) {fDigitsNew->Delete(); delete fDigitsNew;}
if(fClusters) {fClusters->Delete(); delete fClusters;}
- if(GetDebug()) Info("dtor","Stop.");
+ 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.
- if(GetDebug()) Info("Hit2SDigits","Start.");
+ AliDebug(1,"Start.");
for(Int_t iEventN=0;iEventN<GetLoader()->GetRunLoader()->GetAliRun()->GetEventsPerRun();iEventN++){//events loop
GetLoader()->GetRunLoader()->GetEvent(iEventN);//get next event
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())->Glob2Loc(pHit->OutX3());//hit position in the chamber local system
+ TVector2 x2 = C(pHit->C())->Mrs2Pc(pHit->OutX3());//hit position in photocathode plane
Int_t iTotQdc=P()->TotQdc(x2,pHit->Eloss());//total charge produced by hit, 0 if hit in dead zone
if(iTotQdc==0) continue;
TVector area=P()->Loc2Area(x2);//determine affected pads, dead zones analysed inside
- if(GetDebug()) Info("Hits2SDigits","hit(%6.2f,%6.2f)->area(%3.0f,%3.0f)-(%3.0f,%3.0f) QDC=%4i",x2.X(),x2.Y(),area[0],area[1],area[2],area[3],iTotQdc);
+ AliDebug(1,Form("hit(%6.2f,%6.2f)->area(%3.0f,%3.0f)-(%3.0f,%3.0f) QDC=%4i",x2.X(),x2.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(x2,pad);
- if(GetDebug()) Info("Hits2SDigits","current pad(%3.0f,%3.0f) with QDC =%6.2f",pad[0],pad[1],padQdc);
+ AliDebug(1,Form("current pad(%3.0f,%3.0f) with QDC =%6.2f",pad[0],pad[1],padQdc));
if(padQdc>0.1) AddSDigit(pHit->C(),pad,padQdc,GetLoader()->GetRunLoader()->Stack()->Particle(pHit->GetTrack())->GetPdgCode(),pHit->GetTrack());
}//affected pads loop
}//hits loop
}//events loop
GetLoader()->UnloadHits(); GetLoader()->GetRunLoader()->UnloadHeader(); GetLoader()->GetRunLoader()->UnloadKinematics();
GetLoader()->UnloadSDigits();
- if(GetDebug()) Info("Hit2SDigits","Stop.");
+ AliDebug(1,"Stop.");
}//Hits2SDigits()
//__________________________________________________________________________________________________
void AliRICH::BuildGeometry()
{
//Builds a TNode geometry for event display
- if(GetDebug())Info("BuildGeometry","Start.");
+ AliInfo("Start.");
TNode *node, *subnode, *top;
top=gAlice->GetGeometry()->GetNode("alice");
-
- new TBRIK("S_RICH","S_RICH","void",71.09999,11.5,73.15);
- Float_t wid=P()->SectorSizeX();
- Float_t len=P()->SectorSizeY();
- new TBRIK("PHOTO","PHOTO","void",wid/2,0.1,len/2);
-
+ Float_t widx=P()->SectorSizeX();
+ Float_t leny=P()->SectorSizeY();
+ Float_t dz = P()->Zfreon()+P()->Zwin()+P()->Pc2Win();
+ Float_t dead = P()->DeadZone();
+
+ new TBRIK("RICH","RICH","void",widx+dead/2,leny+leny/2+dead,dz+0.1); //RICH chamber
+ new TBRIK("RPC" ,"RPC" ,"void",widx/2,leny/2,0.01); //RICH sector
+
for(int i=1;i<=kNchambers;i++){
top->cd();
- node = new TNode(Form("RICH%i",i),Form("RICH%i",i),"S_RICH",C(i)->X(),C(i)->Y(),C(i)->Z(),C(i)->RotMatrixName());
+ node = new TNode(Form("RICH%i",i),Form("RICH%i",i),"RICH",C(i)->Center().X(),C(i)->Center().Y(),C(i)->Center().Z(),C(i)->RotMatrixName());
node->SetLineColor(kRed);
node->cd();
- subnode = new TNode("PHOTO1","PHOTO1","PHOTO",wid+P()->DeadZone(),5,len/2+P()->DeadZone()/2,"");
+ subnode = new TNode("PHOTO1","PHOTO1","RPC",-widx/2-dead/2,-leny-dead/2,dz,"");
subnode->SetLineColor(kGreen);
fNodes->Add(subnode);
- subnode = new TNode("PHOTO1","PHOTO1","PHOTO",0,5,len/2+P()->DeadZone()/2,"");
+ subnode = new TNode("PHOTO1","PHOTO1","RPC", widx/2+dead/2,-leny-dead/2,dz,"");
subnode->SetLineColor(kGreen);
fNodes->Add(subnode);
- subnode = new TNode("PHOTO1","PHOTO1","PHOTO",-wid-P()->DeadZone(),5,len/2+P()->DeadZone()/2,"");
+ subnode = new TNode("PHOTO1","PHOTO1","RPC",-widx/2-dead/2, 0,dz,"");
subnode->SetLineColor(kGreen);
fNodes->Add(subnode);
- subnode = new TNode("PHOTO1","PHOTO1","PHOTO",wid+P()->DeadZone(),5,-len/2-P()->DeadZone()/2,"");
+ subnode = new TNode("PHOTO1","PHOTO1","RPC", widx/2+dead/2, 0,dz,"");
subnode->SetLineColor(kGreen);
fNodes->Add(subnode);
- subnode = new TNode("PHOTO1","PHOTO1","PHOTO",0,5,-len/2 -P()->DeadZone()/2,"");
+ subnode = new TNode("PHOTO1","PHOTO1","RPC",-widx/2-dead/2, leny+dead/2,dz,"");
subnode->SetLineColor(kGreen);
fNodes->Add(subnode);
- subnode = new TNode("PHOTO1","PHOTO1","PHOTO",-wid-P()->DeadZone(),5,-len/2 - P()->DeadZone()/2,"");
+ subnode = new TNode("PHOTO1","PHOTO1","RPC", widx/2+dead/2, leny+dead/2,dz,"");
subnode->SetLineColor(kGreen);
fNodes->Add(subnode);
fNodes->Add(node);
- }
- if(GetDebug())Info("BuildGeometry","Stop.");
+ }
+
+ AliDebug(1,"Stop.");
}//void AliRICH::BuildGeometry()
//______________________________________________________________________________
void AliRICH::CreateMaterials()
{
// Definition of available RICH materials
-#include "Opticals.h"
- Float_t a=0,z=0,den=0,radl=0,absl=0;
+ Int_t material=0; //tmp material id number
+ Float_t a=0,z=0,den=0,radl=0,absl=0; //tmp material parameters
Float_t tmaxfd=-10.0, deemax=-0.2, stemax=-0.1,epsil=0.001, stmin=-0.001;
Int_t isxfld = gAlice->Field()->Integ();
Float_t sxmgmx = gAlice->Field()->Max();
- Int_t material;
- AliMaterial(material=1, "Air $",a=14.61,z=7.3, den=0.001205,radl=30420.0,absl=67500);//(Air)
- AliMedium(1, "DEFAULT MEDIUM AIR$",material, 0, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
+ Float_t aAir[4]={12.,14.,16.,36.}; Float_t zAir[4]={6.,7.,8.,18.}; Float_t wAir[4]={0.000124,0.755267,0.231781,0.012827};//total 0.9999999
+ AliMixture(++material, "RichAir",aAir,zAir,den=0.00120479,4,wAir); //1 (Air) 0.01% C 75% N 23% O 1% Ar
+ AliMedium(kAir, "RichAir",material, 0, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
- AliMaterial( 6, "HON", a=12.01,z=6.0, den=0.1, radl=18.8, absl=0); //(C)-equivalent radl
- AliMedium(2, "HONEYCOMB$", 6, 0, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
+ AliMaterial(++material, "RichRohacell", a=12.01,z=6.0, den=0.1, radl=18.8, absl=0); //2 Rohacell 51 C-equiv radl rad cover
+ AliMedium(kRoha, "RichRohacell", material, 0, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
- AliMaterial(16, "CSI", a=12.01,z=6.0, den=0.1, radl=18.8, absl=0); //CsI-radl equivalent
- AliMedium(kCSI, "CSI$", 16, 1, isxfld, sxmgmx,tmaxfd, stemax, deemax, epsil, stmin);
+ Float_t aQuartz[2]={28.09,16.0}; Float_t zQuartz[2]={14.00, 8.0}; Float_t wQuartz[2]={1,2};
+ AliMixture(++material, "RichSiO2",aQuartz,zQuartz,den=2.64,-2, wQuartz); //3 Quarz (SiO2) -trasparent rad window
+ AliMedium(kSiO2, "RichSiO2",material, 1, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
- AliMaterial(11, "GRI", a=63.54,z=29.0,den=8.96, radl=1.43, absl=0); //anode grid (Cu)
- AliMedium(7, "GRID$", 11, 0, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
+ Float_t aFreon[2]={12,19}; Float_t zFreon[2]={6,9}; Float_t wmatFreon[2]={6,14};
+ AliMixture(++material, "RichC6F14",aFreon,zFreon,den=1.68,-2,wmatFreon); //4 Freon (C6F14)
+ AliMedium(kC6F14, "RichC6F14",material, 1, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
+
+ Float_t aMethane[2]={12.01,1}; Float_t zMethane[2]={6,1}; Float_t wMethane[2]={1,4};
+ AliMixture (++material, "RichCH4", aMethane, zMethane, den=7.17e-4,-2, wMethane); //5,9 methane (CH4) normal and for Gap
+ AliMedium(kCH4, "RichCH4" , material, 1, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
+ AliMixture (++material, "RichCH4", aMethane, zMethane, den=7.17e-4,-2, wMethane); //5,9 methane (CH4) normal and for Gap
+ AliMedium(kGap, "RichCH4gap$", material, 1, isxfld, sxmgmx, tmaxfd, 0.1 , -deemax, epsil, -stmin);
+
+ AliMaterial(++material, "RichCsI", a=12.01,z=6.0, den=0.1, radl=18.8, absl=0); //6 CsI-radl equivalent
+ AliMedium(kCsI, "RichCsI$", material, 1, isxfld, sxmgmx,tmaxfd, stemax, deemax, epsil, stmin);
- AliMaterial(50, "ALUM", a=26.98,z=13.0,den=2.699, radl=8.9, absl=0); //aluminium sheet (Al)
- AliMedium(10, "ALUMINUM$", 50, 1, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
+ AliMaterial(++material, "GridCu", a=63.54,z=29.0,den=8.96, radl=1.43, absl=0); //7 anode grid (Cu)
+ AliMedium(kGridCu, "GRID$", material, 0, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
- AliMaterial(material=31, "COPPER$", a=63.54,z=29.0,den=8.96, radl=1.4, absl=0); //(Cu)
- AliMedium(12, "PCB_COPPER", material, 0, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
+ AliMixture (++material, "OpSiO2",aQuartz, zQuartz, den=2.64, -2, wQuartz); //8 Quarz (SiO2) - opaque
+ AliMedium(kOpSiO2, "QUARTZO$",material, 1, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
- Float_t aQuartz[2]={28.09,16.0}; Float_t zQuartz[2]={14.00, 8.0}; Float_t wmatQuartz[2]={1,2};
- AliMixture (20, "QUA",aQuartz,zQuartz,den=2.64,-2, wmatQuartz);//Quarz (SiO2) - trasnparent
- AliMedium(3, "QUARTZ$", 20, 1, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
+ AliMaterial(++material, "ALUM", a=26.98,z=13.0,den=2.699, radl=8.9, absl=0); //10 aluminium sheet (Al)
+ AliMedium(kAl, "ALUMINUM$", material, 1, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
- AliMixture (21, "QUAO",aQuartz, zQuartz, den=2.64, -2, wmatQuartz);//Quarz (SiO2) - opaque
- AliMedium(8, "QUARTZO$", 21, 1, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
+ Float_t aGlass[5]={12.01,28.09,16,10.8,23}; Float_t zGlass[5]={6,14,8,5,11}; Float_t wGlass[5]={0.5,0.105,0.355,0.03,0.01};
+ AliMixture(++material,"GLASS",aGlass, zGlass, den=1.74, 5, wGlass); //11 Glass 50%-C 10.5%-Si 35.5%-O 3%-B 1%-Na
+ AliMedium(kGlass, "GLASS", material, 0, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
- Float_t aFreon[2]={12,19}; Float_t zFreon[2]={6,9}; Float_t wmatFreon[2]={6,14};
- AliMixture (material=30, "C6F14",aFreon,zFreon,den=1.68,-2,wmatFreon);//Freon (C6F14)
- AliMedium(4, "C6F14$",material, 1, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
+ AliMaterial(++material, "COPPER$", a=63.54,z=29.0,den=8.96, radl=1.4, absl=0); //12 Cu
+ AliMedium(kCu, "PCB_COPPER", material, 0, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
+
+
+ AliMaterial(++material, "W$", a=183.84,z=74.0,den=19.3, radl=0.35, absl=185.0/den); //13 W - anod wires
+ AliMedium(kW, "W", material, 0, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
- Float_t aMethane[2]={12.01,1}; Float_t zMethane[2]={6,1}; Float_t wmatMethane[2]={1,4};
- AliMixture (material=40, "CH4", aMethane, zMethane, den=7.17e-4,-2, wmatMethane);//methane (CH4)
- AliMedium(kCH4, "CH4$" , material, 1, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
- AliMedium(kGAP, "CH4GAP$", material, 1, isxfld, sxmgmx,tmaxfd, 0.1, -deemax, epsil, -stmin);
if(P()->IsRadioSrc()){
- AliMaterial(material=45, "STEEL$", a=63.54,z=29.0,den=8.96, radl=1.4, absl=0); //Steel
+ AliInfo("Special radioactive source materials");
+ AliMaterial(++material, "STEEL$", a=63.54,z=29.0,den=8.96, radl=1.4, absl=0); //14 Steel
AliMedium(kSteel, "STEEL", material, 0, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
- AliMaterial(material=46, "PERPEX$", a=63.54,z=29.0,den=8.96, radl=1.4, absl=0); //Perpex
+ AliMaterial(++material, "PERPEX$", a=63.54,z=29.0,den=8.96, radl=1.4, absl=0); //15 Perpex
AliMedium(kPerpex, "PERPEX", material, 0, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
- AliMaterial(material=47, "STRONZIUM$", a=87.62,z=38.0,den=2.54, radl=4.24, absl=0); //Sr90
+ AliMaterial(++material, "STRONZIUM$", a=87.62,z=38.0,den=2.54, radl=4.24, absl=0); //16 Sr90
AliMedium(kSr90, "STRONZIUM", material, 0, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
}
- Float_t aGlass[5]={12.01, 28.09, 16., 10.8, 23.};
- Float_t zGlass[5]={ 6., 14., 8., 5., 11.};
- Float_t wGlass[5]={ 0.5, 0.105, 0.355, 0.03, 0.01};
- AliMixture(material=32, "GLASS",aGlass, zGlass, den=1.74, 5, wGlass);//Glass 50%C+10.5%Si+35.5%O+3% + 1%
- AliMedium(11, "GLASS", material, 0, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
-
- Int_t *idtmed = fIdtmed->GetArray()-999;
- gMC->SetCerenkov(idtmed[1000], kNbins, aPckov, aAbsCH4, aQeAll, aIdxCH4);
- gMC->SetCerenkov(idtmed[1001], kNbins, aPckov, aAbsCH4, aQeAll, aIdxCH4);
- gMC->SetCerenkov(idtmed[1002], kNbins, aPckov, aAbsSiO2, aQeAll, aIdxSiO2);
- gMC->SetCerenkov(idtmed[1003], kNbins, aPckov, aAbsC6F14, aQeAll, aIdxC6F14);
- gMC->SetCerenkov(idtmed[1004], kNbins, aPckov, aAbsCH4, aQeAll, aIdxCH4);
- gMC->SetCerenkov(idtmed[1005], kNbins, aPckov, aAbsCsI, aQeCsI, aIdxCH4);
- gMC->SetCerenkov(idtmed[1006], kNbins, aPckov, aAbsGrid, aQeAll, aIdxGrid);
- gMC->SetCerenkov(idtmed[1007], kNbins, aPckov, aAbsOpSiO2, aQeAll, aIdxOpSiO2);
- gMC->SetCerenkov(idtmed[1008], kNbins, aPckov, aAbsCH4, aQeAll, aIdxCH4);
- gMC->SetCerenkov(idtmed[1009], kNbins, aPckov, aAbsGrid, aQeAll, aIdxGrid);
- gMC->SetCerenkov(idtmed[1010], kNbins, aPckov, aAbsOpSiO2, aQeAll, aIdxOpSiO2);
-
+//Optical properties:
+#include "Opticals.h"
+ gMC->SetCerenkov((*fIdtmed)[kAir] , kNbins, aPckov, aAbsCH4, aQeAll, aIdxCH4); //1 Air
+ gMC->SetCerenkov((*fIdtmed)[kRoha] , kNbins, aPckov, aAbsCH4, aQeAll, aIdxCH4); //2 Honeycomb
+ gMC->SetCerenkov((*fIdtmed)[kSiO2] , kNbins, aPckov, aAbsSiO2, aQeAll, aIdxSiO2); //3 Quartz SiO2
+ gMC->SetCerenkov((*fIdtmed)[kC6F14] , kNbins, aPckov, aAbsC6F14, aQeAll, aIdxC6F14); //4 Freon C6F14
+ gMC->SetCerenkov((*fIdtmed)[kCH4] , kNbins, aPckov, aAbsCH4, aQeAll, aIdxCH4); //5 Methane CH4
+ gMC->SetCerenkov((*fIdtmed)[kCsI] , kNbins, aPckov, aAbsCsI, aQeCsI, aIdxCH4); //6 CsI
+ gMC->SetCerenkov((*fIdtmed)[kGridCu] , kNbins, aPckov, aAbsGrid, aQeAll, aIdxGrid); //7 grid Cu
+ gMC->SetCerenkov((*fIdtmed)[kOpSiO2] , kNbins, aPckov, aAbsOpSiO2, aQeAll, aIdxOpSiO2); //8 Opaque quartz SiO2
+ gMC->SetCerenkov((*fIdtmed)[kGap] , kNbins, aPckov, aAbsCH4, aQeAll, aIdxCH4); //9 Special methane gap
+ gMC->SetCerenkov((*fIdtmed)[kAl] , kNbins, aPckov, aAbsGrid, aQeAll, aIdxGrid); //10 Aluminium
+ gMC->SetCerenkov((*fIdtmed)[kGlass] , kNbins, aPckov, aAbsOpSiO2, aQeAll, aIdxOpSiO2); //11 Glass
}//void AliRICH::CreateMaterials()
//__________________________________________________________________________________________________
Float_t AliRICH::Fresnel(Float_t ene,Float_t pdoti, Bool_t pola)const
void AliRICH::MakeBranch(Option_t* option)
{
//Create Tree branches for the RICH.
- if(GetDebug())Info("MakeBranch","Start with option= %s.",option);
+ AliDebug(1,Form("Start with option= %s.",option));
const Int_t kBufferSize = 4000;
for(Int_t i=0;i<kNchambers;i++)
MakeBranchInTree(fLoader->TreeR(),Form("%sClusters%d",GetName(),i+1), &((*fClusters)[i]), kBufferSize, 0);
}//R
- if(GetDebug())Info("MakeBranch","Stop.");
+ AliDebug(1,"Stop.");
}//void AliRICH::MakeBranch(Option_t* option)
//__________________________________________________________________________________________________
void AliRICH::SetTreeAddress()
{
//Set branch address for the Hits and Digits Tree.
- if(GetDebug())Info("SetTreeAddress","Start.");
+ AliDebug(1,"Start.");
TBranch *branch;
if(fLoader->TreeH()){//H
- if(GetDebug())Info("SetTreeAddress","tree H is requested.");
+ AliDebug(1,"tree H is requested.");
CreateHits();//branch map will be in AliDetector::SetTreeAddress
}//H
AliDetector::SetTreeAddress();//this is after TreeH because we need to guarantee that fHits array is created
if(fLoader->TreeS()){//S
- if(GetDebug())Info("SetTreeAddress","tree S is requested.");
+ AliDebug(1,"tree S is requested.");
branch=fLoader->TreeS()->GetBranch(GetName()); if(branch){CreateSDigits(); branch->SetAddress(&fSdigits);}
}//S
if(fLoader->TreeD()){//D
- if(GetDebug())Info("SetTreeAddress","tree D is requested.");
+ AliDebug(1,"tree D is requested.");
for(int i=0;i<kNchambers;i++){
branch=fLoader->TreeD()->GetBranch(Form("%s%d",GetName(),i+1));
if(branch){CreateDigits(); branch->SetAddress(&((*fDigitsNew)[i]));}
}//D
if(fLoader->TreeR()){//R
- if(GetDebug())Info("SetTreeAddress","tree R is requested.");
+ AliDebug(1,"tree R is requested.");
for(int i=0;i<kNchambers;i++){
branch=fLoader->TreeR()->GetBranch(Form("%sClusters%d" ,GetName(),i+1));
if(branch){CreateClusters(); branch->SetAddress(&((*fClusters)[i]));}
}
}//R
- if(GetDebug())Info("SetTreeAddress","Stop.");
+ AliDebug(1,"Stop.");
}//void AliRICH::SetTreeAddress()
//__________________________________________________________________________________________________
void AliRICH::Print(Option_t *option)const
void AliRICH::CreateGeometry()
{
//Creates detailed geometry simulation (currently GEANT volumes tree)
- if(GetDebug())Info("CreateGeometry","Start main.");
-//Opaque quartz thickness
- Float_t oquaThickness = .5;
-//CsI dimensions
- Float_t pcX=P()->PcSizeX();
- Float_t pcY=P()->PcSizeY();
-
- Int_t *idtmed = fIdtmed->GetArray()-999;
-
- Int_t i;
- Float_t zs;
- Int_t idrotm[1099];
+ AliDebug(1,"Start main.");
+ Double_t cm=1,mm=0.1*cm,mkm=0.001*cm;
Float_t par[3];
- par[0]=68.8;par[1]=13 ;par[2]=70.86; gMC->Gsvolu("RICH","BOX ",(*fIdtmed)[kAl], par,3);//External aluminium box
- par[0]=66.3;par[1]=13 ;par[2]=68.35; gMC->Gsvolu("SRIC","BOX ",(*fIdtmed)[kAir],par,3);//Air
- par[0]=66.3;par[1]=0.025;par[2]=68.35; gMC->Gsvolu("ALUM","BOX ",(*fIdtmed)[kAl], par,3);//Aluminium sheet
-//Air 2 (cutting the lower part of the box)
-// par[0]=1.25; par[1] = 3; par[2] = 70.86; gMC->Gsvolu("AIR2", "BOX ", idtmed[1000], par, 3);
-//Air 3 (cutting the lower part of the box)
-// par[0]=66.3; par[1] = 3; par[2] = 1.2505; gMC->Gsvolu("AIR3", "BOX ", idtmed[1000], par, 3);
-//Honeycomb
- par[0]=66.3;par[1]=0.188; par[2] = 68.35; gMC->Gsvolu("HONE", "BOX ", idtmed[1001], par, 3);
-
- //par[0] = 66.5; par[1] = .025; par[2] = 63.1;
-//Quartz
- par[0]=P()->QuartzWidth()/2;par[1]=P()->QuartzThickness()/2;par[2]=P()->QuartzLength()/2;
- gMC->Gsvolu("QUAR", "BOX ", idtmed[1002], par, 3);
-//Spacers (cylinders)
- par[0]=0.;par[1]=.5;par[2]=P()->FreonThickness()/2; gMC->Gsvolu("SPAC", "TUBE", idtmed[1002], par, 3);
-//Feet (freon slabs supports)
- par[0] = .7; par[1] = .3; par[2] = 1.9; gMC->Gsvolu("FOOT", "BOX", idtmed[1009], par, 3);
-//Opaque quartz
- par[0]=P()->QuartzWidth()/2;par[1]= .2;par[2]=P()->QuartzLength()/2;
- gMC->Gsvolu("OQUA", "BOX ", idtmed[1007], par, 3);
-//Frame of opaque quartz
- par[0]=P()->OuterFreonWidth()/2;par[1]=P()->FreonThickness()/2;par[2]=P()->OuterFreonLength()/2;
- gMC->Gsvolu("OQF1", "BOX ", idtmed[1007], par, 3);
- par[0]=P()->InnerFreonWidth()/2;par[1]=P()->FreonThickness()/2;par[2]=P()->InnerFreonLength()/2;
- gMC->Gsvolu("OQF2", "BOX ", idtmed[1007], par, 3);
-//Freon
- par[0]=P()->OuterFreonWidth()/2 - oquaThickness;
- par[1]=P()->FreonThickness()/2;
- par[2]=P()->OuterFreonLength()/2 - 2*oquaThickness;
- gMC->Gsvolu("FRE1", "BOX ", idtmed[1003], par, 3);
-
- par[0]=P()->InnerFreonWidth()/2 - oquaThickness;
- par[1]=P()->FreonThickness()/2;
- par[2]=P()->InnerFreonLength()/2 - 2*oquaThickness;
- gMC->Gsvolu("FRE2", "BOX ", idtmed[1003], par, 3);
-//Methane
- par[0]=pcX/2;par[1]=P()->GapThickness()/2;par[2]=pcY/2; gMC->Gsvolu("META","BOX ",idtmed[1004], par, 3);
-//Methane gap
- par[0]=pcX/2;par[1]=P()->GapAmp()/2;par[2]=pcY/2;gMC->Gsvolu("GAP ","BOX ",(*fIdtmed)[kGAP],par,3);
-//CsI PC
- par[0]=pcX/2;par[1]=.25;par[2]=pcY/2; gMC->Gsvolu("CSI ", "BOX ", (*fIdtmed)[kCSI], par, 3);
-//Anode grid
- par[0] = 0.;par[1] = .001;par[2] = 20.; gMC->Gsvolu("GRID", "TUBE", idtmed[1006], par, 3);
-
-//Wire supports
-//Bar of metal
- par[0]=pcX/2;par[1]=1.05;par[2]=1.05; gMC->Gsvolu("WSMe", "BOX ", idtmed[1009], par, 3);
-//Ceramic pick up (base)
- par[0]=pcX/2;par[1]= .25;par[2]=1.05; gMC->Gsvolu("WSG1", "BOX ", idtmed[1010], par, 3);
-//Ceramic pick up (head)
- par[0] = pcX/2;par[1] = .1;par[2] = .1; gMC->Gsvolu("WSG2", "BOX ", idtmed[1010], par, 3);
-
-//Aluminium supports for methane and CsI
-//Short bar
- par[0]=pcX/2;par[1]=P()->GapThickness()/2 + .25; par[2] = (68.35 - pcY/2)/2;
- gMC->Gsvolu("SMSH", "BOX", idtmed[1009], par, 3);
-//Long bar
- par[0]=(66.3 - pcX/2)/2;par[1]=P()->GapThickness()/2+.25;par[2]=pcY/2+68.35-pcY/2;
- gMC->Gsvolu("SMLG", "BOX", idtmed[1009], par, 3);
-
-//Aluminium supports for freon
-//Short bar
- par[0] = P()->QuartzWidth()/2; par[1] = .3; par[2] = (68.35 - P()->QuartzLength()/2)/2;
- gMC->Gsvolu("SFSH", "BOX", idtmed[1009], par, 3);
-//Long bar
- par[0] = (66.3 - P()->QuartzWidth()/2)/2; par[1] = .3;
- par[2] = P()->QuartzLength()/2 + 68.35 - P()->QuartzLength()/2;
- gMC->Gsvolu("SFLG", "BOX", idtmed[1009], par, 3);
-//PCB backplane
- par[0] = pcX/2;par[1] = .25; par[2] = pcY/4 -.5025; gMC->Gsvolu("PCB ", "BOX", idtmed[1011], par, 3);
-
-//Backplane supports
-//Aluminium slab
- par[0] = 33.15;par[1] = 2;par[2] = 21.65; gMC->Gsvolu("BACK", "BOX", idtmed[1009], par, 3);
-//Big hole
- par[0] = 9.05; par[1] = 2; par[2] = 4.4625; gMC->Gsvolu("BKHL", "BOX", idtmed[1000], par, 3);
-//Small hole
- par[0] = 5.7;par[1] = 2;par[2] = 4.4625; gMC->Gsvolu("BKHS", "BOX", idtmed[1000], par, 3);
-//Place holes inside backplane support
- gMC->Gspos("BKHS", 1, "BACK", .8 + 5.7,0., .6 + 4.4625, 0, "ONLY");
- gMC->Gspos("BKHS", 2, "BACK", -.8 - 5.7,0., .6 + 4.4625, 0, "ONLY");
- gMC->Gspos("BKHS", 3, "BACK", .8 + 5.7,0., -.6 - 4.4625, 0, "ONLY");
- gMC->Gspos("BKHS", 4, "BACK", -.8 - 5.7,0., -.6 - 4.4625, 0, "ONLY");
- gMC->Gspos("BKHS", 5, "BACK", .8 + 5.7,0., .6 + 8.925 + 1.2 + 4.4625, 0, "ONLY");
- gMC->Gspos("BKHS", 6, "BACK", -.8 - 5.7,0., .6 + 8.925 + 1.2 + 4.4625, 0, "ONLY");
- gMC->Gspos("BKHS", 7, "BACK", .8 + 5.7,0., -.6 - 8.925 - 1.2 - 4.4625, 0, "ONLY");
- gMC->Gspos("BKHS", 8, "BACK", -.8 - 5.7,0., -.6 - 8.925 - 1.2 - 4.4625, 0, "ONLY");
- gMC->Gspos("BKHL", 1, "BACK", .8 + 11.4 + 1.6 + 9.05, 0., .6 + 4.4625, 0, "ONLY");
- gMC->Gspos("BKHL", 2, "BACK", -.8 - 11.4 - 1.6 - 9.05, 0., .6 + 4.4625, 0, "ONLY");
- gMC->Gspos("BKHL", 3, "BACK", .8 + 11.4 + 1.6 + 9.05, 0., -.6 - 4.4625, 0, "ONLY");
- gMC->Gspos("BKHL", 4, "BACK", -.8 - 11.4 - 1.6 - 9.05, 0., -.6 - 4.4625, 0, "ONLY");
- gMC->Gspos("BKHL", 5, "BACK", .8 + 11.4+ 1.6 + 9.05, 0., .6 + 8.925 + 1.2 + 4.4625, 0, "ONLY");
- gMC->Gspos("BKHL", 6, "BACK", -.8 - 11.4 - 1.6 - 9.05, 0., .6 + 8.925 + 1.2 + 4.4625, 0, "ONLY");
- gMC->Gspos("BKHL", 7, "BACK", .8 + 11.4 + 1.6 + 9.05, 0., -.6 - 8.925 - 1.2 - 4.4625, 0, "ONLY");
- gMC->Gspos("BKHL", 8, "BACK", -.8 - 11.4 - 1.6 - 9.05, 0., -.6 - 8.925 - 1.2 - 4.4625, 0, "ONLY");
-//Place material inside RICH
- gMC->Gspos("SRIC", 1, "RICH", 0.,0., 0., 0, "ONLY");
-// gMC->Gspos("AIR2", 1, "RICH", 66.3 + 1.2505, 1.276-P()->GapThickness()/2-P()->QuartzThickness()-P()->FreonThickness()- .4 - .6 - .05 - .376 -.5 - 3.35, 0., 0, "ONLY");
-// gMC->Gspos("AIR2", 2, "RICH", -66.3 - 1.2505,1.276-P()->GapThickness()/2-P()->QuartzThickness()-P()->FreonThickness()- .4 - .6 - .05 - .376 -.5 - 3.35, 0., 0, "ONLY");
-// gMC->Gspos("AIR3", 1, "RICH", 0., 1.276-P()->GapThickness()/2 - P()->QuartzThickness() - P()->FreonThickness()- .4 - .6 - .05 - .376 -.5 - 3.35, -68.35 - 1.25, 0, "ONLY");
-// gMC->Gspos("AIR3", 2, "RICH", 0., 1.276 - P()->GapThickness()/2 - P()->QuartzThickness() - P()->FreonThickness()- .4 - .6 - .05 - .376 -.5 - 3.35, 68.35 + 1.25, 0, "ONLY");
- gMC->Gspos("ALUM", 1, "SRIC", 0., 1.276 - P()->GapThickness()/2 - P()->QuartzThickness() - P()->FreonThickness()- .4 - .6 - .05 - .376 -.025, 0., 0, "ONLY");
- gMC->Gspos("HONE", 1, "SRIC", 0., 1.276- P()->GapThickness()/2 - P()->QuartzThickness() - P()->FreonThickness()- .4 - .6 - .05 - .188, 0., 0, "ONLY");
- gMC->Gspos("ALUM", 2, "SRIC", 0., 1.276 - P()->GapThickness()/2 - P()->QuartzThickness() - P()->FreonThickness()- .4 - .6 - .025, 0., 0, "ONLY");
- gMC->Gspos("FOOT", 1, "SRIC", 64.95, 1.276 - P()->GapThickness()/2 - P()->QuartzThickness() - P()->FreonThickness()- .4 - .3, 36.9, 0, "ONLY");
- gMC->Gspos("FOOT", 2, "SRIC", 21.65, 1.276 - P()->GapThickness()/2 - P()->QuartzThickness() - P()->FreonThickness()- .4 - .3 , 36.9, 0, "ONLY");
- gMC->Gspos("FOOT", 3, "SRIC", -21.65, 1.276 - P()->GapThickness()/2 - P()->QuartzThickness() - P()->FreonThickness()- .4 - .3, 36.9, 0, "ONLY");
- gMC->Gspos("FOOT", 4, "SRIC", -64.95, 1.276 - P()->GapThickness()/2 - P()->QuartzThickness() - P()->FreonThickness()- .4 - .3, 36.9, 0, "ONLY");
- gMC->Gspos("FOOT", 5, "SRIC", 64.95, 1.276 - P()->GapThickness()/2 - P()->QuartzThickness() - P()->FreonThickness()- .4 - .3, -36.9, 0, "ONLY");
- gMC->Gspos("FOOT", 6, "SRIC", 21.65, 1.276 - P()->GapThickness()/2 - P()->QuartzThickness() - P()->FreonThickness()- .4 - .3, -36.9, 0, "ONLY");
- gMC->Gspos("FOOT", 7, "SRIC", -21.65, 1.276 - P()->GapThickness()/2 - P()->QuartzThickness() - P()->FreonThickness()- .4 - .3, -36.9, 0, "ONLY");
- gMC->Gspos("FOOT", 8, "SRIC", -64.95, 1.276 - P()->GapThickness()/2 - P()->QuartzThickness() - P()->FreonThickness()- .4 - .3, -36.9, 0, "ONLY");
- gMC->Gspos("OQUA", 1, "SRIC", 0., 1.276 - P()->GapThickness()/2 - P()->QuartzThickness() - P()->FreonThickness()- .2, 0., 0, "ONLY");
-// Methane supports
- gMC->Gspos("SMLG", 1, "SRIC", pcX/2 + (66.3 - pcX/2)/2, 1.276 + .25, 0., 0, "ONLY");
- gMC->Gspos("SMLG", 2, "SRIC", - pcX/2 - (66.3 - pcX/2)/2, 1.276 + .25, 0., 0, "ONLY");
- gMC->Gspos("SMSH", 1, "SRIC", 0., 1.276 + .25, pcY/2 + (68.35 - pcY/2)/2, 0, "ONLY");
- gMC->Gspos("SMSH", 2, "SRIC", 0., 1.276 + .25, - pcY/2 - (68.35 - pcY/2)/2, 0, "ONLY");
-//Freon supports
- Float_t suppY = 1.276 - P()->GapThickness()/2- P()->QuartzThickness() -P()->FreonThickness() - .2 + .3; //y position of freon supports
- gMC->Gspos("SFLG", 1, "SRIC", P()->QuartzWidth()/2 + (66.3 - P()->QuartzWidth()/2)/2, suppY, 0., 0, "ONLY");
- gMC->Gspos("SFLG", 2, "SRIC", - P()->QuartzWidth()/2 - (66.3 - P()->QuartzWidth()/2)/2, suppY, 0., 0, "ONLY");
- gMC->Gspos("SFSH", 1, "SRIC", 0., suppY, P()->QuartzLength()/2 + (68.35 - P()->QuartzLength()/2)/2, 0, "ONLY");
- gMC->Gspos("SFSH", 2, "SRIC", 0., suppY, - P()->QuartzLength()/2 - (68.35 - P()->QuartzLength()/2)/2, 0, "ONLY");
- AliMatrix(idrotm[1019], 0., 0., 90., 0., 90., 90.);
-//Place spacers
- Int_t nspacers = 30;
- for (i = 0; i < nspacers/3; i++) {
- zs = -11.6/2 + (TMath::Abs(nspacers/6) - i) * 12.2;
- gMC->Gspos("SPAC", i, "FRE1", 10.5, 0., zs, idrotm[1019], "ONLY"); //Original settings
- }
- for (i = nspacers/3; i < (nspacers*2)/3; i++) {
- zs = -11.6/2 + (nspacers/3 + TMath::Abs(nspacers/6) - i) * 12.2;
- gMC->Gspos("SPAC", i, "FRE1", 0, 0., zs, idrotm[1019], "ONLY"); //Original settings
- }
- for (i = (nspacers*2)/3; i < nspacers; ++i) {
- zs = -11.6/2 + ((nspacers*2)/3 + TMath::Abs(nspacers/6) - i) * 12.2;
- gMC->Gspos("SPAC", i, "FRE1", -10.5, 0., zs, idrotm[1019], "ONLY"); //Original settings
- }
- for (i = 0; i < nspacers/3; i++) {
- zs = -11.6/2 + (TMath::Abs(nspacers/6) - i) * 12.2;
- gMC->Gspos("SPAC", i, "FRE2", 10.5, 0., zs, idrotm[1019], "ONLY"); //Original settings
- }
- for (i = nspacers/3; i < (nspacers*2)/3; i++) {
- zs = -11.6/2 + (nspacers/3 + TMath::Abs(nspacers/6) - i) * 12.2;
- gMC->Gspos("SPAC", i, "FRE2", 0, 0., zs, idrotm[1019], "ONLY"); //Original settings
- }
- for (i = (nspacers*2)/3; i < nspacers; ++i) {
- zs = -11.6/2 + ((nspacers*2)/3 + TMath::Abs(nspacers/6) - i) * 12.2;
- gMC->Gspos("SPAC", i, "FRE2", -10.5, 0., zs, idrotm[1019], "ONLY"); //Original settings
- }
- gMC->Gspos("FRE1", 1, "OQF1", 0., 0., 0., 0, "ONLY");
- gMC->Gspos("FRE2", 1, "OQF2", 0., 0., 0., 0, "ONLY");
- gMC->Gspos("OQF1", 1, "SRIC", P()->OuterFreonWidth()/2 + P()->InnerFreonWidth()/2 + 2, 1.276 - P()->GapThickness()/2- P()->QuartzThickness() -P()->FreonThickness()/2, 0., 0, "ONLY"); //Original settings (31.3)
- gMC->Gspos("OQF2", 2, "SRIC", 0., 1.276 - P()->GapThickness()/2 - P()->QuartzThickness() - P()->FreonThickness()/2, 0., 0, "ONLY"); //Original settings
- gMC->Gspos("OQF1", 3, "SRIC", - (P()->OuterFreonWidth()/2 + P()->InnerFreonWidth()/2) - 2, 1.276 - P()->GapThickness()/2 - P()->QuartzThickness() - P()->FreonThickness()/2, 0., 0, "ONLY"); //Original settings (-31.3)
- gMC->Gspos("QUAR", 1, "SRIC", 0., 1.276 - P()->GapThickness()/2 - P()->QuartzThickness()/2, 0., 0, "ONLY");
- gMC->Gspos("GAP ", 1, "META", 0., P()->GapThickness()/2 - P()->GapAmp()/2 - 0.0001, 0., 0, "ONLY");
- gMC->Gspos("META", 1, "SRIC", 0., 1.276, 0., 0, "ONLY");
- gMC->Gspos("CSI ", 1, "SRIC", 0., 1.276 + P()->GapThickness()/2 + .25, 0., 0, "ONLY");
-//Wire support placing
- gMC->Gspos("WSG2", 1, "GAP ", 0., P()->GapAmp()/2 - .1, 0., 0, "ONLY");
- gMC->Gspos("WSG1", 1, "CSI ", 0., 0., 0., 0, "ONLY");
- gMC->Gspos("WSMe", 1, "SRIC ", 0., 1.276 + P()->GapThickness()/2 + .5 + 1.05, 0., 0, "ONLY");
-//Backplane placing
- gMC->Gspos("BACK", 1, "SRIC ", -33.15, 1.276 + P()->GapThickness()/2 + .5 + 2.1 + 2, 43.3, 0, "ONLY");
- gMC->Gspos("BACK", 2, "SRIC ", 33.15, 1.276 + P()->GapThickness()/2 + .5 + 2.1 + 2 , 43.3, 0, "ONLY");
- gMC->Gspos("BACK", 3, "SRIC ", -33.15, 1.276 + P()->GapThickness()/2 + .5 + 2.1 + 2, 0., 0, "ONLY");
- gMC->Gspos("BACK", 4, "SRIC ", 33.15, 1.276 + P()->GapThickness()/2 + .5 + 2.1 + 2, 0., 0, "ONLY");
- gMC->Gspos("BACK", 5, "SRIC ", 33.15, 1.276 + P()->GapThickness()/2 + .5 + 2.1 + 2, -43.3, 0, "ONLY");
- gMC->Gspos("BACK", 6, "SRIC ", -33.15, 1.276 + P()->GapThickness()/2 + .5 + 2.1 + 2, -43.3, 0, "ONLY");
-//PCB placing
- gMC->Gspos("PCB ", 1, "SRIC ", 0., 1.276 + P()->GapThickness()/2 + .5 + 1.05, pcX/4 + .5025 + 2.5, 0, "ONLY");
- gMC->Gspos("PCB ", 2, "SRIC ", 0., 1.276 + P()->GapThickness()/2 + .5 + 1.05, -pcX/4 - .5025 - 2.5, 0, "ONLY");
-
+ Int_t id=0;
+
//place chambers into mother volume ALIC
- for(int i=1;i<=kNchambers;i++){
- AliMatrix(idrotm[1000+i],C(i)->ThetaXd(),C(i)->PhiXd(),
- C(i)->ThetaYd(),C(i)->PhiYd(),
- C(i)->ThetaZd(),C(i)->PhiZd());
- gMC->Gspos("RICH",i,"ALIC",C(i)->X(),C(i)->Y(),C(i)->Z(),idrotm[1000+i], "ONLY");
- }
+ par[0]=(6*mm+1681*mm+6*mm)/2;par[1]=(6*mm+1466*mm+6*mm)/2;par[2]=(80*mm+40*mm)*2/2;gMC->Gsvolu("RICH","BOX ",(*fIdtmed)[kCH4],par,3);//2033P1 z p84 TDR
+ if(P()->IsRadioSrc()){
+ AliInfo("Special test beam geometry");
+ gMC->Gspos("RICH",1,"ALIC",0,0,0,0, "ONLY"); //single RICH chamber without rotation test beam geometry
+ }else
+ for(int i=1;i<=kNchambers;i++){ //full RICH with 7 chambers
+ AliMatrix(id,C(i)->ThetaXd(),C(i)->PhiXd(), C(i)->ThetaYd(),C(i)->PhiYd(), C(i)->ThetaZd(),C(i)->PhiZd());
+ gMC->Gspos("RICH",i,"ALIC",C(i)->Center().X(),C(i)->Center().Y(),C(i)->Center().Z(),id, "ONLY");
+ }
+//Pad Panel frame 6 sectors
+ par[0]=648*mm/2;par[1]= 411*mm/2;par[2]=40 *mm/2;gMC->Gsvolu("RPPF","BOX ",(*fIdtmed)[kAl] ,par,3);//PPF 2001P2 inner size of the slab by 1mm more
+ par[0]=181*mm/2;par[1]=89.25*mm/2;par[2]=38.3*mm/2;gMC->Gsvolu("PPFL","BOX ",(*fIdtmed)[kAir] ,par,3);//large whole
+ par[0]=114*mm/2;par[1]=89.25*mm/2;par[2]=38.3*mm/2;gMC->Gsvolu("PPFS","BOX ",(*fIdtmed)[kAir] ,par,3);//small whole
+ par[0]=644*mm/2;par[1]= 407*mm/2;par[2]= 1.7*mm/2;gMC->Gsvolu("RPC ","BOX ",(*fIdtmed)[kCsI] ,par,3);//by 0.2 mm more then actual size (PCB 2006P1)
+
+ gMC->Gspos("RPPF",1,"RICH", -335*mm, -433*mm, 8*cm+20*mm, 0,"ONLY");//F1 2040P1 z p.84 TDR
+ gMC->Gspos("RPPF",2,"RICH", +335*mm, -433*mm, 8*cm+20*mm, 0,"ONLY");
+ gMC->Gspos("RPPF",3,"RICH", -335*mm, 0*mm, 8*cm+20*mm, 0,"ONLY");
+ gMC->Gspos("RPPF",4,"RICH", +335*mm, 0*mm, 8*cm+20*mm, 0,"ONLY");
+ gMC->Gspos("RPPF",5,"RICH", -335*mm, +433*mm, 8*cm+20*mm, 0,"ONLY");
+ gMC->Gspos("RPPF",6,"RICH", +335*mm, +433*mm, 8*cm+20*mm, 0,"ONLY");
+ gMC->Gspos("RPC ",1,"RPPF", 0*mm, 0*mm, -19.15*mm, 0,"ONLY");//PPF 2001P2
+ gMC->Gspos("PPFL",1,"RPPF", -224.5*mm, -151.875*mm, 0.85*mm, 0,"ONLY");
+ gMC->Gspos("PPFL",2,"RPPF", -224.5*mm, - 50.625*mm, 0.85*mm, 0,"ONLY");
+ gMC->Gspos("PPFL",3,"RPPF", -224.5*mm, + 50.625*mm, 0.85*mm, 0,"ONLY");
+ gMC->Gspos("PPFL",4,"RPPF", -224.5*mm, +151.875*mm, 0.85*mm, 0,"ONLY");
+ gMC->Gspos("PPFS",1,"RPPF", - 65.0*mm, -151.875*mm, 0.85*mm, 0,"ONLY");
+ gMC->Gspos("PPFS",2,"RPPF", - 65.0*mm, - 50.625*mm, 0.85*mm, 0,"ONLY");
+ gMC->Gspos("PPFS",3,"RPPF", - 65.0*mm, + 50.625*mm, 0.85*mm, 0,"ONLY");
+ gMC->Gspos("PPFS",4,"RPPF", - 65.0*mm, +151.875*mm, 0.85*mm, 0,"ONLY");
+ gMC->Gspos("PPFS",5,"RPPF", + 65.0*mm, -151.875*mm, 0.85*mm, 0,"ONLY");
+ gMC->Gspos("PPFS",6,"RPPF", + 65.0*mm, - 50.625*mm, 0.85*mm, 0,"ONLY");
+ gMC->Gspos("PPFS",7,"RPPF", + 65.0*mm, + 50.625*mm, 0.85*mm, 0,"ONLY");
+ gMC->Gspos("PPFS",8,"RPPF", + 65.0*mm, +151.875*mm, 0.85*mm, 0,"ONLY");
+ gMC->Gspos("PPFL",5,"RPPF", +224.5*mm, -151.875*mm, 0.85*mm, 0,"ONLY");
+ gMC->Gspos("PPFL",6,"RPPF", +224.5*mm, - 50.625*mm, 0.85*mm, 0,"ONLY");
+ gMC->Gspos("PPFL",7,"RPPF", +224.5*mm, + 50.625*mm, 0.85*mm, 0,"ONLY");
+ gMC->Gspos("PPFL",8,"RPPF", +224.5*mm, +151.875*mm, 0.85*mm, 0,"ONLY");
+//Gap - anod wires 6 copies to RICH
+ par[0]=648*mm/2;par[1]= 411*mm/2 ;par[2]=4.45*mm/2;gMC->Gsvolu("RGAP ","BOX ",(*fIdtmed)[kCH4] ,par,3);//xy as PPF 2001P2 z WP 2099P1
+ par[0]= 0*mm ;par[1]= 20*mkm/2 ;par[2]= 648*mm/2;gMC->Gsvolu("RANO","TUBE",(*fIdtmed)[kW] ,par,3);//WP 2099P1 z = gap x PPF 2001P2
+ AliMatrix(id,180,0, 90,90, 90,0); //wires along x
+
+ gMC->Gspos("RGAP",1,"RICH", -335*mm, -433*mm,8*cm-2.225*mm, 0,"ONLY"); //F1 2040P1 z WP 2099P1
+ gMC->Gspos("RGAP",2,"RICH", +335*mm, -433*mm,8*cm-2.225*mm, 0,"ONLY");
+ gMC->Gspos("RGAP",3,"RICH", -335*mm, 0*mm,8*cm-2.225*mm, 0,"ONLY");
+ gMC->Gspos("RGAP",4,"RICH", +335*mm, 0*mm,8*cm-2.225*mm, 0,"ONLY");
+ gMC->Gspos("RGAP",5,"RICH", -335*mm, +433*mm,8*cm-2.225*mm, 0,"ONLY");
+ gMC->Gspos("RGAP",6,"RICH", +335*mm, +433*mm,8*cm-2.225*mm, 0,"ONLY");
+ for(int i=1;i<=96;i++)
+ gMC->Gspos("RANO",i,"RGAP", 0*mm, -411/2*mm+i*4*mm, 0.185*mm, id,"ONLY"); //WP 2099P1
+//Radiator 3 modules
+ par[0]=1330*mm/2 ;par[1]= 413*mm/2 ;par[2]= 24*mm/2; gMC->Gsvolu("RRAD","BOX ",(*fIdtmed)[kC6F14] ,par,3); // Rad 2011P1
+ par[0]=1330*mm/2 ;par[1]= 413*mm/2 ;par[2]= 4*mm/2; gMC->Gsvolu("RRFR","BOX ",(*fIdtmed)[kRoha] ,par,3); //front
+ par[0]=1330*mm/2 ;par[1]= 413*mm/2 ;par[2]= 5*mm/2; gMC->Gsvolu("RRWI","BOX ",(*fIdtmed)[kSiO2] ,par,3); //window
+ par[0]=1330*mm/2 ;par[1]= 5*mm/2 ;par[2]= 15*mm/2; gMC->Gsvolu("RRLO","BOX ",(*fIdtmed)[kRoha] ,par,3); //long side
+ par[0]= 10*mm/2 ;par[1]= 403*mm/2 ;par[2]= 15*mm/2; gMC->Gsvolu("RRSH","BOX ",(*fIdtmed)[kRoha] ,par,3); //short side
+ par[0]= 0 ;par[1]= 10*mm/2 ;par[2]= 15*mm/2; gMC->Gsvolu("RRSP","TUBE",(*fIdtmed)[kSiO2] ,par,3); //spacer
+
+ gMC->Gspos("RRAD",1,"RICH", 0*mm,-434*mm, -12*mm, 0,"ONLY"); //radiator to RICH
+ gMC->Gspos("RRAD",2,"RICH", 0*mm, 0*mm, -12*mm, 0,"ONLY");
+ gMC->Gspos("RRAD",3,"RICH", 0*mm,+434*mm, -12*mm, 0,"ONLY");
+ gMC->Gspos("RRFR",1,"RRAD", 0*mm, 0*mm, -10.0*mm, 0,"ONLY"); //front cover
+ gMC->Gspos("RRWI",1,"RRAD", 0*mm, 0*mm, 9.5*mm, 0,"ONLY"); //quartz window (back cover)
+ gMC->Gspos("RRLO",1,"RRAD", 0*mm,-204*mm, -0.5*mm, 0,"ONLY"); //long side
+ gMC->Gspos("RRLO",2,"RRAD", 0*mm,+204*mm, -0.5*mm, 0,"ONLY"); //long side
+ gMC->Gspos("RRSH",1,"RRAD",-660*mm, 0*mm, -0.5*mm, 0,"ONLY"); //short side
+ gMC->Gspos("RRSH",2,"RRAD",+660*mm, 0*mm, -0.5*mm, 0,"ONLY"); //short side
+ for(int i=0;i<3;i++)
+ for(int j=0;j<10;j++)
+ gMC->Gspos("RRSP",10*i+j,"RRAD",-1330*mm/2+116*mm+j*122*mm,(i-1)*105*mm,-0.5*mm,0,"ONLY");
-
- if(GetDebug())Info("CreateGeometry","Stop main.");
+//Sandbox
+ par[0]=1419*mm/2 ;par[1]=1378*mm/2;par[2]=50.5*mm/2; gMC->Gsvolu("RSNB","BOX ",(*fIdtmed)[kAir] ,par,3); //2072P1
+ par[0]=1419*mm/2 ;par[1]=1378*mm/2;par[2]= 0.5*mm/2; gMC->Gsvolu("RSCO","BOX ",(*fIdtmed)[kAl] ,par,3); //cover
+ par[0]=1359*mm/2 ;par[1]=1318*mm/2;par[2]=49.5*mm/2; gMC->Gsvolu("RSHO","BOX ",(*fIdtmed)[kRoha] ,par,3); //honeycomb structure
+
+ gMC->Gspos("RSNB",1,"RICH", 0*mm, 0*mm, -73.75*mm, 0,"ONLY"); //p.84 TDR sandbox to rich
+ gMC->Gspos("RSHO",1,"RSNB", 0*mm, 0*mm, 0*mm, 0,"ONLY"); //2072P1 honeycomv to sandbox
+ gMC->Gspos("RSCO",1,"RSNB", 0*mm, 0*mm, +25*mm, 0,"ONLY"); //cover to sandbox
+ gMC->Gspos("RSCO",2,"RSNB", 0*mm, 0*mm, -25*mm, 0,"ONLY"); //cover to sandbox
+
+ AliDebug(1,"Stop main.");
}//void AliRICH::CreateGeometry()
//__________________________________________________________________________________________________
-void AliRICH::Reconstruct()const
-{
- AliRICHClusterFinder finder(const_cast<AliRICH*>(this));
- finder.Exec();
-}
-//__________________________________________________________________________________________________
void AliRICH::ControlPlots()
{
-//Creates a set of hists to control the results of simulation
+// Creates a set of hists to control the results of simulation. Hists are in file $HOME/RCP.root
TH1F *pHxD=0,*pHyD=0,*pNumClusH1=0,
*pQdcH1=0, *pSizeH1=0,
Bool_t isDig =!GetLoader()->LoadDigits();
Bool_t isClus=!GetLoader()->LoadRecPoints();
- if(!isDig && !isClus){Error("ControlPlots","No digits and clusters! Nothing to do.");return;}
+ if(!isDig && !isClus){AliError("No digits and clusters! Nothing to do.");return;}
TStopwatch sw;TDatime time;
TFile *pFile = new TFile("$(HOME)/RCP.root","RECREATE");
if(isDig){
- cout<<"Digits available\n";
+ AliInfo("Digits available");
pHxD=new TH1F("HitDigitDiffX","Hit-Digits diff X all chambers;diff [cm]",100,-10,10);
pHyD=new TH1F("HitDigitDiffY","Hit-Digits diff Y all chambers;diff [cm]",100,-10,10);
}//isDig
if(isDig){
for(Int_t iDigN=0;iDigN<R()->Digits(iChamN)->GetEntries();iDigN++){//digits loop
AliRICHdigit *pDig=(AliRICHdigit*)R()->Digits(iChamN)->At(iDigN);
- AliRICHhit *pHit=Hit(pDig->GetTrack(0));
- TVector2 hitV2=R()->C(iChamN)->Glob2Loc(pHit->OutX3()); TVector2 digV2=R()->P()->Pad2Loc(pDig->Pad());
+ AliRICHhit *pHit=Hit(pDig->GetTrack(0));//get first hit of this digit
+ TVector2 hitV2=R()->C(iChamN)->Mrs2Pc(pHit->OutX3()); TVector2 digV2=R()->P()->Pad2Loc(pDig->Pad());//center of pad for digit
pHxD->Fill(hitV2.X()-digV2.X()); pHyD->Fill(hitV2.Y()-digV2.Y());
}//digits loop
}//isDig
void AliRICH::PrintHits(Int_t iEvtN)
{
//Prints a list of RICH hits for a given event. Default is event number 0.
- Info("PrintHits","List of RICH hits for event %i",iEvtN);
+ AliInfo(Form("List of RICH hits for event %i",iEvtN));
R()->GetLoader()->GetRunLoader()->GetEvent(iEvtN);
if(R()->GetLoader()->LoadHits()) return;
iTotalHits+=R()->Hits()->GetEntries();
}
R()->GetLoader()->UnloadHits();
- Info("PrintHits","totally %i hits",iTotalHits);
+ AliInfo(Form("totally %i hits",iTotalHits));
}
//__________________________________________________________________________________________________
void AliRICH::PrintSDigits(Int_t iEvtN)
Info("PrintClusters","totally %i clusters",iTotalClusters);
}
//__________________________________________________________________________________________________
-void AliRICH::FillESD(AliESD *pESD)const
+void AliRICH::PrintTracks(Int_t iEvtN)
{
-//This methode fills AliESDtrack with information from RICH
- Info("FillESD","Start with %i tracks",pESD->GetNumberOfTracks());
- const Double_t masses[5]={0.000511,0.105658,0.139567,0.493677,0.93828};//electron,muon,pion,kaon,proton
- const Double_t refIndex = 1.29052;
-
- Double_t thetaExp = 0.7;
- Double_t thetaTh[5];
- Double_t sinThetaThNorm;
- Double_t sigmaThetaTh[5];
- Double_t height[5];
- Double_t totalHeight=0;
-
- for(Int_t iTrackN=0;iTrackN<pESD->GetNumberOfTracks();iTrackN++){//ESD tracks loop
- AliESDtrack *pTrack = pESD->GetTrack(iTrackN);
-
- pTrack->Print("");
-// TVector2 x2=P()->HelixCross(pTrack);//returns cross point of track with RICH PC in LRS
- Double_t pmod = pTrack->GetP();
-
- for(Int_t iPart=4;iPart>=0;iPart--){
- Double_t cosThetaTh = TMath::Sqrt(masses[iPart]*masses[iPart]+pmod*pmod)/(refIndex*pmod);
- if(cosThetaTh>=1) {break;}
- thetaTh[iPart] = TMath::ACos(cosThetaTh);
- sinThetaThNorm = TMath::Sin(thetaTh[iPart])/TMath::Sqrt(1-1/(refIndex*refIndex));
- sigmaThetaTh[iPart] = (0.014*(1/sinThetaThNorm-1) + 0.0043)*1.25;
- height[iPart] = TMath::Gaus(thetaExp,thetaTh[iPart],sigmaThetaTh[iPart]);
- totalHeight +=height[iPart];
- }
-
- pTrack->SetRICHsignal(thetaExp);
- Double_t richPID[5];
- for(Int_t iPart=0;iPart<5;iPart++) richPID[iPart] = height[iPart]/totalHeight;
- pTrack->SetRICHpid(richPID);
- }//ESD tracks loop
+//prints a list of tracks (including secondaru) for a given event
+ Info("PrintTracks","List of all tracks for event %i",iEvtN);
+ R()->GetLoader()->GetRunLoader()->GetEvent(iEvtN);
+ if(R()->GetLoader()->GetRunLoader()->LoadHeader()) return;
+ if(R()->GetLoader()->GetRunLoader()->LoadKinematics()) return;
+ Int_t iTracksCounter=0;
+ Info("PrintTracks","totally %i tracks",iTracksCounter);
}
public:
enum ClusterStatus {kEdge,kShape,kSize,kRaw,kResolved,kEmpty=kBad};
AliRICHcluster():TObject(),fCFM(0),fSize(0),fShape(0),fQdc(0),fChamber(0),fX(0),fY(0),fStatus(kEmpty),fDigits(0) {}
- virtual ~AliRICHcluster() {Reset();}
+ virtual ~AliRICHcluster() {AliDebug(1,"Start");/*Reset();*/}
void Reset() {DeleteDigits();fCFM=fSize=fShape=fQdc=fChamber=0;fX=fY=0;fStatus=kEmpty;} //cleans the cluster
- void DeleteDigits() {if(fDigits) delete fDigits; fDigits=0;} //deletes the list of digits
+ void DeleteDigits() {if(fDigits) {delete fDigits;} fDigits=0;} //deletes the list of digits
AliRICHcluster& operator=(const AliRICHcluster&) {return *this;}
Int_t Nlocals() const{return fSize-10000*(fSize/10000);} //number of local maximums
Int_t Size() const{return fSize/10000;} //number of digits in cluster
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
- 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;} //
+ virtual void Print(Option_t *option="")const; //
+ inline void AddDigit(AliRICHdigit *pDig); //
+ inline void CoG(Int_t nLocals); //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;}
+ 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
protected:
Int_t fCFM; //1000000*Ncerenkovs+1000*Nfeedbacks+Nmips
Int_t fSize; //10000*(how many digits belong to this cluster) + nLocalMaxima
virtual void StepManager() =0; //interface from AliMC
virtual void Hits2SDigits(); //interface from AliSimulation
virtual AliDigitizer* CreateDigitizer(AliRunDigitizer* man) const {return new AliRICHDigitizer(man);} //interface from AliSimulation
- virtual void Reconstruct() const; //interface from AliReconstruction
- virtual void FillESD(AliESD *pESD) const; //interface from AliReconstruction
- virtual void Print(Option_t *option) const; //prints current RICH status
+// virtual void Reconstruct() const; //interface from AliReconstruction
+// virtual void FillESD(AliESD *pESD) const; //interface from AliReconstruction
virtual void SetTreeAddress(); //interface from AliLoader
virtual void MakeBranch(Option_t *opt=" "); //interface from AliLoader
virtual void CreateMaterials(); //interface from AliMC
AliRICH* R() {return this;} //provides pointer to RICH main object
TVector Counters() const{return fCounters;} //provides a set of counters
void ControlPlots(); //utility
+ virtual void Print(Option_t *option="") const; //prints current RICH status
void PrintHits (Int_t iEvent=0); //utility
void PrintSDigits (Int_t iEvent=0); //utility
void PrintDigits (Int_t iEvent=0); //utility
void PrintClusters(Int_t iEvent=0); //utility
+ void PrintTracks (Int_t iEvent=0); //utility
void AddHit(Int_t c,Int_t tid,TVector3 i3,TVector3 o3,Double_t eloss=0){TClonesArray &tmp=*fHits;new(tmp[fNhits++])AliRICHhit(c,tid,i3,o3,eloss);}
inline void AddSDigit(Int_t c,TVector pad,Double_t q,Int_t pid,Int_t tid);
{Int_t c=cl.C()-1;TClonesArray &tmp=*((TClonesArray*)fClusters->At(c));new(tmp[fNclusters[c]++])AliRICHcluster(cl);}
AliRICHhit* Hit(Int_t tid); //returns pointer ot RICH hit for a given tid
protected:
- enum {kAir=1,kCSI=6,kGAP=9,kAl=10,kCH4=5,kSteel=15,kPerpex=16,kSr90=17};
+ enum {kAir=1,kRoha,kSiO2,kC6F14,kCH4,kCsI,kGridCu,kOpSiO2,kGap,kAl,kGlass,kCu,kW,kSteel,kPerpex,kSr90};
AliRICHParam *fpParam; //main RICH parametrization
//fHits and fDigits belong to AliDetector
TClonesArray *fSdigits; //! list of sdigits
void AliRICH::CreateHits()
{
if(fHits) return;
- if(GetDebug())Info("CreateHits","creating hits container.");
+ AliDebug(1,"creating hits container.");
fHits=new TClonesArray("AliRICHhit",10000); fNhits=0;
}
//__________________________________________________________________________________________________
void AliRICH::CreateSDigits()
{
if(fSdigits) return;
- if(GetDebug())Info("CreateSDigits","creating sdigits container.");
+ AliDebug(1,"creating sdigits container.");
fSdigits=new TClonesArray("AliRICHdigit",10000); fNsdigits=0;
}
//__________________________________________________________________________________________________
void AliRICH::CreateDigits()
{
if(fDigitsNew) return;
- if(GetDebug())Info("CreateDigits","creating digits containers.");
+ AliDebug(1,"creating digits containers.");
fDigitsNew = new TObjArray(kNchambers);
for(Int_t i=0;i<kNchambers;i++) {fDigitsNew->AddAt(new TClonesArray("AliRICHdigit",10000), i); fNdigitsNew[i]=0;}
}
void AliRICH::CreateClusters()
{
if(fClusters) return;
- if(GetDebug())Info("CreateClusters","creating clusters containers.");
+ AliDebug(1,"creating clusters containers.");
fClusters = new TObjArray(kNchambers);
for(Int_t i=0;i<kNchambers;i++) {fClusters->AddAt(new TClonesArray("AliRICHcluster",10000), i); fNclusters[i]=0;}
}
#include "AliRICHChamber.h"
#include <TRotMatrix.h>
+#include <AliLog.h>
ClassImp(AliRICHChamber)
//______________________________________________________________________________
AliRICHChamber::AliRICHChamber(Int_t iChamber):TNamed()
{
//main ctor. Defines all geometry parameters for the given module.
- SetToZenith();//put to up position
+// 7 6 ^
+// 5 4 3 |
+// 2 1 |
+// <-----z y . x
+// horizontal angle between chambers 19.5 grad
+// vertical angle between chambers 20 grad
+ RotY(90);//rotate around y
+ fCenterX3.SetXYZ(490,0,0);fPcX3.SetXYZ(490+8,0,0);fRadX3.SetXYZ(490-2,0,0); //shift center along x by 490 cm
+
switch(iChamber){
- case 1:
- RotateX(-AliRICHParam::AngleYZ());
- RotateZ(-AliRICHParam::AngleXY());
- fName="RICHc1";fTitle="RICH chamber 1";
+ case 0: //special test beam configuration without rotation.
+ break;
+ case 1:
+ RotY(19.5); RotZ(-20); //right and down
break;
case 2:
- RotateZ(-AliRICHParam::AngleXY());
- fName="RICHc2";fTitle="RICH chamber 2";
+ RotZ(-20); //down
break;
case 3:
- RotateX(-AliRICHParam::AngleYZ());
- fName="RICHc3";fTitle="RICH chamber 3";
+ RotY(19.5); //right
break;
case 4:
- fName="RICHc4";fTitle="RICH chamber 4"; //no turns
- break;
+ break; //no rotation
case 5:
- RotateX(AliRICHParam::AngleYZ());
- fName="RICHc5";fTitle="RICH chamber 5";
+ RotY(-19.5); //left
break;
case 6:
- RotateZ(AliRICHParam::AngleXY());
- fName="RICHc6";fTitle="RICH chamber 6";
+ RotZ(20); //up
break;
case 7:
- RotateX(AliRICHParam::AngleYZ());
- RotateZ(AliRICHParam::AngleXY());
- fName="RICHc7";fTitle="RICH chamber 7";
+ RotY(-19.5); RotZ(20); //left and up
break;
default:
Fatal("named ctor","Wrong chamber number %i, check CreateChamber ctor",iChamber);
}//switch(iModuleN)
- RotateZ(AliRICHParam::AngleRot());//apply common rotation
+ fName=Form("RICHc%i",iChamber);fTitle=Form("RICH chamber %i",iChamber);
+ RotZ(30); //apply common rotation
fpRotMatrix=new TRotMatrix("rot"+fName,"rot"+fName, Rot().ThetaX()*TMath::RadToDeg(), Rot().PhiX()*TMath::RadToDeg(),
Rot().ThetaY()*TMath::RadToDeg(), Rot().PhiY()*TMath::RadToDeg(),
Rot().ThetaZ()*TMath::RadToDeg(), Rot().PhiZ()*TMath::RadToDeg());
void AliRICHChamber::Print(Option_t *opt) const
{
// Debug printout
-// printf("%s r=%8.3f theta=%5.1f phi=%5.1f x=%8.3f y=%8.3f z=%8.3f %6.2f,%6.2f %6.2f,%6.2f %6.2f,%6.2f\n",fName.Data(),
-// Rho(), ThetaD(),PhiD(), X(), Y(), Z(),
-// ThetaXd(),PhiXd(),ThetaYd(),PhiYd(),ThetaZd(),PhiZd());
- fCenterV3.Print(opt);fPcX3.Print(opt);
+ ToAliInfo(fCenterX3.Print(opt));
}//Print()
//__________________________________________________________________________________________________
AliRICHChamber(const AliRICHChamber &chamber):TNamed(chamber) {;}
virtual ~AliRICHChamber() {;}
AliRICHChamber& operator=(const AliRICHChamber&) {return *this;}
-
+
+ static Double_t AlphaFeedback(Int_t ) {return 0.030;} //determines number of feedback photons
+
TRotMatrix* RotMatrix() const{return fpRotMatrix;}
- TString RotMatrixName() const{return "rot"+fName;}
+ TString RotMatrixName() const{return "rot"+fName;}
TRotation Rot() const{return fRot;}
- Double_t Rho() const{return fCenterV3.Mag();} //gives distance to chamber center in MRS
- Double_t ThetaD() const{return fCenterV3.Theta()*TMath::RadToDeg();} //gives polar angle of chamber center in MRS
- Double_t PhiD() const{return fCenterV3.Phi() *TMath::RadToDeg();} //gives azimuthal angle of chamber center in MRS
+ Double_t Rho() const{return fCenterX3.Mag();} //gives distance to chamber center in MRS
+ Double_t ThetaD() const{return fCenterX3.Theta()*TMath::RadToDeg();} //gives polar angle of chamber center in MRS
+ Double_t PhiD() const{return fCenterX3.Phi() *TMath::RadToDeg();} //gives azimuthal angle of chamber center in MRS
Double_t ThetaXd() const{return fRot.ThetaX() *TMath::RadToDeg();}
Double_t PhiXd() const{return fRot.PhiX() *TMath::RadToDeg();}
Double_t ThetaYd() const{return fRot.ThetaY() *TMath::RadToDeg();}
Double_t PhiYd() const{return fRot.PhiY() *TMath::RadToDeg();}
Double_t ThetaZd() const{return fRot.ThetaZ() *TMath::RadToDeg();}
Double_t PhiZd() const{return fRot.PhiZ() *TMath::RadToDeg();}
- void RotateX(Double_t a) {fRot.RotateX(a);fCenterV3.RotateX(a);fPcX3.RotateX(a);} //rotate chamber around X by "a" degrees
- void RotateY(Double_t a) {fRot.RotateY(a);fCenterV3.RotateY(a);fPcX3.RotateY(a);} //rotate chamber around Y by "a" degrees
- void RotateZ(Double_t a) {fRot.RotateZ(a);fCenterV3.RotateZ(a);fPcX3.RotateZ(a);} //rotate chamber around Z by "a" degrees
- Double_t X() const{return fCenterV3.X();}
- Double_t Y() const{return fCenterV3.Y();}
- Double_t Z() const{return fCenterV3.Z();}
- TVector2 Glob2Loc(TVector3 x3)const{x3-=fPcX3;x3.Transform(fRot.Inverse());return TVector2(x3.Z()+0.5*AliRICHParam::PcSizeX(),-x3.X()+0.5*AliRICHParam::PcSizeY());}//Y and Z are misplaced?????
- TVector3 Loc2Glob(TVector2 x2)const{TVector3 x3(-x2.Y()+0.5*AliRICHParam::PcSizeY(),0,x2.X()-0.5*AliRICHParam::PcSizeX());x3.Transform(fRot); x3+=fPcX3;return x3;}
-
- TVector2 Glob2Loc(TLorentzVector x4) const{return Glob2Loc(x4.Vect());}
-
- void Print(Option_t *sOption)const;//virtual
-
-
- inline void SetToZenith();
- TRotMatrix *GetRotMatrix() const{return fpRotMatrix;}
+ void RotX(Double_t a) {a*=TMath::DegToRad();fRot.RotateX(a);fCenterX3.RotateX(a);fRadX3.RotateX(a);fPcX3.RotateX(a);}//degrees around X
+ void RotY(Double_t a) {a*=TMath::DegToRad();fRot.RotateY(a);fCenterX3.RotateY(a);fRadX3.RotateY(a);fPcX3.RotateY(a);}//degrees around Y
+ void RotZ(Double_t a) {a*=TMath::DegToRad();fRot.RotateZ(a);fCenterX3.RotateZ(a);fRadX3.RotateZ(a);fPcX3.RotateZ(a);}//degrees around Z
+ TVector3 Rad() const{return fRadX3;} //provides center of radiator position in MRS, cm
+ TVector3 Pc() const{return fPcX3;} //provides center of photocathond position in MRS, cm
+ TVector3 Center() const{return fCenterX3;} //provides center of chamber position in MRS, cm
+ void Print(Option_t *sOption)const; //virtual interface from TObject
+//Transformations for photcathode plane
+ TVector2 Mrs2Pc(TVector3 x3)const{x3-=fPcX3;x3.Transform(fRot.Inverse());return TVector2(-x3.X()+0.5*AliRICHParam::PcSizeX(),x3.Y()+0.5*AliRICHParam::PcSizeY());}
+ TVector3 Pc2Mrs(TVector2 x2)const{TVector3 x3(-x2.X()+0.5*AliRICHParam::PcSizeX(),x2.Y()-0.5*AliRICHParam::PcSizeY(),0);x3.Transform(fRot); x3+=fPcX3;return x3;}
+ TVector2 Mrs2Pc(TLorentzVector x4) const{return Mrs2Pc(x4.Vect());}
+//Transformations for radiator plane
+ TVector2 Mrs2Rad(TVector3 x3)const{x3-=fRadX3;x3.Transform(fRot.Inverse());return TVector2(-x3.X()+0.5*AliRICHParam::PcSizeX(),x3.Y()+0.5*AliRICHParam::PcSizeY());}
+ TVector3 Rad2Mrs(TVector2 x2)const{TVector3 x3(-x2.X()+0.5*AliRICHParam::PcSizeX(),x2.Y()-0.5*AliRICHParam::PcSizeY(),0);x3.Transform(fRot); x3+=fRadX3;return x3;}
+ TVector3 PMrs2Loc(TVector3 p3)const{TVector3 ploc=Rot().Invert()*p3;ploc.SetXYZ(-ploc.Px(),ploc.Py(),ploc.Pz()); return ploc;}
protected:
- TVector3 fCenterV3; //chamber center position in MRS (cm)
+ TVector3 fCenterX3; //chamber center position in MRS (cm)
+ TVector3 fRadX3; //radiator entrance center position in MRS (cm)
TVector3 fPcX3; //PC center position in MRS (cm)
TRotation fRot; //chamber rotation in MRS
TRotMatrix *fpRotMatrix; //rotation matrix of the chamber with respect to MRS
- ClassDef(AliRICHChamber,6) //single RICH chamber description
+ ClassDef(AliRICHChamber,7) //single RICH chamber description
};//class AliRICHChamber
-//__________________________________________________________________________________________________
-void AliRICHChamber::SetToZenith()
-{
-//Put the chamber to zenith. Position of PC is shifted in X-Z plane since the origin of chamber local system is in
-//left hand down coner.
- fCenterV3.SetXYZ(0,AliRICHParam::Offset()-AliRICHParam::GapThickness()/2 ,0);
- fPcX3.SetXYZ(0,AliRICHParam::Offset()-AliRICHParam::GapThickness()/2+5.276+0.25,0);
-}
-//__________________________________________________________________________________________________
+
#endif //AliRICHChamber_h
//__________________________________________________________________________________________________
AliRICHClusterFinder::AliRICHClusterFinder(AliRICH *pRICH)
{//main ctor
- fRICH = pRICH;//this goes first as GetDebug depends on it.
- if(GetDebug())Info("main ctor","Start.");
+ fRICH = pRICH;
+ AliDebug(1,"main ctor Start.");
fDigitMap = 0;
fRawCluster.Reset();
fResolvedCluster.Reset();
- if(GetDebug())Info("main ctor","Stop.");
+ AliDebug(1,"main ctor Stop.");
}//main ctor
//__________________________________________________________________________________________________
void AliRICHClusterFinder::Exec()
{
//Main method of cluster finder. Loops on events and chambers, everything else is done in FindClusters()
- if(GetDebug()) Info("Exec","Start.");
+ AliDebug(1,"Exec Start.");
R()->GetLoader() ->LoadDigits();
- R()->GetLoader()->GetRunLoader()->LoadHeader();
- R()->GetLoader()->GetRunLoader()->LoadKinematics();
+// R()->GetLoader()->GetRunLoader()->LoadHeader();
+ R()->GetLoader()->GetRunLoader()->LoadKinematics(); //header is already loaded
for(Int_t iEventN=0;iEventN<gAlice->GetEventsPerRun();iEventN++){//events loop
- if(GetDebug()) Info("Exec","Processing event %i...",iEventN);
+ AliDebug(1,Form("Processing event %i...",iEventN));
R()->GetLoader()->GetRunLoader()->GetEvent(iEventN);
R()->GetLoader()->MakeTree("R"); R()->MakeBranch("R");
R()->ResetDigits(); R()->ResetClusters();
R()->GetLoader()->TreeD()->GetEntry(0);
- for(Int_t iChamber=1;iChamber<=kNCH;iChamber++){//chambers loop
+ for(Int_t iChamber=1;iChamber<=kNchambers;iChamber++){//chambers loop
FindClusters(iChamber);
}//chambers loop
R()->GetLoader()->TreeR()->Fill(); R()->GetLoader()->WriteRecPoints("OVERWRITE");//write out clusters for current event
R()->ResetClusters();
R()->GetLoader() ->UnloadDigits();
R()->GetLoader() ->UnloadRecPoints();
- R()->GetLoader()->GetRunLoader()->UnloadHeader();
+// R()->GetLoader()->GetRunLoader()->UnloadHeader();
R()->GetLoader()->GetRunLoader()->UnloadKinematics();
- if(GetDebug()) Info("Exec","Stop.");
+ AliDebug(1,"Stop.");
}//Exec()
//__________________________________________________________________________________________________
void AliRICHClusterFinder::FindClusters(Int_t iChamber)
{
//Loops on digits for a given chamber, forms raw clusters, then tries to resolve them if requested
Int_t iNdigits=R()->Digits(iChamber)->GetEntriesFast();
- if(GetDebug())Info("FindClusters","Start for chamber %i with %i digits.",iChamber,iNdigits);
+ AliDebug(1,Form("Start for chamber %i with %i digits.",iChamber,iNdigits));
if(iNdigits==0)return;//no digits for a given chamber, nothing to do
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
- if(GetDebug()) {Info("FindClusters","After FormRawCluster:");fRawCluster.Print();}
+ AliDebug(1,"After FormRawCluster:");ToAliDebug(1,fRawCluster.Print());
FindLocalMaxima(); //find number of local maxima and initial center of gravity
- if(GetDebug()) {Info("FindClusters","After FindLocalMaxima:");fRawCluster.Print();}
+ AliDebug(1,"After FindLocalMaxima:");ToAliDebug(1,fRawCluster.Print());
if(AliRICHParam::IsResolveClusters()&&fRawCluster.Size()>1&&fRawCluster.Size()<6){
FitCoG(); //serialization of resolved clusters will happen inside
delete fDigitMap;
- if(GetDebug())Info("FindClusters","Stop.");
+ AliDebug(1,"Stop.");
}//FindClusters()
//__________________________________________________________________________________________________
void AliRICHClusterFinder::FindClusterContribs(AliRICHcluster *pCluster)
{
//Finds cerenkov-feedback-mip mixture for a given cluster
- if(GetDebug()) {Info("FindClusterContribs","Start.");pCluster->Print();}
+ AliDebug(1,"Start.");ToAliDebug(1,pCluster->Print());
TObjArray *pDigits = pCluster->Digits();
if(!pDigits) return; //??????????
}//loop on digits of a given cluster
TMath::Sort(contribs.GetSize(),contribs.GetArray(),pindex);
for(Int_t iDigN=0;iDigN<3*pCluster->Size()-1;iDigN++) {//loop on digits to sort tids
- if(GetDebug()) Info("FindClusterContribs","%4i for digit n. %4i",contribs[pindex[iDigN]],iDigN);
+ AliDebug(1,Form("%4i for digit n. %4i",contribs[pindex[iDigN]],iDigN));
if(contribs[pindex[iDigN]]!=contribs[pindex[iDigN+1]]) {
TParticle* particle = R()->GetLoader()->GetRunLoader()->Stack()->Particle(contribs[pindex[iDigN]]);
Int_t code = particle->GetPdgCode();
Double_t charge = 0;
if(particle->GetPDG()) charge=particle->GetPDG()->Charge();
- if(GetDebug()) Info("FindClusterContribs"," charge of particle %f",charge);
+ AliDebug(1,Form(" charge of particle %f",charge));
if(code==50000050) iNckovs++;
if(code==50000051) iNfeeds++;
Int_t code = particle->GetPdgCode();
Double_t charge = 0;
if(particle->GetPDG()) charge=particle->GetPDG()->Charge();
- if(GetDebug()) Info("FindClusterContribs"," charge of particle %f",charge);
+ AliDebug(1,Form(" charge of particle %f",charge));
if(code==50000050) iNckovs++;
if(code==50000051) iNfeeds++;
if(charge!=0) iNmips++;
pCluster->CFM(iNckovs,iNfeeds,iNmips);
//
delete [] pindex;
- if(GetDebug()) {pCluster->Print();Info("FindClusterContribs","Stop.");}
+ ToAliDebug(1,pCluster->Print());
+ AliDebug(1,"Stop.");
}//FindClusterContribs()
//__________________________________________________________________________________________________
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.
- if(GetDebug()) Info("FormRawCluster","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))->Q()));
fRawCluster.AddDigit((AliRICHdigit*) fDigitMap->GetHit(i,j));//take this pad in cluster
fDigitMap->FlagHit(i,j);//flag this pad as taken
void AliRICHClusterFinder::WriteRawCluster()
{
//Add the current raw cluster to the list of clusters
- if(GetDebug()) Info("WriteRawCluster","Start.");
+ AliDebug(1,"Start.");
FindClusterContribs(&fRawCluster);
R()->AddCluster(fRawCluster);
- if(GetDebug()){fRawCluster.Print(); Info("WriteRawCluster","Stop.");}
+ ToAliDebug(1,fRawCluster.Print()); AliDebug(1,"Stop.");
}//WriteRawCluster()
//__________________________________________________________________________________________________
void AliRICHClusterFinder::WriteResolvedCluster()
{
//Add the current resolved cluster to the list of clusters
- if(GetDebug()) Info("WriteResolvedCluster","Start.");
+ AliDebug(1,"Start.");
FindClusterContribs(&fResolvedCluster);
R()->AddCluster(fResolvedCluster);
- if(GetDebug()){fResolvedCluster.Print(); Info("WriteResolvedCluster","Stop.");}
+ ToAliDebug(1,fResolvedCluster.Print()); AliDebug(1,"Stop.");
}//WriteResolvedCluster()
//__________________________________________________________________________________________________
void AliRICHClusterFinder::FitCoG()
{
//Fits cluster of size by the corresponding number of Mathieson shapes.
//This methode is only invoked in case everything is ok to start deconvolution
- if(GetDebug()) {Info("FitCoG","Start with:"); fRawCluster.Print();}
+ AliDebug(1,"Start with:"); ToAliDebug(1,fRawCluster.Print());
Double_t arglist;
Int_t ierflag = 0;
fResolvedCluster.Fill(&fRawCluster,xCoG[i],yCoG[i],qfracCoG[i],fLocalC[i]);
WriteResolvedCluster();
}
- if(GetDebug()) Info("CoG","Stop.");
+ AliDebug(1,"Stop.");
}//FitCoG()
//__________________________________________________________________________________________________
void RICHMinMathieson(Int_t &npar, Double_t *, Double_t &chi2, Double_t *par, Int_t )
AliRICHhit *pHit = (AliRICHhit*)Hits[i].At(j);
if(pHit->C()==iChamber){
TVector3 hitGlobX3= pHit->OutX3();
- TVector2 hitLocX2 = pRich->C(iChamber)->Glob2Loc(hitGlobX3);
+ TVector2 hitLocX2 = pRich->C(iChamber)->Mrs2Pc(hitGlobX3);
pHitsH2->Fill(hitLocX2.X(),hitLocX2.Y(),200);
}//if
}//hits loop
if(!isClusters) {l.SetTextColor(kBlue) ;l.DrawLatex(0.8,0.01,"No CLUSTERS");}
Display->Update();
Display->Modified();
- getchar();
+ gPad->WaitPrimitive();
//deals with digits
if(isDigits){
pRich->GetLoader()->TreeD()->GetEntry(0);
pDigitsH2->Draw("same");
Display->Update();
Display->Modified();
- getchar();
+ gPad->WaitPrimitive();
}//if(isDigits)
//deals with clusters
if(isClusters){
pClustersH2->Draw("same");
Display->Update();
Display->Modified();
- getchar();
+ gPad->WaitPrimitive();
}//if(isClusters)
}//chambers loop
delete [] Hits;
--- /dev/null
+#include "AliRICHHelix.h"
+#include <AliLog.h>
+
+ClassImp(AliRICHHelix)
+
+void AliRICHHelix::Print(Option_t *opt) const
+{
+// Debug printout
+ fX0.Print(opt);fP0.Print(opt);
+ AliInfo("Point of interest:");
+ fX.Print(opt); fP.Print(opt);
+}//Print()
--- /dev/null
+#ifndef AliRICHHelix_h
+#define AliRICHHelix_h
+
+#include <TObject.h>
+#include <TVector3.h>
+#include "AliRICHParam.h"
+#include "AliRICHChamber.h"
+
+class AliRICHHelix: public TObject
+{
+public:
+ AliRICHHelix() {;}
+ AliRICHHelix(TVector3 x0,TVector3 p0,Int_t q,Double_t b=0.4) {fX0=x0;fP0=p0;fQ=q;fBz=b;} //takes position and momentum at parametrised point
+ virtual ~AliRICHHelix() {;}
+ inline void Propagate(Double_t lenght);
+ Bool_t Intersection(TVector3 plane) {return Intersection(plane,plane.Unit());} // special plane given by point only
+ inline Bool_t Intersection(TVector3 planePoint,TVector3 planeNorm); //intersection with plane given by point and normal vector
+ inline Int_t RichIntersect(AliRICHParam *pParam);
+ TVector3 X()const {return fX;}
+ TVector3 P()const {return fP;}
+ TVector3 X0()const {return fX0;}
+ TVector3 P0()const {return fP0;}
+ TVector2 PosPc() const {return fPosPc;} //returns position of intersection with PC
+ TVector2 PosRad()const {return fPosRad;} //returns position of intersection with radiator
+ TVector3 Ploc()const {return fPloc;} //returns momentum at the position of intersection with radiator
+
+ void Print(Option_t *sOption)const; //virtual interface from TObject
+protected:
+ TVector3 fX0; //helix position in parametrised point, cm in MRS
+ TVector3 fP0; //helix momentum in parametrised point, GeV/c in MRS
+ TVector3 fX; //helix position in point of interest, cm in MRS
+ TVector3 fP; //helix momentum in point of interest, GeV/c in MRS
+ Double_t fLength; //helix length in point of interest
+ Int_t fQ; //sign of track charge (value not provided by current ESD)
+ Double_t fBz; //magnetic field along z value in Tesla under assumption of uniformity
+ TVector2 fPosPc; //position on PC in local system
+ TVector2 fPosRad; //position on radiator in local system
+ TVector3 fPloc; //momentum in local system
+ ClassDef(AliRICHHelix,0) //General helix
+};//class AliRICHHelix
+//__________________________________________________________________________________________________
+void AliRICHHelix::Propagate(Double_t length)
+{
+// Propogates the helix to the position of interest defined by helix length s
+// Assumes uniform magnetic field along z direction.
+ const Double_t c = 0.00299792458;//this value provides that coordinates are in cm momentum in GeV/c
+ Double_t a = -c*fBz*fQ;
+
+ Double_t rho = a/fP0.Mag();
+ fX.SetX( fX0.X()+fP0.X()*TMath::Sin(rho*length)/a-fP0.Y()*(1-TMath::Cos(rho*length))/a );
+ fX.SetY( fX0.Y()+fP0.Y()*TMath::Sin(rho*length)/a+fP0.X()*(1-TMath::Cos(rho*length))/a );
+ fX.SetZ( fX0.Z()+fP0.Z()*length/fP0.Mag() );
+ fP.SetX( fP0.X()*TMath::Cos(rho*length)-fP0.Y()*TMath::Sin(rho*length) );
+ fP.SetY( fP0.Y()*TMath::Cos(rho*length)+fP0.X()*TMath::Sin(rho*length) );
+ fP.SetZ( fP0.Z() );
+ fLength=length;
+}
+//__________________________________________________________________________________________________
+Bool_t AliRICHHelix::Intersection(TVector3 planePoint,TVector3 planeNorm)
+{
+// Finds point of intersection (if exists) of the helix to the plane given by point and normal vector.
+// Returns kTrue if helix intersects the plane, kFALSE otherwise.
+// Stores result in current helix fields fX and fP.
+
+ Double_t s=(planePoint-fX0)*planeNorm,dist=99999,distPrev=dist;
+
+ while(TMath::Abs(dist)>0.0001){
+ Propagate(s); //calculates helix at the distance s from x0 ALONG the helix
+ dist=(fX-planePoint)*planeNorm; //distance between current helix position and plane
+ if(TMath::Abs(dist) > TMath::Abs(distPrev)) { return kFALSE;}
+ distPrev=dist;
+ s-=dist;
+ }
+ return kTRUE;
+}
+//__________________________________________________________________________________________________
+Int_t AliRICHHelix::RichIntersect(AliRICHParam *pParam)
+{
+// Searchs for intersection of this helix with all RICH chambers, returns chamber number or 0 if no intersection
+// On exit fPosRad contain position of intersection in Local System with radiator
+// fPosPc contains the same for photocathode
+ for(Int_t iChamberN=1;iChamberN<=kNchambers;iChamberN++){//chamber loop
+ if(Intersection(pParam->C(iChamberN)->Rad())){//there is intersection with radiator plane
+ fPosRad=pParam->C(iChamberN)->Mrs2Rad(fX);//position on radiator plane
+ if(pParam->IsAccepted(fPosRad)){//intersection within radiator (even if in dead zone)
+ if(Intersection(pParam->C(iChamberN)->Pc())){//there is intersection with photocathode
+ fPosPc=pParam->C(iChamberN)->Mrs2Pc(fX);//position on photcathode plane
+ if(pParam->IsAccepted(fPosPc)){//intersection within pc (even if in dead zone)
+ fPloc=pParam->C(iChamberN)->PMrs2Loc(fP);
+ return iChamberN;
+ }//if inside PC
+ }//if for PC
+ }//if inside radiator
+ }//if for radiator
+ }//chamber loop
+ return 0;
+}
+//__________________________________________________________________________________________________
+#endif
//__________________________________________________________________________________________________
void AliRICHParam::Print(Option_t*)
{
- ::Info("","Pads in chamber (%3i,%3i) in sector (%2i,%2i)",NpadsX(),NpadsY(),NpadsXsec(),NpadsYsec());
- fpChambers->Print();
-}
+ AliInfo(Form("Pads in chamber (%3i,%3i) in sector (%2i,%2i)",NpadsX(),NpadsY(),NpadsXsec(),NpadsYsec()));
+ ToAliInfo(fpChambers->Print());
+}//Print()
//__________________________________________________________________________________________________
void AliRICHParam::CreateChambers()
{
//Create all RICH Chambers on each call. Previous chambers deleted.
if(fpChambers) delete fpChambers;
- fpChambers=new TObjArray(kNchambers);
+ if(IsRadioSrc()){
+ fpChambers=new TObjArray(1);//test beam configuration 1 chamber
+ fpChambers->AddAt(new AliRICHChamber(0),0);
+ }else{
+ fpChambers=new TObjArray(kNchambers);//normal configuration 7 chambers
+ for(int iChamberN=0;iChamberN<kNchambers;iChamberN++) fpChambers->AddAt(new AliRICHChamber(iChamberN+1),iChamberN);
+ }
fpChambers->SetOwner();
- for(int iChamberN=0;iChamberN<kNchambers;iChamberN++) fpChambers->AddAt(new AliRICHChamber(iChamberN+1),iChamberN);
-}//void AliRICH::CreateChambers()
+}//CreateChambers()
#include <TMath.h>
#include <TObjArray.h>
#include <TObject.h>
+#include <TMath.h>
#include <TRandom.h>
#include <TVector.h>
#include <TVector2.h>
#include <TVector3.h>
+#include <TRandom.h>
+#include <TError.h>
+#include <TObjArray.h>
+#include <AliLog.h>
+#include <TClass.h>
+
-static const int kNCH=7; //number of RICH chambers ???
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
static const int kBad=-101; //useful static const to mark initial (uninitalised) values
static const int kNsectors=6; //number of sectors per chamber
-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
AliRICHParam():TObject(),fpChambers(0) {CreateChambers();}
virtual ~AliRICHParam() {delete fpChambers;}
void CreateChambers();
- AliRICHChamber* C(Int_t i) {return (AliRICHChamber*)fpChambers->UncheckedAt(i-1);} //returns pointer to chamber i
+ AliRICHChamber* C(Int_t i) {return (AliRICHChamber*)fpChambers->UncheckedAt(i-1);} //returns pointer to chamber i
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 Double_t DeadZone() {return 2.6;} //dead zone size in cm
static Double_t PadSizeX() {return 0.8;} //pad size x in cm
static Double_t PadSizeY() {return 0.84;} //pad size y in cm
+
static Double_t SectorSizeX() {return NpadsX()*PadSizeX()/2;} //sector size x in cm
static Double_t SectorSizeY() {return NpadsY()*PadSizeY()/3;} //sector size y in cm
- static Double_t PcSizeX() {return NpadsX()*PadSizeX()+DeadZone();} //photocathode size x in cm
- static Double_t PcSizeY() {return NpadsY()*PadSizeY()+2*DeadZone();} //photocathode size y in cm
- static Double_t SizeX() {return 132.6;}
- static Double_t SizeY() {return 26;}
- static Double_t SizeZ() {return 136.7;}
- static Double_t Offset() {return 490+1.267;} //distance from IP to center of chamber in cm
- static Double_t AngleYZ() {return 19.5*TMath::DegToRad();} //angle between chambers in YZ plane, rad
- static Double_t AngleXY() {return 20*TMath::DegToRad();} //angle between chambers in XY plane, rad
- static Double_t AngleRot() {return fgAngleRot*TMath::DegToRad();} //RICH rotation around Z, rad
- static Double_t FreonThickness() {return 1.5;}
- static Double_t QuartzThickness() {return 0.5;}
-
- static Double_t GapProx() {return 8.0;} //cm between CsI PC and radiator quartz window
- static Double_t GapColl() {return 7.0;} //cm between CsI PC and third wire grid (collection wires)
- static Double_t GapAnod() {return 0.204;} //cm between CsI PC and first wire grid (anod wires)
- static Double_t GapAmp() {return 0.445;} //cm between CsI PC and second wire grid (cathode wires)
+ 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 collect wires
+ static Double_t PitchColl() {return 0.5;} //cm between collection wires
- static Double_t GapThickness() {return 8.0;}
- static Double_t RadiatorToPads() {return FreonThickness()+QuartzThickness()+GapThickness();}
- static Double_t AnodeCathodeGap() {return 0.2;} //between CsI PC and first wire grid
- static Double_t QuartzLength() {return 133;}
- static Double_t QuartzWidth() {return 127.9;}
- static Double_t OuterFreonLength() {return 133;}
- static Double_t OuterFreonWidth() {return 41.3;}
- static Double_t InnerFreonLength() {return 133;}
- static Double_t InnerFreonWidth() {return 41.3;}
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 MaxQdc() {return 4095;} //QDC number of channels
- static Double_t AlphaFeedback(Int_t ) {return 0.030;} //determines number of feedback photons
static Bool_t IsResolveClusters() {return fgIsResolveClusters;} //go after resolved clusters?
static Bool_t IsWireSag() {return fgIsWireSag;} //take wire sagita in account?
static void SetWireSag(Bool_t status) {fgIsWireSag=status;}
static void SetHV(Int_t sector,Int_t hv){fgHV[sector-1]=hv;}
static void SetAngleRot(Double_t rot) {fgAngleRot =rot;}
+ static Double_t IndOfRefC6F14(Double_t eV) {return eV*0.0172+1.177;} // eV = photon energy in eV
+ static Double_t IndOfRefSiO2(Double_t eV) {Double_t e1=10.666,e2=18.125,f1=46.411,f2= 228.71;
+ return TMath::Sqrt(1.+f1/(e1*e1-TMath::Power(eV,2))+f2/(e2*e2-TMath::Power(eV,2)));}//TDR p.35
+ static Double_t IndOfRefCH4() {return 1.000444;}
inline static TVector Loc2Area(TVector2 x2); //return area affected by hit x2
inline static TVector Loc2Pad(TVector2 x2); //return pad containing given position
static Float_t SigmaThMean() {return fgSigmaThMean;} //QDC electronic noise mean
static Float_t SigmaThSpread() {return fgSigmaThSpread;} //QDC electronic noise width
void Print(const Option_t *opt=""); //virtual
- inline static void PropogateHelix(TVector3 x0,TVector3 p0,Double_t s,TVector3 *x,TVector3 *p);
inline static Int_t Loc2Sec(TVector2 &x2); //return sector, x2->Sector RS
- inline static Int_t Pad2Sec(TVector pad); //return sector
+ inline static Int_t Pad2Sec(const TVector &pad); //return sector
+ static Bool_t IsAccepted(const TVector2 &x2) {return ( x2.X()>=0 && x2.X()<=PcSizeX() && x2.Y()>=0 && x2.Y()<=PcSizeY() )
+? kTRUE:kFALSE;}
protected:
TObjArray *fpChambers; //list of chambers
static Bool_t fgIsWireSag; //wire sagitta ON/OFF flag
Double_t x=v2.X(),y=v2.Y();
if (v2.X() >= x0 && v2.X() <= x1 ) {sector=1;}
else if(v2.X() >= x2 && v2.X() <= x3 ) {sector=2; x=v2.X()-x2;}
- else {sector=kBad; ::Error("Loc2Sec","Position %6.2f,%6.2f is out of chamber in X",v2.X(),v2.Y());return kBad;}
+ else {return kBad;}
if (v2.Y() >= y0 && v2.Y() <= y1 ) {} //sectors 1 or 2
else if(v2.Y() >= y2 && v2.Y() <= y3 ) {sector+=2; y=v2.Y()-y2;} //sectors 3 or 4
else if(v2.Y() >= y4 && v2.Y() <= y5 ) {sector+=4; y=v2.Y()-y4;} //sectors 5 or 6
- else {sector=kBad; ::Error("Loc2Sec","Position %6.2f,%6.2f is out of chamber in Y",v2.X(),v2.Y());return kBad;}
+ else {return kBad;}
v2.Set(x,y);
return sector;
}//Loc2Sec(Double_t x, Double_t y)
return pad;
}
//__________________________________________________________________________________________________
-Int_t AliRICHParam::Pad2Sec(TVector pad)
+Int_t AliRICHParam::Pad2Sec(const TVector &pad)
{
// Determines sector containing the given pad.
Int_t sector=kBad;
if (pad[0] >= 1 && pad[0] <= NpadsXsec() ) {sector=1;}
else if(pad[0] > NpadsXsec() && pad[0] <= NpadsX() ) {sector=2;}
- else ::Error("Pad2Sec","Wrong pad (%3.0f,%3.0f)",pad[0],pad[1]);
+ else AliDebugClass(1,Form("Wrong pad (%3.0f,%3.0f)",pad[0],pad[1]));
if (pad[1] >= 1 && pad[1] <= NpadsYsec() ) {}
else if(pad[1] > NpadsYsec() && pad[1] <= 2*NpadsYsec() ) {sector+=2;}
else if(pad[1] > 2*NpadsYsec() && pad[1] <= NpadsY() ) {sector+=4;}
- else ::Error("Pad2Sec","Wrong pad (%3.0f,%3.0f)",pad[0],pad[1]);
+ else AliDebugClass(1,Form("Wrong pad (%3.0f,%3.0f)",pad[0],pad[1]));
return sector;
}//Pad2Sec()
//__________________________________________________________________________________________________
TVector2 AliRICHParam::Pad2Loc(TVector pad)
{
-// Returns position of the center of the given pad in local system of the chamber
+// Returns position of the center of the given pad in local system of the chamber (cm)
// y ^ 5 6
-// | 3 4 chamber structure
+// | 3 4 sector numbers
// | 1 2
// -------> x
Double_t x=kBad,y=kBad;
else if(pad[0] > NpadsXsec() && pad[0] <= NpadsX())//it's 2 or 4 or 6
x=(pad[0]-0.5)*PadSizeX()+DeadZone();
else
- ::Error("Pad2Loc","Wrong pad (%3.0f,%3.0f)",pad[0],pad[1]);
+ AliDebugClass(1,Form("Wrong pad (%3.0f,%3.0f)",pad[0],pad[1]));
if(pad[1] > 0 && pad[1] <= NpadsYsec())//it's 1 or 2
y=(pad[1]-0.5)*PadSizeY();
else if(pad[1] > 2*NpadsYsec() && pad[1]<= NpadsY())//it's 5 or 6
y=(pad[1]-0.5)*PadSizeY()+2*DeadZone();
else
- ::Error("Pad2Loc","Wrong pad (%3.0f,%3.0f)",pad[0],pad[1]);
+ AliDebugClass(1,Form("Wrong pad (%3.0f,%3.0f)",pad[0],pad[1]));
return TVector2(x,y);
}
Double_t AliRICHParam::GainSag(Double_t x,Int_t sector)
{
// Returns % of gain variation due to wire sagita.
-// All curves are parametrized as per sector basis, so y must be apriory transformed to the Sector RS.
+// All curves are parametrized as per sector basis, so x must be apriory transformed to the Sector RS.
+// Here x is a distance along wires.
x-=SectorSizeX()/2;
if(x>SectorSizeX()) x-=SectorSizeX();
switch(HV(sector)){
{
// Calculates the total charge produced by the eloss in point x2 (Chamber RS).
// Returns this change parametrised in QDC channels, or 0 if the hit in the dead zone.
-// eloss=0 means photons which provided for only 1 electron
-// eloss > 0 for Mip
+// eloss=0 means photon which produces 1 electron only eloss > 0 for Mip
if(Loc2Sec(x2)==kBad) return 0; //hit in the dead zone
Int_t iNelectrons=Int_t(eloss/IonisationPotential()); if(iNelectrons==0) iNelectrons=1;
Double_t qdc=0;
// Calculates the charge fraction induced to given pad by the hit from the given point.
// Integrated Mathieson distribution is used.
TVector2 center2=Pad2Loc(pad);//gives center of requested pad
- Double_t normXmin=(x2.X()-center2.X()-PadSizeX()/2) /AnodeCathodeGap();
- Double_t normXmax=(x2.X()-center2.X()+PadSizeX()/2) /AnodeCathodeGap();
- Double_t normYmin=(x2.Y()-center2.Y()-PadSizeY()/2) /AnodeCathodeGap();
- Double_t normYmax=(x2.Y()-center2.Y()+PadSizeY()/2) /AnodeCathodeGap();
-
- if(Loc2Sec(x2)!=Pad2Sec(pad)) return 0;//requested pad does not belong to the sector of the given hit position
- else return Mathieson(normXmin, normYmin, normXmax, normYmax);
+ Double_t normXmin=(x2.X()-center2.X()-PadSizeX()/2) /Pc2Cath();//parametrise for Mathienson
+ Double_t normXmax=(x2.X()-center2.X()+PadSizeX()/2) /Pc2Cath();
+ Double_t normYmin=(x2.Y()-center2.Y()-PadSizeY()/2) /Pc2Cath();
+ Double_t normYmax=(x2.Y()-center2.Y()+PadSizeY()/2) /Pc2Cath();
+
+//requested pad might not belong to the sector of the given hit position, hence the check:
+ 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)
// 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()));
}
-//__________________________________________________________________________________________________
-void AliRICHParam::PropogateHelix(TVector3 x0,TVector3 p0,Double_t s,TVector3 *x,TVector3 *p)
-{
-// Propogates the helix given by (x0,p0) in MRS to the position of interest defined by helix length s
- const Double_t c = 0.00299792458;
- const Double_t Bz = 0.5; //field in Tesla
- const Double_t q = 1; //charge in electron units
- Double_t a = -c*Bz*q;
-
- Double_t rho = a/p0.Mag();
- p->SetX(p0.X()*TMath::Cos(rho*s)-p0.Y()*TMath::Sin(rho*s));
- p->SetY(p0.Y()*TMath::Cos(rho*s)+p0.X()*TMath::Sin(rho*s));
- p->SetZ(p0.Z());
- x->SetX(x0.X()+p0.X()*TMath::Sin(rho*s)/a-p0.Y()*(1-TMath::Cos(rho*s))/a);
- x->SetY(x0.Y()+p0.Y()*TMath::Sin(rho*s)/a+p0.X()*(1-TMath::Cos(rho*s))/a);
- x->SetZ(x0.Z()+p0.Z()*s/p->Mag());
-}
#endif //AliRICHParam_h
//////////////////////////////////////////////////////////////////////////
#include <Riostream.h>
-#include <TCanvas.h>
-#include <TFile.h>
-#include <TH2.h>
#include <TMath.h>
-#include <TMath.h>
-#include <TMinuit.h>
-#include <TNtuple.h>
-#include <TParticle.h>
-#include <TRandom.h>
#include <TRotation.h>
#include <TVector3.h>
-#include "AliLoader.h"
#include "AliRICH.h"
-#include "AliRICHChamber.h"
#include "AliRICHParam.h"
#include "AliRICHRecon.h"
-#include "AliRun.h"
-#include "AliStack.h"
+#include "AliRICHHelix.h"
+#include <AliLog.h>
#define NPointsOfRing 201
-TMinuit *gAliRICHminuit ;
-
-void fcnrecon(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag);
//__________________________________________________________________________________________________
-AliRICHRecon::AliRICHRecon(const char*name, const char*title)
- :TTask(name,title)
+AliRICHRecon::AliRICHRecon(AliRICHHelix *pHelix,TClonesArray *pClusters,Int_t iMipId)
+ :TTask("RichRec","RichPat")
{
- // main ctor
+// main ctor
+ SetFreonScaleFactor(1);
+ fIsWEIGHT = kFALSE;
fThetaBin=750; fThetaMin = 0.0; fThetaMax = 0.75;
- fDTheta = 0.001; fWindowWidth = 0.060;
- fNpadX = AliRICHParam::NpadsY();
- fNpadY = AliRICHParam::NpadsX();
- fPadSizeX = AliRICHParam::PadSizeY();
- fPadSizeY = AliRICHParam::PadSizeX();
- fRadiatorWidth = AliRICHParam::FreonThickness();
- fQuartzWidth = AliRICHParam::QuartzThickness();
- fGapWidth = AliRICHParam::RadiatorToPads();
- fXmin = -AliRICHParam::PcSizeY()/2.;
- fXmax = AliRICHParam::PcSizeY()/2.;
- fYmin = -AliRICHParam::PcSizeX()/2.;
- fYmax = AliRICHParam::PcSizeX()/2.;
- fRich = (AliRICH*)gAlice->GetDetector("RICH");
- fOutFile=new TFile("Anal.root","RECREATE","My Analysis histos");
- if(fIsDISPLAY) fDisplay = new TCanvas("Display","RICH Display",0,0,1200,750);
- fNtuple=new TNtuple("hn","ntuple",
-"Run:Trig:VertZ:Pmod:Pt:Eta:TrackTheta:TrackPhi:TrackThetaFit:TrackPhiFit:Charge::NPhotons:NPhotonsFit:InRing:MassOfParticle:HoughArea:Multiplicity:TPCLastZ");
+ fDTheta = 0.001; fWindowWidth = 0.060;
+ fRadiatorWidth = AliRICHParam::Zfreon();
+ fQuartzWidth = AliRICHParam::Zwin();
+ fGapWidth = AliRICHParam::Freon2Pc() - fRadiatorWidth - fQuartzWidth;
+ fXmin = -AliRICHParam::PcSizeX()/2.;
+ fXmax = AliRICHParam::PcSizeX()/2.;
+ fYmin = -AliRICHParam::PcSizeY()/2.;
+ fYmax = AliRICHParam::PcSizeY()/2.;
+ SetTrackTheta(pHelix->Ploc().Theta());
+ SetTrackPhi(pHelix->Ploc().Phi());
+ SetMipIndex(iMipId);
+ SetShiftX(pHelix->PosRad().X());
+ SetShiftY(pHelix->PosRad().Y());
+ fpClusters = pClusters;
}
//__________________________________________________________________________________________________
-void AliRICHRecon::StartProcessEvent()
+Double_t AliRICHRecon::ThetaCerenkov()
{
- //start to process for pattern recognition
-
- Float_t trackThetaStored = 0;
- Float_t trackPhiStored = 0;
- Float_t thetaCerenkovStored = 0;
- Int_t houghPhotonsStored = 0;
-
- SetFreonScaleFactor(0.994);
-
- if(fIsDISPLAY)
- {
- DrawEvent(0);
-// Waiting();
- }
-
- AliLoader * richLoader = Rich()->GetLoader();
- AliRunLoader * runLoader = richLoader->GetRunLoader();
-
- if (richLoader->TreeH() == 0x0) richLoader->LoadHits();
- if (richLoader->TreeR() == 0x0) richLoader->LoadRecPoints();
- if (richLoader->TreeD() == 0x0) richLoader->LoadDigits();
- if (runLoader->TreeE() == 0x0) runLoader->LoadHeader();
- if (runLoader->TreeK() == 0x0) runLoader->LoadKinematics();
-
- richLoader->TreeR()->GetEntry(0);
-
- Float_t clusX[7][500],clusY[7][500];
- Int_t clusQ[7][500],clusMul[7][500];
- Int_t nClusters[7];
-
- for (Int_t ich=0;ich<7;ich++) {
- nClusters[ich] = Rich()->Clusters(ich+1)->GetEntries();
- for(Int_t k=0;k<nClusters[ich];k++) {
- AliRICHcluster *pCluster = (AliRICHcluster *)Rich()->Clusters(ich+1)->At(k);
- clusX[ich][k] = pCluster->X();
- clusY[ich][k] = pCluster->Y();
- clusQ[ich][k] = pCluster->Q();
- clusMul[ich][k] = pCluster->Size();
- // pCluster->Print();
- }
- }
-
- Int_t nPrimaries = (Int_t)richLoader->TreeH()->GetEntries();
-
- cout << " N. primaries " << nPrimaries << endl;
-
- for(Int_t i=0;i<nPrimaries;i++){
-
- richLoader->TreeH()->GetEntry(i);
-
- // Rich()->Hits()->Print();
- Int_t iPrim = 0;
-
- AliRICHhit* pHit=0;
-
- for(Int_t j=0;j<Rich()->Hits()->GetEntries();j++) {
-
- pHit = (AliRICHhit*)Rich()->Hits()->At(j);
- if(pHit->GetTrack() < nPrimaries) break;
- iPrim++;
- }
-
- cout << " iPrim " << iPrim << " pHit " << pHit << endl;
-
- if (!pHit) return;
-
- // pHit->Print();
-
- TParticle *pParticle = runLoader->Stack()->Particle(pHit->GetTrack());
- Float_t pmod = pParticle->P();
- Float_t pt = pParticle->Pt();
- Float_t trackEta = pParticle->Eta();
- Int_t q = (Int_t)TMath::Sign(1.,pParticle->GetPDG()->Charge());
-
- // pParticle->Print();
-
- cout << " pmod " << pmod << " pt " << pt << " Eta " << trackEta << " charge " << q << endl;
-
- SetTrackMomentum(pmod);
- SetTrackPt(pt);
- SetTrackEta(trackEta);
- SetTrackCharge(q);
-
- TVector3 pLocal(0,0,0);//?????
-
- TVector2 primLocal =Rich()->C(pHit->C())->Glob2Loc(pHit->InX3());
-
- // Float_t pmodFreo = pLocal.Mag();
- Float_t trackTheta = pLocal.Theta();
- Float_t trackPhi = pLocal.Phi();
-
- // cout << " trackTheta " << trackTheta << " trackPhi " << trackPhi << endl;
-
- SetTrackTheta(trackTheta);
- SetTrackPhi(trackPhi);
-
- Int_t maxInd = 0;
- Float_t minDist = 999.;
-
- // cout << " n Clusters " << nClusters[pHit->Chamber()-1] << " for chamber n. " << pHit->Chamber() << endl;
-
- for(Int_t j=0;j<nClusters[pHit->Chamber()-1];j++)
- {
- Float_t diffx = primLocal.X() - clusX[pHit->Chamber()-1][j];
- Float_t diffy = primLocal.Y() - clusY[pHit->Chamber()-1][j];
-
-
- Float_t diff = sqrt(diffx*diffx + diffy*diffy);
-
- if(diff < minDist)
- {
- minDist = diff;
- maxInd = j;
- }
-
- }
-
- Float_t diffx = primLocal.X() - clusX[pHit->Chamber()-1][maxInd];
- Float_t diffy = primLocal.Y() - clusY[pHit->Chamber()-1][maxInd];
-
- cout << " diffx " << diffx << " diffy " << diffy << endl;
-
-
- SetMipIndex(maxInd);
- SetTrackIndex(i);
-
- Float_t shiftX = 0;//primLocal.X()/primLocal.Z()*(fRadiatorWidth+fQuartzWidth+fGapWidth) + primLocal.X(); ?????
- Float_t shiftY = 0;//primLocal.Y()/primLocal.Z()*(fRadiatorWidth+fQuartzWidth+fGapWidth) + primLocal.Y(); ?????
-
- SetShiftX(shiftX);
- SetShiftY(shiftY);
-
- Float_t *pclusX = &clusX[pHit->Chamber()-1][0];
- Float_t *pclusY = &clusY[pHit->Chamber()-1][0];
-
- SetCandidatePhotonX(pclusX);
- SetCandidatePhotonY(pclusY);
- SetCandidatePhotonsNumber(nClusters[pHit->Chamber()-1]);
-
- Int_t qch = clusQ[pHit->Chamber()-1][maxInd];
-
-
- if(minDist < 3.0 && qch > 120 && maxInd !=0)
- {
-
- if(fIsBACKGROUND)
- {
-
- Float_t xrndm = fXmin + (fXmax-fXmin)*gRandom->Rndm(280964);
- Float_t yrndm = fYmin + (fYmax-fYmin)*gRandom->Rndm(280964);
- SetShiftX(xrndm);
- SetShiftY(yrndm);
-
- }
-
- PatRec();
-
- trackThetaStored = GetTrackTheta();
- trackPhiStored = GetTrackPhi();
- thetaCerenkovStored = GetThetaCerenkov();
- houghPhotonsStored = GetHoughPhotons();
-
- Int_t diffNPhotons = 999;
- Int_t nsteps = 0;
- Float_t diffTrackTheta = 999.;
- Float_t diffTrackPhi = 999.;
-
- while(fIsMINIMIZER && GetHoughPhotons() > 2
- && diffNPhotons !=0
- && diffTrackTheta > 0.0001
- && nsteps < 10)
- {
-
- Int_t houghPhotonsBefore = GetHoughPhotons();
-
- Float_t trackThetaBefore = GetTrackTheta();
- Float_t trackPhiBefore = GetTrackPhi();
-
- Minimization();
-
- PatRec();
-
- diffNPhotons = TMath::Abs(houghPhotonsBefore - GetHoughPhotons());
-
- Float_t trackThetaAfter = GetTrackTheta();
- Float_t trackPhiAfter = GetTrackPhi();
-
- diffTrackTheta = TMath::Abs(trackThetaAfter - trackThetaBefore);
- diffTrackPhi = TMath::Abs(trackPhiAfter - trackPhiBefore);
-
- if(fDebug)
- cout << " houghPhotonsBefore " << houghPhotonsBefore
- << " GetHoughPhotons() " << GetHoughPhotons();
-
- nsteps++;
- }
-
- SetFittedThetaCerenkov(GetThetaCerenkov());
- SetFittedHoughPhotons(GetHoughPhotons());
-
- SetTrackTheta(trackThetaStored);
- SetTrackPhi(trackPhiStored);
- SetThetaCerenkov(thetaCerenkovStored);
- SetHoughPhotons(houghPhotonsStored);
-
- SetMinDist(minDist);
-
- FillHistograms();
-
- if(fIsDISPLAY) DrawEvent(1);
-
- Waiting();
-
- }
- }
- if(fIsDISPLAY) fDisplay->Print("display.ps");
-}//StartProcessEvent()
-//__________________________________________________________________________________________________
-void AliRICHRecon::EndProcessEvent()
-{
- // function called at the end of the event loop
-
- fOutFile->Write();
- fOutFile->Close();
-}
-//__________________________________________________________________________________________________
-void AliRICHRecon::PatRec()
-{
- //pattern recognition method based on Hough transform
-
-
- Float_t trackTheta = GetTrackTheta();
- Float_t trackPhi = GetTrackPhi();
- Float_t pmod = GetTrackMomentum();
- Int_t iMipIndex = GetMipIndex();
+// Pattern recognition method based on Hough transform
+// Return theta Cerenkov for a given track and list of clusters which are set in ctor
+ if(fpClusters->GetEntries()==0) return kBad;//no clusters at all for a given track
Bool_t kPatRec = kFALSE;
Int_t candidatePhotons = 0;
- Float_t shiftX = GetShiftX();
- Float_t shiftY = GetShiftY();
-
- Float_t* candidatePhotonX = GetCandidatePhotonX();
- Float_t* candidatePhotonY = GetCandidatePhotonY();
-
- Int_t candidatePhotonsNumber = GetCandidatePhotonsNumber();
-
- if(fDebug) cout << " n " << candidatePhotonsNumber << endl;
-
SetThetaCerenkov(999.);
SetHoughPhotons(0);
SetHoughPhotonsNorm(0);
- SetHoughRMS(999.);
-
- for (Int_t j=0; j < candidatePhotonsNumber; j++)
- {
-
- SetPhotonIndex(j);
-
- SetPhotonFlag(0);
- SetPhotonEta(-999.);
- SetPhotonWeight(0.);
-
- if (j == iMipIndex) continue;
-
-
- if(candidatePhotonX[j] < -64.) continue; /* avoid artificial clusters from edge uesd by Yale.... */
-
- Float_t xtoentr = candidatePhotonX[j] - shiftX;
- Float_t ytoentr = candidatePhotonY[j] - shiftY;
-
- // Float_t chargehit = fHits_charge[j];
- // if(chargehit > 150) continue;
-
- SetEntranceX(xtoentr);
- SetEntranceY(ytoentr);
-
- FindPhiPoint();
-
- Int_t photonStatus = PhotonInBand();
-
- if(fDebug)
- {
- cout << " Photon n. " << j << " Status " << photonStatus << " accepted " << endl;
- cout << " CandidatePhotonX[j] " << candidatePhotonX[j] << " CandidatePhotonY[j] " << candidatePhotonY[j] << endl;
- }
-
- if(photonStatus == 0) continue;
-
- SetPhotonFlag(1);
-
- FindThetaPhotonCerenkov();
-
- Float_t thetaPhotonCerenkov = GetThetaPhotonCerenkov();
-
- if(fDebug) cout << " theta photon " << thetaPhotonCerenkov << endl;
- SetPhotonEta(thetaPhotonCerenkov);
-
- candidatePhotons++;
-
-
- }
+ for (Int_t j=0; j < fpClusters->GetEntries(); j++){//clusters loop
+ SetPhotonIndex(j);
+ SetPhotonFlag(0);
+ SetPhotonEta(-999.);
+ SetPhotonWeight(0.);
+ if (j == GetMipIndex()) continue; // do not consider MIP cluster as a candidate photon
+ Float_t xtoentr = ((AliRICHcluster*)fpClusters->UncheckedAt(j))->X() - GetShiftX();
+ Float_t ytoentr = ((AliRICHcluster*)fpClusters->UncheckedAt(j))->Y() - GetShiftY();
+ SetEntranceX(xtoentr);
+ SetEntranceY(ytoentr);
+ FindPhiPoint();
+// Int_t photonStatus = PhotonInBand();
+// if(photonStatus == 0) continue;
+ SetPhotonFlag(1);
+ FindThetaPhotonCerenkov();
+ Float_t thetaPhotonCerenkov = GetThetaPhotonCerenkov();
+ AliDebug(1,Form("THETA CERENKOV ---> %i",thetaPhotonCerenkov));
+ SetPhotonEta(thetaPhotonCerenkov);
+ candidatePhotons++;
+ }//clusters loop
if(candidatePhotons >= 1) kPatRec = kTRUE;
- if(!kPatRec) return;
- {
- SetThetaCerenkov(999.);
- SetHoughPhotons(0);
- }
- SetPhotonsNumber(candidatePhotonsNumber);
+ if(!kPatRec) return kBad;
+
+ SetPhotonsNumber(fpClusters->GetEntries());
HoughResponse();
{
SetThetaCerenkov(999.);
SetHoughPhotonsNorm(0.);
- return;
+ return kBad;
}
if(fIsWEIGHT) FindWeightThetaCerenkov();
Float_t thetaCerenkov = GetThetaCerenkov();
+ AliDebug(1,Form("Number of clusters accepted ---> %i",nPhotonHough));
+
SetThetaOfRing(thetaCerenkov);
- FindAreaAndPortionOfRing();
+// FindAreaAndPortionOfRing();
- Float_t nPhotonHoughNorm = ((Float_t)nPhotonHough)/GetPortionOfRing();
- SetHoughPhotonsNorm(nPhotonHoughNorm);
+// Float_t nPhotonHoughNorm = ((Float_t)nPhotonHough)/GetPortionOfRing();
+// SetHoughPhotonsNorm(nPhotonHoughNorm);
// Calculate the area where the photon are accepted...
-
+/*
Float_t thetaInternal = thetaCerenkov - 0.5*fWindowWidth;
SetThetaOfRing(thetaInternal);
FindAreaAndPortionOfRing();
Float_t houghArea = externalArea - internalArea;
SetHoughArea(houghArea);
+*/
+ return GetThetaCerenkov();
- if(fDebug)
- {
- cout << " ----- SUMMARY OF RECONSTRUCTION ----- " << endl;
- cout << " Rings found " << fNrings << " with thetac " << thetaCerenkov << endl;
-
-
- cout << " Nphotons " << GetPhotonsNumber()
- << " Hough " << nPhotonHough
- << " norm " << nPhotonHoughNorm << endl;
-
- cout << " In PatRec:p " << pmod << " theta " << trackTheta << " phi " << trackPhi << endl;
- cout << " ------------------------------------- " << endl;
- }
-
- Int_t nPhotons = GetPhotonsNumber();
-
- Float_t xmean = 0.;
- Float_t x2mean = 0.;
- Int_t nev = 0;
-
- for (Int_t j=0; j < nPhotons;j++)
- {
- SetPhotonIndex(j);
-
- Float_t eta = GetPhotonEta();
-
- if(eta != -999.)
- {
- if(GetPhotonFlag() == 2)
- {
-
-
- xmean += eta;
- x2mean += eta*eta;
- nev++;
- }
- }
- }
-
- if(nev > 0)
- {
- xmean /=(Float_t)nev;
- x2mean /=(Float_t)nev;
- } else {
- xmean = 0.;
- x2mean = 0.;
- }
-
- Float_t vRMS = sqrt(x2mean - xmean*xmean);
-
- SetHoughRMS(vRMS);
-
- if(fDebug) cout << " RMS " << vRMS << endl;
-
-}
-
+}//ThetaCerenkov()
+//__________________________________________________________________________________________________
void AliRICHRecon::FindEmissionPoint()
{
//estimate the emission point in radiator
Float_t emissionPoint = fRadiatorWidth + photonLenghtMin - photonLenghtMax;
SetEmissionPoint(emissionPoint);
+ SetEmissionPoint(fRadiatorWidth/2); // tune the emission point
}
//search band fro photon candidates
// Float_t massOfParticle;
- Float_t beta;
Float_t nfreon;
Float_t thetacer;
Float_t phpad = GetPhiPoint();
- // Float_t pmod = GetTrackMomentum();
- // Float_t trackTheta = GetTrackTheta();
- // Float_t trackPhi = GetTrackPhi();
// inner radius //
SetPhotonEnergy(5.6);
SetEmissionPoint(fRadiatorWidth -0.0001);
- SetMassHypotesis(0.93828);
-
- SetBetaOfParticle();
SetFreonRefractiveIndex();
- beta = GetBetaOfParticle();
nfreon = GetFreonRefractiveIndex();
-
- thetacer = Cerenkovangle(nfreon,beta);
-
thetacer = 0.;
- if(fDebug) cout << " thetacer in photoninband min " << thetacer << endl;
+ AliDebug(1,Form("thetacer in photoninband min %f",thetacer));
FindThetaAtQuartz(thetacer);
SetPhotonEnergy(7.7);
SetEmissionPoint(0.);
// SetMassHypotesis(0.139567);
- SetMassHypotesis(0.);
-
- SetBetaOfParticle();
SetFreonRefractiveIndex();
- beta = GetBetaOfParticle();
nfreon = GetFreonRefractiveIndex();
- thetacer = Cerenkovangle(nfreon,beta);
+ thetacer = Cerenkovangle(nfreon,1);
// thetacer = 0.75;
- if(fDebug) cout << " thetacer in photoninband max " << thetacer << endl;
+ AliDebug(1,Form("thetacer in photoninband max %f",thetacer));
FindThetaAtQuartz(thetacer);
Float_t padradius = sqrt(TMath::Power(xtoentr,2)+TMath::Power(ytoentr,2));
- if(fDebug) printf(" rmin %f r %f rmax %f \n",innerRadius,padradius,outerRadius);
+ AliDebug(1,Form("rmin %f r %f rmax %f",innerRadius,padradius,outerRadius));
if(padradius>=innerRadius && padradius<=outerRadius) return 1;
return 0;
if(trackTheta == 0) {
- if(fDebug) cout << " Theta sol unique " << thetaCerenkov << endl;
-
thetaAtQuartz = thetaCerenkov;
SetThetaAtQuartz(thetaAtQuartz);
return;
Double_t underSqrt = 1 + b*b - c*c;
- if(fDebug)
- {
- cout << " trackTheta " << trackTheta << endl;
- cout << " TrackPhi " << trackPhi << endl;
- cout << " PhiPoint " << phiPoint << endl;
- cout << " ThetaCerenkov " << thetaCerenkov << endl;
- cout << " den b c " << den << " b " << b << " c " << c << endl;
- }
-
if(underSqrt < 0) {
- if(fDebug) cout << " sqrt negative !!!!" << underSqrt << endl;
SetThetaAtQuartz(999.);
return;
}
Double_t thetaSol1 = 2*TMath::ATan(sol1);
Double_t thetaSol2 = 2*TMath::ATan(sol2);
- if(fDebug) cout << " Theta sol 1 " << thetaSol1
- << " Theta sol 2 " << thetaSol2 << endl;
-
if(thetaSol1>0 && thetaSol1 < TMath::Pi()) thetaAtQuartz = (Float_t)thetaSol1;
if(thetaSol2>0 && thetaSol2 < TMath::Pi()) thetaAtQuartz = (Float_t)thetaSol2;
+// AliDebug(1,Form(" Theta @ quartz window %f ",thetaAtQuartz));
+
SetThetaAtQuartz(thetaAtQuartz);
}
const Float_t kTollerance = 0.05;
- // Float_t pmod = GetTrackMomentum();
- // Float_t trackTheta = GetTrackTheta();
- // Float_t trackPhi = GetTrackPhi();
Float_t phiPoint = GetPhiPoint();
Float_t xPoint = GetEntranceX();
Float_t yPoint = GetEntranceY();
- Float_t distPoint = sqrt(xPoint*xPoint + yPoint*yPoint);
+ Float_t distPoint = TMath::Sqrt(xPoint*xPoint + yPoint*yPoint);
- if(fDebug) cout << " DistPoint " << distPoint << endl;
+ AliDebug(1,Form(" DistPoint %f ",distPoint));
// Star minimization...
radiusMean = FromEmissionToCathode();
}
- if(fDebug) cout << " r1 " << radiusMin << " rmean "
- << radiusMean << " r2 " << radiusMax << endl;
+ AliDebug(1,Form(" r1 %f rmean %f r2 %f",radiusMin,radiusMean,radiusMax));
while (TMath::Abs(radiusMean-distPoint) > kTollerance)
{
nIteration++;
if(nIteration>=50) {
- if(fDebug) printf(" max iterations in FindPhotonCerenkov\n");
+ AliDebug(1,Form(" max iterations in FindPhotonCerenkov ",nIteration));
SetThetaPhotonCerenkov(999.);
return;
}
}
+ AliDebug(1,Form(" distpoint %f radius %f ",distPoint,radiusMean));
SetThetaPhotonCerenkov(thetaCerMean);
}
Float_t x0 = xemiss + shiftX;
Float_t y0 = yemiss + shiftY;
- // Float_t pmod = GetTrackMomentum();
- // Float_t trackTheta = GetTrackTheta();
- // Float_t trackPhi = GetTrackPhi();
SetPhotonEnergy(6.85);
SetFreonRefractiveIndex();
return 999;
}
//__________________________________________________________________________________________________
-Float_t AliRICHRecon::PhotonPositionOnCathode()
-{
- // find the photon position on the CsI
- // Float_t massOfParticle;
- Float_t beta;
- Float_t nfreon;
-
-
- SetPhotonEnergy(6.85);
- SetEmissionPoint(fRadiatorWidth/2.);
- SetMassHypotesis(0.139567);
-
- SetBetaOfParticle();
- SetFreonRefractiveIndex();
-
- beta = GetBetaOfParticle();
- nfreon = GetFreonRefractiveIndex();
-
-
- Float_t radius = FromEmissionToCathode();
- if (radius == 999.) return 999.;
-
- return 0;
-}
-
void AliRICHRecon::FindPhotonAnglesInDRS()
{
// Setup the rotation matrix of the track...
Float_t argY = ytoentr - emissionPoint*tan(trackTheta)*sin(trackPhi);
Float_t argX = xtoentr - emissionPoint*tan(trackTheta)*cos(trackPhi);
- Float_t phipad = atan2(argY,argX);
+ Float_t phi = atan2(argY,argX);
- SetPhiPoint(phipad);
+ SetPhiPoint(phi);
}
Int_t nPhotonHough = 0;
Float_t thetaCerenkov = GetThetaCerenkov();
- if(fDebug) cout << " fThetaCerenkov " << thetaCerenkov << endl;
+ AliDebug(1,Form(" fThetaCerenkov ",thetaCerenkov));
Float_t thetaDist= thetaCerenkov - fThetaMin;
Int_t steps = (Int_t)(thetaDist / fDTheta);
tmin = tavg - 0.5*fWindowWidth;
tmax = tavg + 0.5*fWindowWidth;
- if(fDebug) cout << " tmin " << tmin << " tmax " << tmax << endl;
- if(fDebug) cout << " thetac " << thetaCerenkov << endl;
-
// Int_t candidatePhotonsNumber = GetCandidatePhotonsNumber();
Int_t nPhotons = GetPhotonsNumber();
SetHoughPhotons(nPhotonHough);
}
-void AliRICHRecon::DrawEvent(Int_t flag) const
-{
- // draw event with rings
-
- flag=1; // dummy to be removed...
-}
-
-Float_t AliRICHRecon::FindMassOfParticle()
-{
- // find mass of the particle from theta cerenkov
-
- Float_t pmod = GetTrackMomentum();
-
- SetPhotonEnergy(6.85);
- SetFreonRefractiveIndex();
-
- Float_t thetaCerenkov = GetThetaCerenkov();
- FindBetaFromTheta(thetaCerenkov);
-
- Double_t beta = (Double_t)(GetBetaOfParticle());
- Double_t den = 1. - beta*beta;
- if(den<=0.) return 999.;
-
- Double_t gamma = 1./TMath::Sqrt(den);
-
- Float_t mass = pmod/(beta*(Float_t)gamma);
-
- return mass;
-}
-
-
-void AliRICHRecon::FillHistograms()
-{
- // fill histograms..
-
- Float_t fittedTrackTheta, fittedTrackPhi;
-
- Float_t thetaCerenkov = GetThetaCerenkov();
- if(thetaCerenkov == 999.) return;
-
- Float_t vertZ = GetEventVertexZ();
-
- Float_t trackTheta = GetTrackTheta();
- Float_t trackPhi = GetTrackPhi();
- Float_t pmod = GetTrackMomentum();
- Float_t pt = GetTrackPt();
- Float_t trackEta = GetTrackEta();
- Int_t q = GetTrackCharge();
- Float_t tPCLastZ = GetTrackTPCLastZ();
- Float_t minDist = GetMinDist();
-
- fittedTrackTheta = GetFittedTrackTheta();
- fittedTrackPhi = GetFittedTrackPhi();
- Int_t fittednPhotonHough = GetFittedHoughPhotons();
-
- if(fDebug)
- {
- cout << " p " << pmod << " ThetaC " << thetaCerenkov
- << " rings " << fNrings << endl;
- }
-
- Int_t nPhotonHough = GetHoughPhotons();
-// Float_t nPhotonHoughNorm = GetHoughPhotonsNorm();
- Float_t inRing = GetPortionOfRing();
-
- Float_t massOfParticle = FindMassOfParticle();
-
- Float_t houghArea = GetHoughArea();
- Float_t multiplicity = GetEventMultiplicity();
-
-
- Float_t var[20];
-
- var[0] = 0;
- var[1] = 0;
- var[2] = vertZ;
- var[3] = pmod;
- var[4] = pt;
- var[5] = trackEta;
- var[6] = trackTheta;
- var[7] = trackPhi;
- var[8] = fittedTrackTheta;
- var[9] = fittedTrackPhi;
- var[10] = q;
- var[11] = thetaCerenkov;
- var[12] = (Float_t)nPhotonHough;
- var[13] = (Float_t)fittednPhotonHough;
- var[14] = inRing;
- var[15] = massOfParticle;
- var[16] = houghArea;
- var[17] = multiplicity;
- var[18] = tPCLastZ;
- var[19] = minDist;
-
- fNtuple->Fill(var);
-
-
- fittedTrackTheta = GetFittedTrackTheta();
- fittedTrackPhi = GetFittedTrackPhi();
-
-
-
- if(thetaCerenkov > 0.505 && thetaCerenkov < 0.605) {
- SetPhotonEnergy(6.85);
- SetFreonRefractiveIndex();
- }
-
- Int_t nPhotons = GetPhotonsNumber();
-
- for (Int_t j=0; j < nPhotons;j++)
- SetPhotonIndex(j);
-}//FillHistograms()
-//__________________________________________________________________________________________________
-void AliRICHRecon::Minimization()
-{
- // minimization to find the best theta and phi of the track
-
- Double_t arglist;
- Int_t ierflag = 0;
-
- static Double_t vstart[2];
- static Double_t lower[2], upper[2];
- static Double_t step[2]={0.001,0.001};
-
- Double_t trackThetaNew,trackPhiNew;
- TString chname;
- Double_t eps, b1, b2;
- Int_t ierflg;
-
- gAliRICHminuit = new TMinuit(2);
- gAliRICHminuit->SetObjectFit((TObject *)this);
- gAliRICHminuit->SetFCN(fcnrecon);
- gAliRICHminuit->mninit(5,10,7);
-
- vstart[0] = (Double_t)GetTrackTheta();
- vstart[1] = (Double_t)GetTrackPhi();
-
- lower[0] = vstart[0] - 0.03;
- if(lower[0] < 0) lower[0] = 0.;
- upper[0] = vstart[0] + 0.03;
- lower[1] = vstart[1] - 0.03;
- upper[1] = vstart[1] + 0.03;
-
-
- gAliRICHminuit->mnparm(0,"theta",vstart[0],step[0],lower[0],upper[0],ierflag);
- gAliRICHminuit->mnparm(1," phi ",vstart[1],step[1],lower[1],upper[1],ierflag);
-
- arglist = -1;
-
- // gAliRICHminuit->FixParameter(0);
-
- gAliRICHminuit->SetPrintLevel(-1);
-// gAliRICHminuit->mnexcm("SET PRI",&arglist, 1, ierflag);
- gAliRICHminuit->mnexcm("SET NOGR",&arglist, 1, ierflag);
- gAliRICHminuit->mnexcm("SET NOW",&arglist, 1, ierflag);
- arglist = 1;
- gAliRICHminuit->mnexcm("SET ERR", &arglist, 1,ierflg);
- arglist = -1;
-
- // gAliRICHminuit->mnscan();
-
-// gAliRICHminuit->mnexcm("SIMPLEX",&arglist, 0, ierflag);
- gAliRICHminuit->mnexcm("MIGRAD",&arglist, 0, ierflag);
- gAliRICHminuit->mnexcm("EXIT" ,&arglist, 0, ierflag);
-
- gAliRICHminuit->mnpout(0,chname, trackThetaNew, eps , b1, b2, ierflg);
- gAliRICHminuit->mnpout(1,chname, trackPhiNew, eps , b1, b2, ierflg);
-
- //values after the fit...
- SetFittedTrackTheta((Float_t)trackThetaNew);
- SetFittedTrackPhi((Float_t)trackPhiNew);
-
- delete gAliRICHminuit;
-
-}
-
-void AliRICHRecon::EstimationOfTheta()
-{
- // theta estimate
-
- Int_t nPhotons = 0;
-
- Float_t shiftX = GetShiftX();
- Float_t shiftY = GetShiftY();
-
- Float_t *candidatePhotonX = GetCandidatePhotonX();
- Float_t *candidatePhotonY = GetCandidatePhotonY();
-
- Int_t nPhotonsCandidates = GetCandidatePhotonsNumber();
-
- // cout << "MINIM: Nphotons " << nPhotonsCandidates << endl;
-
- for (Int_t j=0; j < nPhotonsCandidates; j++)
- {
-
- SetPhotonIndex(j);
-
- if(!GetPhotonFlag()) continue;
-
- Float_t xtoentr = candidatePhotonX[j] - shiftX;
- Float_t ytoentr = candidatePhotonY[j] - shiftY;
-
- SetEntranceX(xtoentr);
- SetEntranceY(ytoentr);
-
- FindPhiPoint();
-
- FindThetaPhotonCerenkov();
-
- Float_t thetaPhotonCerenkov = GetThetaPhotonCerenkov();
-
- // cout << " ACCEPTED!!! " << thetaPhotonCerenkov << endl;
-
- SetPhotonEta(thetaPhotonCerenkov);
-
- nPhotons++;
-
- }
-
- Float_t xmean = 0.;
- Float_t x2mean = 0.;
- Int_t nev = 0;
-
- for (Int_t j=0; j < nPhotonsCandidates;j++)
- {
- SetPhotonIndex(j);
-
- Float_t eta = GetPhotonEta();
-
- if(eta != -999.)
- {
- if(GetPhotonFlag() == 2)
- {
- xmean += eta;
- x2mean += eta*eta;
- nev++;
- }
- }
- }
-
- if(nev > 0)
- {
- xmean /=(Float_t)nev;
- x2mean /=(Float_t)nev;
- } else {
- xmean = 0.;
- x2mean = 0.;
- }
-
- Float_t vRMS = sqrt(x2mean - xmean*xmean);
-
- // cout << " RMS " << vRMS;
-
- SetEstimationOfTheta(xmean);
- SetEstimationOfThetaRMS(vRMS);
-}
-
-void fcnrecon(Int_t& /*npar*/, Double_t* /*gin*/, Double_t &f, Double_t *par, Int_t)
-{
- // function to be minimized
- AliRICHRecon *gMyRecon = (AliRICHRecon*)gAliRICHminuit->GetObjectFit();
-
- Float_t p0 = (Float_t)par[0];
- Float_t p1 = (Float_t)par[1];
-
- gMyRecon->SetTrackTheta(p0);
- gMyRecon->SetTrackPhi(p1);
-
- gMyRecon->EstimationOfTheta();
- Float_t vRMS = gMyRecon->GetEstimationOfThetaRMS();
-
- Int_t houghPhotons = gMyRecon->GetHoughPhotons();
-
-
- f = (Double_t)(1000*vRMS/(Float_t)houghPhotons);
-
-// if(fDebug) cout << " f " << f
-// << " theta " << par[0] << " phi " << par[1]
-// << " HoughPhotons " << houghPhotons << endl;
-//
-// if(fDebug&&iflag == 3)
-// {
-// cout << " --- end convergence...summary --- " << endl;
-// cout << " theta " << par[0] << endl;
-// cout << " phi " << par[1] << endl;
-// }
-}
-
-void AliRICHRecon::Waiting()
-{
- // wait, wait....
- if(!fIsDISPLAY) return;
- cout << " Press any key to continue...";
-
-// gSystem->ProcessEvents();
- getchar();
-
- cout << endl;
-
- return;
-}
-
#include <TTask.h>
-class AliRICH;
-class TFile;
-class TNtuple;
-class TCanvas;
+class AliRICHHelix;
class AliRICHRecon : public TTask
{
public :
- AliRICHRecon(const char*, const char*);
- ~AliRICHRecon(){EndProcessEvent();}
+ AliRICHRecon(AliRICHHelix *pHelix,TClonesArray *pClusters,Int_t iMipId);
+ virtual ~AliRICHRecon(){;}
- AliRICH* Rich() {return fRich;} //main pointer to RICH
- void StartProcessEvent(); //
- void EndProcessEvent(); //
- //PH void InitRecon(); //
- void PatRec(); //
- void Minimization(); //
- void FillHistograms(); //
+ Double_t ThetaCerenkov(); // it returns reconstructed Theta Cerenkov
void FindThetaPhotonCerenkov(); //
void FindAreaAndPortionOfRing(); //
void FindEmissionPoint(); //
void FindWeightThetaCerenkov(); //
void EstimationOfTheta(); //
void FindIntersectionWithDetector(); //
- void Waiting(); //
- Float_t FindMassOfParticle(); //
Float_t Cerenkovangle(Float_t n, Float_t b);//
- Float_t PhotonPositionOnCathode(); //
Int_t PhotonInBand(); //
- void DrawEvent(Int_t flag) const; //
Int_t CheckDetectorAcceptance() const; //
Int_t GetFittedHoughPhotons() const{ return fFittedHoughPhotons;} //
Int_t GetPhotonFlag() const{ return fPhotonFlag[fPhotonIndex];} //
Int_t GetTrackIndex() const{ return fTrackIndex;} //
Int_t GetCandidatePhotonsNumber() const{ return fCandidatePhotonsNumber;} //
Int_t GetHoughPhotons() const{ return fHoughPhotons;} //
- Float_t GetEventVertexZ() const{ return fEventVertZ;} //
- Float_t GetEventMultiplicity() const{ return fEventMultiplicity;} //
Float_t GetPhotonEnergy() const{ return fPhotonEnergy;} //
Float_t GetFreonRefractiveIndex() const{ return fFreonRefractiveIndex;} //
Float_t GetQuartzRefractiveIndex() const{ return fQuartzRefractiveIndex;} //
Float_t GetGasRefractiveIndex() const{ return fGasRefractiveIndex;} //
- Float_t GetFreonScaleFactor() const{ return fFreonScaleFactor;} //
Float_t GetEmissionPoint() const{ return fLengthEmissionPoint;} //
Float_t GetMassHypotesis() const{ return fMassHypotesis;} //
Float_t GetBetaOfParticle() const{ return fTrackBeta;} //
Float_t GetEntranceY() const{ return fYtoentr;} //
Float_t GetThetaCerenkov() const{ return fThetaCerenkov;} //
Float_t GetThetaPhotonCerenkov() const{ return fThetaPhotonCerenkov;} //
- Float_t GetTrackMomentum() const{ return fTrackMomentum;} //
- Float_t GetTrackEta() const{ return fTrackEta;} //
Float_t GetTrackTheta() const{ return fTrackTheta;} //
Float_t GetTrackPhi() const{ return fTrackPhi;} //
- Float_t GetTrackPt() const{ return fTrackPt;} //
- Float_t GetTrackTPCLastZ() const{ return fTrackTPCLastZ;} //
- Float_t GetMinDist() const{ return fMinDist;} //
Float_t GetXPointOnCathode() const{ return fPhotonLimitX;} //
Float_t GetYPointOnCathode() const{ return fPhotonLimitY;} //
Float_t GetThetaPhotonInDRS() const{ return fThetaPhotonInDRS;} //
Float_t GetPhotonEta() const{ return fPhotonEta[fPhotonIndex];} //
Float_t GetPhotonWeight() const{ return fPhotonWeight[fPhotonIndex];} //
Float_t GetHoughRMS() const{ return fHoughRMS;} //
- Float_t* GetCandidatePhotonX() const{ return fCandidatePhotonX;} //
- Float_t* GetCandidatePhotonY() const{ return fCandidatePhotonY;} //
- Float_t GetHoughPhotonsNorm() const{ return fHoughPhotonsNorm;} //
Float_t GetFittedTrackTheta() const{ return fFittedTrackTheta;} //
Float_t GetFittedTrackPhi() const{ return fFittedTrackPhi;} //
Float_t GetFittedThetaCerenkov() const{ return fFittedThetaCerenkov;} //
Float_t GetEstimationOfTheta() const{ return fEstimationOfTheta;} //
Float_t GetEstimationOfThetaRMS() const{ return fEstimationOfThetaRMS;} //
- void SetEventVertexZ(Float_t EventVertZ) { fEventVertZ = EventVertZ;} //
- void SetEventMultiplicity(Float_t EventMultiplicity) { fEventMultiplicity = EventMultiplicity;} //
void SetPhotonEnergy(Float_t PhotonEnergy) { fPhotonEnergy = PhotonEnergy;} //
void SetFreonRefractiveIndex() {fFreonRefractiveIndex = fFreonScaleFactor*(1.177+0.0172*fPhotonEnergy);}//
void SetQuartzRefractiveIndex() {fQuartzRefractiveIndex = sqrt(1+(46.411/(113.763556-TMath::Power(fPhotonEnergy,2)))+(228.71/(328.51563-TMath::Power(fPhotonEnergy,2))));}//
void SetGasRefractiveIndex() { fGasRefractiveIndex = 1.;} //
void SetFreonScaleFactor(Float_t FreonScaleFactor) {fFreonScaleFactor = FreonScaleFactor;} //
void SetEmissionPoint(Float_t LengthEmissionPoint) { fLengthEmissionPoint = LengthEmissionPoint;} //
- void SetMassHypotesis(Float_t mass) {fMassHypotesis = mass;} //
- void SetBetaOfParticle() { fTrackBeta = fTrackMomentum/sqrt(TMath::Power(fTrackMomentum,2)+TMath::Power(fMassHypotesis,2));}//
void SetEntranceX(Float_t Xtoentr) { fXtoentr = Xtoentr;} //
void SetEntranceY(Float_t Ytoentr) { fYtoentr = Ytoentr;} //
void SetThetaPhotonInTRS(Float_t Theta) {fThetaPhotonInTRS = Theta;} //
void SetRadiusOuterRing(Float_t OuterRadius) {fOuterRadius = OuterRadius;} //
void SetThetaCerenkov(Float_t ThetaCer) {fThetaCerenkov = ThetaCer;} //
void SetThetaPhotonCerenkov(Float_t ThetaPhotCer) {fThetaPhotonCerenkov = ThetaPhotCer;} //
- void SetTrackMomentum(Float_t TrackMomentum) {fTrackMomentum = TrackMomentum;} //
- void SetTrackEta(Float_t TrackEta) {fTrackEta = TrackEta;} //
void SetTrackTheta(Float_t TrackTheta) { fTrackTheta = TrackTheta;} //
void SetTrackPhi(Float_t TrackPhi) { fTrackPhi = TrackPhi;} //
- void SetTrackPt(Float_t TrackPt) { fTrackPt = TrackPt;} //
void SetTrackCharge(Int_t TrackCharge) { fTrackCharge = TrackCharge;} //
- void SetTrackTPCLastZ(Float_t TrackTPCLastZ) { fTrackTPCLastZ = TrackTPCLastZ;} //
- void SetMinDist(Float_t MinDist) { fMinDist = MinDist;} //
void SetShiftX(Float_t ShiftX) { fShiftX = ShiftX;} //
void SetShiftY(Float_t ShiftY) { fShiftY = ShiftY;} //
void SetDetectorWhereX(Float_t Xcoord) { fXcoord = Xcoord;} //
void SetHoughRMS(Float_t HoughRMS) { fHoughRMS = HoughRMS;} //
void SetMipIndex(Int_t MipIndex) { fMipIndex = MipIndex;} //
void SetTrackIndex(Int_t TrackIndex) { fTrackIndex = TrackIndex;} //
- void SetCandidatePhotonX(Float_t *CandidatePhotonX) { fCandidatePhotonX = CandidatePhotonX;} //
- void SetCandidatePhotonY(Float_t *CandidatePhotonY) { fCandidatePhotonY = CandidatePhotonY;} //
- void SetCandidatePhotonsNumber(Int_t CandidatePhotonsNumber) { fCandidatePhotonsNumber = CandidatePhotonsNumber;}//
void SetHoughPhotons(Int_t HoughPhotons) { fHoughPhotons = HoughPhotons;} //
void SetHoughPhotonsNorm(Float_t HoughPhotonsNorm) { fHoughPhotonsNorm = HoughPhotonsNorm;} //
void SetFittedTrackTheta(Float_t FittedTrackTheta) { fFittedTrackTheta = FittedTrackTheta;} //
Float_t FromEmissionToCathode(); //
protected:
- AliRICH* fRich; // main poiter to RICH
+ TClonesArray *fpClusters; // poiter to clusters
Int_t fTrackCharge; // charge track
Int_t fMipIndex; // index for Mip
Int_t fTrackIndex; // index for track
Int_t fCandidatePhotonsNumber; // number of candidate photons
Int_t fHoughPhotons; // n. photons after Hough
Int_t fFittedHoughPhotons; // n. photons after Hough and after minimization
- Float_t fEventVertZ; // z coord. of the primary vertex
- Float_t fEventMultiplicity; // event primary multiplicity
+
Float_t fTrackTheta; // Theta of track at RICH
Float_t fTrackPhi; // Phi of track at RICH
- Float_t fTrackMomentum; // track momentum
- Float_t fTrackEta; // track pseudorapidity
- Float_t fTrackPt; // pt of track
- Float_t fTrackTPCLastZ; // z last point of the TPC
Float_t fMinDist; // min distance between extrapolated track and MIP
Float_t fTrackBeta; // beta of the track
Float_t fXtoentr; // X entrance to RICH
Float_t fPhotonLimitX; // X phys limit for photon
Float_t fPhotonLimitY; // Y phys limit for photon
Float_t fDistanceFromCluster; // distance from cluster
- Float_t fMassHypotesis; // reconstructed mass
Float_t fCerenkovAnglePad; // cherenkov angle of pad
Float_t fThetaPhotonCerenkov; // theta cerenkov for photon
Float_t fShiftX; // x shift to entrance in radiator
Float_t fYcoord; // ..
Float_t fIntersectionX; // ..
Float_t fIntersectionY; // ..
+ Float_t fMassHypotesis; //
Float_t fThetaOfRing; // theta of ring
Float_t fAreaOfRing; // area of the ring
Float_t fPortionOfRing; // fraction of the accepted ring
Int_t fThetaBin; // bin in theta
Float_t fThetaMin,fThetaMax; // min max
Float_t fXmin,fXmax,fYmin,fYmax; // xy min max
- TFile *fOutFile; // histogram output
- TNtuple *fNtuple; // output ntuple
- TCanvas *fDisplay; // for display
Int_t fNrings; //current number of reconstructed rings
- Bool_t fDebug; // flag for debug
- Bool_t fIsDISPLAY; // flag for display
Bool_t fIsWEIGHT; // flag to consider weight procedure
Bool_t fIsBACKGROUND; // flag to simulate bkg
- Bool_t fIsMINIMIZER; // flag to start with (theta,phi) minimization
- Int_t fNpadX; // npadx
- Int_t fNpadY; // npady
- Float_t fPadSizeX; // sizepad x
- Float_t fPadSizeY; // sizepad y
Float_t fRadiatorWidth; // radiator width
Float_t fQuartzWidth; // quartz width
Float_t fGapWidth; // gap width
* provided "as is" without express or implied warranty. *
**************************************************************************/
-/* $Id$ */
///////////////////////////////////////////////////////////////////////////////
// //
#include "AliRICHReconstructor.h"
-#include "AliRunLoader.h"
-#include "AliRun.h"
#include "AliRICHClusterFinder.h"
-
+#include "AliRICHHelix.h"
+#include <AliRunLoader.h>
+#include <AliRun.h>
+#include <AliESD.h>
ClassImp(AliRICHReconstructor)
-
-//_____________________________________________________________________________
-void AliRICHReconstructor::Reconstruct(AliRunLoader* runLoader) const
+//__________________________________________________________________________________________________
+void AliRICHReconstructor::Reconstruct(AliRunLoader* pAL) const
{
-// reconstruct clusters
-
- AliRICH* rich = GetRICH(runLoader);
- if (!rich) return;
- AliRICHClusterFinder clusterer(rich);
- clusterer.Exec();
+//Finds clusters out of digits
+ AliDebug(1,"Start cluster finder.");AliRICHClusterFinder clus(GetRICH(pAL)); clus.Exec();
}
-
-//_____________________________________________________________________________
-void AliRICHReconstructor::FillESD(AliRunLoader* /*runLoader*/,
- AliESD* /*esd*/) const
+//__________________________________________________________________________________________________
+void AliRICHReconstructor::FillESD(AliRunLoader* /*pAL*/, AliESD* pESD) const
{
-// nothing to be done
-
-}
-
+//This methode fills AliESDtrack with information from RICH
+ AliDebug(1,Form("Start with %i tracks",pESD->GetNumberOfTracks()));
+/* const Double_t masses[5]={0.000511,0.105658,0.139567,0.493677,0.93828};//electron,muon,pion,kaon,proton
+ const Double_t refIndex = 1.29052;
-//_____________________________________________________________________________
-AliRICH* AliRICHReconstructor::GetRICH(AliRunLoader* runLoader) const
+ Double_t thetaTh[5];
+ Double_t sinThetaThNorm;
+ Double_t sigmaThetaTh[5];
+ Double_t height[5];
+ Double_t totalHeight=0;
+ Double_t x[3],p[3]; //tmp storage for track parameters
+
+ for(Int_t iTrackN=0;iTrackN<pESD->GetNumberOfTracks();iTrackN++){//ESD tracks loop
+ AliESDtrack *pTrack = pESD->GetTrack(iTrackN);// get next reconstructed track
+// if((pTrack->GetStatus()&AliESDtrack::kTOFout)==0) continue; //ignore tracks not recontructed by TOF
+ pTrack->GetXYZ(x); pTrack->GetPxPyPz(p); Double_t pmod=pTrack->GetP();//get running track parameters
+ continue;
+// AliRICHHelix helix(x[0],x[1],x[2],p[0],p[1],p[2]); //construct helix from running track parameters
+
+// TVector rad(1,5); TVector pc(1,5);
+// helix.RichIntersect(GetRICH(pAL)->P(),rad,pc); //returns cross point of track with RICH PC in LRS
+
+ for(Int_t iPart=4;iPart>=0;iPart--){
+ Double_t cosThetaTh = TMath::Sqrt(masses[iPart]*masses[iPart]+pmod*pmod)/(refIndex*pmod);
+ if(cosThetaTh>=1) {break;}
+ thetaTh[iPart] = TMath::ACos(cosThetaTh);
+ sinThetaThNorm = TMath::Sin(thetaTh[iPart])/TMath::Sqrt(1-1/(refIndex*refIndex));
+ sigmaThetaTh[iPart] = (0.014*(1/sinThetaThNorm-1) + 0.0043)*1.25;
+ height[iPart] = TMath::Gaus(thetaExp,thetaTh[iPart],sigmaThetaTh[iPart]);
+ totalHeight +=height[iPart];
+ }
+
+ Double_t richPID[5];
+ for(Int_t iPart=0;iPart<5;iPart++) richPID[iPart] = height[iPart]/totalHeight;
+ pTrack->SetRICHpid(richPID);
+ }//ESD tracks loop
+ */
+}//FillESD
+//__________________________________________________________________________________________________
+AliRICH* AliRICHReconstructor::GetRICH(AliRunLoader* pAL) const
{
// get the RICH detector
- if (!runLoader->GetAliRun()) runLoader->LoadgAlice();
- if (!runLoader->GetAliRun()) {
- Error("GetRICH", "couldn't get AliRun object");
- return NULL;
- }
- AliRICH* rich = (AliRICH*) runLoader->GetAliRun()->GetDetector("RICH");
- if (!rich) {
- Error("GetRICH", "couldn't get RICH detector");
- return NULL;
- }
- return rich;
+ if (!pAL->GetAliRun()) pAL->LoadgAlice();
+ if (!pAL->GetAliRun()) {AliError("couldn't get AliRun object"); return NULL; }
+ AliRICH* pRich = (AliRICH*) pAL->GetAliRun()->GetDetector("RICH");
+ if (!pRich) {AliError("couldn't get RICH detector"); return NULL; }
+ return pRich;
}
+//__________________________________________________________________________________________________
-#ifndef ALIRICHRECONSTRUCTOR_H
-#define ALIRICHRECONSTRUCTOR_H
+#ifndef AliRICHReconstructor_h
+#define AliRICHReconstructor_h
/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
* See cxx source for full Copyright notice */
-/* $Id$ */
-
-#include "AliReconstructor.h"
+#include <AliReconstructor.h>
+#include <AliLog.h>
+#include "AliRICHTracker.h"
class AliRICH;
-
-class AliRICHReconstructor: public AliReconstructor {
+class AliRICHReconstructor: public AliReconstructor
+{
public:
- AliRICHReconstructor(): AliReconstructor() {};
- virtual ~AliRICHReconstructor() {};
-
- virtual void Reconstruct(AliRunLoader* runLoader) const;
- virtual void FillESD(AliRunLoader* runLoader, AliESD* esd) const;
-
+ AliRICHReconstructor(): AliReconstructor() {AliDebug(1,"Start.");}
+ virtual ~AliRICHReconstructor() {AliDebug(1,"Start.");}
+ virtual AliTracker* CreateTracker(AliRunLoader*)const {AliDebug(1,"Start.");return new AliRICHTracker;} //virtual from AliReconstructor
+ virtual void Reconstruct(AliRunLoader* pAL) const; //virtual from AliReconstructor
+ virtual void FillESD(AliRunLoader* pAL, AliESD* pESD) const; //virtual from AliReconstructor
private:
- AliRICH* GetRICH(AliRunLoader* runLoader) const;
+ AliRICH* GetRICH(AliRunLoader* pAL) const;
- ClassDef(AliRICHReconstructor, 0) // class for the RICH reconstruction
+ ClassDef(AliRICHReconstructor, 0) //class for the RICH reconstruction
};
#endif
--- /dev/null
+#include "AliRICHTracker.h"
+#include <AliESD.h>
+#include <TVector3.h>
+#include <TTree.h>
+#include "AliRICH.h"
+#include "AliRICHHelix.h"
+#include <AliMagF.h>
+#include "AliRICHRecon.h"
+#include <AliStack.h>
+#include <TParticle.h>
+
+ClassImp(AliRICHTracker)
+
+
+Int_t AliRICHTracker::PropagateBack(AliESD *pESD)
+{
+ Int_t iNtracks=pESD->GetNumberOfTracks();
+ Double_t b=GetFieldMap()->SolenoidField()/10;// magnetic field in Tesla
+ AliDebug(1,Form("Start with %i tracks in %f Tesla field",iNtracks,b));
+ Double_t xb[3],pb[3];//tmp storage for track parameters
+ TVector3 x0(0,0,0); TVector3 p0(0,0,0);//tmp storage for AliRICHHelix
+
+ AliRICH *pRich=((AliRICH*)gAlice->GetDetector("RICH"));
+
+// pRich->GetLoader()->GetRunLoader()->LoadHeader();pRich->GetLoader()->GetRunLoader()->LoadKinematics();
+// AliStack *pStack = pRich->GetLoader()->GetRunLoader()->Stack();
+// TParticle *pParticle = pStack->Particle(0);
+// p0.SetMagThetaPhi(pParticle->P(),pParticle->Theta(),pParticle->Phi());
+
+ for(Int_t iTrackN=0;iTrackN<iNtracks;iTrackN++){//ESD tracks loop
+ AliESDtrack *pTrack = pESD->GetTrack(iTrackN);// get next reconstructed track
+// if((pTrack->GetStatus()&AliESDtrack::kTOFout)==0) continue; //ignore tracks not recontructed by TOF
+ pTrack->GetXYZ(xb);
+ pTrack->GetPxPyPz(pb);
+ Int_t status=pTrack->GetStatus()&AliESDtrack::kTOFout;//get running track parameters
+ AliDebug(1,Form("Track %i pmod=%f mass=%f stat=%i",iTrackN,pTrack->GetP(),pTrack->GetMass(),status));
+// x.Print();p.Print();
+ x0.SetXYZ(xb[0],xb[1],xb[2]); p0.SetXYZ(xb[0],xb[1],xb[2]);
+ AliRICHHelix helix(x0,p0,pTrack->GetSign(),b);
+ Int_t iChamber=helix.RichIntersect(pRich->P());
+ AliDebug(1,Form("intersection with %i chamber found",iChamber));
+ if(!iChamber) continue;//intersection with no chamber found
+
+ Double_t distMip=9999; //min distance between clusters and track position on PC
+ Int_t iMipId=0; //index of that min distance cluster
+ for(Int_t iClusN=0;iClusN<pRich->Clusters(iChamber)->GetEntries();iClusN++){//clusters loop for intersected chamber
+ AliRICHcluster *pClus=(AliRICHcluster*)pRich->Clusters(iChamber)->UncheckedAt(iClusN);//get pointer to current cluster
+ Double_t distCurrent=pClus->DistTo(helix.PosPc());//ditance between current cluster and helix intersection with PC
+ if(distCurrent<distMip){distMip=distCurrent;iMipId=iClusN;}//find cluster nearest to the track
+
+ AliDebug(1,Form("Ploc (%f,%f,%f) dist= %f",helix.Ploc().Mag(),helix.Ploc().Theta()*TMath::RadToDeg(),
+ helix.Ploc().Phi()*TMath::RadToDeg(),pClus->DistTo(helix.PosPc())));
+ }////clusters loop for intersected chamber
+
+ AliDebug(1,Form("Min distance cluster: %i dist is %f",iMipId,distMip));
+
+ AliRICHRecon recon(&helix,pRich->Clusters(iChamber),iMipId);
+ Double_t thetaCerenkov=recon.ThetaCerenkov(); //search for mean Cerenkov angle for this track
+ AliDebug(1,Form("FINAL Theta Cerenkov=%f",thetaCerenkov));
+ pTrack->SetRICHsignal(thetaCerenkov);
+
+
+ }//ESD tracks loop
+ AliDebug(1,"Stop.");
+ return 0;
+} //pure virtual from AliTracker
+
+Int_t AliRICHTracker::LoadClusters(TTree *pTree)
+{
+// Load clusters for RICH
+ AliDebug(1,"Start."); pTree->GetEntry(0); AliDebug(1,"Stop."); return 0;
+}
--- /dev/null
+#ifndef AliRICHTracker_h
+#define AliRICHTracker_h
+
+#include <AliTracker.h>
+#include <AliLog.h>
+
+class AliCluster;
+class AliESD;
+class TTree;
+
+class AliRICHTracker : public AliTracker
+{
+public:
+ AliRICHTracker() :AliTracker() {AliDebug(1,"Start.");}
+ virtual ~AliRICHTracker() {AliDebug(1,"Start.");}
+ Int_t Clusters2Tracks(AliESD *) {AliDebug(1,"Start.");return 0;} //pure virtual from AliTracker
+ Int_t RefitInward(AliESD *) {AliDebug(1,"Start.");return 0;} //pure virtual from AliTracker
+ void UnloadClusters() {AliDebug(1,"Start.");} //pure virtual from AliTracker
+ AliCluster *GetCluster(Int_t )const {AliDebug(1,"Start.");return 0;} //pure virtual from AliTracker
+ Int_t PropagateBack(AliESD *); //pure virtual from AliTracker
+ Int_t LoadClusters(TTree *); //pure virtual from AliTracker
+
+protected:
+
+ ClassDef(AliRICHTracker,0)
+};//class AliRICHTracker
+
+#endif//AliRICHTracker_h
void AliRICHv0::StepManager()
{
//This StepManager is a provision for different test-learn activities on the current MC layer
+ static Int_t iStepN;
const char *sParticle;
switch(gMC->TrackPid()){
- case kProton:
- sParticle="proton";break;
- case kNeutron:
- sParticle="neutron";break;
- case kGamma:
- sParticle="gamma";break;
- case kCerenkov:
- sParticle="photon";break;
- case kPi0:
- sParticle="Pi0";break;
- case kElectron:
- sParticle="electron";break;
- default:
- sParticle="not known";break;
+ case kProton: sParticle="proton" ;break;
+ case kNeutron: sParticle="neutron" ;break;
+ case kGamma: sParticle="gamma" ;break;
+ case kCerenkov: sParticle="photon" ;break;
+ case kPi0: sParticle="Pi0" ;break;
+ case kPiPlus: sParticle="Pi+" ;break;
+ case kPiMinus: sParticle="Pi-" ;break;
+ case kElectron: sParticle="electron" ;break;
+ default: sParticle="not known" ;break;
}
- Info("","event=%i hunt=%i tid=%i pid=%i(%s) m=%f q=%3.1f dEdX=%9.3f",
+ Info(Form("Step %i",iStepN),"event=%i hunt=%i tid=%i pid=%i(%s) m=%f q=%3.1f dEdX=%9.3f",
gMC->CurrentEvent(),
fIshunt,
gAlice->GetMCApp()->GetCurrentTrackNumber(),
sParticle,
gMC->TrackMass(),
gMC->TrackCharge(),
- gMC->Edep());
- Info("","Flags:alive(%i) disap(%i) enter(%i) exit(%i) inside(%i) out(%i) stop(%i) new(%i)",
+ gMC->Edep()*1e9);
+ Info("Flags","alive(%i) disap(%i) enter(%i) exit(%i) inside(%i) out(%i) stop(%i) new(%i)",
gMC->IsTrackAlive(),
gMC->IsTrackDisappeared(),
gMC->IsTrackEntering(),
gMC->IsTrackOut(),
gMC->IsTrackStop(),
gMC->IsNewTrack());
- Int_t copy0,copy1,copy2,copy3;
+ Int_t copy0=0,copy1=0,copy2=0,copy3=0;
Int_t vid0=gMC->CurrentVolID(copy0);
Int_t vid1=gMC->CurrentVolOffID(1,copy1);
Int_t vid2=gMC->CurrentVolOffID(2,copy2);
Int_t vid3=gMC->CurrentVolOffID(3,copy3);
- Info("","vid0=%i(%s)c%i vid1=%i(%s)c%i vid2=%i(%s)c%i vid3=%i(%s)c%i %s-%s-%s-%s",
+ Info("Volumes","vid0=%i(%s)c%i vid1=%i(%s)c%i vid2=%i(%s)c%i vid3=%i(%s)c%i %s-%s-%s-%s",
vid0,gMC->VolName(vid0),copy0,
vid1,gMC->VolName(vid1),copy1,
vid2,gMC->VolName(vid2),copy2,
Float_t a,z,den,rad,abs; a=z=den=rad=abs=kBad;
Int_t mid=gMC->CurrentMaterial(a,z,den,rad,abs);
- Info("","mid=%i a=%7.2f z=%7.2f den=%7.2f rad=%7.2f abs=%7.2f",mid,a,z,den,rad,abs);
+ Info("Material","id=%i a=%7.2f z=%7.2f den=%9.4f rad=%9.2f abs=%9.2f",mid,a,z,den,rad,abs);
+
+ Int_t iTmedId=gMC->GetMedium();
+ const char *sTmed;
+ switch(iTmedId){
+ case kAir: sTmed="Air" ;break;
+ case kRoha: sTmed="Rohacell" ;break;
+ case kSiO2: sTmed="SiO2" ;break;
+ case kC6F14: sTmed="C6F14" ;break;
+ case kGridCu: sTmed="GridCu" ;break;
+ case kOpSiO2: sTmed="OpSiO2" ;break;
+ case kGap: sTmed="Gap" ;break;
+ case kGlass: sTmed="Glass" ;break;
+ case kCu: sTmed="Cu" ;break;
+ case kSteel: sTmed="Steel" ;break;
+ case kPerpex: sTmed="Perpex" ;break;
+ case kCH4: sTmed="CH4" ;break;
+ case kCsI: sTmed="CsI" ;break;
+ case kAl: sTmed="Al" ;break;
+ case kW: sTmed="W" ;break;
+ default: sTmed="not known";break;
+ }
+ Info("Medium","id=%i (%s)",iTmedId,sTmed);
TLorentzVector x4;
gMC->TrackPosition(x4);
Float_t glo[3],loc[3];
glo[0]=x4.X();glo[1]=x4.Y();glo[2]=x4.Z();
gMC->Gmtod(glo,loc,1);
- Info("","glo(%+8.3f,%+8.3f,%+8.3f) r=%8.3f theta=%8.3f phi=%8.3f",
+ Info("X4 glo","(x,y,z)=(%+8.3f,%+8.3f,%+8.3f) (r,theta,phi)=(%8.3f,%8.3f,%8.3f)",
glo[0],glo[1],glo[2],x4.Rho(),x4.Theta()*TMath::RadToDeg(),x4.Phi()*TMath::RadToDeg());
- Info("","loc(%+8.3f,%+8.3f,%8.3f) by gMC->Gmtod()", loc[0],loc[1],loc[2]);
- if(gMC->VolId("CSI ")==gMC->CurrentVolID(copy0)){
+ Info("X4 loc","(x,y,z)=(%+8.3f,%+8.3f,%8.3f) by gMC->Gmtod()",loc[0],loc[1],loc[2]);
+ if(gMC->VolId("GAP ")==gMC->CurrentVolID(copy0)){
Int_t iChamber;
gMC->CurrentVolOffID(2,iChamber);
- TVector2 x2=C(iChamber)->Glob2Loc(x4);
- Info("","loc(%+8.3f,%+8.3f) by Glob2Loc", x2.X(),x2.Y());
+ TVector2 x2=C(iChamber)->Mrs2Pc(x4);
+ Info("X4 loc","(%+8.3f,%+8.3f) by Mrs2Pc", x2.X(),x2.Y());
}
- Info("","end of current step\n");
+ Info(Form("Step %i",iStepN++),"end of current step\n");
}//StepManager()
-void AliRICHv0::CreateGeometry()
-{
- if(GetDebug())Info("CreateGeometry","Start v0.");
-
- Double_t cm=1,mm=0.1;
- //place radioactive source compartment if needed
- Double_t containerLen=21.8*mm/2 ,containerR=38*mm/2;
- Double_t screwLen=15*mm/2 ,screwR=5*mm/2;
- Double_t srcLen=10*mm/2 ,srcR=2*mm/2;
- Double_t perpexLen=20*mm/2 ,perpexR=34*mm/2;
- Double_t perpexWholeLen=10*mm/2 ,perpexWholeR=4*mm/2;
- Double_t alBottomLen=containerLen-perpexLen ,alWholeR=5*mm/2;
-//volumes
- Double_t par[3];
-
-// Double_t anodWireD=20*mkm, cathWireD=100*mkm, collWireD=100*mkm;
-
- Double_t pcZ2=0.5*mm;
- par[0]=68.8*cm; par[1]=70.86*cm; par[2]=13*cm; gMC->Gsvolu("RICH","BOX ",(*fIdtmed)[kCH4],par,3); //RICH
- par[0]=P()->PcSizeX()/2; par[1]=P()->PcSizeY()/2; par[2]=pcZ2; gMC->Gsvolu("CSI ","BOX ",(*fIdtmed)[kCSI],par,3); //CSI
- par[0]=P()->PcSizeX()/2; par[1]=P()->PcSizeY()/2; par[2]=P()->GapAmp(); gMC->Gsvolu("GAP ","BOX ",(*fIdtmed)[kCH4],par,3); //GAP
-// par[0]=0; par[1]=; par[2]=P()->PcSizeX()/2;gMC->Gsvolu("WIAN","TUBE",(*fIdtmed)[kW],par,3);//Anod wire
-// par[0]=0; par[1]=srcR; par[2]=P()->PcSizeX()/2;gMC->Gsvolu("WICA","TUBE",(*fIdtmed)[kCu],par,3);//Cathod wire
-// par[0]=0; par[1]=srcR; par[2]=P()->PcSizeX()/2;gMC->Gsvolu("WICO","TUBE",(*fIdtmed)[kCu],par,3);//Collect wire
-
- par[0]=0 ;par[1]=containerR ;par[2]=containerLen; gMC->Gsvolu("CONT","TUBE",(*fIdtmed)[kCH4] ,par,3); //container
- par[0]=perpexR;par[1]=containerR ;par[2]=containerLen; gMC->Gsvolu("ALWA","TUBE",(*fIdtmed)[kAl] ,par,3); //Al cylindric wall
- par[0]=0 ;par[1]=perpexR ;par[2]=alBottomLen; gMC->Gsvolu("ALBO","TUBE",(*fIdtmed)[kAl] ,par,3); //Al bottom
- par[0]=0 ;par[1]=alWholeR ;par[2]=alBottomLen; gMC->Gsvolu("ALWH","TUBE",(*fIdtmed)[kCH4] ,par,3); //Whole in Al bottom
- par[0]=0 ;par[1]=perpexR ;par[2]=perpexLen; gMC->Gsvolu("PEPL","TUBE",(*fIdtmed)[kPerpex],par,3); //Perpex plug
- par[0]=0 ;par[1]=perpexWholeR;par[2]=perpexWholeLen;gMC->Gsvolu("PEWH","TUBE",(*fIdtmed)[kCH4] ,par,3); //Whole in Perpex
- par[0]=0 ;par[1]=screwR ;par[2]=screwLen; gMC->Gsvolu("SCRE","TUBE",(*fIdtmed)[kSteel] ,par,3); //Screw
- par[0]=0 ;par[1]=srcR ;par[2]=srcLen; gMC->Gsvolu("SOUR","TUBE",(*fIdtmed)[kSteel] ,par,3); //Source itself
-//nodes
- gMC->Gspos("RICH",1,"ALIC",0,0,-9 ,0,"ONLY"); //RICH in ALIC
- gMC->Gspos("CSI ",1,"RICH",0,0,-(P()->GapProx()+pcZ2) ,0,"ONLY"); //CsI in ALIC
- gMC->Gspos("GAP ",1,"RICH",0,0,-(P()->GapProx()-P()->GapAmp()/2),0,"ONLY"); //GAP in ALIC
-
- gMC->Gspos("CONT",1,"RICH",0,0,1*cm ,0,"ONLY"); //Sr90 container in RICH
- gMC->Gspos("ALWA",1,"CONT",0,0,0 ,0,"ONLY"); //Al wall
- gMC->Gspos("ALBO",1,"CONT",0,0,-containerLen+alBottomLen ,0,"ONLY"); //Al bottom
- gMC->Gspos("PEPL",1,"CONT",0,0,containerLen-perpexLen ,0,"ONLY"); //Perpex plug
-
- gMC->Gspos("ALWH",1,"ALBO",6*mm,0,0 ,0,"ONLY"); //Whole in Al bottom
-
- gMC->Gspos("PEWH",1,"PEPL",6*mm,0,-perpexLen+perpexWholeLen ,0,"ONLY"); //Whole in Perpex plug
- gMC->Gspos("SOUR",1,"PEPL",6*mm,0, perpexLen-srcLen ,0,"ONLY"); //Source in Perpex plug
- gMC->Gspos("SCRE",1,"PEPL",0,0,perpexLen-screwLen ,0,"ONLY"); //Screw in Perpex plug
- if(GetDebug())Info("CreateGeometry","Stop v0.");
-}//CreateGeometry()
virtual void Init() {;} //interface from AliRICH
virtual Int_t IsVersion() const{return 0;} //interface from AliRICH
virtual void StepManager(); //interface from AliRICH
- virtual void CreateGeometry(); //interface from AliRICH
protected:
ClassDef(AliRICHv0,1) //RICH coarse version for material budget study and debuging
};
void AliRICHv1::StepManager()
{
// Full Step Manager.
-// 3- Ckovs absorbed in Collection electrods
-// 5- Ckovs absorbed in Cathode wires
-// 6- Ckovs absorbed in Anod wires
+// 3- Ckovs absorbed on Collection electrods
+// 5- Ckovs absorbed on Cathode wires
+// 6- Ckovs absorbed on Anod wires
Int_t copy;
static Int_t iCurrentChamber;
//history of Cerenkovs
if(gMC->TrackPid()==kCerenkov){
- if(gMC->IsNewTrack() && (gMC->CurrentVolID(copy)==gMC->VolId("FRE1")||gMC->CurrentVolID(copy)==gMC->VolId("FRE2"))) fCounters(0)++;// 0- Ckovs produced in Freon
- if(!gMC->IsTrackAlive() && (gMC->CurrentVolID(copy)==gMC->VolId("FRE1")||gMC->CurrentVolID(copy)==gMC->VolId("FRE2"))) fCounters(1)++;// 1- Ckovs absorbed in Freon
- if(!gMC->IsTrackAlive() && gMC->CurrentVolID(copy)==gMC->VolId("QUAR")) fCounters(2)++;//2- Ckovs absorbed in Quartz
- if(!gMC->IsTrackAlive() && gMC->CurrentVolID(copy)==gMC->VolId("META")) fCounters(4)++;//4- Ckovs absorbed in CH4
+ if( gMC->IsNewTrack() && gMC->CurrentVolID(copy)==gMC->VolId("RRAD")) fCounters(0)++;// 0- Ckovs produced in radiator
+ if(!gMC->IsTrackAlive() && gMC->CurrentVolID(copy)==gMC->VolId("RRAD")) fCounters(1)++;// 1- Ckovs absorbed in radiator
+ if(!gMC->IsTrackAlive() && gMC->CurrentVolID(copy)==gMC->VolId("RRWI")) fCounters(2)++;// 2- Ckovs absorbed in radiator window
+ if(!gMC->IsTrackAlive() && gMC->CurrentVolID(copy)==gMC->VolId("RICH")) fCounters(4)++;// 4- Ckovs absorbed in CH4
}
//Treat photons
static TLorentzVector cerX4;
- if((gMC->TrackPid()==kCerenkov||gMC->TrackPid()==kFeedback)&&gMC->CurrentVolID(copy)==gMC->VolId("CSI ")){//photon in CSI
- if(gMC->Edep()>0){//CF+CSI+DE
+ if((gMC->TrackPid()==kCerenkov||gMC->TrackPid()==kFeedback)&&gMC->CurrentVolID(copy)==gMC->VolId("RPC ")){//photon in PC
+ if(gMC->Edep()>0){//photon in PC +DE
if(IsLostByFresnel()){
- if(gMC->TrackPid()==kCerenkov) fCounters(7)++;// 7- Ckovs reflected on CsI
+ if(gMC->TrackPid()==kCerenkov) fCounters(7)++;// 7- Ckovs reflected from CsI
gMC->StopTrack();
return;
}
- gMC->TrackPosition(cerX4); gMC->CurrentVolOffID(2,iCurrentChamber);
+ gMC->TrackPosition(cerX4); gMC->CurrentVolOffID(2,iCurrentChamber);//RICH-RPPF-RPC
AddHit(iCurrentChamber,gAlice->GetMCApp()->GetCurrentTrackNumber(),cerX4.Vect(),cerX4.Vect());//HIT for PHOTON in conditions CF+CSI+DE
- fCounters(8)++;//4- Ckovs converted in CsI
+ fCounters(8)++;//4- Ckovs converted to electron on CsI
GenerateFeedbacks(iCurrentChamber);
- }//CF+CSI+DE
- }//CF in CSI
+ }//photon in PC and DE >0
+ }//photon in PC
//Treat charged particles
static Float_t eloss;
static TLorentzVector mipInX4,mipOutX4;
- if(gMC->TrackCharge() && gMC->CurrentVolID(copy)==gMC->VolId("GAP ")){//MIP in GAP
- gMC->CurrentVolOffID(3,iCurrentChamber);
+ if(gMC->TrackCharge() && gMC->CurrentVolID(copy)==gMC->VolId("RGAP")){//MIP in GAP
+ gMC->CurrentVolOffID(1,iCurrentChamber);//RICH-RGAP
if(gMC->IsTrackEntering()||gMC->IsNewTrack()) {//MIP in GAP entering or newly created
eloss=0;
gMC->TrackPosition(mipInX4);
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)){
- if(GetDebug()) Info("IsLostByFresnel","Photon lost");
+ AliDebug(1,"Photon lost");
return kTRUE;
}else
return kFALSE;
TLorentzVector x4;
gMC->TrackPosition(x4);
- TVector2 x2=C(iChamber)->Glob2Loc(x4);
- Int_t sector=P()->Loc2Sec(x2); if(sector==kBad) return; //hit in dead zone, nothing to produce
+ TVector2 x2=C(iChamber)->Mrs2Pc(x4);//hit position on photocathode plane
+ TVector2 xspe=x2;
+ Int_t sector=P()->Loc2Sec(xspe); if(sector==kBad) return; //hit in dead zone, nothing to produce
Int_t iTotQdc=P()->TotQdc(x2,eloss);
- Int_t iNphotons=gMC->GetRandom()->Poisson(P()->AlphaFeedback(sector)*iTotQdc);
- if(GetDebug())Info("GenerateFeedbacks","N photons=%i",iNphotons);
+ Int_t iNphotons=gMC->GetRandom()->Poisson(P()->C(iChamber)->AlphaFeedback(sector)*iTotQdc);
+ AliDebug(1,Form("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
outputNtracksStored,
1.0);
}//feedbacks loop
- if(GetDebug()) Info("GenerateFeedbacks","Stop.");
+ AliDebug(1,"Stop.");
}//GenerateFeedbacks()
kPMD=8192,kDIPO=16384,kEMCAL=32768,kVZERO=65536,kMUON=131072,kZDC=262144,kSHILD=524288};
enum EProcesses {kDCAY=1,kPAIR=2,kCOMP=4,kPHOT=8,kPFIS=16,kDRAY=32,kANNI=64,kBREM=128,kMUNU=256,kCKOV=512,kHADR=1024,kLOSS=2048,kMULS=4096,
kRAYL=8192};
- enum EGenTypes {kGun1,kGun7,kPythia7,kHijing,kHijingPara};
+ enum EGenTypes {kGun1=100,kGun7,kPythia7,kHijing,kHijingPara};
KirConfig(const char*sFileName);
~KirConfig() {Info("ctor","");fMain->Cleanup(); delete fMain;}
//Generator
pVerFrame->AddFrame(pGenGrpFrm=new TGGroupFrame(pHorFrame,"Generator"));
pGenGrpFrm->AddFrame(fGenTypeCombo=new TGComboBox(pGenGrpFrm,100));
- fGenTypeCombo->AddEntry("gun to single chamber",kGun1); fGenTypeCombo->AddEntry("gun to all chambers",kGun7);
+ fGenTypeCombo->AddEntry("gun to single chamber",kGun1);
+ fGenTypeCombo->AddEntry("gun to all chambers",kGun7);
fGenTypeCombo->AddEntry("7 guns on top of Pythia",kPythia7);
- fGenTypeCombo->AddEntry("HIJING",kHijing); fGenTypeCombo->AddEntry("parametrized HIJING",kHijingPara);
+ fGenTypeCombo->AddEntry("HIJING",kHijing);
+ fGenTypeCombo->AddEntry("parametrized HIJING",kHijingPara);
+ fGenTypeCombo->AddEntry("Sr90 source",kSr90);
fGenTypeCombo->Select(kHijingPara);
fGenTypeCombo->Resize(160,20);
// Magnetic Field
pVerFrame->AddFrame(pFldGrpFrm=new TGGroupFrame(pHorFrame,"Magnetic Field"));
pFldGrpFrm->AddFrame(fMagFldChkBtn=new TGCheckButton(pFldGrpFrm,"On/Off"));
- fMagFldChkBtn->SetState(kButtonUp);
+ fMagFldChkBtn->SetState(kButtonDown);
//Detectors
pHorFrame->AddFrame(fDetButGrp=new TGButtonGroup(pHorFrame,"Detectors"));
fDetButGrp->Connect("Pressed(Int_t)","KirConfig",this,"AddDetector(Int_t)");
switch(fGenTypeCombo->GetSelected()){
case kHijingPara:
fprintf(fp," AliGenHIJINGpara *pGen=new AliGenHIJINGpara(%i);\n",(int)fGenNprimEntry->GetNumber());
- fprintf(fp," pGen->SetMomentumRange(0,999); pGen->SetPhiRange(0,360); pGen->SetThetaRange(%f,%f);\n",Eta2Theta(8),Eta2Theta(-8));
+ fprintf(fp," pGen->SetMomentumRange(0,999); pGen->SetThetaRange(%f,%f); pGen->SetPhiRange(0,360);\n",Eta2Theta(8),Eta2Theta(-8));
fprintf(fp," pGen->SetOrigin(0,0,0); pGen->SetSigma(0,0,0);\n");
fprintf(fp," pGen->Init();\n");
break;
case kGun1:
fprintf(fp," AliGenFixed *pGen=new AliGenFixed(1);\n");
- fprintf(fp," pGen->SetPart(%i); pGen->SetMomentum(%3.1f); pGen->SetOrigin(0,0,0);\n",fGenPartIdCombo->GetSelected(),float(fGenMomCombo->GetSelected())/10);
- fprintf(fp," pGen->SetPhiRange(pRICH->C(%i)->PhiD()); pGen->SetThetaRange(pRICH->C(%i)->ThetaD()-2);\n",fGenChamberCombo->GetSelected(),fGenChamberCombo->GetSelected());
+ fprintf(fp," pGen->SetPart(%i); pGen->SetOrigin(0,0,0);\n",fGenPartIdCombo->GetSelected());
+ fprintf(fp," pGen->SetMomentum(%3.1f); pGen->SetTheta(pRICH->C(%i)->ThetaD()-2); pGen->SetPhi(pRICH->C(%i)->PhiD()); \n",
+ float(fGenMomCombo->GetSelected())/10,fGenChamberCombo->GetSelected(),fGenChamberCombo->GetSelected());
fprintf(fp," pGen->Init();\n");
break;
case kGun7:
fprintf(fp," for(int i=1;i<=7;i++){\n");
fprintf(fp," AliGenFixed *pFixed=new AliGenFixed(1);\n");
fprintf(fp," pFixed->SetPart(%i); pFixed->SetMomentum(2.5+i*0.4); pFixed->SetOrigin(0,0,0);\n",fGenPartIdCombo->GetSelected());
- fprintf(fp," pFixed->SetPhiRange(pRICH->C(i)->PhiD()); pFixed->SetThetaRange(pRICH->C(i)->ThetaD()-2);\n");
+ fprintf(fp," pFixed->SetPhi(pRICH->C(i)->PhiD()); pFixed->SetTheta(pRICH->C(i)->ThetaD()-2);\n");
fprintf(fp," pCocktail->AddGenerator(pFixed,Form(\"Fixed %i\",i),1);\n }\n");
fprintf(fp," pCocktail->Init();\n");
break;
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," pGen->SetJetQuenching(0); pGen->SetShadowing(0);\n");
+ fprintf(fp," pGen->KeepFullEvent(); pGen->SetSelectAll(0); \n");
+ fprintf(fp," pGen->SetImpactParameterRange(0., 5.);\n");
+ fprintf(fp," pGen->Init();\n");
break;
}
//central before RICH detectors
Double_t dx,dy,dz,r1,r2,a=0,z=0,den=0,radlen=0,intlen=0;
Double_t cm=1,m=100*cm,mm=0.1*cm,mkm=0.001*cm;
//Medium
- pAir=new TGeoMedium("Air" ,1, g->GetMaterial("Air"));
- new TGeoMaterial("CH4" ,a=26.98 ,z=13,den=0.4224 ); pCH4=new TGeoMedium("CH4" ,2, g->GetMaterial("CH4"));
- new TGeoMaterial("CsI" ,a=26.98 ,z=13,den=2.7 ,radlen=24.01,intlen=70.6); pCsI=new TGeoMedium("CsI" ,3, g->GetMaterial("CsI"));
- new TGeoMaterial("Al" ,a=26.98 ,z=13,den=2.7 ,radlen=24.01,intlen=70.6); pAl =new TGeoMedium("Al" ,4, g->GetMaterial("Al"));
- new TGeoMaterial("Fe" ,a=55.845,z=26,den=7.87,radlen=13.84,intlen=82.8); pFe =new TGeoMedium("Fe" ,5, g->GetMaterial("Fe"));
- new TGeoMaterial("W" ,a=183.84,z=27,den=19.3,radlen= 9.59*cm,intlen=0.35*cm);pW =new TGeoMedium("W" ,6, g->GetMaterial("W"));
- new TGeoMaterial("Cu" ,a=55.845,z=26,den=7.87,radlen=13.84,intlen=82.8); pCu =new TGeoMedium("Cu" ,7, g->GetMaterial("Cu"));
- new TGeoMaterial("Perpex" ,a=55.845,z=26,den=7.87,radlen=13.84,intlen=82.8); new TGeoMedium("Perpex" ,8, g->GetMaterial("Perpex"));
- new TGeoMaterial("Rohacell" ,a=12.01,z=6,den=0.1,radlen=18.8,intlen=0); new TGeoMedium("Rohacell" ,9, g->GetMaterial("Rohacell"));
- new TGeoMaterial("SiO2" ,a=0,z=0,den=0); new TGeoMedium("SiO2" ,10, g->GetMaterial("SiO2"));
- new TGeoMaterial("C6F14" ,a=0,z=0,den=0); new TGeoMedium("C6F14" ,11, g->GetMaterial("C6F14"));
+ pAir =new TGeoMedium("Air" ,1, g->GetMaterial("Air"));
+ new TGeoMaterial("CH4" ,a=26.98 ,z=13,den=0.4224 ); pCH4 =new TGeoMedium("CH4" ,2, g->GetMaterial("CH4"));
+ new TGeoMaterial("CsI" ,a=26.98 ,z=13,den=2.7 ,radlen=24.01*cm,intlen=70.6*cm); pCsI =new TGeoMedium("CsI" ,3, g->GetMaterial("CsI"));
+ new TGeoMaterial("Al" ,a=26.98 ,z=13,den=2.7 ,radlen=24.01*cm,intlen=70.6*cm); pAl =new TGeoMedium("Al" ,4, g->GetMaterial("Al"));
+ new TGeoMaterial("W" ,a=183.84,z=27,den=19.3,radlen= 9.59*cm,intlen=0.35*cm); pW =new TGeoMedium("W" ,5, g->GetMaterial("W"));
+ new TGeoMaterial("Cu" ,a=55.845,z=26,den=7.87,radlen=13.84*cm,intlen=82.8*cm); pCu =new TGeoMedium("Cu" ,6, g->GetMaterial("Cu"));
+ new TGeoMaterial("Rohacell",a=12.01 ,z=6 ,den=0.1,radlen=18.8,intlen=0); pRoha =new TGeoMedium("Rohacell" ,7, g->GetMaterial("Rohacell"));
+ new TGeoMaterial("SiO2" ,a=0 ,z=0 ,den=0); pSiO2= new TGeoMedium("SiO2" ,8,g->GetMaterial("SiO2"));
+ new TGeoMaterial("C6F14" ,a=0 ,z=0 ,den=0); pC6F14=new TGeoMedium("C6F14" ,9,g->GetMaterial("C6F14"));
+//for Sr90 source
+ new TGeoMaterial("Perpex" ,a=55.845,z=26,den=7.87,radlen=13.84*cm,intlen=82.8*cm); pPerpex=new TGeoMedium("Perpex",10, g->GetMaterial("Perpex"));
+ new TGeoMaterial("Steel" ,a=55.845,z=26,den=7.87,radlen=13.84*cm,intlen=82.8*cm); pSteel=new TGeoMedium("Steel" ,11, g->GetMaterial("Steel"));
+ new TGeoMaterial("Mylar" ,a=55.845,z=26,den=7.87,radlen=13.84*cm,intlen=82.8*cm); pMylar=new TGeoMedium("Mylar" ,12, g->GetMaterial("Mylar"));
+ new TGeoMaterial("Sr" ,a=87.62 ,z=38,den=7.87,radlen=13.84*cm,intlen=82.8*cm); pSr =new TGeoMedium("Sr" ,13, g->GetMaterial("Sr"));
//Mother
TGeoVolume *pMother=g->MakeBox("Mother",pAir,dx=10*m/2,dy=10*m/2,dz=30*m/2); //arbitrary values
//Rich
pRich->AddNode(pPpf,copy=5,new TGeoTranslation(-335*mm,+433*mm,8*cm+20*mm));
pRich->AddNode(pPpf,copy=6,new TGeoTranslation(+335*mm,+433*mm,8*cm+20*mm));
pPpf->AddNode( pPc ,copy=1,new TGeoTranslation( 0*mm, 0*mm,-19.15*mm));//PPF 2001P2
- pPpf->AddNode(g->GetVolume("PPFlarge"),copy=1,new TGeoTranslation(-224.5*mm,-151.875*mm, 0.85*mm));
+ pPpf->AddNode(pPpfLarge,copy=1,new TGeoTranslation(-224.5*mm,-151.875*mm, 0.85*mm));
pPpf->AddNode(g->GetVolume("PPFlarge"),copy=2,new TGeoTranslation(-224.5*mm,- 50.625*mm, 0.85*mm));
pPpf->AddNode(g->GetVolume("PPFlarge"),copy=3,new TGeoTranslation(-224.5*mm,+ 50.625*mm, 0.85*mm));
pPpf->AddNode(g->GetVolume("PPFlarge"),copy=4,new TGeoTranslation(-224.5*mm,+151.875*mm, 0.85*mm));
TGeoVolume *pRadShort =g->MakeBox( "RadShort" ,g->GetMedium("Neoceram") ,dx= 10*mm/2 ,dy= 403*mm/2 ,dz= 15*mm/2);
TGeoVolume *pRadSpacer=g->MakeTube("RadSpacer",g->GetMedium("SiO2") ,r1= 0 ,r2=10*mm/2 ,dz= 15*mm/2);
- pRich->AddNode(pRad,copy=1,new TGeoTranslation( 0*mm,-434*mm, -12*mm));
- pRich->AddNode(pRad,copy=2,new TGeoTranslation( 0*mm, 0*mm, -12*mm));
- pRich->AddNode(pRad,copy=3,new TGeoTranslation( 0*mm,+434*mm, -12*mm));
- pRad->AddNode(g->GetVolume("RadFront") ,copy=1,new TGeoTranslation( 0*mm, 0*mm, -10.0*mm));
- pRad->AddNode(g->GetVolume("RadWin") ,copy=1,new TGeoTranslation( 0*mm, 0*mm, 9.5*mm));
- pRad->AddNode(g->GetVolume("RadLong") ,copy=1,new TGeoTranslation( 0*mm,-204*mm, -0.5*mm));
- pRad->AddNode(g->GetVolume("RadLong") ,copy=2,new TGeoTranslation( 0*mm,+204*mm, -0.5*mm));
- pRad->AddNode(g->GetVolume("RadShort") ,copy=1,new TGeoTranslation(-660*mm, 0*mm, -0.5*mm));
- pRad->AddNode(g->GetVolume("RadShort") ,copy=2,new TGeoTranslation(+660*mm, 0*mm, -0.5*mm));
+ pRich->AddNode(pRad ,copy=1,new TGeoTranslation( 0*mm,-434*mm, -12*mm));
+ pRich->AddNode(pRad ,copy=2,new TGeoTranslation( 0*mm, 0*mm, -12*mm));
+ pRich->AddNode(pRad ,copy=3,new TGeoTranslation( 0*mm,+434*mm, -12*mm));
+
+ pRad ->AddNode(pRadFront ,copy=1,new TGeoTranslation( 0*mm, 0*mm, -10.0*mm));
+ pRad ->AddNode(pRadWin ,copy=1,new TGeoTranslation( 0*mm, 0*mm, 9.5*mm));
+ pRad ->AddNode(pRadLong ,copy=1,new TGeoTranslation( 0*mm,-204*mm, -0.5*mm));
+ pRad ->AddNode(pRadLong ,copy=2,new TGeoTranslation( 0*mm,+204*mm, -0.5*mm));
+ pRad ->AddNode(pRadShort ,copy=1,new TGeoTranslation(-660*mm, 0*mm, -0.5*mm));
+ pRad ->AddNode(pRadShort ,copy=2,new TGeoTranslation(+660*mm, 0*mm, -0.5*mm));
for(int i=0;i<3;i++)
for(int j=0;j<10;j++)
pRad->AddNode(pRadSpacer,copy=10*i+j,new TGeoTranslation(-1330*mm/2+116*mm+j*122*mm,(i-1)*105*mm,-0.5*mm));
}else{//Sr90
- Double_t containerZ=21.8*mm/2,containerR=38*mm/2;
- Double_t screwLen=15*mm/2 ,screwR=5*mm/2;
- Double_t srcLen=10*mm/2 ,srcR=2*mm/2;
- Double_t perpexZ=20*mm/2 ,perpexR=34*mm/2;
- Double_t perpexWholeLen=10*mm/2,perpexWholeR=4*mm/2;
- Double_t alBottomLen=containerZ-perpexZ,alWholeR=5*mm/2;
- g->MakeTube("Sr90" ,pAir ,0 ,containerR ,containerZ);
- g->MakeTube("Sr90AlWall" ,pAl ,perpexR,containerR ,containerZ);
- g->MakeTube("Sr90AlBottom" ,pAl ,0 ,perpexR ,alBottomLen);
- g->MakeTube("Sr90AlWhole" ,pAir ,0 ,alWholeR ,alBottomLen);
- g->MakeTube("Sr90PerpexPlug" ,g->GetMedium("Perpex"),0 ,perpexR ,perpexZ);
- g->MakeTube("Sr90PerpexWhole",pAir ,0 ,perpexWholeR,perpexWholeLen);
- g->MakeTube("Sr90Screw" ,g->GetMedium("Steel") ,0 ,screwR ,screwLen);
- g->MakeTube("Sr90Source" ,g->GetMedium("Steel") ,0 ,srcR ,srcLen);
+ pSrc =g->MakeTube("Src" ,pCH4 , 0 , 70*mm/2 , 30*mm/2); //top container
+ pAlGlass =g->MakeTube("SrcAlGlass" ,pAl , 0 , 38*mm/2 ,21.8*mm/2); //Al glass wall
+ pPerpexPlug =g->MakeTube("SrcPerpex" ,pPerpex , 0 , 34*mm/2 , 20*mm/2); //Perpex plug
+ pSteelScrew =g->MakeTube("SrcSteelScrew" ,pSteel , 0 , 5*mm/2 , 15*mm/2); //Steel screw in the center
+ pSr90Screw =g->MakeTube("SrcSr90Screw" ,pSteel , 0 , 2*mm/2 , 10*mm/2); //Steel screw to support Sr90
+ pSr90 =g->MakeTube("SrcSr90" ,pSr , 0 , 1*mm/2 , 1*mm/2); //Sr90 source
+ pHolePerpex =g->MakeTube("SrcHolePerpex" ,pAir , 0 , 4*mm/2 , 10*mm/2); //Air hole in perpex plug
+ pHoleAl =g->MakeTube("SrcHoleAl" ,pAir , 0 , 4*mm/2 , 1.8*mm/2); //Air hole in Al glass bottom
+ pMylarFoil =g->MakeTube("SrcMylarFoil" ,pMylar , 0 , 30*mm/2 , 50*mkm/2); //Mylar foil
- pRich->AddNode(pContainer,1,new TGeoTranslation(30*cm,0,containerZ));
- pContainer ->AddNode(pAlWall,1);
- pContainer ->AddNode(pAlBottom,1,new TGeoTranslation(0,0,-containerZ+alBottomLen));
- pContainer ->AddNode(pPerpexPlug,1,new TGeoTranslation(0,0,containerZ-perpexZ));
- pAlBottom ->AddNode(pAlWhole,1,new TGeoTranslation(6*mm,0,0));
- pPerpexPlug->AddNode(pPerpexWhole,1,new TGeoTranslation(6*mm,0,-perpexZ+perpexWholeLen));
- pPerpexPlug->AddNode(pSource,1 ,new TGeoTranslation(6*mm,0, perpexZ-srcLen));
- pPerpexPlug->AddNode(pScrew,1,new TGeoTranslation(0,0,perpexZ-screwLen));
+ pRich->AddNode(pSrc,1,new TGeoTranslation(30*cm,0,alVesselZ));
+ pSrc ->AddNode(pMylarFoil,1,new TGeoTranslation(0,0,21.8*mm/2+50*mkm/2));
+ pSrc ->AddNode(pAlGlass,1,new TGeoTranslation(0,0,0));//Al glass to fake Src volume
+ pAlGlass->AddNode(pHoleAl,1,new TGeoTranslation(6*mm, 0*mm, 10.9*mm));
+ pAlGlass->AddNode(pPerpexPlug,1,new TGeoTranslation(0*mm,0*mm,-0.9*mm));
+ pPerpexPlug->AddNode(pSteelScrew,1,new TGeoTranslation(shift,0, perpexPlugZ-steelScrewZ));
+ pPerpexPlug->AddNode(pHolePerpex,1,new TGeoTranslation(shift,0,-perpexPlugZ+airWholeZ));
+ pPerpexPlug->AddNode(pSr90Screw ,1,new TGeoTranslation(shift,0, perpexPlugZ-sr90ScrewZ));
+ pSr90Screw->AddNode(pSr90 ,1,new TGeoTranslation(0,0, -(sr90ScrewZ-sr90Z)));
}//if(!isSr90)
//Sandbox
- TGeoVolume *pSandBox=g->MakeBox( "SandBox" ,g->GetMedium("Air") ,dx=1419*mm/2 ,dy=1378*mm/2 ,dz=50.5*mm/2); //2072P1
- g->MakeBox( "SandCover",g->GetMedium("Al") ,dx=1419*mm/2 ,dy=1378*mm/2 ,dz= 0.5*mm/2);
- g->MakeBox( "SandComb" ,g->GetMedium("Rohacell") ,dx=1359*mm/2 ,dy=1318*mm/2 ,dz=49.5*mm/2);
+ TGeoVolume *pSandBox =g->MakeBox( "SandBox" ,pAir ,dx=1419*mm/2 ,dy=1378*mm/2 ,dz=50.5*mm/2); //2072P1
+ TGeoVolume *pSandCover=g->MakeBox( "SandCover",pAl ,dx=1419*mm/2 ,dy=1378*mm/2 ,dz= 0.5*mm/2);
+ TGeoVolume *pSandComb =g->MakeBox( "SandComb" ,pRoha ,dx=1359*mm/2 ,dy=1318*mm/2 ,dz=49.5*mm/2);
pRich->AddNode(pSandBox,copy=1,new TGeoTranslation( 0*mm,0*mm, -73.75*mm)); //p.84 TDR
- pSandBox->AddNode(g->GetVolume("SandComb") ,copy=1,new TGeoTranslation( 0*mm,0*mm, 0*mm)); //2072P1
- pSandBox->AddNode(g->GetVolume("SandCover") ,copy=1,new TGeoTranslation( 0*mm,0*mm, +25*mm));
- pSandBox->AddNode(g->GetVolume("SandCover") ,copy=2,new TGeoTranslation( 0*mm,0*mm, -25*mm));
+ pSandBox->AddNode(pSandComb ,copy=1,new TGeoTranslation( 0*mm,0*mm, 0*mm)); //2072P1
+ pSandBox->AddNode(pSandCover ,copy=1,new TGeoTranslation( 0*mm,0*mm, +25*mm));
+ pSandBox->AddNode(pSandCover ,copy=2,new TGeoTranslation( 0*mm,0*mm, -25*mm));
//Colors
g->GetVolume("Pc") ->SetLineColor(kGreen);
g->GetVolume("PPF") ->SetLineColor(33);
Float_t aIdxSiO2[kNbins];
Float_t aIdxOpSiO2[kNbins];
- Float_t e1=10.666,e2=18.125,f1=46.411,f2= 228.71; //RICH TDR page 35
for (i=0;i<kNbins;i++){
- aIdxC6F14[i] = aPckov[i] *0.0172*1e9+1.177;
- aIdxSiO2[i] = TMath::Sqrt(1.+f1/(e1*e1-TMath::Power(aPckov[i]*1e9,2))+f2/(e2*e2-TMath::Power(aPckov[i]*1e9,2)));//TDR p.35
- aIdxCH4[i] =1.000444;
- aIdxGrid[i] =1;
+ aIdxC6F14[i] = (Float_t)AliRICHParam::IndOfRefC6F14(aPckov[i]*1e9);
+ aIdxSiO2[i] = (Float_t)AliRICHParam::IndOfRefSiO2(aPckov[i]*1e9);
+ aIdxCH4[i] = (Float_t)AliRICHParam::IndOfRefCH4();
+ aIdxGrid[i] =1;
aIdxOpSiO2[i] =1;
}
#pragma link C++ class AliRICHDigitizer+;
#pragma link C++ class AliRICHDisplFast+;
#pragma link C++ class AliRICHReconstructor+;
+#pragma link C++ class AliRICHHelix+;
+#pragma link C++ class AliGenRadioactive+;
+#pragma link C++ class AliRICHTracker+;
+#pragma link C++ enum RadioNuclei_t;
#endif
void RichBatch(const Int_t iNevents,const Bool_t isDebug,const char *sConfigFileName)
{
- if(isDebug) gAlice->SetDebug(1);
+ gSystem->Exec("rm -rf *.root hlt hough ZZZ*");
+ if(isDebug) AliLog::SetGlobalDebugLevel(AliLog::kDebug);
- Info("my/RichBatch.C","%i event(s) requested, debug %i,config file %s",iNevents,isDebug,sConfigFileName);
TStopwatch sw;TDatime time;
-
+
+
AliSimulation *pSim=new AliSimulation; pSim->Run(iNevents); delete pSim;
- AliReconstruction *pRec=new AliReconstruction; pRec->SetRunLocalReconstruction("RICH"); pRec->SetFillESD("RICH"); pRec->Run(); delete pRec;
+ AliReconstruction *pRec=new AliReconstruction;
+ pRec->SetRunLocalReconstruction("ITS TPC TRD TOF RICH");
+ pRec->SetFillESD("ITS TPC TRD TOF RICH");
+ pRec->Run(); delete pRec;
+
+ //pRec->SetRunLocalReconstruction("RICH"); pRec->SetFillESD("RICH");
+
+ cout<<"\n!!!!!!!!!!!!Info in <my/RichBatch.C>: Start creating Control Plots \n";
- AliRunLoader* pAL = AliRunLoader::Open();
- pAL->LoadgAlice();
+ AliRunLoader* pAL = AliRunLoader::Open(); pAL->LoadgAlice();
AliRICH *pRICH=(AliRICH*)pAL->GetAliRun()->GetDetector("RICH");
if(pRICH) pRICH->ControlPlots();
- cout<<"\nInfo in <my/RichBatch.C>: Start time: ";time.Print();
- cout<<"Info in <my/RichBatch.C>: Stop time: ";time.Set(); time.Print();
- cout<<"Info in <my/RichBatch.C>: Time used: ";sw.Print();
+ cout<<"\n!!!!!!!!!!!!Info in <my/RichBatch.C>: Start time: ";time.Print();
+ cout<<"!!!!!!!!!!!!Info in <my/RichBatch.C>: Stop time: ";time.Set(); time.Print();
+ cout<<"!!!!!!!!!!!!Info in <my/RichBatch.C>: Time used: ";sw.Print();
gSystem->Exec("touch ZZZ______finished_______ZZZ");
+// gSystem->Exec("playwave .kde/my/end.wav");
}
How to retrive pointer to alice run loader:
use pRICH->GetLoader()->GetRunLoader() (all detector classes inherit from AliDetector wich has GetLoader())
use methode AliRun::GetRunLoader for gAlice (depricated)
-How to get pointers to deifferent root trees:
+How to get pointers to different root trees:
TreeE belongs to AliRunLoader, available after AliRunLoader::LoadHeader()
TreeK belongs to AliRunLoader, available after AliRunLoader::LoadKinematics()
TreeH belongs to AliLoader , available after AliLoader::LoadHits()
How to work with the stack of particles?
- - pointer to the stack is returned by gAlice->Stack() (global gAlice of type AliRun) or AliRunLoader::Stack() but
- before one needs to load event header by AliRunLoader::LoadHeader() otherwise both methods return 0.
- Moreover loading header gives the information about number of particles only.
- To retrive the list of particle one also needs to load kinematics by AliRunLoader::LoadKinematics()
- - total amount of particles in stack for a given event: AliRunLoader::Stack()->GetNtrack() or AliRun::GetEvent() (after LoadHeader())
- - total amount of primiry particles in stack for a given event: AliRunLoader::Stack()->GetNprimary() or AliLoader::TreeH()->GetEntries() (after LoadHeader())
+ - first of all, the stack includes primary as well as secondary particles
+ - pointer to the stack is returned by AliRun::Stack() (global gAlice of type AliRun - depricated way to do)
+ or by AliRunLoader::Stack() but before one needs to load event header by AliRunLoader::LoadHeader() otherwise both methods return 0.
+ Moreover loading header gives the information about number of particles only.
+ To retrive the list of particle one also needs to load kinematics by AliRunLoader::LoadKinematics()
+ - total amount of particles in stack for a given event: AliStack::GetNtrack() or AliRun::GetEvent() (after LoadHeader())
+ - total amount of primiry particles in stack for a given event: AliStack::GetNprimary() or AliLoader::TreeH()->GetEntries() (after LoadHeader())
How to retrive hits:
- Hits a stored on primiry by primiry basis. To retrieve all hits one needs to do:
+ Hits a stored on primary by primary basis. To retrieve all hits one needs to do:
initialise the root tree and containers: AliLoader::LoadHits()
- read number of primiries in current event:
+ read number of primaries in current event:
loop on the list of primiries:
One needs to say:
pRich->GetLoader()->LoadSDigits(); this one open file, get the tree and invoke AliRICH::SetTreeAddress()
-
-
-gAlice->GetMCApp()->GetCurrentTrackNumber()
-
-
What are the debug methodes avail:
- AliModule::GetDebug()
+ AliModule::GetDebug() now depricated, use AliDebug printout instead
AliModule::SetDebug()
AliRun::GetDebug()
AliRun::SetDebug()
How to avoid using gAlice?
Rich()->GetLoader()->GetRunLoader()->GetAliRun() returns gAlice global pointer.
-
-
AliRICHClusterFinder.cxx AliRICHMap.cxx\
AliRICHChamber.cxx AliRICHRecon.cxx\
AliRICHDisplFast.cxx\
- AliRICHDigitizer.cxx AliRICHReconstructor.cxx
+ AliRICHDigitizer.cxx AliGenRadioactive.cxx AliRICHHelix.cxx AliRICHTracker.cxx AliRICHReconstructor.cxx
HDRS = $(SRCS:.cxx=.h)
DHDR= RICHLinkDef.h
void ps(Int_t event=0) {R()->PrintSDigits(event);} //utility print sdigits
void pd(Int_t event=0) {R()->PrintDigits(event);} //utility print digits
void pc(Int_t event=0) {R()->PrintClusters(event);}//utility print clusters
+void pt(Int_t event=0) {R()->PrintTracks(event);} //utility print tracks
//__________________________________________________________________________________________________
void pp(int tid)
//__________________________________________________________________________________________________
void PrintParticleInfo(int tid)
{
+// Prints particle info for a given TID
TParticle *p=al->Stack()->Particle(tid);
cout<<p->GetName()<<"("<<tid<<")";
if(p->GetMother(0)!=-1){cout<<" from "; PrintParticleInfo(p->GetMother(0));}
//__________________________________________________________________________________________________
Int_t prim(Int_t tid)
{
+// Provides mother TID for the given TID
R()->GetLoader()->GetRunLoader()->LoadHeader(); R()->GetLoader()->GetRunLoader()->LoadKinematics();
if(tid<0||tid>=R()->GetLoader()->GetRunLoader()->Stack()->GetNtrack())
Bool_t ReadAlice()
{
Info("ReadAlice","Tring to read ALICE from SIMULATED FILE.");
- AliLoader::SetDebug(0);
if(gAlice) delete gAlice;
if(!(al=AliRunLoader::Open("galice.root","AlicE","update"))){
gSystem->Exec("rm -rf *.root *.dat");
Error("menu.C::ReadAlice","galice.root broken, removing all this garbage then init new one");
new AliRun("gAlice","Alice experiment system");
- gAlice->SetDebug(-1);
gAlice->Init("Config.C");
r=(AliRICH*)gAlice->GetDetector("RICH");
return kFALSE;
a->SetDebug(0);
//RICH
if(!(r=(AliRICH*)gAlice->GetDetector("RICH"))) Warning("RICH/menu.C::ReadAlice","No RICH in file");
- r->SetDebug(0);
if(!(rl=al->GetLoader("RICHLoader"))) Warning("RICH/menu.C::ReadAlice","No RICH loader in file");
Info("ReadAlice","Run contains %i event(s)",gAlice->GetEventsPerRun());
TPolyLine3D *pXaxis=new TPolyLine3D(2,X);pXaxis->SetLineColor(kRed); pXaxis->Draw();
TPolyLine3D *pYaxis=new TPolyLine3D(2,Y);pYaxis->SetLineColor(kGreen); pYaxis->Draw();
TPolyLine3D *pZaxis=new TPolyLine3D(2,Z);pZaxis->SetLineColor(kBlue); pZaxis->Draw();
-//chambers
+//Draw PC for all chambers by trasfering Pc plane using Pc2Mrs methode
Int_t iNpointsX=50,iNpointsY=50;
for(Int_t iChamberN=1;iChamberN<=7;iChamberN++){//chamber loop
TPolyMarker3D *pChamber=new TPolyMarker3D(iNpointsX*iNpointsY);
Int_t i=0;
for(Double_t x=0;x<p.PcSizeX();x+=p.PcSizeX()/iNpointsX)
for(Double_t y=0;y<p.PcSizeY();y+=p.PcSizeY()/iNpointsY){//step loop
- TVector3 v3=p.C(iChamberN)->Loc2Glob(TVector2(x,y));
- pChamber->SetPoint(i++,v3.X(),v3.Y(),v3.Z());
- pChamber->SetMarkerSize(1);
- pChamber->SetMarkerColor(iChamberN);
+ TVector3 v3=p.C(iChamberN)->Pc2Mrs(TVector2(x,y));//from regular grid of local PC points to MRS presentation
+ pChamber->SetPoint(i++,v3.X(),v3.Y(),v3.Z());//Pc plane poing in MRS
}//step loop
+ pChamber->SetMarkerSize(1);
+ pChamber->SetMarkerColor(iChamberN);
pChamber->Draw();
t.SetNDC();t.SetTextColor(iChamberN); t.DrawText(0.1,iChamberN*0.1,Form("Chamber %i",iChamberN));
}//chamber loop
- gPad->GetView()->RotateView(94,45);
+// gPad->GetView()->RotateView(94,45);
}//void TestSeg()
//__________________________________________________________________________________________________
void TestMenu()
}
fclose(pRawFile);
}
+//__________________________________________________________________________________________________
+void FillContribs(Int_t part,Int_t C, Bool_t print)
+{
+ static Int_t contrib[10][7][2];
+ static Bool_t kEnter=kFALSE;
+
+ C--; // chamber numbering from 1 to 7
+ if(!kEnter) {
+ for(Int_t i=0;i<10;i++) {
+ for(Int_t j=0;j<7;j++) {
+ for(Int_t k=0;k<2;k++) contrib[i][j][k]=0;
+ }
+ }
+ kEnter=kTRUE;
+ }
+
+ if(print) {
+ for(Int_t k=0;k<2;k++) {cout << "----------------Positives to RICH ---------------" << endl;
+ cout << " chamber 1 2 3 4 5 6 7 " << endl;
+ cout << " -------------------------------------------------" << endl;
+ for(Int_t i=0;i<4;i++) {
+ if(i==0) cout << " e";
+ if(i==1) cout << "pi";
+ if(i==2) cout << " K";
+ if(i==3) cout << " p";
+ if(k==0) cout << "+ ";
+ if(k==1) cout << "- ";
+ for(Int_t j=0;j<7;j++) {
+ cout << contrib[i][j][k] << " ";
+ }
+ cout << endl;
+ }
+ }
+ }
+
+ // +ves
+ if(part==kPositron) contrib[0][C][0]++;
+ if(part==kPiPlus) contrib[1][C][0]++;
+ if(part==kKPlus) contrib[2][C][0]++;
+ if(part==kProton) contrib[3][C][0]++;
+
+ // -ves
+ if(part==kElectron) contrib[0][C][1]++;
+ if(part==kPiMinus) contrib[1][C][1]++;
+ if(part==kKMinus) contrib[2][C][1]++;
+ if(part==kProtonBar) contrib[3][C][1]++;
+}
+//__________________________________________________________________________________________________
+void ParticleContribs()
+{
+
+ TH1F *distH1 = new TH1F("distH1","distH1",100,0.,5.);
+ AliRICH *pRich = (AliRICH*)gAlice->GetDetector("RICH");
+ Bool_t isHits =!pRich->GetLoader()->LoadHits();
+
+ pRich->GetLoader()->GetRunLoader()->LoadHeader(); pRich->GetLoader()->GetRunLoader()->LoadKinematics();
+
+
+ if(!isHits){Error("Exec","No hits. No contribs to calculate.");return;}
+
+ for(Int_t iEventN=0;iEventN<gAlice->GetEventsPerRun();iEventN++){//events Loop
+ pRich->GetLoader()->GetRunLoader()->GetEvent(iEventN);
+ cout << " event n. " << iEventN << endl;
+ Int_t nPrimaries = (Int_t)pRich->GetLoader()->TreeH()->GetEntries();
+ TObjArray * Hits = new TObjArray[nPrimaries];
+
+ for(Int_t i=0;i<nPrimaries;i++) {
+ pRich->GetLoader()->TreeH()->GetEntry(i);
+ Int_t nHits = pRich->Hits()->GetEntries();
+ for(Int_t k=0;k<nHits;k++) Hits[i].Add(pRich->Hits()->At(k));
+
+ }
+//deals with hits
+ for(Int_t i=0;i<nPrimaries;i++){//prims loop
+ pRich->GetLoader()->TreeH()->GetEntry(i);
+ Int_t nHits = pRich->Hits()->GetEntries();
+ for(Int_t j=0;j<nHits;j++){//hits loop
+ AliRICHhit *pHit = (AliRICHhit*)Hits[i].At(j);
+ TParticle *pParticle = pRich->GetLoader()->GetRunLoader()->GetAliRun()->Stack()->Particle(pHit->GetTrack());
+// if(pParticle->P()>1) FillContribs(pParticle->GetPdgCode(),pHit->C(),kFALSE);
+ FillContribs(pParticle->GetPdgCode(),pHit->C(),kFALSE);
+ if(pParticle->GetPDG()->Charge()!=0) distH1->Fill(pHit->Length());
+ }//hits loop
+ }//prims loop
+ }// event loop
+ FillContribs(0,0,kTRUE);
+ distH1->Draw();
+
+}//ParticleContribs()
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff