-///////////////////////////////////////////////////////////////////////////////
-// //
-// Ring Imaging Cherenkov //
-// This class contains the basic functions for the Ring Imaging Cherenkov //
-// detector. Functions specific to one particular geometry are //
-// contained in the derived classes //
-// //
-//Begin_Html
-/*
-<img src="picts/AliRICHClass.gif">
-*/
-//End_Html
-// //
-// //
-///////////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////
+// Manager and hits classes for set:RICH //
+////////////////////////////////////////////////
-#include <TBRIK.h>
+#include <TBRIK.h>
+#include <TTUBE.h>
#include <TNode.h>
#include <TRandom.h>
+#include <TObject.h>
+#include <TVector.h>
+#include <TObjArray.h>
-#include <TClass.h>
#include "AliRICH.h"
+#include "AliRICHHitMap.h"
#include "AliRun.h"
-#include "TGeant3.h"
-
-
+#include "AliMC.h"
+#include "iostream.h"
+#include "AliCallf77.h"
+
+// Static variables for the pad-hit iterator routines
+static Int_t sMaxIterPad=0;
+static Int_t sCurIterPad=0;
+static TClonesArray *fClusters2;
+static TClonesArray *fHits2;
+
ClassImp(AliRICH)
-
-//_____________________________________________________________________________
+
+//___________________________________________
AliRICH::AliRICH()
{
- //
- // Default constructor for RICH
- //
- fIshunt = 0;
- fHits = 0;
- fMips = 0;
- fCkovs = 0;
- fPadhits = 0;
- fNmips = 0;
- fNckovs = 0;
- fNpadhits = 0;
-
- fChslope = 0;
- fAlphaFeed= 0;
- fSxcharge = 0;
- fIritri = 0;
+ fIshunt = 0;
+ fHits = 0;
+ fClusters = 0;
+ fNclusters = 0;
+ fDchambers = 0;
+ fRecClusters= 0;
+ fCerenkovs = 0;
+ fNdch = 0;
}
-
-//_____________________________________________________________________________
+
+//___________________________________________
AliRICH::AliRICH(const char *name, const char *title)
- : AliDetector(name,title)
+ : AliDetector(name,title)
{
- //
- // Standard constructor for RICH
- //
- fIshunt = 0;
- fNmips = 0;
- fNckovs = 0;
- fNpadhits = 0;
- //
- // Allocate space for different components
- fHits = new TClonesArray("AliRICHhit", 100);
- fMips = new TClonesArray("AliRICHmip", 100);
- fCkovs = new TClonesArray("AliRICHckov", 100);
- fPadhits = new TClonesArray("AliRICHpadhit", 100);
- //
- // Set parameters to default value
- fChslope = 40;
- fAlphaFeed= 0.04;
- fSxcharge = 0.18;
- fIritri = 0;
- //
- SetMarkerColor(6);
- SetMarkerStyle(20);
- SetMarkerSize(0.5);
+//Begin_Html
+/*
+ <img src="gif/alirich.gif">
+*/
+//End_Html
+
+ fHits = new TClonesArray("AliRICHhit",1000 );
+ fClusters = new TClonesArray("AliRICHcluster",10000);
+ fCerenkovs = new TClonesArray("AliRICHCerenkov",20000);
+ fNclusters = 0;
+ fIshunt = 0;
+
+ fNdch = new Int_t[7];
+
+ fDchambers = new TObjArray(7);
+
+ Int_t i;
+
+ for (i=0; i<7 ;i++) {
+ (*fDchambers)[i] = new TClonesArray("AliRICHdigit",10000);
+ fNdch[i]=0;
+ }
+
+ fRecClusters=new TObjArray(7);
+ for (i=0; i<7;i++)
+ (*fRecClusters)[i] = new TObjArray(1000);
+
+//
+// Transport angular cut
+ fAccCut=0;
+ fAccMin=2;
+ fAccMax=9;
+
+ SetMarkerColor(kRed);
}
-//_____________________________________________________________________________
+//___________________________________________
AliRICH::~AliRICH()
{
- //
- // Destructor for RICH
- //
- fIshunt = 0;
- delete fHits;
- fMips->Delete(); delete fMips;
- fCkovs->Delete(); delete fCkovs;
- fPadhits->Delete(); delete fPadhits;
+ fIshunt = 0;
+ delete fHits;
+ delete fClusters;
+ delete fCerenkovs;
+ for (Int_t i=0;i<7;i++)
+ delete (*fRecClusters)[i];
+ delete fRecClusters;
+
}
-//_____________________________________________________________________________
+//___________________________________________
void AliRICH::AddHit(Int_t track, Int_t *vol, Float_t *hits)
{
- //
- // Add a RICH hit
- //
- switch ((int) hits[0]) {
- case 0:
- {
- //
- // Simple hit
- TClonesArray &lhits = *fHits;
- new(lhits[fNhits++]) AliHit(fIshunt,track);
- AliHit *lhit = (AliHit*) fHits->AddrAt(fNhits-1);
- lhit->fX = hits[19];
- lhit->fY = hits[20];
- lhit->fZ = hits[21];
- } break;
- case 1:
- //
- // MIP hit
- AddMipHit(track,vol,hits);
- break;
- case 2:
- //
- // Cherenkov hit
- AddCkovHit(track,vol,hits);
- break;
- case 3:
- //
- // Pad hit
- AddPadHit(track,vol,hits);
- break;
- case 4:
- // Update a mip hit
- UpdateMipHit(hits);
- break;
- default:
- printf("Error: AliRICH::AddHit flag %d not defined./n",(int) hits[0]);
- return;
- }
+ TClonesArray &lhits = *fHits;
+ new(lhits[fNhits++]) AliRICHhit(fIshunt,track,vol,hits);
}
-
//_____________________________________________________________________________
-void AliRICH::AddMipHit(Int_t track, Int_t *vol, Float_t *hits)
-{
- // Adds a mip hit in the RICH.
- TClonesArray &lhits = *fMips;
- new(lhits[fNmips++]) AliRICHmip(fIshunt,track,vol,hits,
- fNckovs,fNpadhits);
+void AliRICH::AddCerenkov(Int_t track, Int_t *vol, Float_t *cerenkovs)
+{
+ TClonesArray &lcerenkovs = *fCerenkovs;
+ new(lcerenkovs[fNcerenkovs++]) AliRICHCerenkov(fIshunt,track,vol,cerenkovs);
}
-
-//_____________________________________________________________________________
-void AliRICH::AddCkovHit(Int_t track, Int_t *vol, Float_t *hits)
-{
- //
- // Adds a cerenkov hit in the RICH.
- //
- TClonesArray &lhits = *fCkovs;
- AliRICHmip *lmip = (AliRICHmip*) fMips->AddrAt(fNmips-1);
- //
- // If this ckov come from a mip update the mip lastckov entry.
- Int_t fmipslocal=-1;
- if (lmip->GetZ() != -999.0) {
- fmipslocal = fNmips-1;
- lmip->SetLastCkov(fNckovs);
- }
- new(lhits[fNckovs++]) AliRICHckov(fIshunt,track,vol,hits,
- fmipslocal,fNpadhits);
+//___________________________________________
+void AliRICH::AddCluster(Int_t *clhits)
+{
+ TClonesArray &lclusters = *fClusters;
+ new(lclusters[fNclusters++]) AliRICHcluster(clhits);
}
-
//_____________________________________________________________________________
-void AliRICH::AddPadHit(Int_t track, Int_t *vol, Float_t *hits)
-{
- //
- // Adds pad hits in the RICH
- //
- TClonesArray &lhits = *fPadhits;
- // Update the last padhit of the respective particle:
- if ((int) hits[1]==50) { // a ckov
- ((AliRICHckov *) fCkovs->AddrAt(fNckovs-1))->SetLastpad(fNpadhits);
- new(lhits[fNpadhits++]) AliRICHpadhit(fIshunt,track,vol,hits,-1,fNckovs-1);
- }else { // a mip
- ((AliRICHmip *) fMips->AddrAt(fNmips-1))->SetLastpad(fNpadhits);
- new(lhits[fNpadhits++]) AliRICHpadhit(fIshunt,track,vol,hits,fNmips-1,-1);
- }
+void AliRICH::AddDigits(Int_t id, Int_t *tracks, Int_t *charges, Int_t *digits)
+{
+ //
+ // Add a RICH digit to the list
+ //
+
+ TClonesArray &ldigits = *((TClonesArray*)(*fDchambers)[id]);
+ new(ldigits[fNdch[id]++]) AliRICHdigit(tracks,charges,digits);
}
+
//_____________________________________________________________________________
-void AliRICH::BuildGeometry()
+void AliRICH::AddRecCluster(Int_t iCh, Int_t iCat, AliRICHRecCluster* Cluster)
{
- //
- // Builds a TNode geometry for event display
- //
- TNode *Node, *Top;
-
- const int kColorRICH = kGreen;
- //
- Top=gAlice->GetGeometry()->GetNode("alice");
-
- new TRotMatrix("rot993","rot993",90,0,70.69,90,19.30999,-90);
- new TRotMatrix("rot994","rot994",90,-20,90,70,0,0);
- new TRotMatrix("rot995","rot995",90,0,90,90,0,0);
- new TRotMatrix("rot996","rot996",90,20,90,110,0,0);
- new TRotMatrix("rot997","rot997",90,340,108.1999,70,18.2,70);
- new TRotMatrix("rot998","rot998",90,0,109.3099,90,19.30999,90);
- new TRotMatrix("rot999","rot999",90,20,108.1999,110,18.2,110);
- new TBRIK("S_RICH","S_RICH","void",71.09999,11.5,73.15);
- Top->cd();
- Node = new TNode("RICH1","RICH1","S_RICH",0,471.8999,165.2599,"rot993");
- Node->SetLineColor(kColorRICH);
- fNodes->Add(Node);
- Top->cd();
- Node = new TNode("RICH2","RICH2","S_RICH",171,470,0,"rot994");
- Node->SetLineColor(kColorRICH);
- fNodes->Add(Node);
- Top->cd();
- Node = new TNode("RICH3","RICH3","S_RICH",0,500,0,"rot995");
- Node->SetLineColor(kColorRICH);
- fNodes->Add(Node);
- Top->cd();
- Node = new TNode("RICH4","RICH4","S_RICH",-171,470,0,"rot996");
- Node->SetLineColor(kColorRICH);
- fNodes->Add(Node);
- Top->cd();
- Node = new TNode("RICH5","RICH5","S_RICH",161.3999,443.3999,-165.3,"rot997");
- Node->SetLineColor(kColorRICH);
- fNodes->Add(Node);
- Top->cd();
- Node = new TNode("RICH6","RICH6","S_RICH",0,471.8999,-165.3,"rot998");
- Node->SetLineColor(kColorRICH);
- fNodes->Add(Node);
- Top->cd();
- Node = new TNode("RICH7","RICH7","S_RICH",-161.399,443.3999,-165.3,"rot999");
- Node->SetLineColor(kColorRICH);
- fNodes->Add(Node);
+ //
+ // Add a RICH reconstructed cluster to the list
+ //
+ TObjArray* ClusterList = RecClusters(iCh,iCat);
+ ClusterList->Add(Cluster);
}
-//*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
-//
-// Section for Rich Step
-#define MAXPH 1000
-
-static Float_t sVloc[3];
-static Float_t sVectIn[3];
-static Float_t sDetot;
-static Float_t sYanode[MAXPH];
-static Float_t sXpad[MAXPH];
-static Float_t sYpad[MAXPH];
-static Int_t sNpx;
-static Int_t sNpy;
-static Float_t sDxp;
-static Float_t sDyp;
-static Float_t sDlx;
-static Float_t sDly;
-static Float_t sDpad;
-
-
-static Int_t sNsecon;
-//static Float_t sQint[2];
-
-#define MAXSEC 2000
-
-static Int_t sNpads;
-static Int_t sIpx[MAXSEC];
-static Int_t sIpy[MAXSEC];
-static Int_t sIqpad[MAXSEC];
-static Int_t sNphlink[MAXSEC];
-static Int_t sNphoton;
-
-static Int_t sNfeed, sNfeedd, sNdir;
-
-static Float_t sPparent;
-static Float_t sThincParent;
-static Int_t sIloss[MAXPH];
-static Int_t sIprod[MAXPH];
-static Float_t sXphit[MAXPH];
-static Float_t sYphit[MAXPH];
-
-static Float_t sSycharge;
-static Float_t sXsox, sYsox, sZsox;
-
-static Int_t sMckov[MAXPH];
-static Int_t idpartgx;
-static Float_t phisx;
-static Int_t nmodsx;
-static Float_t psx;
-static Float_t xsx;
-static Float_t ysx;
-static Float_t thetasx;
-//_____________________________________________________________________________
-void AliRICH::Init()
+//___________________________________________
+void AliRICH::BuildGeometry()
+
{
- //
- // Initialise RICH after that it has been built
- //
- const Float_t sconv=2*TMath::Sqrt(2*TMath::Log(2));
- const Float_t yK3=1.20;
- Float_t ansp;
- Int_t i;
- //
- sNpx=162;
- sNpy=162;
- sDxp=0.80;
- sDyp=0.80;
- ansp=sDyp/2;
- sDlx=sNpx*sDxp/2;
- sDly=sNpy*sDyp/2;
- sDpad=0.2;
- //
- for(i=0;i<sNpx;i++) {
- sXpad[i]=i*sDxp;
- sYpad[i]=i*sDyp;
- }
- for(i=0;i<2*sNpy+1;i++) sYanode[i]=ansp/2+i*ansp;
- //
-
- sSycharge=4*TMath::ATanH(1/TMath::Sqrt(2*yK3))/TMath::Pi()
- /(1-0.5*TMath::Sqrt(yK3))/sDpad/sconv;
- //
+ //
+ // Builds a TNode geometry for event display
+ //
+ TNode *Node, *Top;
+
+ const int kColorRICH = kGreen;
+ //
+ Top=gAlice->GetGeometry()->GetNode("alice");
+
+
+ new TBRIK("S_RICH","S_RICH","void",71.09999,11.5,73.15);
+
+ Top->cd();
+ Float_t pos1[3]={0,471.8999,165.2599};
+ Chamber(0).SetChamberTransform(pos1[0],pos1[1],pos1[2],new TRotMatrix("rot993","rot993",90,0,70.69,90,19.30999,-90));
+ Node = new TNode("RICH1","RICH1","S_RICH",pos1[0],pos1[1],pos1[2],"rot993");
+
+
+ Node->SetLineColor(kColorRICH);
+ fNodes->Add(Node);
+ Top->cd();
+
+ Float_t pos2[3]={171,470,0};
+ Chamber(1).SetChamberTransform(pos2[0],pos2[1],pos2[2],new TRotMatrix("rot994","rot994",90,-20,90,70,0,0));
+ Node = new TNode("RICH2","RICH2","S_RICH",pos2[0],pos2[1],pos2[2],"rot994");
+
+
+ Node->SetLineColor(kColorRICH);
+ fNodes->Add(Node);
+ Top->cd();
+ Float_t pos3[3]={0,500,0};
+ Chamber(2).SetChamberTransform(pos3[0],pos3[1],pos3[2],new TRotMatrix("rot995","rot995",90,0,90,90,0,0));
+ Node = new TNode("RICH3","RICH3","S_RICH",pos3[0],pos3[1],pos3[2],"rot995");
+
+
+ Node->SetLineColor(kColorRICH);
+ fNodes->Add(Node);
+ Top->cd();
+ Float_t pos4[3]={-171,470,0};
+ Chamber(3).SetChamberTransform(pos4[0],pos4[1],pos4[2], new TRotMatrix("rot996","rot996",90,20,90,110,0,0));
+ Node = new TNode("RICH4","RICH4","S_RICH",pos4[0],pos4[1],pos4[2],"rot996");
+
+
+ Node->SetLineColor(kColorRICH);
+ fNodes->Add(Node);
+ Top->cd();
+ Float_t pos5[3]={161.3999,443.3999,-165.3};
+ Chamber(4).SetChamberTransform(pos5[0],pos5[1],pos5[2],new TRotMatrix("rot997","rot997",90,340,108.1999,70,18.2,70));
+ Node = new TNode("RICH5","RICH5","S_RICH",pos5[0],pos5[1],pos5[2],"rot997");
+
+ Node->SetLineColor(kColorRICH);
+ fNodes->Add(Node);
+ Top->cd();
+ Float_t pos6[3]={0., 471.9, -165.3,};
+ Chamber(5).SetChamberTransform(pos6[0],pos6[1],pos6[2],new TRotMatrix("rot998","rot998",90,0,109.3099,90,19.30999,90));
+ Node = new TNode("RICH6","RICH6","S_RICH",pos6[0],pos6[1],pos6[2],"rot998");
+
+
+ Node->SetLineColor(kColorRICH);
+ fNodes->Add(Node);
+ Top->cd();
+ Float_t pos7[3]={-161.399,443.3999,-165.3};
+ Chamber(6).SetChamberTransform(pos7[0],pos7[1],pos7[2],new TRotMatrix("rot999","rot999",90,20,108.1999,110,18.2,110));
+ Node = new TNode("RICH7","RICH7","S_RICH",pos7[0],pos7[1],pos7[2],"rot999");
+ Node->SetLineColor(kColorRICH);
+ fNodes->Add(Node);
+
}
+//___________________________________________
+Int_t AliRICH::DistancetoPrimitive(Int_t , Int_t )
+{
+ return 9999;
+}
-//___________________________________________________________________________
-void AliRICH::StepManager()
+//___________________________________________
+void AliRICH::MakeBranch(Option_t* option)
{
- //
- // Called at every step in the RICH
- //
-
- TGeant3 *geant3 = (TGeant3*) gMC;
-
- const Float_t xshift[3] = { 41.3, 0, -41.3 };
- const Int_t nrooth = 25;
-
- static Int_t ixold=-1, iyold=-1;
-
- // System generated locals
- Int_t j, i1;
- Float_t r1, r2;
-
- // Local variables
- Float_t ranf[2], rrhh[nrooth], phiangle, cost, vmod;
- //Int_t idpartsx;
- Int_t i;
- Float_t t, vxloc[3];
- Int_t ll, irivol[2];
- Int_t lcase;
- //Int_t iprimx;
- Int_t ix, iy;
- Float_t stwght;
- Float_t cophi;
- Float_t dir[3];
- Int_t ihitrak;
- Int_t medprod;
-
- Int_t nmult=0;
- //Float_t xtrig[200], ytrig[200];
- //Int_t itrig[200];
-
-
-
- // ILOSS = 0 NOT LOST
- // 1 REFLECTED FREON-QUARZ
- // 2 REFLECTED QUARZ-METHANE
- // 3 REFLECTED METHANE-CSI
- // 11 ABSORBED IN FREON
- // 12 ABSORBED IN QUARZ
- // 13 ABSORBED IN METHANE
- // 21 CSI QE
- // IPROD = 1 PRODUCED IN FREON
- // IPROD = 2 PRODUCED IN QUARZ
-
- // new (changed NROOTH from 10 to 25!!!!!!!!!!!!!)
-
-
- Int_t *idtmed = fIdtmed->GetArray()-999;
-
- //--------------------------------------------------------------------------
-
- // MIP inside CsI
-
- if (geant3->Gckine()->charge) {
+ // Create Tree branches for the RICH.
- // Charged particles treatment
- if (fIritri && !geant3->Gctrak()->upwght) {
- if (geant3->Gctmed()->numed == idtmed[fIritri-1]) {
- if (geant3->Gcking()->ngkine > 0) {
- strncpy((char *)&lcase,"HADR",4);
- if (geant3->Gcking()->kcase == lcase) {
- i1 = geant3->Gcking()->ngkine;
- for (i = 1; i <= i1; ++i) {
- gMC->Gmtod(geant3->Gckin3()->gpos[i-1], vxloc, 1);
- gMC->Gmtod(geant3->Gcking()->gkin[i-1], dir, 2);
- if (geant3->Gcking()->gkin[i-1][4] == 8. ||
- geant3->Gcking()->gkin[i-1][4] == 9.) {
- ++nmult;
- // Computing 2nd power
- r1 = dir[0];
- // Computing 2nd power
- r2 = dir[2];
- //theta = TMath::ATan2(TMath::Sqrt(r1*r1+r2*r2),dir[1]);
- //xtrig[nmult - 1] = theta;
- //ytrig[nmult - 1] = vxloc[1] + .25;
- //itrig[nmult - 1] = (Int_t) geant3->Gcking()->gkin[i-1][4];
- }
- }
- }
- }
- }
- if ((geant3->Gctmed()->numed == idtmed[1006-1] &&
- geant3->Gctrak()->inwvol == 2) ||
- geant3->Gctrak()->istop) {
- if (!nmult) {
- printf("NOT TRIGGERED\n");
- sDetot = 0.;
- sNsecon = 0;
- sNpads = 0;
- sNphoton = 0;
- sNfeed = 0;
- sNfeedd = 0;
- sNdir = 0;
- geant3->Gctrak()->istory = 0;
- geant3->Gctrak()->upwght = 0.;
- geant3->Gcflag()->ieotri = 1;
- nmult = 0;
- //sQint[0] = 0.;
- //sQint[1] = 0.;
- } else {
- printf("TRIGGERED %d\n",nmult);
- }
- }
- }
- // MIP inside Methane
- if (geant3->Gctmed()->numed == idtmed[1009-1]) {
-
- // new
- // If particle produced already Cerenkov Photons (istory=1)
- // update the impact point only
- if (geant3->Gctrak()->istory == 1) {
- // Direction of incidence and where did it hit ?
- gMC->Gmtod(geant3->Gctrak()->vect, sVloc, 1);
- gMC->Gmtod(&geant3->Gctrak()->vect[3], dir, 2);
- phiangle = TMath::ATan2(dir[2], dir[0]);
- if (phiangle < 0.) phiangle += 2*TMath::Pi();
- i1 = nrooth;
- for (ll = 0; ll < i1; ++ll) rrhh[ll] = 0;
- irivol[0] = geant3->Gcvolu()->number[geant3->Gcvolu()->nlevel - 4];
- irivol[1] = geant3->Gcvolu()->number[geant3->Gcvolu()->nlevel - 4];
- // NMODSX
- ihitrak = gAlice->CurrentTrack();
- rrhh[0] = 4.;
- // flag to say this is update
- rrhh[1] = sVloc[0] + sDlx;
- // XSX
- rrhh[2] = sVloc[2] + sDly;
- // YSX
- rrhh[3] = (Float_t) geant3->Gcvolu()->number[geant3->Gcvolu()->nlevel - 4];
- // NMODSX
- // Computing 2nd power
- r1 = dir[0];
- // Computing 2nd power
- r2 = dir[2];
- rrhh[4] = TMath::ATan2(TMath::Sqrt(r1 * r1 + r2 * r2), dir[1]);
- // theta
- rrhh[5] = phiangle;
- // PRINT *, 'ITRA = ',ITRA,'ISTAK = ',ISTAK
- AddHit(ihitrak,irivol,rrhh);
- }
- // enew
- // Record particle properties
- // If particle produced already Cerenkov Photons (istory=1)
-
- // update the impact point only
- if (geant3->Gctrak()->istory != 2) {
- if (!geant3->Gctrak()->istory) {
- ++sNsecon;
-
- // Is this a primary particle ?
- //iprimx = 1;
- //if (geant3->Gctrak()->upwght) iprimx = 0;
-
- // Where did it come from ?
- sXsox = geant3->Gckine()->vert[0];
- sYsox = geant3->Gckine()->vert[1];
- sZsox = geant3->Gckine()->vert[2];
-
- // Momentum
- psx = geant3->Gctrak()->vect[6];
-
- // Particle type and parent
- //idpartsx = geant3->Gckine()->ipart;
- r1 = geant3->Gctrak()->upwght / 100.;
- idpartgx = Int_t(r1+0.5);
- if (!geant3->Gctrak()->upwght) {
- sPparent = geant3->Gctrak()->vect[6];
- sThincParent = thetasx;
- }
-
- // Direction of incidence and where did it hit ?
- gMC->Gmtod(geant3->Gctrak()->vect, sVloc, 1);
- gMC->Gmtod(&geant3->Gctrak()->vect[3], dir, 2);
- // Computing 2nd power
- r1 = dir[0];
- // Computing 2nd power
- r2 = dir[2];
- thetasx = TMath::ATan2(TMath::Sqrt(r1 * r1 + r2 * r2), dir[1]);
- phisx = TMath::ATan2(dir[2], dir[0]);
- if (phisx < 0.) phisx += 2*TMath::Pi();
- ysx = sVloc[2] + sDly;
- xsx = sVloc[0] + sDlx;
- nmodsx = geant3->Gcvolu()->number[geant3->Gcvolu()->nlevel - 4];
- // new
- i1 = nrooth;
- for (ll = 0; ll < i1; ++ll) rrhh[ll] = 0;
- irivol[0] = geant3->Gcvolu()->number[geant3->Gcvolu()->nlevel - 4];
- irivol[1] = nmodsx;
- ihitrak = gAlice->CurrentTrack();
- rrhh[0] = 1.;
- // Flag to say this is MIP
- rrhh[1] = (Float_t) geant3->Gckine()->ipart;
- rrhh[2] = xsx;
- rrhh[3] = ysx;
- rrhh[4] = (Float_t) nmodsx;
- // Module Number
- rrhh[5] = thetasx;
- rrhh[6] = geant3->Gctrak()->tofg;
- // in seconds
- rrhh[7] = (Float_t) idpartgx;
- // mips specific
- rrhh[8] = phisx;
- rrhh[9] = psx;
- // charge of current particle in electron charge unit;
- rrhh[10] = geant3->Gckine()->charge;
- rrhh[11] = -999.;
- // Zo of ckov generation
- rrhh[12] = 0.;
- // no ckov !!!
- AddHit(ihitrak, irivol, rrhh);
- // end of new
-
- // Earmark track as being recorded in methane gap
-
- geant3->Gctrak()->istory = 2;
- }
- }
-
- // Signal generation in methane gap
- gMC->Gmtod(geant3->Gctrak()->vect, sVloc, 1);
- gMC->Gmtod(geant3->Gckine()->vert, vxloc, 1);
- ix = (Int_t) ((sVloc[0] + sDlx) / sDxp);
- iy = (Int_t) ((sVloc[2] + sDly) / sDyp);
-
- // Is this the first step?
- if (ixold == -1 && iyold == -1) {
- ixold = ix;
- iyold = iy;
- for(j=0;j<3;j++) sVectIn[j]=geant3->Gctrak()->vect[j];
- }
-
- // Mip left gap
- if (geant3->Gctrak()->inwvol == 2 || geant3->Gctrak()->istop) {
- sDetot += geant3->Gctrak()->destep;
- if (sDetot > 0.) RichIntegration();
- sDetot = 0.;
- ixold = -1;
- iyold = -1;
-
- // Mip left current pad
- } else if (ixold != ix || iyold != iy) {
- if (sDetot > 0.) RichIntegration();
- for(j=0;j<3;j++) sVectIn[j]=geant3->Gctrak()->vect[j];
- sDetot = geant3->Gctrak()->destep;
- ixold = ix;
- iyold = iy;
- } else {
- sDetot += geant3->Gctrak()->destep;
- }
- }
- }
-
- // End charged particles treatment
- //
- // Treat photons produced in Freon and Quartz
- if (geant3->Gckin2()->ngphot > 0 &&
- (geant3->Gctmed()->numed == idtmed[1004-1] ||
- geant3->Gctmed()->numed == idtmed[1003-1])) {
- if (!geant3->Gctrak()->upwght) {
-
- // If it is a primary, save all generated photons
- i1 = geant3->Gckin2()->ngphot;
- for (i = 1; i <= i1; ++i) {
- ++sNphoton;
- if (sNphoton > MAXPH) {
- sNphoton = MAXPH;
- printf("ATTENTION NPHOTON %d\n",sNphoton);
- continue;
- }
-
- // Production medium
- medprod = 1;
- if (geant3->Gctmed()->numed == idtmed[1003-1]) medprod = 2;
- //
- // Production angle and energy
- vmod=0;
- cost=0;
- for(j=0;j<3;j++) {
- cost+=geant3->Gckin2()->xphot[i-1][3+j]*geant3->Gctrak()->vect[3+j];
- vmod+=geant3->Gckin2()->xphot[i-1][3+j]*
- geant3->Gckin2()->xphot[i-1][3+j];
- }
- cost/=sqrt(vmod);
- sIloss[sNphoton - 1] = 22;
- sIprod[sNphoton - 1] = medprod;
- sXphit[sNphoton - 1] = 0.;
- sYphit[sNphoton - 1] = 0.;
- stwght = geant3->Gctrak()->upwght;
- geant3->Gctrak()->upwght = (Float_t) sNphoton;
- // geant3->Gskpho(i);
- sMckov[sNphoton - 1] = gAlice->CurrentTrack()+i;
- geant3->Gctrak()->upwght = stwght;
- }
- } else {
- stwght = geant3->Gctrak()->upwght;
- geant3->Gctrak()->upwght = 0.;
- // geant3->Gskpho(0);
- geant3->Gctrak()->upwght = stwght;
- }
+ const Int_t buffersize = 4000;
+ char branchname[20];
- // Particle did not yet pass the methane gap
- if (geant3->Gctrak()->istory == 0) {
- geant3->Gctrak()->istory = 1;
- ++sNsecon;
- // Is this a primary particle ?
- //iprimx = 1;
- //if (geant3->Gctrak()->upwght) iprimx = 0;
-
- // Where did it come from ?
- sXsox = geant3->Gckine()->vert[0];
- sYsox = geant3->Gckine()->vert[1];
- sZsox = geant3->Gckine()->vert[2];
-
- // Where did it hit ?
- gMC->Gmtod(geant3->Gctrak()->vect, sVloc, 1);
- gMC->Gmtod(&geant3->Gctrak()->vect[3], dir, 2);
- ysx = sVloc[2] + sDly;
- if (geant3->Gctmed()->numed == idtmed[1004-1]) {
- nmodsx = geant3->Gcvolu()->number[geant3->Gcvolu()->nlevel - 4];
- xsx = sVloc[0] + sDlx +
- xshift[geant3->Gcvolu()->number[geant3->Gcvolu()->nlevel - 2] - 1];
- } else if (geant3->Gctmed()->numed == idtmed[1003-1]) {
- nmodsx = geant3->Gcvolu()->number[geant3->Gcvolu()->nlevel - 3];
- xsx = sVloc[0] + sDlx;
- } else {
- nmodsx = 0;
- }
-
- // Momentum and direction of incidence
- psx = geant3->Gctrak()->vect[6];
- // Computing 2nd power
- r1 = dir[0];
- // Computing 2nd power
- r2 = dir[2];
- thetasx = TMath::ATan2(TMath::Sqrt(r1 * r1 + r2 * r2), dir[1]);
- phisx = TMath::ATan2(dir[2], dir[0]);
- if (phisx < 0.) phisx += 2*TMath::Pi();
-
- // Particle type and parent
- //idpartsx = geant3->Gckine()->ipart;
- r1 = geant3->Gctrak()->upwght / 100.;
- idpartgx = Int_t(r1+0.5);
- if (!geant3->Gctrak()->upwght) {
- sPparent = geant3->Gctrak()->vect[6];
- sThincParent = thetasx;
- }
- // new
- for (ll = 0; ll <nrooth; ++ll) rrhh[ll] = 0;
- irivol[0] = geant3->Gcvolu()->number[geant3->Gcvolu()->nlevel - 4];
- irivol[1] = nmodsx;
- ihitrak = gAlice->CurrentTrack();
- rrhh[0] = 1.;
- // Flag to say that this is MIP
- rrhh[1] = (Float_t) geant3->Gckine()->ipart;
- rrhh[2] = xsx;
- rrhh[3] = ysx;
- rrhh[4] = (Float_t) nmodsx;
- // Module Number
- rrhh[5] = thetasx;
- rrhh[6] = geant3->Gctrak()->tofg;
- // in seconds
- rrhh[7] = (Float_t) idpartgx;
- // mips specific
- rrhh[8] = phisx;
- rrhh[9] = psx;
- // charge of current particle in electron charge unit;
- rrhh[10] = geant3->Gckine()->charge;
- rrhh[11] = sVloc[1];
- // Zo of ckov generation
- rrhh[12] = 1.;
- // ckov generation
- AddHit(ihitrak, irivol, rrhh);
- // enew
+
+ AliDetector::MakeBranch(option);
+ sprintf(branchname,"%sCerenkov",GetName());
+ if (fCerenkovs && gAlice->TreeH()) {
+ gAlice->TreeH()->Branch(branchname,&fCerenkovs, buffersize);
+ printf("Making Branch %s for Cerenkov Hits\n",branchname);
}
- }
-
- // Current particle is cherenkov photon
- if (geant3->Gckine()->ipart == 50) {
- gMC->Gmtod(geant3->Gctrak()->vect, sVloc, 1);
- // WRITE(6,* ) UPWGHT, VLOC(2), NUMED, DESTEP
- // Photon crosses ch4-csi boundary
- // take into account fresnel losses with complex refraction index
- if (geant3->Gctrak()->inwvol == 1 && geant3->Gctmed()->numed == idtmed[1006-1]) {
-
- // fresnel losses commented out for the moment
- // make sure that qe correction for fresnel losses is also switched off !
- // CALL FRESNELCSI
- // IF (ISTOP .EQ. 2) RETURN
- // Put transmission of electrodes in by hand
- gMC->Gmtod(&geant3->Gctrak()->vect[3], dir, 2);
- cophi = TMath::Cos(TMath::ATan2(dir[0], dir[1]));
- t = (1. - .025 / cophi) * (1. - .05 / cophi);
- gMC->Rndm(ranf, 1);
- if (ranf[0] > t) {
- if (geant3->Gctrak()->upwght && Int_t(geant3->Gctrak()->upwght+0.5)<MAXPH)
- sIloss[Int_t(geant3->Gctrak()->upwght+0.5) - 1] = 15;
- geant3->Gctrak()->istop = 2;
- return;
- }
+
+ sprintf(branchname,"%sCluster",GetName());
+ if (fClusters && gAlice->TreeH()) {
+ gAlice->TreeH()->Branch(branchname,&fClusters, buffersize);
+ printf("Making Branch %s for clusters\n",branchname);
}
- // Photon lost energy in CsI
- if (geant3->Gctrak()->destep > 0. && geant3->Gctmed()->numed == idtmed[1006-1]) {
- geant3->Gctrak()->istop = 2;
- r1 = geant3->Gctrak()->upwght / 100.;
- if (Int_t(r1+0.5) > 50) {
- ++sNfeedd;
- } else {
- ++sNdir;
- }
- // WRITE(6,*) 'PHOTON',UPWGHT, MAXPH
- if (geant3->Gctrak()->upwght && Int_t(geant3->Gctrak()->upwght+0.5) < MAXPH)
- sIloss[Int_t(geant3->Gctrak()->upwght+0.5) - 1] = 0;
- // write(6,*) 'photon detected'
- for(j=0;j<3;j++) sVectIn[j]=geant3->Gctrak()->vect[j];
- // new
- // copied from miphit in Freon or Quartz
- // Where did it hit ?
- gMC->Gmtod(&geant3->Gctrak()->vect[3], dir, 2);
-
- // Momentum and direction of incidence
- for (ll = 0; ll < nrooth; ++ll) rrhh[ll]=0;
- irivol[0] = geant3->Gcvolu()->number[geant3->Gcvolu()->nlevel - 4];
- // ???
- irivol[1] = geant3->Gcvolu()->number[geant3->Gcvolu()->nlevel - 4];
- // ???
- ihitrak = gAlice->CurrentTrack();
- rrhh[0] = 2.;
- // Flag to say that this is CK
- rrhh[1] = (Float_t) geant3->Gckine()->ipart;
- rrhh[2] = sVloc[0] + sDlx;
- rrhh[3] = sVloc[2] + sDly;
- rrhh[4] = (Float_t) geant3->Gcvolu()->number[geant3->Gcvolu()->nlevel - 3];
- // ??? Module Number
- // Computing 2nd power
- r1 = dir[0];
- // Computing 2nd power
- r2 = dir[2];
- rrhh[5] = TMath::ATan2(TMath::Sqrt(r1 * r1 + r2 * r2), dir[1]);
- // THETASX
- rrhh[6] = geant3->Gctrak()->tofg;
- // in seconds
- r1 = geant3->Gctrak()->upwght / 100.;
- rrhh[7] = (Float_t) Int_t(r1+0.5);
- // ckov specific
- // Feedback ???
- rrhh[8] = geant3->Gctrak()->getot;
- rrhh[9] = 0.;
- // Stop in CsI
- AddHit(ihitrak, irivol, rrhh);
- // end of new
- RichIntegration();
- return;
+// one branch for digits per chamber
+ Int_t i;
+
+ for (i=0; i<7 ;i++) {
+ sprintf(branchname,"%sDigits%d",GetName(),i+1);
+
+ if (fDchambers && gAlice->TreeD()) {
+ gAlice->TreeD()->Branch(branchname,&((*fDchambers)[i]), buffersize);
+ printf("Making Branch %s for digits in chamber %d\n",branchname,i+1);
+ }
}
- if (geant3->Gctrak()->upwght && Int_t(geant3->Gctrak()->upwght+0.5) < MAXPH)
- {
- // Losses by reflection and absorption and leaving detector
- if (sIloss[Int_t(geant3->Gctrak()->upwght+0.5) - 1] == 22) {
- i1 = geant3->Gctrak()->nmec;
- for (i = 0; i < i1; ++i) {
- // Reflection loss
- if (geant3->Gctrak()->lmec[i] == 106) {
- sIloss[Int_t(geant3->Gctrak()->upwght+0.5)-1]=10;
- if (geant3->Gctmed()->numed == idtmed[1004-1])
- sIloss[Int_t(geant3->Gctrak()->upwght+0.5)-1]=1;
-
- if (geant3->Gctmed()->numed == idtmed[1003-1])
- sIloss[Int_t(geant3->Gctrak()->upwght+0.5)-1]=2;
-
- // Absorption loss
- } else if (geant3->Gctrak()->lmec[i] == 101) {
- geant3->Gctrak()->istop = 2;
- sIloss[Int_t(geant3->Gctrak()->upwght+0.5)-1]=20;
- if (geant3->Gctmed()->numed == idtmed[1004-1])
- sIloss[Int_t(geant3->Gctrak()->upwght+0.5)-1]=11;
-
- if (geant3->Gctmed()->numed == idtmed[1003-1])
- sIloss[Int_t(geant3->Gctrak()->upwght+0.5)-1]=12;
-
- if (geant3->Gctmed()->numed == idtmed[1005-1])
- sIloss[Int_t(geant3->Gctrak()->upwght+0.5)-1]=13;
-
- if (geant3->Gctmed()->numed == idtmed[1009-1])
- sIloss[Int_t(geant3->Gctrak()->upwght+0.5)-1]=13;
-
- if (geant3->Gctmed()->numed == idtmed[1001-1])
- sIloss[Int_t(geant3->Gctrak()->upwght+0.5)-1]=14;
-
- // CsI inefficiency
- if (geant3->Gctmed()->numed == idtmed[1006-1])
- sIloss[Int_t(geant3->Gctrak()->upwght+0.5)-1]=16;
-
-
- // Photon goes out of tracking scope
- } else if (geant3->Gctrak()->lmec[i] == 30)
- sIloss[Int_t(geant3->Gctrak()->upwght+0.5)-1]=21;
-
- }
+// one branch for rec clusters
+ for (i=0; i<7; i++) {
+ sprintf(branchname,"%sRecClus%d",GetName(),i+1);
+ if (fRecClusters && gAlice->TreeD()) {
+ gAlice->TreeR()
+ ->Branch(branchname,"TObjArray",
+ &((*fRecClusters)[i]), buffersize,0);
+ printf("Making Branch %s for clusters in chamber %d\n",
+ branchname,i+1);
}
- }
- }
+ }
}
-
-//_____________________________________________________________________________
-void AliRICH::RichIntegration()
+//___________________________________________
+void AliRICH::SetTreeAddress()
{
- //
- // Integrates over RICH pads
- //
-
- TGeant3 *geant3 = (TGeant3*) gMC;
-
- Int_t i1, i2;
- Float_t r1;
-
- // Local variables
- Float_t rrhh[25];
- Float_t qtot;
- Int_t ifeed;
- Float_t x0;
- Int_t ixmin, ixmax, iymin, iymax;
- Int_t ll, ix, iy;
- Float_t qp;
- Float_t source[3];
- Int_t irivol[2];
- Float_t y0a, qtot_check, xi1, xi2, yi1, yi2;
- Int_t nph = 0, iqp;
- Int_t ihitrak, modulen;
- //Int_t isec[MAXSEC];
- // ILOSS = 0 NOT LOST
- // 1 REFLECTED FREON-QUARZ
- // 2 REFLECTED QUARZ-METHANE
- // 3 REFLECTED METHANE-CSI
- // 11 ABSORBED IN FREON
- // 12 ABSORBED IN QUARZ
- // 13 ABSORBED IN METHANE
- // 21 CSI QE
- // IPROD = 1 PRODUCED IN FREON
- // IPROD = 2 PRODUCED IN QUARZ
- // get charge for MIP or cherenkov photon
-
- if (geant3->Gckine()->ipart == 50) {
- GetCharge(qtot);
- modulen = geant3->Gcvolu()->number[geant3->Gcvolu()->nlevel - 3];
- //sQint[1] = qtot;
- } else {
- GetChargeMip(qtot);
- modulen = geant3->Gcvolu()->number[geant3->Gcvolu()->nlevel - 4];
- //sQint[0] = qtot;
- }
- //
- gMC->Gmtod(sVectIn, sVloc, 1);
- if (TMath::Abs(sVloc[0]) >= sDlx) return;
-
- if (TMath::Abs(sVloc[2]) >= sDly) return;
-
- sVloc[0] += sDlx;
- x0 = sVloc[0];
- sVloc[2] += sDly;
- //y0 = sVloc[2];
- AnodicWires(y0a);
- ixmin = (Int_t) ((x0 - fSxcharge * 5.) / sDxp) + 1;
- if (ixmin < 1) ixmin = 1;
- ixmax = (Int_t) ((x0 + fSxcharge * 5.) / sDxp) + 1;
- if (ixmax > sNpx) ixmax = sNpx;
- iymin = (Int_t) ((y0a - sSycharge * 5.) / sDyp) + 1;
- if (iymin < 1) iymin = 1;
- iymax = (Int_t) ((y0a + sSycharge * 5.) / sDyp) + 1;
- if (iymax > sNpy) iymax = sNpy;
- qtot_check = 0.;
- i1 = ixmax;
- for (ix = ixmin; ix <= i1; ++ix) {
- i2 = iymax;
- for (iy = iymin; iy <= i2; ++iy) {
- xi1 = (sXpad[ix - 1] - x0) / sDpad;
- xi2 = xi1 + sDxp / sDpad;
- yi1 = (sYpad[iy - 1] - y0a) / sDpad;
- yi2 = yi1 + sDyp / sDpad;
- qp = qtot * FMathieson(xi1, xi2) * FMathieson(yi1, yi2);
- iqp = (Int_t) TMath::Abs(qp);
- qtot_check += TMath::Abs(qp);
-
- // FILL COMMON FOR PADS
-
- if (iqp) {
- if (iqp > 4095) {
- iqp = 4096;
- }
- ++sNpads;
- if (sNpads > MAXSEC) {
- // write(6,*) 'attention npads',npads
- sNpads = MAXSEC;
+ // Set branch address for the Hits and Digits Tree.
+ char branchname[20];
+ AliDetector::SetTreeAddress();
+
+ TBranch *branch;
+ TTree *treeH = gAlice->TreeH();
+ TTree *treeD = gAlice->TreeD();
+
+ if (treeH) {
+ if (fClusters) {
+ branch = treeH->GetBranch("RICHCluster");
+ if (branch) branch->SetAddress(&fClusters);
}
- sIpx[sNpads - 1] = ix;
- sIpy[sNpads - 1] = iy;
- sIqpad[sNpads - 1] = iqp;
- // TAG THE ORIGIN OF THE CHARGE DEPOSITION
- r1 = geant3->Gctrak()->upwght / 100.;
- ifeed = Int_t(r1+0.5);
- sNphlink[sNpads - 1] = 0;
- if (geant3->Gckine()->ipart != 50) {
- // MIP
- //isec[sNpads - 1] = sNsecon;
- } else {
- if (ifeed < 50) {
- // CERENKOV PHOTON
-
- nph = Int_t(geant3->Gctrak()->upwght+0.5);
- sNphlink[sNpads - 1] = nph;
- sXphit[nph - 1] = sVloc[0];
- sYphit[nph - 1] = sVloc[2];
- //isec[sNpads - 1] = sNsecon + 100;
- } else if (ifeed == 51) {
- // FEEDBACK FROM CERENKOV
-
- //isec[sNpads - 1] = sNsecon + 300;
- } else if (ifeed == 52) {
- // FEEDBACK FROM MIP
-
- //isec[sNpads - 1] = sNsecon + 200;
- }
+ if (fCerenkovs) {
+ branch = treeH->GetBranch("RICHCerenkov");
+ if (branch) branch->SetAddress(&fCerenkovs);
}
- // Generate the hit for Root IO
- for (ll = 0; ll < 25; ++ll) rrhh[ll] = 0;
- irivol[0] = modulen;
- irivol[1] = nmodsx;
- rrhh[0] = 0.;
- // Flag to say this is a hit
- rrhh[1] = xsx;
- rrhh[2] = ysx;
- rrhh[3] = psx;
- rrhh[4] = thetasx;
- rrhh[5] = phisx;
- rrhh[6] = (Float_t) idpartgx;
- rrhh[7] = sXsox;
- rrhh[8] = sYsox;
- rrhh[9] = sZsox;
- rrhh[10] = (Float_t) sIpx[sNpads - 1];
- rrhh[11] = (Float_t) sIpy[sNpads - 1];
- rrhh[12] = (Float_t) sIqpad[sNpads - 1];
- ihitrak = gAlice->CurrentTrack();
- if (sNphlink[sNpads - 1] > 0) {
- rrhh[13] = sPparent;
- rrhh[14] = sThincParent;
- rrhh[15] = (Float_t) sIloss[nph - 1];
- rrhh[16] = (Float_t) sIprod[nph - 1];
- rrhh[17] = sXphit[nph - 1];
- rrhh[18] = sYphit[nph - 1];
- ihitrak = sMckov[nph - 1];
+
+ }
+
+ if (treeD) {
+ for (int i=0; i<7; i++) {
+ sprintf(branchname,"%sDigits%d",GetName(),i+1);
+ if (fDchambers) {
+ branch = treeD->GetBranch(branchname);
+ if (branch) branch->SetAddress(&((*fDchambers)[i]));
+ }
}
- // This is the current position of the hit
- rrhh[19] = geant3->Gctrak()->vect[0];
- rrhh[20] = geant3->Gctrak()->vect[1];
- rrhh[21] = geant3->Gctrak()->vect[2];
- // PRINT *,' '
- // PRINT *,'RXAHIT(oldhit)=[',RRHH(20),',',RRHH(21),',',
- // + RRHH(22),']'
- // PRINT *, 'ITRA = ',ITRA,'ISTAK = ',ISTAK
- AddHit(ihitrak, irivol, rrhh);
- // new Padhits
- for (ll = 0; ll < 25; ++ll) rrhh[ll] = 0;
- rrhh[0] = 3.;
- rrhh[1] = (Float_t) geant3->Gckine()->ipart;
- rrhh[2] = (Float_t) sIpx[sNpads - 1];
- rrhh[3] = (Float_t) sIpy[sNpads - 1];
- rrhh[4] = (Float_t) modulen;
- rrhh[5] = -1.;
- // !!!Prod
- rrhh[6] = (Float_t) sIqpad[sNpads - 1];
- AddHit(ihitrak, irivol, rrhh);
- // enew
- }
}
- }
- // IF (IPART .NE. 50) WRITE(6,*) 'CHECK',QTOT,QTOT_CHECK
-
- // GENERATE FEEDBACK PHOTONS
-
- source[0] = x0 - sDlx;
- source[1] = sVloc[1] - .2;
- source[2] = y0a - sDly;
- gMC->Gdtom(source, source, 1);
- FeedBack(source, qtot);
- return;
+}
+//___________________________________________
+void AliRICH::ResetHits()
+{
+ // Reset number of clusters and the cluster array for this detector
+ AliDetector::ResetHits();
+ fNclusters = 0;
+ if (fClusters) fClusters->Clear();
+ if (fCerenkovs) fCerenkovs->Clear();
}
-//_____________________________________________________________________________
-void AliRICH::AnodicWires(Float_t &y0a)
+//____________________________________________
+void AliRICH::ResetDigits()
+{
+ //
+ // Reset number of digits and the digits array for this detector
+ //
+ for ( int i=0;i<7;i++ ) {
+ if ((*fDchambers)[i]) (*fDchambers)[i]->Clear();
+ if (fNdch) fNdch[i]=0;
+ }
+}
+//____________________________________________
+void AliRICH::ResetRecClusters()
{
- //
- // Position of anodic wires
- //
- Int_t i1;
-
- // Local variables
- Int_t i;
- Float_t ass_i, y0, ass_i1;
-
- y0 = sVloc[2];
- y0a = -1.;
- i1 = (sNpy << 1) + 1;
- for (i = 1; i <= i1; ++i) {
- if (y0 > sYanode[i - 1] && y0 <= sYanode[i]) {
- ass_i1 = TMath::Abs(sYanode[i] - y0);
- ass_i = TMath::Abs(sYanode[i - 1] - y0);
- if (ass_i1 <= ass_i) y0a = sYanode[i];
- if (ass_i1 > ass_i) y0a = sYanode[i - 1];
- return;
+ //
+ // Reset the rec clusters
+ //
+ for ( int i=0;i<7;i++ ) {
+ if ((*fRecClusters)[i]) (*fRecClusters)[i]->Clear();
}
- }
-}
+}
+//___________________________________________
-//_____________________________________________________________________________
-void AliRICH::GetChargeMip(Float_t &qtot)
+void AliRICH::SetPADSIZ(Int_t id, Int_t isec, Float_t p1, Float_t p2)
{
- //
- // Calculates the charge deposited by a MIP
- //
-
-
- Int_t i1;
-
- // Local variables
- Float_t ranf[2];
- Int_t i;
- Int_t nel;
-
- // NUMBER OF ELECTRONS
-
- nel = (Int_t) (sDetot * 1e9 / 26.);
- qtot = 0.;
- if (!nel) nel = 1;
- i1 = nel;
- for (i = 1; i <= i1; ++i) {
- gMC->Rndm(ranf, 1);
- qtot -= fChslope * TMath::Log(ranf[0]);
- }
-}
+ Int_t i=2*(id-1);
+ ((AliRICHchamber*) (*fChambers)[i]) ->SetPADSIZ(isec,p1,p2);
+ ((AliRICHchamber*) (*fChambers)[i+1])->SetPADSIZ(isec,p1,p2);
+}
-//_____________________________________________________________________________
-void AliRICH::GetCharge(Float_t &qtot)
+//___________________________________________
+void AliRICH::SetMUCHSP(Int_t id, Float_t p1)
{
- //
- // Charge deposited
- //
+ Int_t i=2*(id-1);
+ ((AliRICHchamber*) (*fChambers)[i])->SetMUCHSP(p1);
+ ((AliRICHchamber*) (*fChambers)[i+1])->SetMUCHSP(p1);
+}
+//___________________________________________
+void AliRICH::SetMUSIGM(Int_t id, Float_t p1, Float_t p2)
+{
+ Int_t i=2*(id-1);
+ ((AliRICHchamber*) (*fChambers)[i])->SetMUSIGM(p1,p2);
+ ((AliRICHchamber*) (*fChambers)[i+1])->SetMUSIGM(p1,p2);
+}
- Float_t ranf[1];
-
- gMC->Rndm(ranf, 1);
- qtot = -fChslope * TMath::Log(ranf[0]);
-}
+//___________________________________________
+void AliRICH::SetRSIGM(Int_t id, Float_t p1)
+{
+ Int_t i=2*(id-1);
+ ((AliRICHchamber*) (*fChambers)[i])->SetRSIGM(p1);
+ ((AliRICHchamber*) (*fChambers)[i+1])->SetRSIGM(p1);
+}
-//_____________________________________________________________________________
-void AliRICH::FeedBack(Float_t *source, Float_t qtot)
+//___________________________________________
+void AliRICH::SetMAXADC(Int_t id, Float_t p1)
{
- //
- // Generate FeedBack photons
- //
-
- TGeant3 *geant3 = (TGeant3*) gMC;
-
- Int_t i1, j;
- Float_t r1, r2;
-
- // Local variables
- Float_t cthf, ranf[2], phif, enfp = 0, sthf;
- Int_t i, ifeed;
- Float_t e1[3], e2[3], e3[3];
- Float_t vmod, uswop;
- Float_t fp, random;
- Float_t dir[3], phi;
- Int_t nfp;
- Float_t pol[3], supwght;
-
- // DETERMINE NUMBER OF FEEDBACK PHOTONS
-
- // Function Body
- fp = fAlphaFeed * qtot;
- nfp = gRandom->Poisson(fp);
-
- // GENERATE PHOTONS
-
- geant3->Gckin2()->ngphot = 0;
- i1 = nfp;
- for (i = 1; i <= i1; ++i) {
-
- // DIRECTION
- gMC->Rndm(ranf, 2);
- cthf = ranf[0] * 2. - 1.;
- if (cthf < 0.) continue;
- sthf = TMath::Sqrt((1. - cthf) * (cthf + 1.));
- phif = ranf[1] * 2 * TMath::Pi();
- ++geant3->Gckin2()->ngphot;
- if (geant3->Gckin2()->ngphot > 800) {
- printf("ATTENTION NGPHOT ! %d %f %d\n",
- geant3->Gckin2()->ngphot,fp,nfp);
- break;
- }
-
- // POSITION
- for(j=0;j<3;j++)
- geant3->Gckin2()->xphot[geant3->Gckin2()->ngphot-1][j]=source[j];
-
- // ENERGY
- gMC->Rndm(&random, 1);
- if (random <= .57) {
- enfp = 7.5e-9;
- } else if (random > .57 && random <= .7) {
- enfp = 6.4e-9;
- } else if (random > .7) {
- enfp = 7.9e-9;
- }
- geant3->Gckin2()->xphot[geant3->Gckin2()->ngphot-1][6] = enfp;
- dir[0] = sthf * TMath::Sin(phif);
- dir[1] = cthf;
- dir[2] = sthf * TMath::Cos(phif);
- gMC->Gdtom(dir, &geant3->Gckin2()->xphot[geant3->Gckin2()->ngphot-1][3], 2);
-
- // POLARISATION
- e1[0] = 0.;
- e1[1] = -dir[2];
- e1[2] = dir[1];
-
- e2[0] = -dir[1];
- e2[1] = dir[0];
- e2[2] = 0.;
-
- e3[0] = dir[1];
- e3[1] = 0.;
- e3[2] = -dir[0];
-
- vmod=0;
- for(j=0;j<3;j++) vmod+=e1[j]*e1[j];
- if (!vmod) for(j=0;j<3;j++) {
- uswop=e1[j];
- e1[j]=e3[j];
- e3[j]=uswop;
- }
- vmod=0;
- for(j=0;j<3;j++) vmod+=e2[j]*e2[j];
- if (!vmod) for(j=0;j<3;j++) {
- uswop=e2[j];
- e2[j]=e3[j];
- e3[j]=uswop;
- }
-
- vmod=0;
- for(j=0;j<3;j++) vmod+=e1[j]*e1[j];
- vmod=TMath::Sqrt(1/vmod);
- for(j=0;j<3;j++) e1[j]*=vmod;
-
- vmod=0;
- for(j=0;j<3;j++) vmod+=e2[j]*e2[j];
- vmod=TMath::Sqrt(1/vmod);
- for(j=0;j<3;j++) e2[j]*=vmod;
-
- gMC->Rndm(ranf, 1);
- phi = ranf[0] * 2 * TMath::Pi();
- r1 = TMath::Sin(phi);
- r2 = TMath::Cos(phi);
- for(j=0;j<3;j++) pol[j]=e1[j]*r1+e2[j]*r2;
- gMC->Gdtom(pol, &geant3->Gckin2()->xphot[geant3->Gckin2()->ngphot-1][7], 2);
-
- // TIME OF FLIGHT
- geant3->Gckin2()->xphot[geant3->Gckin2()->ngphot-1][10] = 0.;
-
- // PUT PHOTON ON THE STACK AND LABEL IT AS FEEDBACK (51,52)
- supwght = geant3->Gctrak()->upwght;
- r1 = geant3->Gctrak()->upwght / 100.;
- ifeed = Int_t(r1+0.5);
- ++sNfeed;
- if (geant3->Gckine()->ipart == 50 && ifeed != 52) {
- geant3->Gctrak()->upwght = 5100.;
- } else {
- geant3->Gctrak()->upwght = 5200.;
- }
- // geant3->Gskpho(geant3->Gckin2()->ngphot);
- geant3->Gctrak()->upwght = supwght;
-
- }
- geant3->Gckin2()->ngphot = 0;
+ Int_t i=2*(id-1);
+ ((AliRICHchamber*) (*fChambers)[i])->SetMAXADC(p1);
+ ((AliRICHchamber*) (*fChambers)[i+1])->SetMAXADC(p1);
}
-//_____________________________________________________________________________
-Float_t AliRICH::FMathieson(Float_t lambda1, Float_t lambda2)
+//___________________________________________
+void AliRICH::SetSMAXAR(Float_t p1)
{
- //
- // Mathieson charge distribution
- //
- // CALCULATES INTEGRAL OF GATTI/MATHIESON CHARGE DISTRIBUTION
- // FOR CPC AND CSC CHAMBERS
- // INTEGRATION LIMITS LAMBDA1,LAMBDA2= X/D WHERE D IS THE ANODE CATHODE
- // SEPARATION
- // K3: GEOMETRY PARAMETER
- //
- Float_t u1, u2;
- // const Float_t k1=0.2828;
- const Float_t k2=0.962;
- const Float_t k3=0.6;
- const Float_t k4=0.379;
- const Float_t sqrtk3=TMath::Sqrt(k3);
-
-
- u1 = tanh(lambda1 * k2) * sqrtk3;
- u2 = tanh(lambda2 * k2) * sqrtk3;
-
- return (atan(u2) - atan(u1)) * 2 * k4;
-
+ fMaxStepGas=p1;
}
+//___________________________________________
+void AliRICH::SetSMAXAL(Float_t p1)
+{
+ fMaxStepAlu=p1;
+}
-//_____________________________________________________________________________
-void AliRICH::UpdateMipHit(Float_t *hits)
+//___________________________________________
+void AliRICH::SetDMAXAR(Float_t p1)
{
- //
- // Update some parameters when the current mip hits the detector.
- //
- TClonesArray &lhits = *fMips;
- AliRICHmip *lmip = (AliRICHmip *) lhits.AddrAt(fNmips-1);
- lmip->SetX ( hits[1]);
- lmip->SetY ( hits[2]);
- lmip->SetModule((int) hits[3]);
- lmip->SetTheta ( hits[4]);
- lmip->SetPhi ( hits[5]);
+ fMaxDestepGas=p1;
}
-//*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
+//___________________________________________
+void AliRICH::SetDMAXAL(Float_t p1)
+{
+ fMaxDestepAlu=p1;
+}
+//___________________________________________
+void AliRICH::SetRICHACC(Bool_t acc, Float_t angmin, Float_t angmax)
+{
+ fAccCut=acc;
+ fAccMin=angmin;
+ fAccMax=angmax;
+}
+//___________________________________________
+void AliRICH::SetSegmentationModel(Int_t id, Int_t isec, AliRICHsegmentation *segmentation)
+{
+ ((AliRICHchamber*) (*fChambers)[id])->SegmentationModel(isec, segmentation);
-//_____________________________________________________________________________
-void AliRICH::MakeBranch(Option_t *option)
+}
+//___________________________________________
+void AliRICH::SetResponseModel(Int_t id, Response_t res, AliRICHresponse *response)
{
- //
- // Create a new branch in the current Root Tree.
- // The branch of fHits is automatically split
- //
- AliDetector::MakeBranch(option);
-
- Int_t buffersize = 4000;
- if (gAlice->TreeH()) {
- if(fMips){
- gAlice->TreeH()->Branch("RICHmip",&fMips, buffersize);
- printf("Making Branch %s for mips\n","RICHmip");
- }
- if(fCkovs){
- gAlice->TreeH()->Branch("RICHckov",&fCkovs, buffersize);
- printf("Making Branch %s for ckovs\n","RICHckov");
- }
- if(fPadhits){
- gAlice->TreeH()->Branch("RICHpadhit",&fPadhits, buffersize);
- printf("Making Branch %s for padhits\n","RICHpadhit");
- }
- }
+ ((AliRICHchamber*) (*fChambers)[id])->ResponseModel(res, response);
}
-//_____________________________________________________________________________
-void AliRICH::SetTreeAddress()
+void AliRICH::SetNsec(Int_t id, Int_t nsec)
+{
+ ((AliRICHchamber*) (*fChambers)[id])->SetNsec(nsec);
+}
+
+
+//___________________________________________
+
+void AliRICH::StepManager()
{
- //
- // Set branch address for the Hits and Digits Tree.
- //
- AliDetector::SetTreeAddress();
- TBranch *branch;
- TTree *treeH = gAlice->TreeH();
- if (treeH) {
- if (fMips) {
- branch = treeH->GetBranch("RICHmip");
- if (branch) branch->SetAddress(&fMips);
+ printf("Dummy version of RICH step -- it should never happen!!\n");
+ const Float_t kRaddeg = 180/TMath::Pi();
+ AliMC* pMC = AliMC::GetMC();
+ Int_t nsec, ipart;
+ Float_t x[4], p[4];
+ Float_t pt, th0, th1;
+ char proc[5];
+ if(fAccCut) {
+ if((nsec=pMC->NSecondaries())>0) {
+ pMC->ProdProcess(proc);
+ if((pMC->TrackPid()==113 || pMC->TrackPid()==114) && !strcmp(proc,"DCAY")) {
+
+ // Check angular acceptance
+ //* --- and have muons from resonance decays in the wanted window ---
+ if(nsec != 2) {
+ printf(" AliRICH::StepManager: Strange resonance Decay into %d particles\n",nsec);
+ pMC->StopEvent();
+ } else {
+ pMC->GetSecondary(0,ipart,x,p);
+ pt = TMath::Sqrt(p[0]*p[0]+p[1]*p[1]);
+ th0 = TMath::ATan2(pt,p[2])*kRaddeg;
+ pMC->GetSecondary(1,ipart,x,p);
+ pt = TMath::Sqrt(p[0]*p[0]+p[1]*p[1]);
+ th1 = TMath::ATan2(pt,p[2])*kRaddeg;
+ if(!(fAccMin < th0 && th0 < fAccMax) ||
+ !(fAccMin < th1 && th1 < fAccMax))
+ pMC->StopEvent();
+ }
+ }
+ }
}
- if (fCkovs) {
- branch = treeH->GetBranch("RICHckov");
- if (branch) branch->SetAddress(&fCkovs);
+}
+void AliRICH::ReconstructClusters()
+{
+//
+// Initialize the necessary correspondance table
+//
+ static const Int_t kMaxNpadx = 600;
+ static const Int_t kMaxNpady = 600;
+ Int_t elem[kMaxNpadx*2][kMaxNpady*2];
+//
+// Loop on chambers and on cathode planes
+//
+ for (Int_t ich=0;ich<7;ich++)
+ for (Int_t icat=0;icat<1;icat++) {
+ //
+ // Get ready the current chamber stuff
+ //
+
+ printf ("Olarilole");
+ AliRICHchamber* iChamber= &(this->Chamber(ich));
+ AliRICHsegmentation* segmentation =
+ iChamber->GetSegmentationModel(icat+1);
+ if (!segmentation)
+ continue;
+ TClonesArray *RICHdigits = this->DigitsAddress(ich);
+ if (RICHdigits == 0)
+ continue;
+ cout << "Npx " << segmentation->Npx()
+ << " Npy " << segmentation->Npy() << endl;
+ //
+ // Ready the digits
+ //
+ gAlice->ResetDigits();
+ gAlice->TreeD()->GetEvent(icat+1); // spurious +1 ...
+ Int_t ndigits = RICHdigits->GetEntriesFast();
+ if (ndigits == 0)
+ continue;
+ printf("Found %d digits for cathode %d in chamber %d \n",
+ ndigits,icat,ich+1);
+ AliRICHdigit *mdig;
+ AliRICHRecCluster *Cluster;
+ //
+ // Build the correspondance table
+ //
+ memset(elem,0,sizeof(Int_t)*kMaxNpadx*kMaxNpady*4);
+ Int_t digit;
+ for (digit=0; digit<ndigits; digit++)
+ {
+ mdig = (AliRICHdigit*)RICHdigits->UncheckedAt(digit);
+ elem[kMaxNpadx+mdig->fPadX][kMaxNpady+mdig->fPadY] = digit+1;
+ // because default is 0
+ }
+ //
+ // Declare some useful variables
+ //
+ Int_t Xlist[10];
+ Int_t Ylist[10];
+ Int_t Nlist;
+ Int_t nclust=0;
+ //
+ // loop over digits
+ //
+ for (digit=0;digit<ndigits;digit++) {
+ mdig = (AliRICHdigit*)RICHdigits->UncheckedAt(digit);
+ //
+ // if digit still available, start clustering
+ //
+ if (elem[kMaxNpadx+mdig->fPadX][kMaxNpady+mdig->fPadY]) {
+ Cluster = new AliRICHRecCluster(digit, ich, icat);
+ elem[kMaxNpadx+mdig->fPadX][kMaxNpady+mdig->fPadY]=0;
+ //
+ // loop over the current list of digits
+ // and look for neighbours
+ //
+ for(Int_t clusDigit=Cluster->FirstDigitIndex();
+ clusDigit!=Cluster->InvalidDigitIndex();
+ clusDigit=Cluster->NextDigitIndex()) {
+ AliRICHdigit* pDigit=(AliRICHdigit*)RICHdigits
+ ->UncheckedAt(clusDigit);
+ segmentation->Neighbours(pDigit->fPadX,pDigit->fPadY,
+ &Nlist, Xlist, Ylist);
+ for (Int_t Ilist=0;Ilist<Nlist;Ilist++) {
+ if (elem[kMaxNpadx+Xlist[Ilist]][kMaxNpady
+ +Ylist[Ilist]]) {
+ //
+ // Add the digit at the end and reset the table
+ //
+ Cluster->AddDigit(elem[kMaxNpadx+Xlist[Ilist]][kMaxNpady+Ylist[Ilist]]-1);
+ elem[kMaxNpadx+Xlist[Ilist]][kMaxNpady
+ +Ylist[Ilist]] =0;
+ } // if elem
+ } // for Ilist
+ } // for pDigit
+ //
+ // Store the cluster (good time to do Cluster polishing)
+ //
+ segmentation->FitXY(Cluster,RICHdigits);
+ nclust ++;
+ AddRecCluster(ich,icat,Cluster);
+ }
+ }
+ printf("===> %d Clusters\n",nclust);
+ } // for icat
+}
+
+
+//______________________________________________________________________________
+void AliRICH::Streamer(TBuffer &R__b)
+{
+ // Stream an object of class AliRICH.
+ AliRICHchamber *iChamber;
+ AliRICHsegmentation *segmentation;
+ AliRICHresponse *response;
+ TClonesArray *digitsaddress;
+
+ if (R__b.IsReading()) {
+ Version_t R__v = R__b.ReadVersion(); if (R__v) { }
+ AliDetector::Streamer(R__b);
+ R__b >> fNclusters;
+ R__b >> fClusters; // diff
+ R__b >> fDchambers;
+ R__b.ReadArray(fNdch);
+ //
+ R__b >> fAccCut;
+ R__b >> fAccMin;
+ R__b >> fAccMax;
+ //
+ R__b >> fChambers;
+// Stream chamber related information
+ for (Int_t i =0; i<7; i++) {
+ iChamber=(AliRICHchamber*) (*fChambers)[i];
+ iChamber->Streamer(R__b);
+ if (iChamber->Nsec()==1) {
+ segmentation=iChamber->GetSegmentationModel(1);
+ segmentation->Streamer(R__b);
+ } else {
+ segmentation=iChamber->GetSegmentationModel(1);
+ segmentation->Streamer(R__b);
+ segmentation=iChamber->GetSegmentationModel(2);
+ segmentation->Streamer(R__b);
+ }
+ response=iChamber->GetResponseModel(mip);
+ response->Streamer(R__b);
+ response=iChamber->GetResponseModel(cerenkov);
+ response->Streamer(R__b);
+
+ digitsaddress=(TClonesArray*) (*fDchambers)[i];
+ digitsaddress->Streamer(R__b);
+ }
+
+ } else {
+ R__b.WriteVersion(AliRICH::IsA());
+ AliDetector::Streamer(R__b);
+ R__b << fNclusters;
+ R__b << fClusters; // diff
+ R__b << fDchambers;
+ R__b.WriteArray(fNdch, 7);
+ //
+ R__b << fAccCut;
+ R__b << fAccMin;
+ R__b << fAccMax;
+ //
+ R__b << fChambers;
+// Stream chamber related information
+ for (Int_t i =0; i<7; i++) {
+ iChamber=(AliRICHchamber*) (*fChambers)[i];
+ iChamber->Streamer(R__b);
+ if (iChamber->Nsec()==1) {
+ segmentation=iChamber->GetSegmentationModel(1);
+ segmentation->Streamer(R__b);
+ } else {
+ segmentation=iChamber->GetSegmentationModel(1);
+ segmentation->Streamer(R__b);
+ segmentation=iChamber->GetSegmentationModel(2);
+ segmentation->Streamer(R__b);
+ }
+ response=iChamber->GetResponseModel(mip);
+ response->Streamer(R__b);
+ response=iChamber->GetResponseModel(cerenkov);
+ response->Streamer(R__b);
+
+ digitsaddress=(TClonesArray*) (*fDchambers)[i];
+ digitsaddress->Streamer(R__b);
+ }
}
- if (fPadhits) {
- branch = treeH->GetBranch("RICHpadhit");
- if (branch) branch->SetAddress(&fPadhits);
+}
+AliRICHcluster* AliRICH::FirstPad(AliRICHhit* hit,TClonesArray *clusters )
+{
+//
+ // Initialise the pad iterator
+ // Return the address of the first padhit for hit
+ TClonesArray *theClusters = clusters;
+ Int_t nclust = theClusters->GetEntriesFast();
+ if (nclust && hit->fPHlast > 0) {
+ sMaxIterPad=hit->fPHlast;
+ sCurIterPad=hit->fPHfirst;
+ return (AliRICHcluster*) clusters->UncheckedAt(sCurIterPad-1);
+ } else {
+ return 0;
}
- }
+
}
-//_____________________________________________________________________________
-void AliRICH::ResetHits()
+AliRICHcluster* AliRICH::NextPad(TClonesArray *clusters)
{
- //
- // Reset number of mips and the mips array for RICH
- //
- AliDetector::ResetHits();
- fNmips = 0;
- if (fMips) fMips->Clear();
- // Reset number of Ckovs and the Ckovs array for RICH
- fNckovs = 0;
- if (fCkovs) fCkovs->Clear();
- // Reset number of Padhits and the Padhits array for RICH
- fNpadhits = 0;
- if (fPadhits) fPadhits->Clear();
+ sCurIterPad++;
+ if (sCurIterPad <= sMaxIterPad) {
+ return (AliRICHcluster*) clusters->UncheckedAt(sCurIterPad-1);
+ } else {
+ return 0;
+ }
}
+ClassImp(AliRICHcluster)
+
+//___________________________________________
+AliRICHcluster::AliRICHcluster(Int_t *clhits)
+{
+ fHitNumber=clhits[0];
+ fCathode=clhits[1];
+ fQ=clhits[2];
+ fPadX=clhits[3];
+ fPadY=clhits[4];
+ fQpad=clhits[5];
+ fRSec=clhits[6];
+}
+ClassImp(AliRICHdigit)
//_____________________________________________________________________________
-Int_t AliRICH::DistancetoPrimitive(Int_t , Int_t )
+AliRICHdigit::AliRICHdigit(Int_t *digits)
{
- //
- // Distance from mouse RICH on the display
- // Dummy routine
- //
- return 9999;
+ //
+ // Creates a RICH digit object to be updated
+ //
+ fPadX = digits[0];
+ fPadY = digits[1];
+ fSignal = digits[2];
+
}
-
//_____________________________________________________________________________
-void AliRICH::PreTrack()
+AliRICHdigit::AliRICHdigit(Int_t *tracks, Int_t *charges, Int_t *digits)
{
- //
- // Called before transporting a track
- //
- sDetot=0;
- sNsecon=0;
- sNpads=0;
- sNphoton=0;
- sNfeed=0;
- sNfeedd=0;
- sNdir=0;
+ //
+ // Creates a RICH digit object
+ //
+ fPadX = digits[0];
+ fPadY = digits[1];
+ fSignal = digits[2];
+ for(Int_t i=0; i<10; i++) {
+ fTcharges[i] = charges[i];
+ fTracks[i] = tracks[i];
+ }
}
-//_____________________________________________________________________________
-void AliRICH::Streamer(TBuffer &R__b)
+ClassImp(AliRICHlist)
+
+//____________________________________________________________________________
+AliRICHlist::AliRICHlist(Int_t ich, Int_t *digits):
+ AliRICHdigit(digits)
{
- //
- // Stream an object of class AliRICH.
- //
- if (R__b.IsReading()) {
- Version_t R__v = R__b.ReadVersion(); if (R__v) { }
- AliDetector::Streamer(R__b);
- R__b >> fNmips;
- R__b >> fNckovs;
- R__b >> fNpadhits;
- R__b >> fMips; //diff
- R__b >> fCkovs; //diff
- R__b >> fPadhits; //diff
- R__b >> fChslope;
- R__b >> fAlphaFeed;
- R__b >> fSxcharge;
- R__b >> fIritri;
- } else {
- R__b.WriteVersion(AliRICH::IsA());
- AliDetector::Streamer(R__b);
- R__b << fNmips;
- R__b << fNckovs;
- R__b << fNpadhits;
- R__b << fMips; //diff
- R__b << fCkovs; //diff
- R__b << fPadhits; //diff
- R__b << fChslope;
- R__b << fAlphaFeed;
- R__b << fSxcharge;
- R__b << fIritri;
- }
+ //
+ // Creates a RICH digit list object
+ //
+
+ fChamber = ich;
+ fTrackList = new TObjArray;
+
}
+//_____________________________________________________________________________
-ClassImp(AliRICHv1)
-///////////////////////////////////////////////////////////////////////////////
-// //
-// Ring Imaging Cherenkov //
-// This class contains version 1 of the Ring Imaging Cherenkov //
-// //
-//Begin_Html
-/*
-<img src="picts/AliRICHv1Class.gif">
-*/
-//End_Html
-// //
-///////////////////////////////////////////////////////////////////////////////
-
-//_____________________________________________________________________________
-AliRICHv1::AliRICHv1() : AliRICH()
+ClassImp(AliRICHhit)
+
+//___________________________________________
+AliRICHhit::AliRICHhit(Int_t shunt, Int_t track, Int_t *vol, Float_t *hits):
+ AliHit(shunt, track)
{
- //
- // Default Constructor RICH for version 1
- //
+ fChamber=vol[0];
+ fParticle=hits[0];
+ fX=hits[1];
+ fY=hits[2];
+ fZ=hits[3];
+ fTheta=hits[4];
+ fPhi=hits[5];
+ fTlength=hits[6];
+ fEloss=hits[7];
+ fPHfirst=(Int_t) hits[8];
+ fPHlast=(Int_t) hits[9];
}
-
-//_____________________________________________________________________________
-AliRICHv1::AliRICHv1(const char *name, const char *title)
- : AliRICH(name,title)
+ClassImp(AliRICHreccluster)
+
+ClassImp(AliRICHCerenkov)
+//___________________________________________
+AliRICHCerenkov::AliRICHCerenkov(Int_t shunt, Int_t track, Int_t *vol, Float_t *hits):
+ AliHit(shunt, track)
{
- //
- // Standard constructor for RICH version 1
- //
+ fChamber=vol[0];
+ fX=hits[1];
+ fY=hits[2];
+ fZ=hits[3];
+ fTheta=hits[4];
+ fPhi=hits[5];
+ fTlength=hits[6];
+ fPHfirst=(Int_t) hits[8];
+ fPHlast=(Int_t) hits[9];
}
-//_____________________________________________________________________________
-AliRICHv1::~AliRICHv1()
+ClassImp(AliRICHRecCluster)
+
+//_____________________________________________________________________
+AliRICHRecCluster::AliRICHRecCluster()
{
- //
- // Standard destructor for RICH version 1
- //
+ fDigits=0;
+ fNdigit=-1;
}
-
-//_____________________________________________________________________________
-void AliRICHv1::CreateGeometry()
+
+AliRICHRecCluster::AliRICHRecCluster(Int_t FirstDigit,Int_t Ichamber, Int_t Icathod)
{
- //
- // Create the geometry for RICH version 1
- //
- // Modified by: N. Colonna (INFN - BARI, Nicola.Colonna@ba.infn.it)
- // R.A. Fini (INFN - BARI, Rosanna.Fini@ba.infn.it)
- // R.A. Loconsole (Bari University, loco@riscom.ba.infn.it)
- //
- //Begin_Html
- /*
- <img src="picts/AliRICHv1.gif">
- */
- //End_Html
- //Begin_Html
- /*
- <img src="picts/AliRICHv1Tree.gif">
- */
- //End_Html
-
-
- Int_t *idtmed = fIdtmed->GetArray()-999;
-
- Int_t i;
- Float_t zs;
- Int_t idrotm[1099];
- Float_t par[3];
-
- // --- Define the RICH detector
- // External aluminium box
- par[0] = 71.1;
- par[1] = 11.5;
- par[2] = 73.15;
- gMC->Gsvolu("RICH", "BOX ", idtmed[1009], par, 3);
-
- // Sensitive part of the whole RICH
- par[0] = 64.8;
- par[1] = 11.5;
- par[2] = 66.55;
- gMC->Gsvolu("SRIC", "BOX ", idtmed[1000], par, 3);
-
- // Honeycomb
- par[0] = 63.1;
- par[1] = .188;
- par[2] = 66.55;
- gMC->Gsvolu("HONE", "BOX ", idtmed[1001], par, 3);
-
- // Aluminium sheet
- par[0] = 63.1;
- par[1] = .025;
- par[2] = 66.55;
- gMC->Gsvolu("ALUM", "BOX ", idtmed[1009], par, 3);
-
- // Quartz
- par[0] = 63.1;
- par[1] = .25;
- par[2] = 65.5;
- gMC->Gsvolu("QUAR", "BOX ", idtmed[1002], par, 3);
-
- // Spacers (cylinders)
- par[0] = 0.;
- par[1] = .5;
- par[2] = .5;
- gMC->Gsvolu("SPAC", "TUBE", idtmed[1002], par, 3);
-
- // Opaque quartz
- par[0] = 61.95;
- par[1] = .2;
- par[2] = 66.5;
- gMC->Gsvolu("OQUA", "BOX ", idtmed[1007], par, 3);
-
- // Frame of opaque quartz
- par[0] = 20.65;
- par[1] = .5;
- par[2] = 66.5;
- gMC->Gsvolu("OQUF", "BOX ", idtmed[1007], par, 3);
-
- // Little bar of opaque quartz
- par[0] = 63.1;
- par[1] = .25;
- par[2] = .275;
- gMC->Gsvolu("BARR", "BOX ", idtmed[1007], par, 3);
-
- // Freon
- par[0] = 20.15;
- par[1] = .5;
- par[2] = 65.5;
- gMC->Gsvolu("FREO", "BOX ", idtmed[1003], par, 3);
-
- // Methane
- par[0] = 64.8;
- par[1] = 5.;
- par[2] = 64.8;
- gMC->Gsvolu("META", "BOX ", idtmed[1004], par, 3);
-
- // Methane gap
- par[0] = 64.8;
- par[1] = .2;
- par[2] = 64.8;
- gMC->Gsvolu("GAP ", "BOX ", idtmed[1008], par, 3);
-
- // CsI photocathode
- par[0] = 64.8;
- par[1] = .25;
- par[2] = 64.8;
- gMC->Gsvolu("CSI ", "BOX ", idtmed[1005], par, 3);
-
- // Anode grid
- par[0] = 0.;
- par[1] = .0025;
- par[2] = 20.;
- gMC->Gsvolu("GRID", "TUBE", idtmed[1006], par, 3);
-
- // --- Places the detectors defined with GSVOLU
- // Place material inside RICH
- gMC->Gspos("SRIC", 1, "RICH", 0., 0., 0., 0, "ONLY");
-
- gMC->Gspos("ALUM", 1, "SRIC", 0., -6.075, 0., 0, "ONLY");
- gMC->Gspos("HONE", 1, "SRIC", 0., -5.862, 0., 0, "ONLY");
- gMC->Gspos("ALUM", 2, "SRIC", 0., -5.649, 0., 0, "ONLY");
- gMC->Gspos("OQUA", 1, "SRIC", 0., -5.424, 0., 0, "ONLY");
-
- AliMatrix(idrotm[1019], 0., 0., 90., 0., 90., 90.);
-
- for (i = 1; i <= 9; ++i) {
- zs = (5 - i) * 14.4;
- gMC->Gspos("SPAC", i, "FREO", 6.7, 0., zs, idrotm[1019], "ONLY");
- }
- for (i = 10; i <= 18; ++i) {
- zs = (14 - i) * 14.4;
- gMC->Gspos("SPAC", i, "FREO", -6.7, 0., zs, idrotm[1019], "ONLY");
- }
-
- gMC->Gspos("FREO", 1, "OQUF", 0., 0., 0., 0, "ONLY");
- gMC->Gspos("OQUF", 1, "SRIC", 41.3, -4.724, 0., 0, "ONLY");
- gMC->Gspos("OQUF", 2, "SRIC", 0., -4.724, 0., 0, "ONLY");
- gMC->Gspos("OQUF", 3, "SRIC", -41.3, -4.724, 0., 0, "ONLY");
- gMC->Gspos("BARR", 1, "QUAR", 0., 0., -21.65, 0, "ONLY");
- gMC->Gspos("BARR", 2, "QUAR", 0., 0., 21.65, 0, "ONLY");
- gMC->Gspos("QUAR", 1, "SRIC", 0., -3.974, 0., 0, "ONLY");
- gMC->Gspos("GAP ", 1, "META", 0., 4.8, 0., 0, "ONLY");
- gMC->Gspos("META", 1, "SRIC", 0., 1.276, 0., 0, "ONLY");
- gMC->Gspos("CSI ", 1, "SRIC", 0., 6.526, 0., 0, "ONLY");
-
- // Place RICH inside ALICE apparatus
-
- AliMatrix(idrotm[1000], 90., 0., 70.69, 90., 19.31, -90.);
- AliMatrix(idrotm[1001], 90., -20., 90., 70., 0., 0.);
- AliMatrix(idrotm[1002], 90., 0., 90., 90., 0., 0.);
- AliMatrix(idrotm[1003], 90., 20., 90., 110., 0., 0.);
- AliMatrix(idrotm[1004], 90., 340., 108.2, 70., 18.2, 70.);
- AliMatrix(idrotm[1005], 90., 0., 109.31, 90., 19.31, 90.);
- AliMatrix(idrotm[1006], 90., 20., 108.2, 110., 18.2, 110.);
-
- gMC->Gspos("RICH", 1, "ALIC", 0., 471.9, 165.26, idrotm[1000], "ONLY");
- gMC->Gspos("RICH", 2, "ALIC", 171., 470., 0., idrotm[1001], "ONLY");
- gMC->Gspos("RICH", 3, "ALIC", 0., 500., 0., idrotm[1002], "ONLY");
- gMC->Gspos("RICH", 4, "ALIC", -171., 470., 0., idrotm[1003], "ONLY");
- gMC->Gspos("RICH", 5, "ALIC", 161.4, 443.4, -165.3, idrotm[1004], "ONLY");
- gMC->Gspos("RICH", 6, "ALIC", 0., 471.9, -165.3, idrotm[1005], "ONLY");
- gMC->Gspos("RICH", 7, "ALIC", -161.4, 443.4, -165.3, idrotm[1006], "ONLY");
-
+ fX = 0.;
+ fY = 0.;
+ fDigits = new TArrayI(10);
+ fNdigit=0;
+ AddDigit(FirstDigit);
+ fChamber=Ichamber;
+ fCathod=Icathod;
}
-//_____________________________________________________________________________
-void AliRICHv1::DrawModule()
+void AliRICHRecCluster::AddDigit(Int_t Digit)
{
- //
- // Draw a shaded view of the Ring Imaging Cherenkov
- //
-
- TGeant3 *geant3 = (TGeant3*) gMC;
-
- // Set everything unseen
- gMC->Gsatt("*", "seen", -1);
- //
- // Set ALIC mother transparent
- gMC->Gsatt("ALIC","SEEN",0);
- //
- // Set the volumes visible
-
- gMC->Gsatt("RICH","seen",0);
- gMC->Gsatt("SRIC","seen",0);
- gMC->Gsatt("HONE","seen",1);
- gMC->Gsatt("ALUM","seen",1);
- gMC->Gsatt("QUAR","seen",1);
- gMC->Gsatt("SPAC","seen",1);
- gMC->Gsatt("OQUA","seen",1);
- gMC->Gsatt("OQUF","seen",1);
- gMC->Gsatt("BARR","seen",1);
- gMC->Gsatt("FREO","seen",1);
- gMC->Gsatt("META","seen",1);
- gMC->Gsatt("GAP ","seen",1);
- gMC->Gsatt("CSI ","seen",1);
- gMC->Gsatt("GRID","seen",1);
-
- //
- geant3->Gdopt("hide", "on");
- geant3->Gdopt("shad", "on");
- geant3->Gsatt("*", "fill", 7);
- geant3->SetClipBox(".");
- geant3->Gdopt("hide", "on");
- geant3->Gdopt("shad", "on");
- geant3->Gsatt("*", "fill", 7);
- geant3->SetClipBox(".");
- // geant3->SetClipBox("*", 0, 2000, -2000, 2000, -2000, 2000);
- geant3->DefaultRange();
- gMC->Gdraw("alic", 60, 50, 0, 10, 0, .03, .03);
- geant3->Gdhead(1111, "Ring Imaging Cherenkov version 1");
- geant3->Gdman(16, 6, "MAN");
- geant3->Gdopt("hide", "off");
+ if (fNdigit==fDigits->GetSize()) {
+ //enlarge the list by hand!
+ Int_t *array= new Int_t[fNdigit*2];
+ for(Int_t i=0;i<fNdigit;i++)
+ array[i] = fDigits->At(i);
+ fDigits->Adopt(fNdigit*2,array);
+ }
+ fDigits->AddAt(Digit,fNdigit);
+ fNdigit++;
}
-//_____________________________________________________________________________
-void AliRICHv1::CreateMaterials()
+
+AliRICHRecCluster::~AliRICHRecCluster()
{
- //
- // *** DEFINITION OF AVAILABLE RICH MATERIALS ***
- // ORIGIN : NICK VAN EIJNDHOVEN
- // Modified by: N. Colonna (INFN - BARI, Nicola.Colonna@ba.infn.it)
- // R.A. Fini (INFN - BARI, Rosanna.Fini@ba.infn.it)
- // R.A. Loconsole (Bari University, loco@riscom.ba.infn.it)
- //
- Int_t ISXFLD = gAlice->Field()->Integ();
- Float_t SXMGMX = gAlice->Field()->Max();
-
- Float_t ppckov[14] = { 5.63e-9,5.77e-9,5.9e-9,6.05e-9,6.2e-9,6.36e-9,6.52e-9,
- 6.7e-9,6.88e-9,7.08e-9,7.3e-9,7.51e-9,7.74e-9,8e-9 };
- Float_t rindex_quarz[14] = { 1.528309,1.533333,
- 1.538243,1.544223,1.550568,1.55777,
- 1.565463,1.574765,1.584831,1.597027,
- 1.611858,1.6277,1.6472,1.6724 };
- Float_t rindex_quarzo[14] = { 1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1. };
- Float_t rindex_methane[14] = { 1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1. };
- Float_t rindex_gri[14] = { 1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1. };
- Float_t absco_freon[14] = { 179.0987,179.0987,
- 179.0987,179.0987,179.0987,35.7,12.54,5.92,4.92,3.86,1.42,.336,.134,0. };
- Float_t absco_quarz[14] = { 20.126,16.27,13.49,11.728,9.224,8.38,7.44,7.17,
- 6.324,4.483,1.6,.323,.073,0. };
- Float_t absco_quarzo[14] = { 1e-5,1e-5,1e-5,1e-5,1e-5,1e-5,1e-5,1e-5,1e-5,
- 1e-5,1e-5,1e-5,1e-5,1e-5 };
- Float_t absco_csi[14] = { 1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,
- 1e-4,1e-4,1e-4,1e-4 };
- Float_t absco_methane[14] = { 1e6,1e6,1e6,1e6,1e6,1e6,1e6,1e6,1e6,1e6,1e6,
- 1e6,1e6,1e6 };
- Float_t absco_gri[14] = { 1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,
- 1e-4,1e-4,1e-4,1e-4 };
- Float_t effic_all[14] = { 1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1. };
- Float_t effic_csi[14] = { 4.74e-4,.00438,.009,.0182,.0282,.0653,.1141,.163,
- .2101,.2554,.293,.376,.3861,.418 };
- Float_t effic_gri[14] = { 1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1. };
-
- Float_t afre[2], agri, amet[2], aqua[2], ahon, zfre[2], zgri, zhon,
- zmet[2], zqua[2];
- Int_t nlmatfre;
- Float_t densquao;
- Int_t nlmatmet, nlmatqua;
- Float_t wmatquao[2], rindex_freon[14];
- Int_t i;
- Float_t aquao[2], epsil, stmin, zquao[2];
- Int_t nlmatquao;
- Float_t radlal, densal, tmaxfd, deemax, stemax;
- Float_t aal, zal, radlgri, densfre, radlhon, densgri, denshon,densqua, densmet, wmatfre[2], wmatmet[2], wmatqua[2];
-
- Int_t *idtmed = fIdtmed->GetArray()-999;
-
- TGeant3 *geant3 = (TGeant3*) gMC;
-
- // --- Photon energy (GeV)
- // --- Refraction indexes
- for (i = 0; i < 14; ++i) {
- rindex_freon[i] = ppckov[i] * .01095 * 1e9 + 1.2177;
- }
- // need to be changed
-
- // --- Absorbtion lenghts (in cm)
- // DATA ABSCO_QUARZ /
- // & 5 * 1000000.,
- // & 29.85, 7.34, 4.134, 1.273, 0.722,
- // & 0.365, 0.365, 0.365, 0. /
- // need to be changed
-
- // --- Detection efficiencies (quantum efficiency for CsI)
- // --- Define parameters for honeycomb.
- // Used carbon of equivalent rad. lenght
-
- ahon = 12.01;
- zhon = 6.;
- denshon = 2.265;
- radlhon = 18.8;
-
- // --- Parameters to include in GSMIXT, relative to Quarz (SiO2)
-
- aqua[0] = 28.09;
- aqua[1] = 16.;
- zqua[0] = 14.;
- zqua[1] = 8.;
- densqua = 2.64;
- nlmatqua = -2;
- wmatqua[0] = 1.;
- wmatqua[1] = 2.;
-
- // --- Parameters to include in GSMIXT, relative to opaque Quarz (SiO2)
-
- aquao[0] = 28.09;
- aquao[1] = 16.;
- zquao[0] = 14.;
- zquao[1] = 8.;
- densquao = 2.64;
- nlmatquao = -2;
- wmatquao[0] = 1.;
- wmatquao[1] = 2.;
-
- // --- Parameters to include in GSMIXT, relative to Freon (C6F14)
-
- afre[0] = 12.;
- afre[1] = 19.;
- zfre[0] = 6.;
- zfre[1] = 9.;
- densfre = 1.7;
- nlmatfre = -2;
- wmatfre[0] = 6.;
- wmatfre[1] = 14.;
-
- // --- Parameters to include in GSMIXT, relative to methane (CH4)
-
- amet[0] = 12.01;
- amet[1] = 1.;
- zmet[0] = 6.;
- zmet[1] = 1.;
- densmet = 7.17e-4;
- nlmatmet = -2;
- wmatmet[0] = 1.;
- wmatmet[1] = 4.;
-
- // --- Parameters to include in GSMIXT, relative to anode grid (Cu)
-
- agri = 63.54;
- zgri = 29.;
- densgri = 8.96;
- radlgri = 1.43;
-
- // --- Parameters to include in GSMATE related to aluminium sheet
-
- aal = 26.98;
- zal = 13.;
- densal = 2.7;
- radlal = 8.9;
-
- AliMaterial(1, "Air $", 14.61, 7.3, .001205, 30420., 67500);
- AliMaterial(6, "HON", ahon, zhon, denshon, radlhon, 0);
- AliMaterial(16, "CSI", ahon, zhon, denshon, radlhon, 0);
- AliMixture(20, "QUA", aqua, zqua, densqua, nlmatqua, wmatqua);
- AliMixture(21, "QUAO", aquao, zquao, densquao, nlmatquao, wmatquao);
- AliMixture(30, "FRE", afre, zfre, densfre, nlmatfre, wmatfre);
- AliMixture(40, "MET", amet, zmet, densmet, nlmatmet, wmatmet);
- AliMixture(41, "METG", amet, zmet, densmet, nlmatmet, wmatmet);
- AliMaterial(11, "GRI", agri, zgri, densgri, radlgri, 0);
- AliMaterial(50, "ALUM", aal, zal, densal, radlal, 0);
-
- tmaxfd = -10.;
- stemax = -.1;
- deemax = -.2;
- epsil = .001;
- stmin = -.001;
-
- AliMedium(1, "DEFAULT MEDIUM AIR$", 1, 0, ISXFLD, SXMGMX, tmaxfd, stemax, deemax, epsil, stmin);
- AliMedium(2, "HONEYCOMB$", 6, 0, ISXFLD, SXMGMX, tmaxfd, stemax, deemax, epsil, stmin);
- AliMedium(3, "QUARZO$", 20, 1, ISXFLD, SXMGMX, tmaxfd, stemax, deemax, epsil, stmin);
- AliMedium(4, "FREON$", 30, 1, ISXFLD, SXMGMX, tmaxfd, stemax, deemax, epsil, stmin);
- AliMedium(5, "METANO$", 40, 1, ISXFLD, SXMGMX, tmaxfd, stemax, deemax, epsil, stmin);
- AliMedium(6, "CSI$", 16, 1, ISXFLD, SXMGMX,tmaxfd, stemax, deemax, epsil, stmin);
- AliMedium(7, "GRIGLIA$", 11, 0, ISXFLD, SXMGMX, tmaxfd, stemax, deemax, epsil, stmin);
- AliMedium(8, "QUARZOO$", 21, 1, ISXFLD, SXMGMX, tmaxfd, stemax, deemax, epsil, stmin);
- AliMedium(9, "GAP$", 41, 1, ISXFLD, SXMGMX,tmaxfd, .1, -deemax, epsil, -stmin);
- AliMedium(10, "ALUMINUM$", 50, 1, ISXFLD, SXMGMX, tmaxfd, stemax, deemax, epsil, stmin);
-
-
- // Switch on delta-ray production in the methane and freon gaps
-
- gMC->Gstpar(idtmed[1002], "LOSS", 1.);
- gMC->Gstpar(idtmed[1003], "LOSS", 1.);
- gMC->Gstpar(idtmed[1004], "LOSS", 1.);
- gMC->Gstpar(idtmed[1008], "LOSS", 1.);
- gMC->Gstpar(idtmed[1005], "LOSS", 1.);
- gMC->Gstpar(idtmed[1002], "HADR", 1.);
- gMC->Gstpar(idtmed[1003], "HADR", 1.);
- gMC->Gstpar(idtmed[1004], "HADR", 1.);
- gMC->Gstpar(idtmed[1008], "HADR", 1.);
- gMC->Gstpar(idtmed[1005], "HADR", 1.);
- gMC->Gstpar(idtmed[1002], "DCAY", 1.);
- gMC->Gstpar(idtmed[1003], "DCAY", 1.);
- gMC->Gstpar(idtmed[1004], "DCAY", 1.);
- gMC->Gstpar(idtmed[1008], "DCAY", 1.);
- gMC->Gstpar(idtmed[1005], "DCAY", 1.);
- geant3->Gsckov(idtmed[1000], 14, ppckov, absco_methane, effic_all, rindex_methane);
- geant3->Gsckov(idtmed[1001], 14, ppckov, absco_methane, effic_all, rindex_methane);
- geant3->Gsckov(idtmed[1002], 14, ppckov, absco_quarz, effic_all,rindex_quarz);
- geant3->Gsckov(idtmed[1003], 14, ppckov, absco_freon, effic_all,rindex_freon);
- geant3->Gsckov(idtmed[1004], 14, ppckov, absco_methane, effic_all, rindex_methane);
- geant3->Gsckov(idtmed[1005], 14, ppckov, absco_csi, effic_csi, rindex_methane);
- geant3->Gsckov(idtmed[1006], 14, ppckov, absco_gri, effic_gri, rindex_gri);
- geant3->Gsckov(idtmed[1007], 14, ppckov, absco_quarzo, effic_all, rindex_quarzo);
- geant3->Gsckov(idtmed[1008], 14, ppckov, absco_methane, effic_all, rindex_methane);
- geant3->Gsckov(idtmed[1009], 14, ppckov, absco_gri, effic_gri, rindex_gri);
+ if (fDigits)
+ delete fDigits;
}
-ClassImp(AliRICHhit)
-
-//_____________________________________________________________________________
-AliRICHhit::AliRICHhit(Int_t shunt, Int_t track, Int_t *vol,
- Float_t *hits, Int_t fNpadhits) :
- AliHit(shunt,track)
+Int_t AliRICHRecCluster::FirstDigitIndex()
{
- //
- // Standard constructor for RICH hit
- //
- Int_t i;
- for (i=0;i<2;i++) fVolume[i] = vol[i];
- fTrack = (int) track;
- //Hit information
- fPart = (int) hits[ 1];
- //AliHit information, position of the hit
- fX = hits[ 2];
- fY = hits[ 3];
- //Pad information
- fFirstpad = fNpadhits;
- fLastpad = fNpadhits-1;
-
- //Hit information
- fModule = (int) hits[ 4];
- fTheta = hits[ 5];
- fArrivaltime= hits[ 6];
- fFeed = (int) hits[ 7];
+ fCurrentDigit=0;
+ return fDigits->At(fCurrentDigit);
}
-ClassImp(AliRICHmip)
+Int_t AliRICHRecCluster::NextDigitIndex()
+{
+ fCurrentDigit++;
+ if (fCurrentDigit<fNdigit)
+ return fDigits->At(fCurrentDigit);
+ else
+ return InvalidDigitIndex();
+}
-//_____________________________________________________________________________
-AliRICHmip::AliRICHmip(Int_t shunt, Int_t track, Int_t *vol,
- Float_t *hits, Int_t fNckovs, Int_t fNpadhits) :
- AliRICHhit(shunt,track,vol,hits,fNpadhits)
+Int_t AliRICHRecCluster::NDigits()
{
- //
- // Standard constructor for RICH Mip hit
- //
- fPhi = hits[ 8];
- fPs = hits[ 9];
- fQ = hits[10];
- fZ = hits[11];
- //Ckov information
- if ((int) hits[12]){
- fFirstCkov = fNckovs;
- fLastCkov = fNckovs-1;
- } else {
- fFirstCkov = -1;
- fLastCkov = -1;
- }
+ return fNdigit;
}
-ClassImp(AliRICHckov)
-
-//_____________________________________________________________________________
-AliRICHckov::AliRICHckov(Int_t shunt, Int_t track, Int_t *vol,
- Float_t *hits, Int_t fNmips, Int_t fNpadhits) :
- AliRICHhit(shunt,track,vol,hits,fNpadhits)
+void AliRICHRecCluster::Finish()
{
- //
- // Standard creator for RICH Cherenkov information
- //
- fEnergy = hits[8];
- fStop = (int) hits[9];
-
- //Parent info
- fParent = fNmips;
+ // In order to reconstruct coordinates, one has to
+ // get back to the digits which is not trivial here,
+ // because we don't know where digits are stored!
+ // Center Of Gravity, or other method should be
+ // a property of AliRICH class!
}
-ClassImp(AliRICHpadhit)
-//_____________________________________________________________________________
-AliRICHpadhit::AliRICHpadhit(Int_t shunt, Int_t track, Int_t *vol,
- Float_t *hits, Int_t fNmips, Int_t fNckovs):
- AliHit(shunt,track)
+void AliRICH::Digitise(Int_t nev,Option_t *option,Text_t *filename)
{
- //
- // Standard constructor for RICH pad hit
- //
- Int_t i;
- for (i=0;i<2;i++) fVolume[i] = vol[i];
-
- // Hit information
- fX = (int) hits[ 2];
- fY = (int) hits[ 3];
- fModule = (int) hits[ 4];
- fParentMip = fNmips;
- fParentCkov = fNckovs;
- fProd = (int) hits[ 5];
- fCharge = hits[ 6];
- fZ = -999.0; // Not implemented
+ // keep galice.root for signal and name differently the file for
+ // background when add! otherwise the track info for signal will be lost !
+
+ static Bool_t first=true;
+ static TTree *TH1;
+ static TFile *File;
+ char *Add = strstr(option,"Add");
+
+ AliRICHchamber* iChamber;
+ AliRICHsegmentation* segmentation;
+
+
+ Int_t trk[50];
+ Int_t chtrk[50];
+ TObjArray *list=new TObjArray;
+ Int_t digits[3];
+
+ AliRICH *RICH = (AliRICH *) gAlice->GetDetector("RICH");
+ AliRICHHitMap* HitMap[10];
+ for (Int_t i=0; i<10; i++) {HitMap[i]=0;}
+ if (Add ) {
+ if(first) {
+ fFileName=filename;
+ cout<<"filename"<<fFileName<<endl;
+ File=new TFile(fFileName);
+ cout<<"I have opened "<<fFileName<<" file "<<endl;
+ fHits2 = new TClonesArray("AliRICHhit",1000 );
+ fClusters2 = new TClonesArray("AliRICHcluster",10000);
+ first=false;
+ }
+ File->cd();
+ File->ls();
+ // Get Hits Tree header from file
+ if(fHits2) fHits2->Clear();
+ if(fClusters2) fClusters2->Clear();
+ if(TH1) delete TH1;
+ TH1=0;
+ //
+ char treeName[20];
+ sprintf(treeName,"TreeH%d",nev);
+ TH1 = (TTree*)gDirectory->Get(treeName);
+ if (!TH1) {
+ printf("ERROR: cannot find Hits Tree for event:%d\n",nev);
+ }
+ // Set branch addresses
+ TBranch *branch;
+ char branchname[20];
+ sprintf(branchname,"%s",GetName());
+ if (TH1 && fHits2) {
+ branch = TH1->GetBranch(branchname);
+ if (branch) branch->SetAddress(&fHits2);
+ }
+ if (TH1 && fClusters2) {
+ branch = TH1->GetBranch("RICHCluster");
+ if (branch) branch->SetAddress(&fClusters2);
+ }
+ }
+ //
+ // loop over cathodes
+ //
+ AliRICHHitMap* hm;
+
+ for (int icat=0; icat<1; icat++) {
+ for (Int_t i=0; i<7; i++) {
+ if (HitMap[i]) {
+ hm=HitMap[i];
+ delete hm;
+ HitMap[i]=0;
+ }
+ }
+ Int_t counter=0;
+ for (Int_t i =0; i<7; i++) {
+ iChamber=(AliRICHchamber*) (*fChambers)[i];
+ if (iChamber->Nsec()==1 && icat==1) {
+ continue;
+ } else {
+ segmentation=iChamber->GetSegmentationModel(icat+1);
+ }
+ HitMap[i] = new AliRICHHitMapA1(segmentation, list);
+ }
+ printf("Start loop over tracks \n");
+//
+// Loop over tracks
+//
+
+ TTree *TH = gAlice->TreeH();
+ Int_t ntracks =(Int_t) TH->GetEntries();
+ for (Int_t track=0; track<ntracks; track++) {
+ gAlice->ResetHits();
+ TH->GetEvent(track);
+//
+// Loop over hits
+ for(AliRICHhit* mHit=(AliRICHhit*)RICH->FirstHit(-1);
+ mHit;
+ mHit=(AliRICHhit*)RICH->NextHit())
+ {
+ Int_t nch = mHit->fChamber-1; // chamber number
+ if (nch >7) continue;
+ iChamber = &(RICH->Chamber(nch));
+
+//
+// Loop over pad hits
+ for (AliRICHcluster* mPad=
+ (AliRICHcluster*)RICH->FirstPad(mHit,fClusters);
+ mPad;
+ mPad=(AliRICHcluster*)RICH->NextPad(fClusters))
+ {
+ Int_t cathode = mPad->fCathode; // cathode number
+ Int_t ipx = mPad->fPadX; // pad number on X
+ Int_t ipy = mPad->fPadY; // pad number on Y
+ Int_t iqpad = mPad->fQpad; // charge per pad
+//
+//
+
+ if (cathode != (icat+1)) continue;
+ // fill the info array
+ Float_t thex, they;
+ segmentation=iChamber->GetSegmentationModel(cathode);
+ segmentation->GetPadCxy(ipx,ipy,thex,they);
+ TVector *trinfo_p= new TVector(2);
+ TVector &trinfo = *trinfo_p;
+ trinfo(0)=(Float_t)track;
+ trinfo(1)=(Float_t)iqpad;
+
+ digits[0]=ipx;
+ digits[1]=ipy;
+ digits[2]=iqpad;
+
+ AliRICHlist* pdigit;
+ // build the list of fired pads and update the info
+ if (!HitMap[nch]->TestHit(ipx, ipy)) {
+ list->AddAtAndExpand(
+ new AliRICHlist(nch,digits),counter);
+ HitMap[nch]->SetHit(ipx, ipy, counter);
+ counter++;
+ pdigit=(AliRICHlist*)list->At(list->GetLast());
+ // list of tracks
+ TObjArray *trlist=(TObjArray*)pdigit->TrackList();
+ trlist->Add(&trinfo);
+ } else {
+ pdigit=(AliRICHlist*) HitMap[nch]->GetHit(ipx, ipy);
+ // update charge
+ (*pdigit).fSignal+=iqpad;
+ // update list of tracks
+ TObjArray* trlist=(TObjArray*)pdigit->TrackList();
+ Int_t last_entry=trlist->GetLast();
+ TVector *ptrk_p=(TVector*)trlist->At(last_entry);
+ TVector &ptrk=*ptrk_p;
+ Int_t last_track=Int_t(ptrk(0));
+ Int_t last_charge=Int_t(ptrk(1));
+ if (last_track==track) {
+ last_charge+=iqpad;
+ trlist->RemoveAt(last_entry);
+ trinfo(0)=last_track;
+ trinfo(1)=last_charge;
+ trlist->AddAt(&trinfo,last_entry);
+ } else {
+ trlist->Add(&trinfo);
+ }
+ // check the track list
+ Int_t nptracks=trlist->GetEntriesFast();
+ if (nptracks > 2) {
+ printf("Attention - nptracks > 2 %d \n",nptracks);
+ printf("cat,nch,ix,iy %d %d %d %d \n",icat+1,nch,ipx,ipy);
+ for (Int_t tr=0;tr<nptracks;tr++) {
+ TVector *pptrk_p=(TVector*)trlist->At(tr);
+ TVector &pptrk=*pptrk_p;
+ trk[tr]=Int_t(pptrk(0));
+ chtrk[tr]=Int_t(pptrk(1));
+ }
+ } // end if nptracks
+ } // end if pdigit
+ } //end loop over clusters
+ } // hit loop
+ } // track loop
+
+ Int_t nentr1=list->GetEntriesFast();
+ printf(" \n counter, nentr1 %d %d\n",counter,nentr1);
+
+ // open the file with background
+
+ if (Add ) {
+ ntracks =(Int_t)TH1->GetEntries();
+ printf("background - icat,ntracks1 %d %d\n",icat,ntracks);
+ printf("background - Start loop over tracks \n");
+//
+// Loop over tracks
+//
+ for (Int_t track=0; track<ntracks; track++) {
+
+ if (fHits2) fHits2->Clear();
+ if (fClusters2) fClusters2->Clear();
+
+ TH1->GetEvent(track);
+//
+// Loop over hits
+ AliRICHhit* mHit;
+ for(int i=0;i<fHits2->GetEntriesFast();++i)
+ {
+ mHit=(AliRICHhit*) (*fHits2)[i];
+ Int_t nch = mHit->fChamber-1; // chamber number
+ if (nch >9) continue;
+ iChamber = &(RICH->Chamber(nch));
+ Int_t rmin = (Int_t)iChamber->RInner();
+ Int_t rmax = (Int_t)iChamber->ROuter();
+//
+// Loop over pad hits
+ for (AliRICHcluster* mPad=
+ (AliRICHcluster*)RICH->FirstPad(mHit,fClusters2);
+ mPad;
+ mPad=(AliRICHcluster*)RICH->NextPad(fClusters2))
+ {
+ Int_t cathode = mPad->fCathode; // cathode number
+ Int_t ipx = mPad->fPadX; // pad number on X
+ Int_t ipy = mPad->fPadY; // pad number on Y
+ Int_t iqpad = mPad->fQpad; // charge per pad
+ if (track==3 && nch==0 && icat==0) printf("bgr - track,iqpad,ipx,ipy %d %d %d %d\n",track,iqpad,ipx,ipy);
+//
+//
+ if (cathode != (icat+1)) continue;
+ // fill the info array
+ Float_t thex, they;
+ segmentation=iChamber->GetSegmentationModel(cathode);
+ segmentation->GetPadCxy(ipx,ipy,thex,they);
+ Float_t rpad=TMath::Sqrt(thex*thex+they*they);
+ if (rpad < rmin || iqpad ==0 || rpad > rmax) continue;
+
+ TVector *trinfo_p;
+ trinfo_p = new TVector(2);
+ TVector &trinfo = *trinfo_p;
+ trinfo(0)=-1; // tag background
+ trinfo(1)=-1;
+
+ digits[0]=ipx;
+ digits[1]=ipy;
+ digits[2]=iqpad;
+
+
+ if (track <4 && icat==0 && nch==0)
+ printf("bgr - HitMap[nch]->TestHit(ipx, ipy),track %d %d\n",
+ HitMap[nch]->TestHit(ipx, ipy),track);
+ AliRICHlist* pdigit;
+ // build the list of fired pads and update the info
+ if (!HitMap[nch]->TestHit(ipx, ipy)) {
+ list->AddAtAndExpand(new AliRICHlist(nch,digits),counter);
+
+ HitMap[nch]->SetHit(ipx, ipy, counter);
+ counter++;
+ printf("bgr new elem in list - counter %d\n",counter);
+
+ pdigit=(AliRICHlist*)list->At(list->GetLast());
+ // list of tracks
+ TObjArray *trlist=(TObjArray*)pdigit->TrackList();
+ trlist->Add(&trinfo);
+ } else {
+ pdigit=(AliRICHlist*) HitMap[nch]->GetHit(ipx, ipy);
+ // update charge
+ (*pdigit).fSignal+=iqpad;
+ // update list of tracks
+ TObjArray* trlist=(TObjArray*)pdigit->TrackList();
+ Int_t last_entry=trlist->GetLast();
+ TVector *ptrk_p=(TVector*)trlist->At(last_entry);
+ TVector &ptrk=*ptrk_p;
+ Int_t last_track=Int_t(ptrk(0));
+ if (last_track==-1) {
+ continue;
+ } else {
+ trlist->Add(&trinfo);
+ }
+ // check the track list
+ Int_t nptracks=trlist->GetEntriesFast();
+ if (nptracks > 0) {
+ for (Int_t tr=0;tr<nptracks;tr++) {
+ TVector *pptrk_p=(TVector*)trlist->At(tr);
+ TVector &pptrk=*pptrk_p;
+ trk[tr]=Int_t(pptrk(0));
+ chtrk[tr]=Int_t(pptrk(1));
+ }
+ } // end if nptracks
+ } // end if pdigit
+ } //end loop over clusters
+ } // hit loop
+ } // track loop
+ Int_t nentr2=list->GetEntriesFast();
+ printf(" \n counter2, nentr2 %d %d \n",counter,nentr2);
+ TTree *fAli=gAlice->TreeK();
+ if (fAli) File =fAli->GetCurrentFile();
+ File->cd();
+ } // if Add
+
+ Int_t tracks[10];
+ Int_t charges[10];
+ cout<<"start filling digits \n "<<endl;
+ Int_t nentries=list->GetEntriesFast();
+ printf(" \n \n nentries %d \n",nentries);
+
+ // start filling the digits
+
+ for (Int_t nent=0;nent<nentries;nent++) {
+ AliRICHlist *address=(AliRICHlist*)list->At(nent);
+ if (address==0) continue;
+ Int_t ich=address->fChamber;
+ Int_t q=address->fSignal;
+ iChamber=(AliRICHchamber*) (*fChambers)[ich];
+ // add white noise and do zero-suppression and signal truncation
+ Float_t MeanNoise = gRandom->Gaus(1, 0.2);
+ Float_t ZeroSupp=5*MeanNoise;
+ Float_t Noise = gRandom->Gaus(0, MeanNoise);
+ q+=(Int_t)Noise;
+ if ( q <= ZeroSupp) continue;
+ digits[0]=address->fPadX;
+ digits[1]=address->fPadY;
+ digits[2]=q;
+
+ TObjArray* trlist=(TObjArray*)address->TrackList();
+ Int_t nptracks=trlist->GetEntriesFast();
+
+ // this was changed to accomodate the real number of tracks
+ if (nptracks > 10) {
+ cout<<"Attention - nptracks > 10 "<<nptracks<<endl;
+ nptracks=10;
+ }
+ if (nptracks > 2) {
+ printf("Attention - nptracks > 2 %d \n",nptracks);
+ printf("cat,ich,ix,iy,q %d %d %d %d %d \n",icat,ich,digits[0],digits[1],q);
+ }
+ for (Int_t tr=0;tr<nptracks;tr++) {
+ TVector *pp_p=(TVector*)trlist->At(tr);
+ TVector &pp =*pp_p;
+ tracks[tr]=Int_t(pp(0));
+ charges[tr]=Int_t(pp(1));
+ } //end loop over list of tracks for one pad
+ if (nptracks < 10 ) {
+ for (Int_t i=nptracks; i<10; i++) {
+ tracks[i]=0;
+ charges[i]=0;
+ }
+ }
+ // fill digits
+ RICH->AddDigits(ich,tracks,charges,digits);
+
+ delete address;
+ }
+ cout<<"I'm out of the loops for digitisation"<<endl;
+ gAlice->TreeD()->Fill();
+ TTree *TD=gAlice->TreeD();
+ Stat_t ndig=TD->GetEntries();
+ cout<<"number of digits "<<ndig<<endl;
+ TClonesArray *fDch;
+ for (int i=0;i<7;i++) {
+ fDch= RICH->DigitsAddress(i);
+ int ndig=fDch->GetEntriesFast();
+ printf (" i, ndig %d %d \n",i,ndig);
+ }
+ RICH->ResetDigits();
+
+ list->Clear();
+
+ } //end loop over cathodes
+ char hname[30];
+ sprintf(hname,"TreeD%d",nev);
+ gAlice->TreeD()->Write(hname);
}
+
+
+
+
+
+
+
////////////////////////////////////////////////
// Manager and hits classes for set:RICH //
////////////////////////////////////////////////
-
#include "AliDetector.h"
#include "AliHit.h"
+#include "AliRICHConst.h"
+#include "AliDigit.h"
+#include <TVector.h>
+#include <TObjArray.h>
+#include <TFile.h>
+#include <TTree.h>
+#include <TRotMatrix.h>
-class AliRICH : public AliDetector {
-
-protected:
- Int_t fNmips; //Number of mips in RICH
- Int_t fNckovs; //Number of cerenkovs in RICH
- Int_t fNpadhits; //Number of pad hits in RICH
-
- TClonesArray *fMips; //List of mips
- TClonesArray *fCkovs; //List of cerenkovs
- TClonesArray *fPadhits; //List of Padhits
-
- Float_t fChslope; //Charge slope
- Float_t fAlphaFeed; //Feed-back coefficient
- Float_t fSxcharge; //Charge slope along x
- Int_t fIritri; //Trigger flag
+static const int NCH=7;
+typedef enum {mip, cerenkov} Response_t;
-public:
- AliRICH();
- AliRICH(const char *name, const char *title);
- virtual ~AliRICH();
- virtual void AddHit(Int_t, Int_t*, Float_t*);
- virtual void AddMipHit(Int_t, Int_t*, Float_t*);
- virtual void AddCkovHit(Int_t, Int_t*, Float_t*);
- virtual void AddPadHit(Int_t, Int_t*, Float_t*);
- virtual void BuildGeometry();
- virtual void CreateGeometry() {}
- virtual void CreateMaterials() {}
- Int_t DistancetoPrimitive(Int_t px, Int_t py);
- inline virtual int GetNmips() {return fNmips;}
- inline virtual int GetNckovs() {return fNckovs;}
- inline virtual int GetNpadhits() {return fNpadhits;}
- virtual Int_t IsVersion() const =0;
- virtual void Init();
- inline TClonesArray *Mips() {return fMips;}
- inline TClonesArray *Ckovs() {return fCkovs;}
- inline TClonesArray *Padhits() {return fPadhits;}
- void FinishEvent(void){;}
- virtual void MakeBranch(Option_t *);
- void SetTreeAddress(void);
- virtual void StepManager();
- virtual void PreTrack();
-
- virtual void SetSP(Float_t chslope){ fChslope=chslope;}
- virtual void SetFEED(Float_t alphafeed){fAlphaFeed=alphafeed;}
- virtual void SetSIGM(Float_t sxcharge){fSxcharge=sxcharge;}
- virtual void SetTRIG(Int_t iritri) {fIritri=iritri;}
- virtual void ResetHits();
- virtual void UpdateMipHit(Float_t*);
- virtual void RichIntegration();
- virtual void AnodicWires(Float_t &);
- virtual void GetChargeMip(Float_t &);
- virtual void GetCharge(Float_t &);
- virtual void FeedBack(Float_t *, Float_t );
- virtual Float_t FMathieson(Float_t , Float_t );
-
- ClassDef(AliRICH,1) // Base class for RICH
-};
+class AliRICHcluster;
+class AliRICHchamber;
+class AliRICHRecCluster;
+class AliRICHCerenkov;
-class AliRICHv1 : public AliRICH {
-
-public:
- AliRICHv1();
- AliRICHv1(const char *name, const char *title);
- virtual ~AliRICHv1();
- virtual void CreateGeometry();
- virtual void CreateMaterials();
- virtual Int_t IsVersion() const {return 1;}
- virtual void DrawModule();
-
-
- ClassDef(AliRICHv1,1) // RICH version 1
+//----------------------------------------------
+//----------------------------------------------
+//
+// Chamber segmentation virtual base class
+//
+class AliRICHsegmentation :
+public TObject {
+
+ public:
+
+ // Set Chamber Segmentation Parameters
+ virtual void SetPADSIZ(Float_t p1, Float_t p2) =0;
+ virtual void SetDAnod(Float_t D) =0;
+ // Transform from pad (wire) to real coordinates and vice versa
+ virtual Float_t GetAnod(Float_t xhit) =0;
+ virtual void GetPadIxy(Float_t x ,Float_t y ,Int_t &ix,Int_t &iy)=0;
+ virtual void GetPadCxy(Int_t ix,Int_t iy,Float_t &x ,Float_t &y )=0;
+ //
+ // Initialisation
+ virtual void Init(AliRICHchamber*) =0;
+ //
+ // Get member data
+ virtual Float_t Dpx() =0;
+ virtual Float_t Dpy() =0;
+ virtual Int_t Npx() =0;
+ virtual Int_t Npy() =0;
+ //
+ // Iterate over pads
+ virtual void FirstPad(Float_t xhit, Float_t yhit, Float_t dx, Float_t dy) =0;
+ virtual void NextPad()=0;
+ virtual Int_t MorePads() =0;
+ // Get next neighbours
+ virtual void Neighbours
+ (Int_t iX, Int_t iY, Int_t* Nlist, Int_t Xlist[10], Int_t Ylist[10]) =0;
+ // Provisory RecCluster coordinates reconstructor
+ virtual void FitXY(AliRICHRecCluster* Cluster,TClonesArray* RICHdigits) =0;
+ //
+ // Current pad cursor during disintegration
+ virtual Int_t Ix() =0;
+ virtual Int_t Iy() =0;
+ virtual Int_t ISector() =0;
+ //
+ // Signal Generation Condition during Stepping
+ virtual Int_t SigGenCond(Float_t x, Float_t y, Float_t z) = 0;
+ virtual void SigGenInit(Float_t x, Float_t y, Float_t z) = 0;
+ virtual void IntegrationLimits
+ (Float_t& x1, Float_t& x2, Float_t& y1, Float_t& y2) = 0;
+ //
+ // Identification
+ virtual char* YourName() =0;
+ ClassDef(AliRICHsegmentation,1)
+ };
+//----------------------------------------------
+//
+// Chamber response virtual base class
+//
+class AliRICHresponse :
+public TObject {
+ public:
+ //
+ // Configuration methods
+ virtual void SetRSIGM(Float_t p1) =0;
+ virtual void SetMUCHSP(Float_t p1) =0;
+ virtual void SetMUSIGM(Float_t p1, Float_t p2) =0;
+ virtual void SetMAXADC(Float_t p1) =0;
+ //
+ // Get member data
+ virtual Float_t Chslope() =0;
+ virtual Float_t ChwX() =0;
+ virtual Float_t ChwY() =0;
+ virtual Float_t Nsigma() =0;
+ virtual Float_t adc_satm() =0;
+ //
+ // Chamber response methods
+ // Pulse height from scored quantity (eloss)
+ virtual Float_t IntPH(Float_t eloss) =0;
+ virtual Float_t IntPH() =0;
+ virtual Int_t FeedBackPhotons(Float_t *source, Float_t qtot) =0;
+ // Charge disintegration
+ virtual Float_t IntXY(AliRICHsegmentation *) =0;
+ //
+ // Identification
+ virtual char* YourName() =0;
+ // Mathieson parameters
+ virtual void SetSqrtKx3(Float_t p1) =0;
+ virtual void SetKx2(Float_t p1) =0;
+ virtual void SetKx4(Float_t p1) =0;
+ virtual void SetSqrtKy3(Float_t p1) =0;
+ virtual void SetKy2(Float_t p1) =0;
+ virtual void SetKy4(Float_t p1) =0;
+ virtual void SetPitch(Float_t p1) =0;
+ ClassDef(AliRICHresponse,1)
+ };
+
+//----------------------------------------------
+class AliRICHchamber :
+public TObject
+{
+
+ public:
+
+//Rotation matrices for each chamber
+
+ TRotMatrix *fChamberMatrix;
+ Float_t fChamberTrans[3];
+
+ public:
+ AliRICHchamber();
+ ~AliRICHchamber(){}
+//
+// Set and get GEANT id
+ Int_t GetGid() {return fGid;}
+ void SetGid(Int_t id) {fGid=id;}
+//
+// Initialisation and z-Position
+ void Init();
+ void SetZPOS(Float_t p1) {fzPos=p1;}
+ Float_t ZPosition() {return fzPos;}
+//
+
+// Set inner radius of sensitive volume
+ void SetRInner(Float_t rmin) {frMin=rmin;}
+// Set outer radius of sensitive volum
+ void SetROuter(Float_t rmax) {frMax=rmax;}
+
+// Return inner radius of sensitive volume
+ Float_t RInner() {return frMin;}
+// Return outer radius of sensitive volum
+ Float_t ROuter() {return frMax;}
+
+//Transformation from Global to local coordinates, chamber-dependant
+ void LocaltoGlobal(Float_t pos[3],Float_t Localpos[3]);
+
+//Setting chamber specific rotation matrices
+
+ void SetChamberTransform(Float_t Trans1,Float_t Trans2,Float_t Trans3,TRotMatrix *Matrix)
+
+ {
+ fChamberMatrix=Matrix;
+ fChamberTrans[0]=Trans1;
+ fChamberTrans[1]=Trans2;
+ fChamberTrans[2]=Trans3;
+ }
+
+// Configure response model
+ void ResponseModel(Response_t res, AliRICHresponse* thisResponse);
+
+ //
+// Configure segmentation model
+ void SegmentationModel(Int_t i, AliRICHsegmentation* thisSegmentation) {
+ (*fSegmentation)[i-1] = thisSegmentation;
+ }
+//
+// Get reference to response model
+ AliRICHresponse* GetResponseModel(Response_t res);
+//
+// Get reference to segmentation model
+ AliRICHsegmentation* GetSegmentationModel(Int_t isec) {
+ return (AliRICHsegmentation *) (*fSegmentation)[isec-1];
+ }
+ Int_t Nsec() {return fnsec;}
+ void SetNsec(Int_t nsec) {fnsec=nsec;}
+//
+// Member function forwarding to the segmentation and response models
+//
+// Calculate pulse height from energy loss
+ Float_t IntPH(Float_t eloss) {return ((AliRICHresponse*) (*fResponse)[0])->IntPH(eloss);}
+ Float_t IntPH() {return ((AliRICHresponse*) (*fResponse)[1])->IntPH(); }
+//
+// Ask segmentation if signal should be generated
+ Int_t SigGenCond(Float_t x, Float_t y, Float_t z)
+ {
+ if (fnsec==1) {
+ return ((AliRICHsegmentation*) (*fSegmentation)[0])
+ ->SigGenCond(x, y, z) ;
+ } else {
+ return (((AliRICHsegmentation*) (*fSegmentation)[0])
+ ->SigGenCond(x, y, z)) ||
+ (((AliRICHsegmentation*) (*fSegmentation)[1])
+ ->SigGenCond(x, y, z)) ;
+ }
+ }
+//
+// Initialisation of segmentation for hit
+ void SigGenInit(Float_t x, Float_t y, Float_t z)
+ {
+
+ if (fnsec==1) {
+ ((AliRICHsegmentation*) (*fSegmentation)[0])->SigGenInit(x, y, z) ;
+ } else {
+ ((AliRICHsegmentation*) (*fSegmentation)[0])->SigGenInit(x, y, z) ;
+ ((AliRICHsegmentation*) (*fSegmentation)[1])->SigGenInit(x, y, z) ;
+ }
+ }
+
+// Configuration forwarding
+//
+ void SetRSIGM(Float_t p)
+ {
+ ((AliRICHresponse*) (*fResponse)[0])->SetRSIGM(p);
+ ((AliRICHresponse*) (*fResponse)[1])->SetRSIGM(p);
+ }
+ void SetMUCHSP(Float_t p)
+ {
+ ((AliRICHresponse*) (*fResponse)[0])->SetMUCHSP(p);
+ ((AliRICHresponse*) (*fResponse)[1])->SetMUCHSP(p);
+ }
+ void SetMUSIGM(Float_t p1, Float_t p2)
+ {
+ ((AliRICHresponse*) (*fResponse)[0])->SetMUSIGM(p1,p2);
+ ((AliRICHresponse*) (*fResponse)[1])->SetMUSIGM(p1,p2);
+ }
+ void SetMAXADC(Float_t p)
+ {
+ ((AliRICHresponse*) (*fResponse)[0])->SetMAXADC(p);
+ ((AliRICHresponse*) (*fResponse)[1])->SetMAXADC(p);
+ }
+ void SetSqrtKx3(Float_t p)
+ {
+ ((AliRICHresponse*) (*fResponse)[0])->SetSqrtKx3(p);
+ ((AliRICHresponse*) (*fResponse)[1])->SetSqrtKx3(p);
+ }
+ void SetKx2(Float_t p)
+ {
+ ((AliRICHresponse*) (*fResponse)[0])->SetKx2(p);
+ ((AliRICHresponse*) (*fResponse)[1])->SetKx2(p);
+ }
+ void SetKx4(Float_t p)
+ {
+ ((AliRICHresponse*) (*fResponse)[0])->SetKx4(p);
+ ((AliRICHresponse*) (*fResponse)[1])->SetKx4(p);
+ }
+ void SetSqrtKy3(Float_t p)
+ {
+ ((AliRICHresponse*) (*fResponse)[0])->SetSqrtKy3(p);
+ ((AliRICHresponse*) (*fResponse)[1])->SetSqrtKy3(p);
+ }
+ void SetKy2(Float_t p)
+ {
+ ((AliRICHresponse*) (*fResponse)[0])->SetKy2(p);
+ ((AliRICHresponse*) (*fResponse)[1])->SetKy2(p);
+ }
+ void SetKy4(Float_t p)
+ {
+ ((AliRICHresponse*) (*fResponse)[0])->SetKy4(p);
+ ((AliRICHresponse*) (*fResponse)[1])->SetKy4(p);
+ }
+
+ void SetPitch(Float_t p)
+ {
+ ((AliRICHresponse*) (*fResponse)[0])->SetPitch(p);
+ ((AliRICHresponse*) (*fResponse)[1])->SetPitch(p);
+ }
+
+ void SetPADSIZ(Int_t isec, Float_t p1, Float_t p2) {
+ ((AliRICHsegmentation*) (*fSegmentation)[isec-1])->SetPADSIZ(p1,p2);
+ }
+//
+// Cluster formation method
+ void DisIntegration(Float_t, Float_t, Float_t, Int_t&x, Float_t newclust[6][500], Response_t res);
+ ClassDef(AliRICHchamber,1)
+
+ private:
+
+// Maximum and Minimum Chamber size
+ Float_t frMin;
+ Float_t frMax;
+// GEANT volume if for sensitive volume of this chamber
+ Int_t fGid;
+// z-position of this chamber
+ Float_t fzPos;
+// The segmentation models for the cathode planes
+// fnsec=1: one plane segmented, fnsec=2: both planes are segmented.
+ Int_t fnsec;
+ TObjArray *fSegmentation;
+ TObjArray *fResponse;
+
};
-//_____________________________________________________________________________
-class AliRICHhit: public AliHit {
-public:
- Int_t fVolume[2]; //array of volumes
-
- //Pad informations
- Int_t fFirstpad; //First index in padhits
- Int_t fLastpad; //Last index in padhits
-
- //Hit information
- Int_t fModule; //Module number
- Float_t fTheta; //Theta of the particle generating the hit
-
- Float_t fArrivaltime;// Time of hit.
- Int_t fPart; //Particle type
-
- // we don't know what is this for :
- Int_t fFeed; //Type of feedback (origin of charge deposition)
-
-public:
- AliRICHhit() {}
- AliRICHhit(Int_t shunt, Int_t track, Int_t *vol, Float_t *hits,
- Int_t fNpadhits);
- virtual ~AliRICHhit(){}
-
- void SetLastpad(Int_t lastpad){fLastpad=lastpad;}
-
- ClassDef(AliRICHhit,1) // Hits for set RICH
-};
-//_____________________________________________________________________________
-class AliRICHmip: public AliRICHhit
-{
-public:
- // Hit information keep
- Float_t fPhi; //Phi of the particle generating the hit
- Float_t fPs; //Momentum of the particle generating the hit
- Float_t fQ; //Charge of the particle
-
- // Generated cerenkov information (Z of generation stored in fZ of AliHit)
- Int_t fFirstCkov; //Index in the ckov TcloneArray of the first generated
- Int_t fLastCkov; //Here the last.
-
+
+class AliRICHcluster : public TObject {
public:
- AliRICHmip() {}
- AliRICHmip(Int_t shunt, Int_t track, Int_t *vol, Float_t *hits,
- Int_t fNckovs, Int_t fNpadhits);
- virtual ~AliRICHmip() {}
-
- Float_t GetZ() { return fZ;}
- void SetX(Float_t x) { fX = x; }
- void SetY(Float_t y) { fY = y; }
- void SetZ(Float_t z) { fZ = z; }
- void SetLastCkov(Int_t last){ fLastCkov = last; }
- void SetModule(Int_t module){ fModule = module;}
- void SetTheta(Float_t theta){ fTheta = theta; }
- void SetPhi(Float_t phi) { fPhi = phi; }
-
- ClassDef(AliRICHmip,1) //Mip hits for RICH
-};
-
-//_____________________________________________________________________________
-class AliRICHckov: public AliRICHhit
-{
+
+ Int_t fHitNumber; // Hit number
+ Int_t fCathode; // Cathode number
+ Int_t fQ ; // Total charge
+ Int_t fPadX ; // Pad number along X
+ Int_t fPadY ; // Pad number along Y
+ Int_t fQpad ; // Charge per pad
+ Int_t fRSec ; // R -sector of pad
+
public:
- // Hit information keep
- Float_t fEnergy; //Photon energy
- Int_t fStop; //Stop mechanism (cut, threshold, ...)
-
- //Parent info
- Int_t fParent; //Index in array of mips of parent which generatethis
+ AliRICHcluster() {
+ fHitNumber=fQ=fPadX=fPadY=fQpad=fRSec=0;
+ }
+ AliRICHcluster(Int_t *clhits);
+ virtual ~AliRICHcluster() {;}
+
+ ClassDef(AliRICHcluster,1) //Cluster object for set:RICH
+ };
+
+
+ class AliRICHreccluster : public TObject {
public:
- AliRICHckov() {}
- AliRICHckov(Int_t shunt, Int_t track, Int_t *vol, Float_t *hits,
- Int_t fNmips, Int_t fNpadhits);
- virtual ~AliRICHckov() {}
-
- ClassDef(AliRICHckov,1) //Cerenkov hits for RICH
-};
+
+ Int_t fTracks[3]; //labels of overlapped tracks
+
+ Int_t fQ ; // Q of cluster (in ADC counts)
+ Float_t fX ; // X of cluster
+ Float_t fY ; // Y of cluster
+
+ public:
+ AliRICHreccluster() {
+ fTracks[0]=fTracks[1]=fTracks[2]=0;
+ fQ=0; fX=fY=0;
+ }
+ virtual ~AliRICHreccluster() {;}
+
+ ClassDef(AliRICHreccluster,1) //Cluster object for set:RICH
+ };
+
+//_____________________________________________________________________________
+class AliRICHdigit : public TObject {
+ public:
+ Int_t fPadX; // Pad number along x
+ Int_t fPadY ; // Pad number along y
+ Int_t fSignal; // Signal amplitude
+
+
+ Int_t fTcharges[10]; // charge per track making this digit (up to 10)
+ Int_t fTracks[10]; // tracks making this digit (up to 10)
+
+
+
+ public:
+ AliRICHdigit() {}
+ AliRICHdigit(Int_t *digits);
+ AliRICHdigit(Int_t *tracks, Int_t *charges, Int_t *digits);
+ virtual ~AliRICHdigit() {}
+
+
+ ClassDef(AliRICHdigit,1) //Digits for set:RICH
+ };
//_____________________________________________________________________________
-class AliRICHpadhit: public AliHit
-{
-public:
- Int_t fVolume[2]; //array of volumes
-
- // Hit information
- Int_t fX; //Integer x position in pad
- Int_t fY; //Integer y position in pad
- Int_t fModule; //Module number
- // Particle info
- Int_t fParentMip; //Parent particle
- Int_t fParentCkov; //Parent CKOV
- // physics info
- Int_t fProd; //Production mechanism
- Float_t fCharge; //Charge deposited
-public:
- AliRICHpadhit(){}
- AliRICHpadhit(Int_t shunt, Int_t track, Int_t *vol, Float_t *hits,
- Int_t fNmips,Int_t fNckovs);
- virtual ~AliRICHpadhit() {}
-
- ClassDef(AliRICHpadhit,1) //Pad hits for RICH
-};
+class AliRICHlist : public AliRICHdigit {
+ public:
+
+ Int_t fRpad; // r_pos of pad
+ Int_t fChamber; // chamber number of pad
+ TObjArray *fTrackList;
+
+
+ public:
+ AliRICHlist() {fTrackList=0;}
+ AliRICHlist(Int_t ich, Int_t *digits);
+ virtual ~AliRICHlist() {}
+
+ TObjArray *TrackList() {return fTrackList;}
+
+ ClassDef(AliRICHlist,1) //Digits for set:RICH
+ };
+//___________________________________________
+
+
+//___________________________________________
+
+class AliRICHhit : public AliHit {
+ public:
+ Int_t fChamber; // Chamber number
+ Float_t fParticle; // Geant3 particle type
+ Float_t fTheta ; // Incident theta angle in degrees
+ Float_t fPhi ; // Incident phi angle in degrees
+ Float_t fTlength; // Track length inside the chamber
+ Float_t fEloss; // ionisation energy loss in gas
+ Int_t fPHfirst; // first padhit
+ Int_t fPHlast; // last padhit
+ public:
+ AliRICHhit() {}
+ AliRICHhit(Int_t fIshunt, Int_t track, Int_t *vol, Float_t *hits);
+ virtual ~AliRICHhit() {}
+
+ ClassDef(AliRICHhit,1) //Hits object for set:RICH
+ };
+
+//------------------------------------------------
+// Cerenkov photon object
+//------------------------------------------------
+class AliRICHCerenkov: public AliHit {
+ public:
+ Int_t fChamber; // Chamber number
+ Float_t fTheta ; // Incident theta angle in degrees
+ Float_t fPhi ; // Incident phi angle in degrees
+ Float_t fTlength; // Track length inside the chamber
+ Int_t fPHfirst; // first padhit
+ Int_t fPHlast; // last padhit
+ public:
+ AliRICHCerenkov() {}
+ AliRICHCerenkov(Int_t fIshunt, Int_t track, Int_t *vol, Float_t *Cerenkovs);
+ virtual ~AliRICHCerenkov() {}
+
+ ClassDef(AliRICHCerenkov,1) //Cerenkovs object for set:RICH
+ };
+
+//--------------------------------------------------
+
+class AliRICH : public AliDetector {
+ public:
+ AliRICH();
+ AliRICH(const char *name, const char *title);
+ virtual ~AliRICH();
+ virtual void AddHit(Int_t, Int_t*, Float_t*);
+ virtual void AddCerenkov(Int_t, Int_t*, Float_t*);
+ virtual void AddCluster(Int_t*);
+ virtual void AddDigits(Int_t, Int_t*, Int_t*, Int_t*);
+ virtual void AddRecCluster(Int_t iCh, Int_t iCat,
+ AliRICHRecCluster* Cluster);
+ virtual void BuildGeometry();
+ virtual void CreateGeometry() {}
+ virtual void CreateMaterials() {}
+ virtual void StepManager();
+ Int_t DistancetoPrimitive(Int_t px, Int_t py);
+ virtual Int_t IsVersion() const =0;
+//
+ TClonesArray *Clusters() {return fClusters;}
+ TClonesArray *Cerenkovs() {return fCerenkovs;}
+ virtual void MakeBranch(Option_t *opt=" ");
+ void SetTreeAddress();
+ virtual void ResetHits();
+ virtual void ResetDigits();
+ virtual void ResetRecClusters();
+ virtual void ReconstructClusters();
+ virtual void Digitise(Int_t,Option_t *opt=" ",Text_t *name=" ");
+//
+// Configuration Methods (per station id)
+//
+// Set Chamber Segmentation Parameters
+// id refers to the station and isec to the cathode plane
+ virtual void SetPADSIZ(Int_t id, Int_t isec, Float_t p1, Float_t p2);
+
+// Set Signal Generation Parameters
+ virtual void SetRSIGM(Int_t id, Float_t p1);
+ virtual void SetMUCHSP(Int_t id, Float_t p1);
+ virtual void SetMUSIGM(Int_t id, Float_t p1, Float_t p2);
+ virtual void SetMAXADC(Int_t id, Float_t p1);
+// Set Segmentation and Response Model
+ virtual void SetSegmentationModel(Int_t id, Int_t isec, AliRICHsegmentation *segmentation);
+ virtual void SetResponseModel(Int_t id, Response_t res, AliRICHresponse *response);
+ virtual void SetNsec(Int_t id, Int_t nsec);
+// Set Stepping Parameters
+ virtual void SetSMAXAR(Float_t p1);
+ virtual void SetSMAXAL(Float_t p1);
+ virtual void SetDMAXAR(Float_t p1);
+ virtual void SetDMAXAL(Float_t p1);
+ virtual void SetRICHACC(Bool_t acc=0, Float_t angmin=2, Float_t angmax=9);
+// Response Simulation
+ virtual void MakePadHits(Float_t xhit,Float_t yhit,Float_t eloss,Int_t id,Response_t res);
+// Return reference to Chamber #id
+ virtual AliRICHchamber& Chamber(Int_t id) {return *((AliRICHchamber *) (*fChambers)[id]);}
+// Retrieve pad hits for a given Hit
+ virtual AliRICHcluster* FirstPad(AliRICHhit *, TClonesArray *);
+ virtual AliRICHcluster* NextPad(TClonesArray *);
+// Return pointers to digits
+ TObjArray *Dchambers() {return fDchambers;}
+ Int_t *Ndch() {return fNdch;}
+ virtual TClonesArray *DigitsAddress(Int_t id) {return ((TClonesArray *) (*fDchambers)[id]);}
+// Return pointers to reconstructed clusters
+ virtual TObjArray *RecClusters(Int_t iCh, Int_t iCat)
+ {return ( (TObjArray*) (*fRecClusters)[iCh+iCat*10]);}
+
+
+ protected:
+ TObjArray *fChambers; // List of Tracking Chambers
+ Int_t fNclusters; // Number of clusters
+ Int_t fNcerenkovs; // Number of cerenkovs
+ TClonesArray *fClusters; // List of clusters
+ TObjArray *fDchambers; // List of digits
+ TObjArray *fRecClusters; // List of clusters
+ TClonesArray *fCerenkovs; // List of cerenkovs
+ Int_t *fNdch; // Number of digits
+ Text_t *fFileName; // Filename for event mixing
+
+
+ TObjArray *fRawClusters; // List of raw clusters
+ Int_t *fNrawch; // Number of raw clusters
+ TObjArray *fCathCorrel; // List of correlated clusters
+ Int_t *fNcorch; // Number of correl clusters
+ TTree *fTreeC; // Cathode correl index tree
+
+//
+ Bool_t fAccCut; //Transport acceptance cut
+ Float_t fAccMin; //Minimum acceptance cut used during transport
+ Float_t fAccMax; //Minimum acceptance cut used during transport
+//
+
+// Stepping Parameters
+ Float_t fMaxStepGas; // Maximum step size inside the chamber gas
+ Float_t fMaxStepAlu; // Maximum step size inside the chamber aluminum
+ Float_t fMaxDestepGas; // Maximum relative energy loss in gas
+ Float_t fMaxDestepAlu; // Maximum relative energy loss in aluminum
+
+ protected:
+
+ ClassDef(AliRICH,1) //Hits manager for set:RICH
+ };
+//___________________________________________
+class AliRICHRecCluster : public TObject {
+ public:
+ AliRICHRecCluster() ;
+ AliRICHRecCluster(Int_t FirstDigit,Int_t Ichamber, Int_t Icathod) ;
+ virtual ~AliRICHRecCluster();
+ virtual void AddDigit(Int_t Digit);
+ virtual Int_t FirstDigitIndex();
+ virtual Int_t NextDigitIndex();
+ virtual Int_t InvalidDigitIndex() {return -1;}
+
+ virtual Int_t NDigits();
+ virtual void Finish(); // Nothing yet ...
+ virtual Int_t GetCathod() {return fCathod;}
+ virtual Int_t GetChamber() {return fChamber;}
+
+ public:
+ Float_t fX; // reconstructed x
+ Float_t fY; // reconstructed y
+
+ protected:
+ TArrayI *fDigits; // List of digits indexes for that cluster
+ Int_t fNdigit; // Number of digits indexes stored;
+ Int_t fCathod; // Number of the cathod to be used;
+ Int_t fChamber; // Number of the chamber to be used;
+ Int_t fCurrentDigit; // Current Digit inside an iteration
+
+ ClassDef(AliRICHRecCluster,1) //Cluster object for set:RICH
+ };
+//___________________________________________
#endif
+
+
+
+
+
+
+
+
+
+
+
+
--- /dev/null
+#ifndef RICHConst_H
+#define RICHConst_H
+/////////////////////////////////////////////////////////////////////////////
+//
+//---------------------------------------------------------------------
+// ALICE RICH chambers geometry
+//--------------------------------------------------------------------
+//
+const Float_t zend = 511.+0.15-2*0.001; // z-out position of first chamber
+
+
+/////////////////////////////////////////////////////////////////////////////
+//
+//---------------------------------------------------------------------
+// ALICE RICH Electronics Parameters
+//--------------------------------------------------------------------
+//
+//
+
+const Float_t adc_satm = 1024; // dynamic range (10 bits)
+const Float_t zero_supm = 6.; // zero suppression
+const Float_t sig_noise = 500.; // electronics noise (no. of electrons)
+
+#endif
--- /dev/null
+#include "AliRICHHitMap.h"
+
+ClassImp(AliRICHHitMap)
+ClassImp(AliRICHHitMapA1)
+
+AliRICHHitMapA1::AliRICHHitMapA1(AliRICHsegmentation *seg, TObjArray *dig)
+{
+ fSegmentation = seg;
+ fNpx = fSegmentation->Npx();
+ fNpy = fSegmentation->Npy();
+ fMaxIndex=2*(fNpx+1)*2*(fNpy+1)+2*fNpy;
+
+ fHitMap = new Int_t[fMaxIndex];
+ fDigits = dig;
+ fNdigits = fDigits->GetEntriesFast();
+ Clear();
+}
+
+
+AliRICHHitMapA1::~AliRICHHitMapA1()
+{
+// if (fDigits) delete fDigits;
+ if (fHitMap) delete[] fHitMap;
+}
+
+void AliRICHHitMapA1::Clear()
+{
+ memset(fHitMap,0,sizeof(int)*fMaxIndex);
+}
+
+Int_t AliRICHHitMapA1::CheckedIndex(Int_t ix, Int_t iy)
+{
+ Int_t index=2*fNpy*(ix+fNpx)+(iy+fNpy);
+ if (index > fMaxIndex) {
+ printf("\n \n \n Try to read/write outside array !!!! \n \n %d %d %d %d %d %d",ix,iy, fMaxIndex, index, fNpx, fNpy);
+ return fMaxIndex-1;
+ } else {
+ return index;
+ }
+}
+
+
+void AliRICHHitMapA1::FillHits()
+{
+ Int_t ndigits = fDigits->GetEntriesFast();
+ printf("\n Filling hits into HitMap\n");
+ printf("FindRawClusters -- ndigits %d \n",ndigits);
+ if (!ndigits) return;
+ AliRICHdigit *dig;
+ for (Int_t ndig=0; ndig<fNdigits; ndig++) {
+ dig = (AliRICHdigit*)fDigits->UncheckedAt(ndig);
+ SetHit(dig->fPadX,dig->fPadY,ndig);
+ }
+}
+
+
+void AliRICHHitMapA1::SetHit(Int_t ix, Int_t iy, Int_t idigit)
+{
+// fHitMap[kMaxNpady*(ix+fNpx)+(iy+fNpy)]=idigit+1;
+ fHitMap[CheckedIndex(ix, iy)]=idigit+1;
+}
+
+void AliRICHHitMapA1::DeleteHit(Int_t ix, Int_t iy)
+{
+// fHitMap[kMaxNpady*(ix+fNpx)+(iy+fNpy)]=0;
+ fHitMap[CheckedIndex(ix, iy)]=0;
+}
+
+void AliRICHHitMapA1::FlagHit(Int_t ix, Int_t iy)
+{
+ fHitMap[CheckedIndex(ix, iy)]=
+ -TMath::Abs(fHitMap[CheckedIndex(ix, iy)]);
+}
+
+Int_t AliRICHHitMapA1::GetHitIndex(Int_t ix, Int_t iy)
+{
+ return TMath::Abs(fHitMap[CheckedIndex(ix, iy)])-1;
+}
+
+TObject* AliRICHHitMapA1::GetHit(Int_t ix, Int_t iy)
+{
+ Int_t index=GetHitIndex(ix,iy);
+ // Force crash if index does not exist ! (Manu)
+ return (index <0) ? 0 : fDigits->UncheckedAt(GetHitIndex(ix,iy));
+}
+
+Flag_t AliRICHHitMapA1::TestHit(Int_t ix, Int_t iy)
+{
+ Int_t inf=fHitMap[CheckedIndex(ix, iy)];
+ if (inf < 0) {
+ return used;
+ } else if (inf == 0) {
+ return empty;
+ } else {
+ return unused;
+ }
+}
+
+
--- /dev/null
+#ifndef AliRICHHitMap_H
+#define AliRICHHitMap_H
+
+#include "AliRICH.h"
+#include "TArrayI.h"
+typedef enum {empty, used, unused} Flag_t;
+const Int_t kMaxNpadx=1200, kMaxNpady=1200;
+
+class AliRICHHitMap :
+public TObject {
+ public:
+ virtual void FillHits() =0;
+ virtual void Clear() =0;
+ virtual void SetHit(Int_t ix, Int_t iy, Int_t idigit) =0;
+ virtual void DeleteHit(Int_t ix, Int_t iy) =0;
+ virtual Int_t GetHitIndex(Int_t ix, Int_t iy) =0;
+ virtual TObject * GetHit(Int_t ix, Int_t iy) =0;
+ virtual void FlagHit(Int_t ix, Int_t iy) =0;
+ virtual Flag_t TestHit(Int_t ix, Int_t iy) =0;
+
+ ClassDef(AliRICHHitMap,1) //virtual base class for muon HitMap
+};
+
+class AliRICHHitMapA1 :
+public AliRICHHitMap
+{
+ private:
+ AliRICHsegmentation *fSegmentation;
+ Int_t fNpx;
+ Int_t fNpy;
+ TObjArray *fDigits;
+ Int_t fNdigits;
+ Int_t *fHitMap;
+ Int_t fMaxIndex;
+
+ public:
+ AliRICHHitMapA1(AliRICHsegmentation *seg, TObjArray *dig);
+ virtual ~AliRICHHitMapA1();
+ virtual void FillHits();
+ virtual void Clear();
+ virtual void SetHit(Int_t ix, Int_t iy, Int_t idigit);
+ virtual void DeleteHit(Int_t ix, Int_t iy);
+ virtual Int_t GetHitIndex(Int_t ix, Int_t iy);
+ virtual TObject* GetHit(Int_t ix, Int_t);
+ virtual void FlagHit(Int_t ix, Int_t iy);
+ virtual Flag_t TestHit(Int_t ix, Int_t iy);
+ private:
+ Int_t CheckedIndex(Int_t ix, Int_t iy);
+ ClassDef(AliRICHHitMapA1,1) // Implements HitMap as a 2-dim array
+};
+#endif
+
+
--- /dev/null
+#include "AliRICHSegResCkv.h"
+#include "TMath.h"
+#include "TRandom.h"
+
+
+//___________________________________________
+ClassImp(AliRICHresponseCkv)
+//___________________________________________
+Float_t AliRICHresponseCkv::IntPH()
+{
+
+ Float_t charge = -fChslope*TMath::Log(gRandom->Rndm());
+ return charge;
+}
+
--- /dev/null
+#ifndef RICHSegResCkv_H
+#define RICHSegResCkv_H
+
+#include "AliRICHSegResV0.h"
+
+class AliRICHresponseCkv : public AliRICHresponseV0 {
+
+ public:
+ AliRICHresponseCkv(){}
+ virtual ~AliRICHresponseCkv(){}
+
+ virtual Float_t IntPH();
+
+ ClassDef(AliRICHresponseCkv,1)
+
+
+ };
+
+#endif
--- /dev/null
+#include "AliRICHSegResV0.h"
+#include "AliRun.h"
+#include "TParticle.h"
+#include "TMath.h"
+#include "TRandom.h"
+
+
+ClassImp(AliRICHsegmentation)
+ClassImp(AliRICHresponse)
+//___________________________________________
+ClassImp(AliRICHsegmentationV0)
+
+void AliRICHsegmentationV0::Init(AliRICHchamber* Chamber)
+{
+ fNpx=(Int_t) (Chamber->ROuter()/fDpx+1);
+ fNpy=(Int_t) (Chamber->ROuter()/fDpy+1);
+}
+
+
+Float_t AliRICHsegmentationV0::GetAnod(Float_t xhit)
+{
+ Float_t wire= (xhit>0)? Int_t(xhit/fWireD)+0.5:Int_t(xhit/fWireD)-0.5;
+ return fWireD*wire;
+}
+
+void AliRICHsegmentationV0::SetPADSIZ(Float_t p1, Float_t p2)
+{
+ fDpx=p1;
+ fDpy=p2;
+}
+void AliRICHsegmentationV0::GetPadIxy(Float_t x, Float_t y, Int_t &ix, Int_t &iy)
+{
+// returns pad coordinates (ix,iy) for given real coordinates (x,y)
+//
+ ix = (x>0)? Int_t(x/fDpx)+1 : Int_t(x/fDpx);
+ iy = (y>0)? Int_t(y/fDpy)+1 : Int_t(y/fDpy);
+ if (iy > fNpy) iy= fNpy;
+ if (iy < -fNpy) iy=-fNpy;
+ if (ix > fNpx) ix= fNpx;
+ if (ix < -fNpx) ix=-fNpx;
+}
+void AliRICHsegmentationV0::
+GetPadCxy(Int_t ix, Int_t iy, Float_t &x, Float_t &y)
+{
+// returns real coordinates (x,y) for given pad coordinates (ix,iy)
+//
+ x = (ix>0) ? Float_t(ix*fDpx)-fDpx/2. : Float_t(ix*fDpx)-fDpx/2.;
+ y = (iy>0) ? Float_t(iy*fDpy)-fDpy/2. : Float_t(iy*fDpy)-fDpy/2.;
+}
+
+void AliRICHsegmentationV0::
+FirstPad(Float_t xhit, Float_t yhit, Float_t dx, Float_t dy)
+{
+ //
+ // Find the wire position (center of charge distribution)
+ Float_t x0a=GetAnod(xhit);
+ //
+ // and take fNsigma*sigma around this center
+ Float_t x01=x0a - dx;
+ Float_t x02=x0a + dx;
+ Float_t y01=yhit - dy;
+ Float_t y02=yhit + dy;
+ //
+ // find the pads over which the charge distributes
+ GetPadIxy(x01,y01,fixmin,fiymin);
+ GetPadIxy(x02,y02,fixmax,fiymax);
+ //
+ // Set current pad to lower left corner
+ fix=fixmin;
+ fiy=fiymin;
+ GetPadCxy(fix,fiy,fx,fy);
+}
+
+void AliRICHsegmentationV0::NextPad()
+{
+ //
+ // Step to next pad in integration region
+ if (fix != fixmax) {
+ fix++;
+ } else if (fiy != fiymax) {
+ fix=fixmin;
+ fiy++;
+ } else {
+ printf("\n Error: Stepping outside integration region\n ");
+ }
+ GetPadCxy(fix,fiy,fx,fy);
+}
+
+Int_t AliRICHsegmentationV0::MorePads()
+
+//
+// Are there more pads in the integration region
+{
+ if (fix == fixmax && fiy == fiymax) {
+ return 0;
+ } else {
+ return 1;
+
+ }
+}
+
+void AliRICHsegmentationV0::SigGenInit(Float_t x,Float_t y,Float_t)
+{
+//
+// Initialises pad and wire position during stepping
+ fxt =x;
+ fyt =y;
+ GetPadIxy(x,y,fixt,fiyt);
+ fiwt=Int_t(x/fWireD)+1;
+
+}
+
+Int_t AliRICHsegmentationV0::SigGenCond(Float_t x,Float_t y,Float_t)
+{
+//
+// Signal will be generated if particle crosses pad boundary or
+// boundary between two wires.
+ Int_t ixt, iyt;
+ GetPadIxy(x,y,ixt,iyt);
+ Int_t iwt=Int_t(x/fWireD)+1;
+
+ if ((ixt != fixt) || (iyt !=fiyt) || (iwt != fiwt)) {
+ return 1;
+ } else {
+ return 0;
+ }
+}
+void AliRICHsegmentationV0::
+IntegrationLimits(Float_t& x1,Float_t& x2,Float_t& y1, Float_t& y2)
+{
+ x1=fxt-fx-fDpx/2.;
+ x2=x1+fDpx;
+ y1=fyt-fy-fDpy/2.;
+ y2=y1+fDpy;
+
+}
+
+void AliRICHsegmentationV0::
+Neighbours(Int_t iX, Int_t iY, Int_t* Nlist, Int_t Xlist[7], Int_t Ylist[7])
+{
+//Is used for the cluster finder, include diagonal elements
+
+ *Nlist=4;Xlist[0]=Xlist[1]=iX;Xlist[2]=iX-1;Xlist[3]=iX+1;
+ Ylist[0]=iY-1;Ylist[1]=iY+1;Ylist[2]=Ylist[3]=iY;
+}
+
+void AliRICHsegmentationV0::
+FitXY(AliRICHRecCluster* ,TClonesArray* )
+ // Default : Centre of gravity method
+{
+ ;
+}
+
+
+//___________________________________________
+ClassImp(AliRICHresponseV0)
+ Float_t AliRICHresponseV0::IntPH(Float_t eloss)
+{
+ // Get number of electrons and return charge
+
+ Int_t nel;
+//E9/26=magic number should parameter
+ nel= Int_t(eloss*1.e9/26.);
+ Float_t charge=0;
+ if (nel == 0) nel=1;
+ for (Int_t i=1;i<=nel;i++) {
+ charge -= fChslope*TMath::Log(gRandom->Rndm());
+ }
+ return charge;
+}
+// -------------------------------------------
+Float_t AliRICHresponseV0::IntXY(AliRICHsegmentation * segmentation)
+{
+
+ const Float_t invpitch = 1/fPitch;
+ Float_t response;
+//
+// Integration limits defined by segmentation model
+//
+
+ Float_t xi1, xi2, yi1, yi2;
+ segmentation->IntegrationLimits(xi1,xi2,yi1,yi2);
+ xi1=xi1*invpitch;
+ xi2=xi2*invpitch;
+ yi1=yi1*invpitch;
+ yi2=yi2*invpitch;
+
+ //
+// The Mathieson function
+ Double_t ux1=fSqrtKx3*TMath::TanH(fKx2*xi1);
+ Double_t ux2=fSqrtKx3*TMath::TanH(fKx2*xi2);
+
+ Double_t uy1=fSqrtKy3*TMath::TanH(fKy2*yi1);
+ Double_t uy2=fSqrtKy3*TMath::TanH(fKy2*yi2);
+
+ response=4.*fKx4*(TMath::ATan(ux2)-TMath::ATan(ux1))*fKy4*(TMath::ATan(uy2)-TMath::ATan(uy1));
+
+ return response;
+
+}
+//___________________________________________
+Int_t AliRICHresponseV0::FeedBackPhotons(Float_t source[3], Float_t qtot)
+{
+ //
+ // Generate FeedBack photons
+ //
+ Int_t j, ipart, nt;
+
+ //Probability of feedback
+ Float_t fAlphaFeed=0.05;
+
+ Int_t sNfeed=0;
+
+ // Local variables
+ Float_t cthf, ranf[2], phif, enfp = 0, sthf, weight;
+ Int_t i, ifeed;
+ Float_t e1[3], e2[3], e3[3];
+ Float_t vmod, uswop;
+ Float_t fp, random;
+ Float_t dir[3], phi;
+ Int_t nfp;
+ Float_t pol[3], mom[3];
+ TLorentzVector position;
+ //
+ // Determine number of feedback photons
+
+ // Get weight of current particle
+ TParticle *current = (TParticle*)
+ (*gAlice->Particles())[gAlice->CurrentTrack()];
+
+ ifeed = Int_t(current->GetWeight()/100+0.5);
+ ipart = gMC->TrackPid();
+ fp = fAlphaFeed * qtot;
+ nfp = gRandom->Poisson(fp);
+
+ // This call to fill the time of flight
+ gMC->TrackPosition(position);
+ //
+ // Generate photons
+ for (i = 0; i <nfp; ++i) {
+
+ // Direction
+ gMC->Rndm(ranf, 2);
+ cthf = ranf[0] * 2 - 1.;
+ if (cthf < 0) continue;
+ sthf = TMath::Sqrt((1 - cthf) * (1 + cthf));
+ phif = ranf[1] * 2 * TMath::Pi();
+ //
+ gMC->Rndm(&random, 1);
+ if (random <= .57) {
+ enfp = 7.5e-9;
+ } else if (random <= .7) {
+ enfp = 6.4e-9;
+ } else {
+ enfp = 7.9e-9;
+ }
+
+ dir[0] = sthf * TMath::Sin(phif);
+ dir[1] = cthf;
+ dir[2] = sthf * TMath::Cos(phif);
+ gMC->Gdtom(dir, mom, 2);
+ mom[0]*=enfp;
+ mom[1]*=enfp;
+ mom[2]*=enfp;
+
+ // Polarisation
+ e1[0] = 0;
+ e1[1] = -dir[2];
+ e1[2] = dir[1];
+
+ e2[0] = -dir[1];
+ e2[1] = dir[0];
+ e2[2] = 0;
+
+ e3[0] = dir[1];
+ e3[1] = 0;
+ e3[2] = -dir[0];
+
+ vmod=0;
+ for(j=0;j<3;j++) vmod+=e1[j]*e1[j];
+ if (!vmod) for(j=0;j<3;j++) {
+ uswop=e1[j];
+ e1[j]=e3[j];
+ e3[j]=uswop;
+ }
+ vmod=0;
+ for(j=0;j<3;j++) vmod+=e2[j]*e2[j];
+ if (!vmod) for(j=0;j<3;j++) {
+ uswop=e2[j];
+ e2[j]=e3[j];
+ e3[j]=uswop;
+ }
+
+ vmod=0;
+ for(j=0;j<3;j++) vmod+=e1[j]*e1[j];
+ vmod=TMath::Sqrt(1/vmod);
+ for(j=0;j<3;j++) e1[j]*=vmod;
+
+ vmod=0;
+ for(j=0;j<3;j++) vmod+=e2[j]*e2[j];
+ vmod=TMath::Sqrt(1/vmod);
+ for(j=0;j<3;j++) e2[j]*=vmod;
+
+ gMC->Rndm(ranf, 1);
+ phi = ranf[0] * 2 * TMath::Pi();
+ for(j=0;j<3;j++) pol[j]=e1[j]*TMath::Sin(phi)+e2[j]*TMath::Cos(phi);
+ gMC->Gdtom(pol, pol, 2);
+
+ // Put photon on the stack and label it as feedback (51, 52)
+ ++sNfeed;
+ if (ipart == 50000050 && ifeed != 50000052) {
+ weight = 5000;
+ } else {
+ weight = 5000;
+ }
+ gAlice->SetTrack(1, gAlice->CurrentTrack(), gMC->PDGFromId(50),
+ mom,source,pol,position[3],
+ "Cherenkov", nt, weight);
+ }
+ return(nfp);
+}
--- /dev/null
+#ifndef RICHSegResV0_H
+#define RICHSegResV0_H
+
+#include "AliRICH.h"
+#include "AliRICHv0.h"
+
+class AliRICHsegmentationV0 :
+public AliRICHsegmentation {
+ public:
+ AliRICHsegmentationV0(){}
+ virtual ~AliRICHsegmentationV0(){}
+ //
+ // Set Chamber Segmentation Parameters
+ virtual void SetPADSIZ(Float_t p1, Float_t p2);
+ virtual void SetDAnod(Float_t D) {fWireD = D;};
+ //
+ // Transform from pad (wire) to real coordinates and vice versa
+ virtual Float_t GetAnod(Float_t xhit);
+ virtual void GetPadIxy(Float_t x ,Float_t y ,Int_t &ix,Int_t &iy);
+ virtual void GetPadCxy(Int_t ix,Int_t iy,Float_t &x ,Float_t &y );
+ //
+ // Initialisation
+ virtual void Init(AliRICHchamber*);
+ //
+ // Get member data
+ virtual Float_t Dpx(){return fDpx;}
+ virtual Float_t Dpy(){return fDpy;}
+ virtual Int_t Npx(){return fNpx;}
+ virtual Int_t Npy(){return fNpy;}
+ //
+ // Iterate over pads
+ virtual void FirstPad(Float_t xhit, Float_t yhit, Float_t dx, Float_t dy);
+ virtual void NextPad();
+ virtual Int_t MorePads();
+ // Get next neighbours
+ virtual void Neighbours
+ (Int_t iX, Int_t iY, Int_t* Nlist, Int_t Xlist[10], Int_t Ylist[10]);
+ // Provisory RecCluster coordinates reconstructor
+ virtual void FitXY(AliRICHRecCluster* Cluster,TClonesArray* RICHdigits);
+ //
+ // Current Pad during Integration
+ virtual Int_t Ix(){return fix;}
+ virtual Int_t Iy(){return fiy;}
+ virtual Int_t ISector(){return 1;}
+ //
+ // Signal Generation Condition during Stepping
+ virtual Int_t SigGenCond(Float_t x, Float_t y, Float_t z);
+ virtual void SigGenInit(Float_t x, Float_t y, Float_t z);
+ virtual void IntegrationLimits
+ (Float_t& x1, Float_t& x2, Float_t& y1, Float_t& y2);
+ //
+ // Identification
+ virtual char* YourName(){return fName;}
+ ClassDef(AliRICHsegmentationV0,1)
+ protected:
+ //
+ // Implementation of the segmentation data
+ // Version 0 models rectangular pads with the same dimensions all
+ // over the cathode plane
+ //
+ // geometry
+ //
+ Float_t fDpx; // x pad width per sector
+ Float_t fDpy; // y pad base width
+ Int_t fNpx;
+ Int_t fNpy; // Number of pads in y
+ Float_t fWireD; // wire pitch
+
+ // Chamber region consideres during disintegration (lower left and upper right corner)
+ //
+ Int_t fixmin;
+ Int_t fixmax;
+ Int_t fiymin;
+ Int_t fiymax;
+ //
+ // Current pad during integration (cursor for disintegration)
+ Int_t fix;
+ Int_t fiy;
+ Float_t fx;
+ Float_t fy;
+ //
+ // Current pad and wire during tracking (cursor at hit centre)
+ Int_t fixt;
+ Int_t fiyt;
+ Int_t fiwt;
+ Float_t fxt;
+ Float_t fyt;
+ //
+ // Version Identifier
+ char *fName;
+};
+
+class AliRICHresponseV0 : //Mathieson response
+public AliRICHresponse {
+ public:
+ AliRICHresponseV0(){}
+ virtual ~AliRICHresponseV0(){}
+
+
+
+ //
+ // Configuration methods
+ //
+ virtual void SetRSIGM(Float_t p1) {fNsigma=p1;}
+ virtual void SetMUCHSP(Float_t p1) {fChslope=p1;}
+ virtual void SetMUSIGM(Float_t p1, Float_t p2) {fChwX=p1; fChwY=p2;}
+ virtual void SetMAXADC(Float_t p1) {fadc_satm=p1;}
+ // Mathieson parameters
+ virtual void SetSqrtKx3(Float_t p1) {fSqrtKx3=p1;};
+ virtual void SetKx2(Float_t p1) {fKx2=p1;};
+ virtual void SetKx4(Float_t p1) {fKx4=p1;};
+ virtual void SetSqrtKy3(Float_t p1) {fSqrtKy3=p1;};
+ virtual void SetKy2(Float_t p1) {fKy2=p1;};
+ virtual void SetKy4(Float_t p1) {fKy4=p1;};
+ virtual void SetPitch(Float_t p1) {fPitch=p1;};
+
+ //
+ // Get member data
+ virtual Float_t Chslope() {return fChslope;}
+ virtual Float_t ChwX() {return fChwX;}
+ virtual Float_t ChwY() {return fChwY;}
+ virtual Float_t Nsigma() {return fNsigma;}
+ virtual Float_t adc_satm() {return fadc_satm;}
+ //
+ // Chamber response methods
+ // Pulse height from scored quantity (eloss)
+ virtual Float_t IntPH(Float_t eloss);
+ virtual Float_t IntPH(){return 0;}
+ virtual Int_t FeedBackPhotons(Float_t *source, Float_t qtot);
+
+ // Charge disintegration
+ virtual Float_t IntXY(AliRICHsegmentation * segmentation);
+ // Identification
+ //
+ virtual char* YourName() {return fName;}
+
+ ClassDef(AliRICHresponseV0,1)
+ protected:
+ Float_t fChslope; // Slope of the charge distribution
+ Float_t fChwX; // Width of the charge distribution in x
+ Float_t fChwY; // Width of the charge distribution in y
+ Float_t fNsigma; // Number of sigma's used for charge distribution
+ Float_t fadc_satm; // Maximum ADC channel
+ Float_t fSqrtKx3; // Mathieson parameters for x
+ Float_t fKx2;
+ Float_t fKx4;
+ Float_t fSqrtKy3; // Mathieson parameters for y
+ Float_t fKy2;
+ Float_t fKy4;
+ Float_t fPitch; //anode-cathode pitch
+ char *fName; // Version Identifier
+};
+
+#endif
--- /dev/null
+
+//////////////////////////////////////////////////////////////////////////
+// //
+// AliDisplay //
+// //
+// Utility class to display ALICE outline, tracks, hits,.. //
+// //
+//////////////////////////////////////////////////////////////////////////
+
+#include <TROOT.h>
+#include <TTree.h>
+#include <TButton.h>
+#include <TColor.h>
+#include <TCanvas.h>
+#include <TView.h>
+#include <TText.h>
+#include <TPolyMarker3D.h>
+#include <TPolyMarker.h>
+#include <TPaveLabel.h>
+#include <TPaveText.h>
+#include <TList.h>
+#include <TDiamond.h>
+#include <TNode.h>
+#include <TArc.h>
+#include <TTUBE.h>
+#include <TSlider.h>
+#include <TSliderBox.h>
+#include <TGaxis.h>
+#include <TGXW.h>
+#include <TMath.h>
+#include <X3DBuffer.h>
+
+#include "AliRun.h"
+#include "AliDetector.h"
+#include "AliRICH.h"
+#include "AliRICHConst.h"
+#include "AliRICHdisplay.h"
+#include "AliRICHpoints.h"
+#include "TParticle.h"
+
+
+ClassImp(AliRICHdisplay)
+
+
+//____________________________________________________________________________
+AliRICHdisplay::AliRICHdisplay()
+{
+ fPoints = 0;
+ fPhits = 0;
+ fPCerenkovs = 0;
+ fCanvas = 0;
+}
+
+//_____________________________________________________________________________
+AliRICHdisplay::AliRICHdisplay(Int_t size)
+{
+// Create an event display object.
+// A canvas named "edisplay" is created with a vertical size in pixels
+//
+// A QUICK Overview of the Event Display functions
+// ===============================================
+//
+// The event display can ve invoked by executing the macro "display.C"
+// A canvas like in the picture below will appear.
+//
+// On the left side of the canvas, the following buttons appear:
+// *Next* to move to the next event
+// *Previous* to move to the previous event
+
+// *Pick* Select this option to be able to point on a track with the
+// mouse. Once on the track, use the right button to select
+// an action. For example, select SetMarkerAttributes to
+// change the marker type/color/size for the track.
+// *Zoom* Select this option (default) if you want to zoom.
+// To zoom, simply select the selected area with the left button.
+// *UnZoom* To revert to the previous picture size.
+//
+// slider R On the left side, the vertical slider can be used to
+// set the default picture size.
+//
+// When you are in Zoom mode, you can click on the black part of the canvas
+// to select special options with the right mouse button.
+
+//
+// When you are in pick mode, you can "Inspect" the object pointed by the mouse.
+// When you are on a track, select the menu item "InspectParticle"
+// to display the current particle attributes.
+//
+// You can activate the Root browser by selecting the Inspect menu
+// in the canvas tool bar menu. Then select "Start Browser"
+// This will open a new canvas with the browser. At this point, you may want
+// to display some histograms (from the Trees). Go to the "File" menu
+// of the browser and click on "New canvas".
+// In the browser, click on item "ROOT files" in the left pane.
+// Click on galice.root.
+// Click on TH
+// Click on TPC for example
+// Click on any variable (eg TPC.fX) to histogram the variable.
+//
+// If you are lost, you can click on HELP in any Root canvas or browser.
+//Begin_Html
+/*
+ <img src="gif/aliRICHdisplay.gif">
+*/
+//End_Html
+
+
+ fPad = 0;
+
+ gAlice->SetDisplay(this);
+
+ // Initialize display default parameters
+ SetRange();
+
+ // Set front view by default
+ fTheta = 90;
+ fPhi = 90;
+ fPsi = 0;
+ fChamber = 1;
+ fCathode = 1;
+ // fRzone = 1.e10;
+ fDrawClusters = kTRUE;
+ fZoomMode = 1;
+ fZooms = 0;
+ fClustersCuts = 0;
+ fPoints = 0;
+ fPCerenkovs = 0;
+ fPhits = 0;
+ // Create colors
+ CreateColors();
+ // Create display canvas
+ Int_t ysize = size;
+ if (ysize < 100) ysize = 750;
+ Int_t xsize = Int_t(size*830./ysize);
+ fCanvas = new TCanvas("Canvas", "RICH Clusters Display",14,47,xsize,ysize);
+ fCanvas->SetEditable(kIsNotEditable);
+ fCanvas->ToggleEventStatus();
+
+ // Create main display pad
+ fPad = new TPad("viewpad", "RICH display",0.15,0,0.9,1);
+ fPad->Draw();
+ fPad->Modified();
+ fPad->SetFillColor(1);
+ fPad->SetBorderSize(2);
+
+ fCanvas->cd();
+
+ // Create colors pad
+ fColPad = new TPad("colpad", "Colors pad",0.9,0,1,1);
+ fColPad->Draw();
+ fColPad->Modified();
+ fColPad->SetFillColor(19);
+ fColPad->SetBorderSize(2);
+ fColPad->cd();
+ DisplayColorScale();
+
+ fCanvas->cd();
+
+ // Create user interface control pad
+ DisplayButtons();
+ fCanvas->cd();
+
+ // Create Range and mode pad
+ Float_t dxtr = 0.15;
+ Float_t dytr = 0.45;
+ fTrigPad = new TPad("trigger", "range and mode pad",0,0,dxtr,dytr);
+ fTrigPad->Draw();
+ fTrigPad->cd();
+ fTrigPad->SetFillColor(22);
+ fTrigPad->SetBorderSize(2);
+ fRangeSlider = new TSlider("range","range",0.7,0.42,0.9,0.98);
+ fRangeSlider->SetObject(this);
+ char pickmode[] = "gAlice->Display()->SetPickMode()";
+ Float_t db = 0.09;
+ fPickButton = new TButton("Pick",pickmode,0.05,0.32,0.65,0.32+db);
+ fPickButton->SetFillColor(38);
+ fPickButton->Draw();
+ char zoommode[] = "gAlice->Display()->SetZoomMode()";
+ fZoomButton = new TButton("Zoom",zoommode,0.05,0.21,0.65,0.21+db);
+ fZoomButton->SetFillColor(38);
+ fZoomButton->Draw();
+ fArcButton = new TArc(.8,fZoomButton->GetYlowNDC()+0.5*db,0.33*db);
+ fArcButton->SetFillColor(kGreen);
+ fArcButton->Draw();
+ char butUnzoom[] = "gAlice->Display()->UnZoom()";
+ TButton *button = new TButton("UnZoom",butUnzoom,0.05,0.05,0.95,0.15);
+ button->SetFillColor(38);
+ button->Draw();
+ AppendPad(); // append display object as last object to force selection
+
+ fCanvas->cd();
+ fCanvas->Update();
+}
+
+
+//_____________________________________________________________________________
+AliRICHdisplay::~AliRICHdisplay()
+{
+ // Delete space point structure
+ if (fPoints) fPoints->Delete();
+ delete fPoints;
+ fPoints = 0;
+ //
+ if (fPhits) fPhits->Delete();
+ delete fPhits;
+ fPhits = 0;
+ //
+ if (fPCerenkovs) fPCerenkovs->Delete();
+ delete fPCerenkovs;
+ fPCerenkovs = 0;
+}
+
+//_____________________________________________________________________________
+void AliRICHdisplay::Clear(Option_t *)
+{
+// Delete graphics temporary objects
+}
+
+//_____________________________________________________________________________
+void AliRICHdisplay::DisplayButtons()
+{
+// Create the user interface buttons
+
+
+ fButtons = new TPad("buttons", "newpad",0,0.45,0.15,1);
+ fButtons->Draw();
+ fButtons->SetFillColor(38);
+ fButtons->SetBorderSize(2);
+ fButtons->cd();
+
+ // Int_t butcolor = 33;
+ Float_t dbutton = 0.08;
+ Float_t y = 0.96;
+ Float_t dy = 0.014;
+ Float_t x0 = 0.05;
+ Float_t x1 = 0.95;
+
+ TButton *button;
+ char but1[] = "gAlice->Display()->ShowNextEvent(1)";
+ button = new TButton("Next",but1,x0,y-dbutton,x1,y);
+ button->SetFillColor(38);
+ button->Draw();
+
+ y -= dbutton +dy;
+ char but2[] = "gAlice->Display()->ShowNextEvent(-1)";
+ button = new TButton("Previous",but2,x0,y-dbutton,x1,y);
+ button->SetFillColor(38);
+ button->Draw();
+
+ // display logo
+ TDiamond *diamond = new TDiamond(0.05,0.015,0.95,0.22);
+ diamond->SetFillColor(50);
+ diamond->SetTextAlign(22);
+ diamond->SetTextColor(5);
+ diamond->SetTextSize(0.11);
+ diamond->Draw();
+ diamond->AddText(".. ");
+ diamond->AddText("ROOT");
+ diamond->AddText("RICH");
+ diamond->AddText("... ");
+ diamond->AddText(" ");
+}
+
+//_____________________________________________________________________________
+void AliRICHdisplay::CreateColors()
+{
+// Create the colors palette used to display clusters
+
+ Int_t k,i;
+ Int_t color;
+ Float_t r,g,b;
+
+ for (k=1;k<=5;k++) {
+ switch(k) {
+ case 1:
+ for (i=1;i<=5;i++) {
+ r=1.;
+ g=i*0.2;
+ b=0.;
+ color=i;
+ color=50+23-color;
+ new TColor(color,r,g,b);
+ }
+ break;
+ case 2:
+ for (i=1;i<=4;i++) {
+ r=1.1-i*0.2;
+ g=1.;
+ b=0.;
+ color=i+5;
+ color=50+23-color;
+ new TColor(color,r,g,b);
+ }
+ break;
+ case 3:
+ for (i=1;i<=4;i++) {
+ r=0.;
+ g=1.;
+ b=i*0.2+0.2;
+ color=i+9;
+ color=50+23-color;
+ new TColor(color,r,g,b);
+ }
+ break;
+ case 4:
+ for (i=1;i<=4;i++) {
+ r=0.;
+ g=1.1-i*0.2;
+ b=1.;
+ color=i+13;
+ color=50+23-color;
+ new TColor(color,r,g,b);
+ }
+ break;
+ case 5:
+ for (i=1;i<=5;i++) {
+ r=i*0.2;
+ g=0.;
+ b=1.;
+ color=i+17;
+ color=50+23-color;
+ new TColor(color,r,g,b);
+ }
+ break;
+ }
+
+ }
+
+}
+
+//_____________________________________________________________________________
+void AliRICHdisplay::DisplayColorScale()
+{
+
+ Int_t i;
+ Int_t color;
+ Float_t xlow, ylow, xup, yup, hs;
+ Float_t x1, y1, x2, y2;
+ x1 = y1 = 0;
+ x2 = y2 = 20;
+
+ gPad->SetFillColor(0);
+ gPad->Clear();
+ gPad->Range(x1,y1,x2,y2);
+ TText *text = new TText(0,0,"");
+ text->SetTextFont(61);
+ text->SetTextSize(0.03);
+ text->SetTextAlign(22);
+
+ TBox *box;
+ char label[8];
+//*-* draw colortable boxes
+ hs = (y2-y1)/Float_t(22);
+ xlow=x1+1;
+ xup=x2-9;
+ for (i=0;i<22;i++) {
+ ylow = y1 + hs*(Float_t(i));
+ yup = y1 + hs*(Float_t(i+1));
+ color = 51+i;
+ Double_t logscale=Double_t(i+1)*(TMath::Log(adc_satm)/22);
+ Int_t scale=(Int_t)TMath::Exp(logscale);
+ sprintf(label,"%d",scale);
+ box = new TBox(xlow, ylow, xup, yup);
+ box->SetFillColor(color);
+ box->Draw();
+ text->DrawText(xup+4, 0.5*(ylow+yup),label);
+ }
+}
+
+//______________________________________________________________________________
+Int_t AliRICHdisplay::DistancetoPrimitive(Int_t px, Int_t)
+{
+// Compute distance from point px,py to objects in event
+
+ gPad->SetCursor(kCross);
+
+ if (gPad == fTrigPad) return 9999;
+
+ const Int_t big = 9999;
+ Int_t dist = big;
+ Float_t xmin = gPad->GetX1();
+ Float_t xmax = gPad->GetX2();
+ Float_t dx = 0.02*(xmax - xmin);
+ Float_t x = gPad->AbsPixeltoX(px);
+ if (x < xmin+dx || x > xmax-dx) return dist;
+
+ if (fZoomMode) return 0;
+ else return 7;
+}
+
+//_____________________________________________________________________________
+void AliRICHdisplay::Draw(Option_t *)
+{
+// Display current event
+
+ fPad->cd();
+
+ DrawView(fTheta, fPhi, fPsi); // see how to draw PGON+inner frames
+
+ // Display the event number and title
+ fPad->cd();
+ DrawTitle();
+}
+
+
+//_____________________________________________________________________________
+
+void AliRICHdisplay::DrawCerenkovs()
+{
+// Draw cerenkovs hits for RICH chambers
+
+ LoadCerenkovs(fChamber);
+ printf("\nDrawCerenkovs\n");
+
+ Int_t ntracks, track;
+ TObjArray *cpoints;
+ AliRICHpoints *pm;
+
+ fHitsCuts = 0;
+ cpoints = fPCerenkovs;
+ printf ("Cpoints: %p",cpoints);
+ if (!cpoints) return;
+ ntracks = cpoints->GetEntriesFast();
+ printf("DrawCerenkovs - ntracks %d \n",ntracks);
+ for (track=0;track<ntracks;track++) {
+ pm = (AliRICHpoints*)cpoints->UncheckedAt(track);
+ if (!pm) continue;
+ pm->Draw();
+ fHitsCuts += pm->GetN();
+ }
+}
+
+//_____________________________________________________________________________
+void AliRICHdisplay::DrawClusters()
+{
+// Draw clusterss for RICH chambers
+
+ Int_t ndigits, digit;
+ TObjArray *points;
+ AliRICHpoints *pm;
+
+ fClustersCuts = 0;
+ points = fPoints;
+ if (!points) return;
+ ndigits = points->GetEntriesFast();
+ printf("DrawClusters - ndigits %d \n",ndigits);
+ for (digit=0;digit<ndigits;digit++){
+ pm = (AliRICHpoints*)points->UncheckedAt(digit);
+ if (!pm) continue;
+ pm->Draw();
+ fClustersCuts +=pm->GetN();
+ }
+}
+
+//_____________________________________________________________________________
+void AliRICHdisplay::DrawHits()
+{
+// Draw hits for RICH chambers
+
+ LoadHits(fChamber);
+
+ Int_t ntracks, track;
+ TObjArray *points;
+ AliRICHpoints *pm;
+
+ fHitsCuts = 0;
+ points = Phits();
+ if (!points) return;
+ ntracks = points->GetEntriesFast();
+ printf("DrawHits - ntracks %d \n",ntracks);
+ for (track=0;track<ntracks;track++) {
+ pm = (AliRICHpoints*)points->UncheckedAt(track);
+ if (!pm) continue;
+ pm->Draw();
+ fHitsCuts += pm->GetN();
+ }
+}
+
+
+//_____________________________________________________________________________
+void AliRICHdisplay::DrawTitle(Option_t *option)
+{
+// Draw the event title
+
+ Float_t xmin = gPad->GetX1();
+ Float_t xmax = gPad->GetX2();
+ Float_t ymin = gPad->GetY1();
+ Float_t ymax = gPad->GetY2();
+ Float_t dx = xmax-xmin;
+ Float_t dy = ymax-ymin;
+
+ if (strlen(option) == 0) {
+ TPaveText *title = new TPaveText(xmin +0.01*dx, ymax-0.09*dy, xmin +0.5*dx, ymax-0.01*dy);
+ title->SetBit(kCanDelete);
+ title->SetFillColor(42);
+ title->Draw();
+ char ptitle[100];
+ sprintf(ptitle,"Alice event: %d, Run:%d",gAlice->GetHeader()->GetEvent(), gAlice->GetHeader()->GetRun());
+ title->AddText(ptitle);
+ Int_t nparticles = gAlice->Particles()->GetEntriesFast();
+ sprintf(ptitle,"Nparticles = %d Nhits = %d Npads fired = %d",nparticles, fHitsCuts,fClustersCuts);
+ title->AddText(ptitle);
+ } else {
+ TPaveLabel *label = new TPaveLabel(xmin +0.01*dx, ymax-0.07*dy, xmin +0.2*dx, ymax-0.01*dy,option);
+ label->SetBit(kCanDelete);
+ label->SetFillColor(42);
+ label->Draw();
+ }
+}
+
+//_____________________________________________________________________________
+void AliRICHdisplay::DrawView(Float_t theta, Float_t phi, Float_t psi)
+{
+// Draw a view of RICH clusters
+
+ gPad->SetCursor(kWatch);
+ gPad->SetFillColor(1);
+ gPad->Clear();
+
+ Int_t iret;
+ TView *view = new TView(1);
+ Float_t range = fRrange*fRangeSlider->GetMaximum();
+ view->SetRange(-range,-range,-range,range, range, range);
+ fZoomX0[0] = -1;
+ fZoomY0[0] = -1;
+ fZoomX1[0] = 1;
+ fZoomY1[0] = 1;
+ fZooms = 0;
+
+ //Display RICH Chamber Geometry
+ gAlice->GetGeometry()->Draw("same");
+
+ //add clusters to the pad
+ DrawClusters();
+ DrawHits();
+ DrawCerenkovs();
+
+ // add itself to the list (must be last)
+ AppendPad();
+
+ view->SetView(phi, theta, psi, iret);
+}
+
+
+//______________________________________________________________________________
+void AliRICHdisplay::ExecuteEvent(Int_t event, Int_t px, Int_t py)
+{
+// Execute action corresponding to the mouse event
+
+ static Float_t x0, y0, x1, y1;
+
+ static Int_t pxold, pyold;
+ static Int_t px0, py0;
+ static Int_t linedrawn;
+ Float_t temp;
+
+ if (px == 0 && py == 0) { //when called by sliders
+ if (event == kButton1Up) {
+ Draw();
+ }
+ return;
+ }
+ if (!fZoomMode && gPad->GetView()) {
+ gPad->GetView()->ExecuteRotateView(event, px, py);
+ return;
+ }
+
+ // something to zoom ?
+ gPad->SetCursor(kCross);
+
+ switch (event) {
+
+ case kButton1Down:
+ gGXW->SetLineColor(-1);
+ gPad->TAttLine::Modify(); //Change line attributes only if necessary
+ x0 = gPad->AbsPixeltoX(px);
+ y0 = gPad->AbsPixeltoY(py);
+ px0 = px; py0 = py;
+ pxold = px; pyold = py;
+ linedrawn = 0;
+ return;
+
+ case kButton1Motion:
+ if (linedrawn) gGXW->DrawBox(px0, py0, pxold, pyold, TGXW::kHollow);
+ pxold = px;
+ pyold = py;
+ linedrawn = 1;
+ gGXW->DrawBox(px0, py0, pxold, pyold, TGXW::kHollow);
+ return;
+
+ case kButton1Up:
+ gPad->GetCanvas()->FeedbackMode(kFALSE);
+ if (px == px0) return;
+ if (py == py0) return;
+ x1 = gPad->AbsPixeltoX(px);
+ y1 = gPad->AbsPixeltoY(py);
+
+ if (x1 < x0) {temp = x0; x0 = x1; x1 = temp;}
+ if (y1 < y0) {temp = y0; y0 = y1; y1 = temp;}
+ gPad->Range(x0,y0,x1,y1);
+ if (fZooms < kMAXZOOM-1) {
+ fZooms++;
+ fZoomX0[fZooms] = x0;
+ fZoomY0[fZooms] = y0;
+ fZoomX1[fZooms] = x1;
+ fZoomY1[fZooms] = y1;
+ }
+ gPad->Modified(kTRUE);
+ return;
+ }
+
+}
+
+//___________________________________________
+void AliRICHdisplay::LoadDigits()
+{
+// Read digits info and store x,y,z info in arrays fPoints
+// Loop on all detectors
+
+ printf("Entering Load-digits");
+
+
+ ResetPoints();
+ AliRICH *RICH = (AliRICH*)gAlice->GetDetector("RICH");
+ AliRICHchamber* iChamber;
+ AliRICHsegmentation* segmentation;
+ Int_t NallDigits=0;
+
+ for (Int_t ich=0; ich<7; ich++) {
+ TClonesArray *RICHdigits = RICH->DigitsAddress(ich);
+ if (RICHdigits == 0) continue;
+ gAlice->ResetDigits();
+ gAlice->TreeD()->GetEvent(1);
+ Int_t ndigits = RICHdigits->GetEntriesFast();
+ printf("Found %d digits in chamber %d \n",ndigits,ich);
+ NallDigits+=ndigits;
+ }
+ if (fPoints == 0) fPoints = new TObjArray(NallDigits);
+ Int_t counter=0;
+ for (Int_t ich=0; ich<7; ich++) {
+ TClonesArray *RICHdigits = RICH->DigitsAddress(ich);
+ if (RICHdigits == 0) continue;
+ gAlice->ResetDigits();
+ gAlice->TreeD()->GetEvent(1);
+ Int_t ndigits = RICHdigits->GetEntriesFast();
+ printf("Found %d digits in chamber %d \n",ndigits,ich);
+ if (ndigits == 0) continue;
+
+
+
+ iChamber = &(RICH->Chamber(ich));
+ printf("LoadPoints - chamber %d \n",ich);
+ segmentation=iChamber->GetSegmentationModel(1);
+ Float_t dpx = segmentation->Dpx();
+ Float_t dpy = segmentation->Dpy();
+ printf("LoadPoints - dpx, dpy %f %f \n",dpx,dpy);
+ AliRICHdigit *mdig;
+ AliRICHpoints *points = 0;
+ //
+ //loop over all digits and store their position
+ Int_t npoints=1;
+
+ for (Int_t digit=0;digit<ndigits;digit++) {
+ mdig = (AliRICHdigit*)RICHdigits->UncheckedAt(digit);
+ points = new AliRICHpoints(npoints);
+ fPoints->AddAt(points,counter);
+ counter++;
+ Int_t charge=mdig->fSignal;
+ Int_t index=Int_t(TMath::Log(charge)/(TMath::Log(adc_satm)/22));
+ Int_t color=51+index;
+ if (color>72) color=72;
+ points->SetMarkerColor(color);
+ points->SetMarkerStyle(21);
+ points->SetMarkerSize(0.5);
+ // get the center of the pad - add on x and y half of pad size
+ Float_t xpad, ypad;
+ segmentation->GetPadCxy(mdig->fPadX, mdig->fPadY,xpad, ypad);
+ printf("\n chamber x,y, %d %f %f ", ich, xpad, ypad);
+
+ Float_t VecLoc[3]={xpad,0,ypad};
+ Float_t VecGlob[3];
+ iChamber->LocaltoGlobal(VecLoc,VecGlob);
+ points->SetParticle(-1);
+ points->SetHitIndex(-1);
+ points->SetTrackIndex(-1);
+ points->SetDigitIndex(digit);
+ points->SetPoint(0,VecGlob[0],VecGlob[1],VecGlob[2]);
+ } // loop over digits
+ } // loop over chambers
+}
+
+
+//___________________________________________
+void AliRICHdisplay::LoadHits(Int_t chamber)
+{
+// Read hits info and store x,y,z info in arrays fPhits
+// Loop on all detectors
+
+ printf("Entering Load-hits");
+
+ fChamber=chamber;
+ ResetPhits();
+
+ AliRICH *RICH = (AliRICH*)gAlice->GetDetector("RICH");
+ AliRICHchamber* iChamber;
+
+ iChamber = &(RICH->Chamber(chamber-1));
+ Float_t zpos=iChamber->ZPosition();
+ printf("LoadHits - zpos %f \n",zpos);
+
+ Int_t ntracks = (Int_t)gAlice->TreeH()->GetEntries();
+ printf("ntracks %d\n",ntracks);
+ Int_t ntrks = gAlice->GetNtrack();
+ printf("ntrks %d\n",ntrks);
+
+ if (fPhits == 0) fPhits = new TObjArray(ntracks);
+ TVector *xp = new TVector(1000);
+ TVector *yp = new TVector(1000);
+ TVector *zp = new TVector(1000);
+ TVector *ptrk = new TVector(1000);
+ TVector *phit = new TVector(1000);
+ for (Int_t track=0; track<ntracks;track++) {
+ gAlice->ResetHits();
+ gAlice->TreeH()->GetEvent(track);
+ TClonesArray *RICHhits = RICH->Hits();
+ if (RICHhits == 0) return;
+ Int_t nhits = RICHhits->GetEntriesFast();
+ if (nhits == 0) continue;
+ AliRICHhit *mHit;
+ AliRICHpoints *points = 0;
+ Int_t npoints=0;
+ for (Int_t hit=0;hit<nhits;hit++) {
+ mHit = (AliRICHhit*)RICHhits->UncheckedAt(hit);
+ (*xp)(npoints)=mHit->fX;
+ (*yp)(npoints)=mHit->fY;
+ (*zp)(npoints)=mHit->fZ;
+ (*ptrk)(npoints)=Float_t(mHit->GetTrack());
+ (*phit)(npoints)=Float_t(hit);
+ npoints++;
+ }
+
+ if (npoints == 0) continue;
+ points = new AliRICHpoints(npoints);
+ for (Int_t p=0;p<npoints;p++) {
+ points->SetMarkerColor(kRed);
+ points->SetMarkerStyle(5);
+ points->SetMarkerSize(1.);
+ points->SetParticle(Int_t((*ptrk)(p)));
+ points->SetHitIndex(Int_t((*phit)(p)));
+ points->SetTrackIndex(track);
+ points->SetDigitIndex(-1);
+ points->SetPoint(p,(*xp)(p),(*yp)(p),(*zp)(p));
+ }
+ xp->Zero();
+ yp->Zero();
+ ptrk->Zero();
+ phit->Zero();
+ fPhits->AddAt(points,track);
+ }
+
+}
+
+//_____________________________________________________________________________
+
+void AliRICHdisplay::LoadCerenkovs(Int_t chamber)
+{
+// Read cerenkov hits info and store x,y,z info in array fPCerenkovs
+// Loop on all detectors
+
+ printf("Entering Load-Cerenkovs");
+
+ fChamber=chamber;
+ ResetPCerenkovs();
+
+ AliRICH *RICH = (AliRICH*)gAlice->GetDetector("RICH");
+ AliRICHchamber* iChamber;
+
+ iChamber = &(RICH->Chamber(chamber-1));
+ Float_t zpos=iChamber->ZPosition();
+ printf("LoadCerenkovs - zpos %f \n",zpos);
+
+ RICH->SetTreeAddress();
+ Int_t ntracks = (Int_t)gAlice->TreeH()->GetEntries();
+ printf("ntracks %d\n",ntracks);
+ Int_t ntrks = gAlice->GetNtrack();
+ printf("ntrks %d\n",ntrks);
+
+ if (fPCerenkovs == 0) fPCerenkovs = new TObjArray(ntracks);
+ TVector *xp = new TVector(1000);
+ TVector *yp = new TVector(1000);
+ TVector *zp = new TVector(1000);
+ TVector *ptrk = new TVector(1000);
+ TVector *phit = new TVector(1000);
+ for (Int_t track=0; track<ntracks;track++) {
+ gAlice->ResetHits();
+ gAlice->TreeH()->GetEvent(track);
+ TClonesArray *RICHCerenkovs = RICH->Cerenkovs();
+ printf("RICHCerenkovs %p\n",RICHCerenkovs);
+ if (RICHCerenkovs == 0) return;
+ Int_t nhits = RICHCerenkovs->GetEntriesFast();
+ if (nhits == 0) continue;
+ printf("nhits %d \n",nhits);
+ AliRICHCerenkov *mCerenkov;
+ AliRICHpoints *cpoints = 0;
+ Int_t npoints=0;
+
+//Display Cerenkov hits in blue
+
+ for (Int_t hit=0;hit<nhits;hit++) {
+ mCerenkov = (AliRICHCerenkov*)RICHCerenkovs->UncheckedAt(hit);
+ (*xp)(npoints)=mCerenkov->fX;
+ (*yp)(npoints)=mCerenkov->fY;
+ (*zp)(npoints)=mCerenkov->fZ;
+ (*ptrk)(npoints)=Float_t(mCerenkov->GetTrack());
+ (*phit)(npoints)=Float_t(hit);
+ printf("track, trk %d %d\n",track,mCerenkov->GetTrack());
+ npoints++;
+ }
+ if (npoints == 0) continue;
+ printf("npoints %d \n",npoints);
+ cpoints = new AliRICHpoints(npoints);
+ for (Int_t p=0;p<npoints;p++) {
+ cpoints->SetMarkerColor(kBlue);
+ cpoints->SetMarkerStyle(3);
+ cpoints->SetMarkerSize(1.);
+ cpoints->SetParticle(Int_t((*ptrk)(p)));
+ Int_t index=cpoints->GetIndex();
+ printf("index %d \n",index);
+ cpoints->SetHitIndex(Int_t((*phit)(p)));
+ cpoints->SetTrackIndex(track);
+ cpoints->SetDigitIndex(-1);
+ cpoints->SetPoint(p,(*xp)(p),(*yp)(p),(*zp)(p));
+ }
+ xp->Zero();
+ yp->Zero();
+ ptrk->Zero();
+ phit->Zero();
+ fPCerenkovs->AddAt(cpoints,track);
+ }
+}
+
+//_____________________________________________________________________________
+void AliRICHdisplay::Paint(Option_t *)
+{
+// Paint miscellaneous items
+
+}
+
+//_____________________________________________________________________________
+void AliRICHdisplay::SetPickMode()
+{
+ fZoomMode = 0;
+
+ fArcButton->SetY1(fPickButton->GetYlowNDC()+0.5*fPickButton->GetHNDC());
+ fTrigPad->Modified();
+}
+
+//_____________________________________________________________________________
+void AliRICHdisplay::SetZoomMode()
+{
+ fZoomMode = 1;
+
+ fArcButton->SetY1(fZoomButton->GetYlowNDC()+0.5*fZoomButton->GetHNDC());
+ fTrigPad->Modified();
+}
+
+//_____________________________________________________________________________
+void AliRICHdisplay::SetChamberAndCathode(Int_t chamber, Int_t cathode)
+{
+// Set chamber and cathode number
+ fChamber = chamber;
+ fCathode = cathode;
+
+ printf("SetChamberAndCathode - fChamber fCathode %d %d\n",fChamber,fCathode);
+ if (!fPad) return;
+ fPad->Clear();
+ LoadDigits();
+ Draw();
+}
+
+//_____________________________________________________________________________
+void AliRICHdisplay::SetRange(Float_t rrange, Float_t zrange)
+{
+// Set view range along R and Z
+ fRrange = rrange;
+ fZrange = zrange;
+
+ if (!fPad) return;
+ fPad->Clear();
+ Draw();
+}
+
+//_____________________________________________________________________________
+void AliRICHdisplay::SetView(Float_t theta, Float_t phi, Float_t psi)
+{
+// change viewing angles for current event
+
+ fPad->cd();
+ fPhi = phi;
+ fTheta = theta;
+ fPsi = psi;
+ Int_t iret = 0;
+
+ TView *view = gPad->GetView();
+ if (view) view->SetView(fPhi, fTheta, fPsi, iret);
+ else Draw();
+
+ gPad->Modified();
+}
+
+//_____________________________________________________________________________
+void AliRICHdisplay::ShowNextEvent(Int_t delta)
+{
+// Display (current event_number+delta)
+// delta = 1 shown next event
+// delta = -1 show previous event
+
+ if (delta) {
+ gAlice->Clear();
+ Int_t current_event = gAlice->GetHeader()->GetEvent();
+ Int_t new_event = current_event + delta;
+ gAlice->GetEvent(new_event);
+ if (!gAlice->TreeD()) return;
+ }
+ LoadDigits();
+ DrawClusters();
+ fPad->cd();
+ Draw();
+
+
+}
+
+//______________________________________________________________________________
+void AliRICHdisplay::UnZoom()
+{
+ if (fZooms <= 0) return;
+ fZooms--;
+ TPad *pad = (TPad*)gPad->GetPadSave();
+ pad->Range(fZoomX0[fZooms],fZoomY0[fZooms], fZoomX1[fZooms],fZoomY1[fZooms]);
+ pad->Modified();
+}
+
+//_____________________________________________________________________________
+void AliRICHdisplay::ResetPoints()
+{
+ //
+ // Reset array of points
+ //
+ if (fPoints) {
+ fPoints->Delete();
+ delete fPoints;
+ fPoints = 0;
+ }
+}
+//_____________________________________________________________________________
+void AliRICHdisplay::ResetPhits()
+{
+ //
+ // Reset array of points
+ //
+ if (fPhits) {
+ fPhits->Delete();
+ delete fPhits;
+ fPhits = 0;
+ }
+}
+//_____________________________________________________________________________
+void AliRICHdisplay::ResetPCerenkovs()
+{
+ //
+ // Reset array of points
+ //
+ if (fPCerenkovs) {
+ fPCerenkovs->Delete();
+ delete fPCerenkovs;
+ fPCerenkovs = 0;
+ }
+}
--- /dev/null
+#ifndef AliRICHdisplay_H
+#define AliRICHdisplay_H
+
+//////////////////////////////////////////////////////////////////////////
+// //
+// AliDisplay //
+// //
+// Utility class to display ALice outline, tracks, hits,.. //
+// //
+//////////////////////////////////////////////////////////////////////////
+
+//#ifndef ROOT_TObject
+#include <TObject.h>
+//#endif
+#include "AliDisplay.h"
+
+class TCanvas;
+class TPad;
+class TList;
+class TSlider;
+class TButton;
+class TArc;
+
+const Int_t kMAXZOOM = 20;
+
+class AliRICHdisplay : /*splaypublic TObject,*/ public AliDisplay {
+
+ private:
+ Int_t fChamber;
+ Int_t fCathode;
+ Int_t fZoomMode; //=1 if in zoom mode
+
+ Bool_t fDrawClusters; //Flag True if Clusters to be drawn
+ Float_t fTheta; //Viewing angle theta
+ Float_t fPhi; //Viewing angle phi
+ Float_t fPsi; //Viewving angle psi (rotation on display)
+ Float_t fRrange; //Size of view in R
+ Float_t fZrange; //Size of view along Z
+ Float_t fZoomX0[20]; //Low x range of zoom number i
+ Float_t fZoomY0[20]; //Low y range of zoom number i
+ Float_t fZoomX1[20]; //High x range of zoom number i
+ Float_t fZoomY1[20]; //High y range of zoom number i
+ Int_t fZooms; //Number of zooms
+ Int_t fHitsCuts; //Number of hits surviving cuts
+ Int_t fClustersCuts; //Number of clusters surviving cuts
+ TCanvas *fCanvas; //Pointer to the display canvas
+ TPad *fTrigPad; //Pointer to the trigger pad
+ TPad *fColPad; //Pointer to the colors pad
+ TPad *fButtons; //Pointer to the buttons pad
+ TPad *fPad; //Pointer to the event display main pad
+ TSlider *fRangeSlider; //Range slider
+ TButton *fPickButton; //Button to activate Pick mode
+ TButton *fZoomButton; //Button to activate Zoom mode
+ TArc *fArcButton; //Gren/Red button to show Pick/Zoom mode
+ TObjArray *fPoints; //Array of points for each cathode
+ TObjArray *fPhits; //Array of hit points for each chamber
+ TObjArray *fPCerenkovs; //Array of cerenkov hits for each chamber
+ public:
+ AliRICHdisplay();
+ AliRICHdisplay(Int_t size);
+ virtual ~AliRICHdisplay();
+ virtual void Clear(Option_t *option="");
+ virtual void DisplayButtons();
+ virtual void CreateColors();
+ virtual void DisplayColorScale();
+ virtual Int_t DistancetoPrimitive(Int_t px, Int_t py);
+ virtual void Draw(Option_t *option="");
+ virtual void DrawCerenkovs();
+ virtual void DrawClusters();
+ virtual void DrawHits();
+ virtual void DrawTitle(Option_t *option="");
+ virtual void DrawView(Float_t theta, Float_t phi, Float_t psi=0);
+ virtual void ExecuteEvent(Int_t event, Int_t px, Int_t py);
+ Int_t GetZoomMode() {return fZoomMode;}
+ Int_t GetChamber() {return fChamber;}
+ Int_t GetCathode() {return fCathode;}
+ virtual void LoadDigits();
+ virtual void LoadCerenkovs(Int_t chamber);
+ virtual void LoadHits(Int_t chamber);
+ TPad *Pad() {return fPad;}
+ TObjArray *Points() {return fPoints;}
+ TObjArray *Phits() {return fPhits;}
+ TObjArray *PCerenkovs() {return fPCerenkovs;}
+ virtual void Paint(Option_t *option="");
+ virtual void SetDrawClusters(Bool_t draw=kTRUE) {fDrawClusters=draw;} // *MENU*
+ virtual void SetChamberAndCathode(Int_t chamber=1, Int_t cathode=1); // *MENU*
+ virtual void SetRange(Float_t rrange=250, Float_t zrange=1050); // *MENU*
+ virtual void SetView(Float_t theta, Float_t phi, Float_t psi=0);
+ virtual void SetPickMode();
+ virtual void SetZoomMode();
+ virtual void ShowNextEvent(Int_t delta=1);
+ virtual void UnZoom(); // *MENU*
+ virtual void ResetPoints();
+ virtual void ResetPhits();
+ virtual void ResetPCerenkovs();
+
+
+ ClassDef(AliRICHdisplay, 0) //Utility class to display RICH clusters...
+ };
+
+#endif
+
--- /dev/null
+///////////////////////////////////////////////////////////////////////////////
+// //
+// This class contains the points for the ALICE event display //
+// //
+//Begin_Html
+/*
+<img src="gif/AliRICHpointsClass.gif">
+*/
+//End_Html
+// //
+// //
+///////////////////////////////////////////////////////////////////////////////
+#include "AliRICHdisplay.h"
+#include "AliRICHpoints.h"
+#include "AliRun.h"
+#include "TPad.h"
+#include "TView.h"
+#include "TMath.h"
+
+const Int_t MAX_Nipx=400, MAX_Nipy=800;
+
+ClassImp(AliRICHpoints)
+
+//_____________________________________________________________________________
+AliRICHpoints::AliRICHpoints()
+{
+ //
+ // Default constructor
+ //
+ fHitIndex = 0;
+ fTrackIndex = 0;
+ fDigitIndex = 0;
+}
+
+//_____________________________________________________________________________
+AliRICHpoints::AliRICHpoints(Int_t npoints)
+ :AliPoints(npoints)
+{
+ //
+ // Standard constructor
+ //
+ fHitIndex = 0;
+ fTrackIndex = 0;
+ fDigitIndex = 0;
+}
+
+//_____________________________________________________________________________
+AliRICHpoints::~AliRICHpoints()
+{
+ //
+ // Default destructor
+ //
+ fHitIndex = 0;
+ fTrackIndex = 0;
+ fDigitIndex = 0;
+}
+
+//_____________________________________________________________________________
+void AliRICHpoints::DumpHit()
+{
+ //
+ // Dump hit corresponding to this point
+ //
+ AliRICHhit *hit = GetHit();
+ if (hit) hit->Dump();
+}
+
+//_____________________________________________________________________________
+void AliRICHpoints::DumpDigit()
+{
+ //
+ // Dump digit corresponding to this point
+ //
+ AliRICHdigit *digit = GetDigit();
+ if (digit) digit->Dump();
+}
+
+//_____________________________________________________________________________
+void AliRICHpoints::InspectHit()
+{
+ //
+ // Inspect hit corresponding to this point
+ //
+ AliRICHhit *hit = GetHit();
+ if (hit) hit->Inspect();
+}
+
+//_____________________________________________________________________________
+void AliRICHpoints::InspectDigit()
+{
+ //
+ // Inspect digit corresponding to this point
+ //
+ AliRICHdigit *digit = GetDigit();
+ if (digit) digit->Inspect();
+}
+
+//_____________________________________________________________________________
+Int_t AliRICHpoints::GetTrackIndex()
+{
+ //
+ // Dump digit corresponding to this point
+ //
+ printf("GetTrackIndex - fTrackIndex %d \n",fTrackIndex);
+ this->Inspect();
+ return fTrackIndex;
+}
+
+//_____________________________________________________________________________
+AliRICHhit *AliRICHpoints::GetHit() const
+{
+ //
+ // Returns pointer to hit index in AliRun::fParticles
+ //
+ AliRICH *RICH = (AliRICH*)gAlice->GetDetector("RICH");
+ gAlice->TreeH()->GetEvent(fTrackIndex);
+ TClonesArray *RICHhits = RICH->Hits();
+ Int_t nhits = RICHhits->GetEntriesFast();
+ if (fHitIndex < 0 || fHitIndex >= nhits) return 0;
+ return (AliRICHhit*)RICHhits->UncheckedAt(fHitIndex);
+}
+
+//_____________________________________________________________________________
+AliRICHdigit *AliRICHpoints::GetDigit() const
+{
+ //
+ // Returns pointer to digit index in AliRun::fParticles
+ //
+
+ AliRICHdisplay *display=(AliRICHdisplay*)gAlice->Display();
+ Int_t chamber=display->GetChamber();
+ Int_t cathode=display->GetCathode();
+
+ AliRICH *RICH = (AliRICH*)gAlice->GetDetector("RICH");
+ TClonesArray *RICHdigits = RICH->DigitsAddress(chamber-1);
+ gAlice->TreeD()->GetEvent(cathode);
+ Int_t ndigits = RICHdigits->GetEntriesFast();
+ if (fDigitIndex < 0 || fDigitIndex >= ndigits) return 0;
+ return (AliRICHdigit*)RICHdigits->UncheckedAt(fDigitIndex);
+}
+//_____________________________________________________________________________
+struct Bin {
+ const AliRICHdigit *dig;
+ int idx;
+ Bin() {dig=0; idx=-1;}
+};
+
+struct PreCluster : public AliRICHreccluster {
+ const AliRICHdigit* summit;
+ int idx;
+ int cut;
+ int npeaks;
+ PreCluster() : AliRICHreccluster() {cut=npeaks=0;}
+};
+//_____________________________________________________________________________
+
+static void FindCluster(AliRICHchamber *iChamber, AliRICHsegmentation *segmentation, int i, int j, Bin bins[MAX_Nipx][MAX_Nipy], PreCluster &c)
+
+{
+
+ //
+ // Find clusters
+ //
+
+ printf("I'm in FindCluster \n");
+
+ Bin& b=bins[i][j];
+ Int_t q=b.dig->fSignal;
+
+ printf("FindCluster - i j q %d %d %d\n",i,j,q);
+
+ if (q<0) {
+ q=-q;
+ c.cut=1;
+ }
+ if (b.idx >= 0 && b.idx != c.idx) {
+ c.idx=b.idx;
+ c.npeaks++;
+ }
+
+ if (q > TMath::Abs(c.summit->fSignal)) c.summit=b.dig;
+
+ Int_t npx = segmentation->Npx();
+ Int_t npy = segmentation->Npy();
+ Float_t x,y;
+ segmentation->GetPadCxy(i-npx, j-npy, x,y);
+ printf("FindCluster - x y %f %f \n",x,y);
+
+
+ c.fX += q*x;
+ c.fY += q*y;
+ c.fQ += q;
+
+ b.dig = 0; b.idx = c.idx;
+
+ if (bins[i-1][j].dig) FindCluster(iChamber,segmentation,i-1,j,bins,c);
+ if (bins[i][j-1].dig) FindCluster(iChamber,segmentation,i,j-1,bins,c);
+ if (bins[i+1][j].dig) FindCluster(iChamber,segmentation,i+1,j,bins,c);
+ if (bins[i][j+1].dig) FindCluster(iChamber,segmentation,i,j+1,bins,c);
+
+}
+
+//_____________________________________________________________________________
+
+void AliRICHpoints::GetCenterOfGravity()
+{
+ //
+ // simple RICH cluster finder from digits -- finds neighbours and
+ // calculates center of gravity for the cluster
+ //
+ const Int_t MAX_Nipx=400, MAX_Nipy=800;
+ printf("\n Hallo world");
+ AliRICHdisplay *display=(AliRICHdisplay*)gAlice->Display();
+ Int_t chamber=display->GetChamber();
+ Int_t cathode=display->GetCathode();
+
+ AliRICH *RICH = (AliRICH*)gAlice->GetDetector("RICH");
+ AliRICHchamber *iChamber;
+ AliRICHsegmentation *segmentation;
+ iChamber =&(RICH->Chamber(chamber-1));
+ segmentation=iChamber->GetSegmentationModel(cathode);
+ Int_t npx = segmentation->Npx();
+ Int_t npy = segmentation->Npy();
+ Float_t zpos=iChamber->ZPosition();
+
+ TClonesArray *RICHdigits = RICH->DigitsAddress(chamber-1);
+ gAlice->TreeD()->GetEvent(cathode);
+ Int_t ndigits = RICHdigits->GetEntriesFast();
+ if (fDigitIndex < 0 || fDigitIndex >= ndigits) return;
+
+ AliRICHdigit *dig;
+ dig=(AliRICHdigit*)RICHdigits->UncheckedAt(fDigitIndex);
+ Int_t ipx=dig->fPadX;
+ Int_t ipy=dig->fPadY;
+ Bin bins[MAX_Nipx][MAX_Nipy];
+ bins[ipx+npx][ipy+npy].dig=dig;
+
+ int ndig;
+ int ncls=0;
+ for (ndig=0; ndig<ndigits; ndig++) {
+ dig = (AliRICHdigit*)RICHdigits->UncheckedAt(ndig);
+ int i=dig->fPadX, j=dig->fPadY;
+ bins[i+npx][j+npy].dig=dig;
+ }
+
+ PreCluster c; c.summit=bins[ipx+npx][ipy+npy].dig; c.idx=ncls;
+ FindCluster(iChamber,segmentation,ipx+npx, ipy+npy, bins, c);
+ if (c.npeaks>1) {
+ printf("GetCenterOfGravity -- more than one peak");
+ }
+ c.fX /= c.fQ;
+ c.fY /= c.fQ;
+ printf("GetCenterOfGravity - c.fX c.fY c.fQ %f %f %d \n",c.fX,c.fY,c.fQ);
+
+ c.fTracks[0]=c.summit->fTracks[0];
+ c.fTracks[1]=c.summit->fTracks[1];
+ c.fTracks[2]=c.summit->fTracks[2];
+ ncls++;
+ AliRICHpoints *points = 0;
+ points = new AliRICHpoints(1);
+ points->SetMarkerColor(kYellow);
+ points->SetMarkerStyle(5);
+ points->SetMarkerSize(1.);
+ points->SetPoint(0,c.fX,c.fY,zpos);
+ points->Draw();
+
+ printf("GetCenterOfGravity -- ncls %d \n",ncls);
+
+}
+
+
+
--- /dev/null
+#ifndef AliRICHpoints_H
+#define AliRICHpoints_H
+
+#include "TPolyMarker3D.h"
+#include "AliRICH.h"
+#include "AliPoints.h"
+
+class AliRICHpoints : public AliPoints {
+protected:
+ Int_t fHitIndex; // Link to hit number
+ Int_t fTrackIndex; // Link to track number
+ Int_t fDigitIndex; // Link to digit
+
+public:
+ AliRICHpoints();
+ AliRICHpoints(Int_t npoints);
+ virtual ~AliRICHpoints();
+
+ Int_t GetHitIndex() {return fHitIndex;}
+ Int_t GetTrackIndex(); // *MENU*
+ Int_t GetDigitIndex() {return fDigitIndex;}
+ AliRICHhit *GetHit() const;
+ AliRICHdigit *GetDigit() const;
+ virtual void InspectHit(); // *MENU*
+ virtual void DumpHit(); // *MENU*
+ virtual void InspectDigit(); // *MENU*
+ virtual void DumpDigit(); // *MENU*
+ virtual void GetCenterOfGravity(); // *MENU*
+ virtual void SetHitIndex(Int_t hitindex) {fHitIndex = hitindex;}
+ virtual void SetTrackIndex(Int_t trackindex) {fTrackIndex = trackindex;}
+ virtual void SetDigitIndex(Int_t digitindex) {fDigitIndex = digitindex;}
+
+ ClassDef(AliRICHpoints,1) //Class to draw detector clusters (is PolyMarker3D)
+};
+#endif
--- /dev/null
+/////////////////////////////////////////////////////////
+// Manager and hits classes for set:RICH version 0 //
+/////////////////////////////////////////////////////////
+
+#include <TTUBE.h>
+#include <TNode.h>
+#include <TRandom.h>
+
+#include "AliRICHv0.h"
+#include "AliRun.h"
+#include "AliMC.h"
+#include "iostream.h"
+#include "AliCallf77.h"
+#include "AliConst.h"
+#include "TGeant3.h"
+
+ClassImp(AliRICHv0)
+
+//___________________________________________
+AliRICHv0::AliRICHv0() : AliRICH()
+{
+ fChambers = 0;
+}
+
+//___________________________________________
+AliRICHv0::AliRICHv0(const char *name, const char *title)
+ : AliRICH(name,title)
+{
+
+ fChambers = new TObjArray(7);
+ for (Int_t i=0; i<7; i++) {
+
+ (*fChambers)[i] = new AliRICHchamber();
+
+ }
+}
+
+
+//___________________________________________
+void AliRICHv0::CreateGeometry()
+{
+ //
+ // Create the geometry for RICH version 1
+ //
+ // Modified by: N. Colonna (INFN - BARI, Nicola.Colonna@ba.infn.it)
+ // R.A. Fini (INFN - BARI, Rosanna.Fini@ba.infn.it)
+ // R.A. Loconsole (Bari University, loco@riscom.ba.infn.it)
+ //
+ //Begin_Html
+ /*
+ <img src="picts/AliRICHv1.gif">
+ */
+ //End_Html
+ //Begin_Html
+ /*
+ <img src="picts/AliRICHv1Tree.gif">
+ */
+ //End_Html
+
+
+ Int_t *idtmed = fIdtmed->GetArray()-999;
+
+ Int_t i;
+ Float_t zs;
+ Int_t idrotm[1099];
+ Float_t par[3];
+
+ // --- Define the RICH detector
+ // External aluminium box
+ par[0] = 71.1;
+ par[1] = 11.5;
+ par[2] = 73.15;
+ gMC->Gsvolu("RICH", "BOX ", idtmed[1009], par, 3);
+
+ // Sensitive part of the whole RICH
+ par[0] = 64.8;
+ par[1] = 11.5;
+ par[2] = 66.55;
+ gMC->Gsvolu("SRIC", "BOX ", idtmed[1000], par, 3);
+
+ // Honeycomb
+ par[0] = 63.1;
+ par[1] = .188;
+ par[2] = 66.55;
+ gMC->Gsvolu("HONE", "BOX ", idtmed[1001], par, 3);
+
+ // Aluminium sheet
+ par[0] = 63.1;
+ par[1] = .025;
+ par[2] = 66.55;
+ gMC->Gsvolu("ALUM", "BOX ", idtmed[1009], par, 3);
+
+ // Quartz
+ par[0] = 63.1;
+ par[1] = .25;
+ par[2] = 65.5;
+ gMC->Gsvolu("QUAR", "BOX ", idtmed[1002], par, 3);
+
+ // Spacers (cylinders)
+ par[0] = 0.;
+ par[1] = .5;
+ par[2] = .5;
+ gMC->Gsvolu("SPAC", "TUBE", idtmed[1002], par, 3);
+
+ // Opaque quartz
+ par[0] = 61.95;
+ par[1] = .2;
+ par[2] = 66.5;
+ gMC->Gsvolu("OQUA", "BOX ", idtmed[1007], par, 3);
+
+ // Frame of opaque quartz
+ par[0] = 20.65;
+ par[1] = .5;
+ par[2] = 66.5;
+ gMC->Gsvolu("OQUF", "BOX ", idtmed[1007], par, 3);
+
+ // Little bar of opaque quartz
+ par[0] = 63.1;
+ par[1] = .25;
+ par[2] = .275;
+ gMC->Gsvolu("BARR", "BOX ", idtmed[1007], par, 3);
+
+ // Freon
+ par[0] = 20.15;
+ par[1] = .5;
+ par[2] = 65.5;
+ gMC->Gsvolu("FREO", "BOX ", idtmed[1003], par, 3);
+
+ // Methane
+ par[0] = 64.8;
+ par[1] = 5.;
+ par[2] = 64.8;
+ gMC->Gsvolu("META", "BOX ", idtmed[1004], par, 3);
+
+ // Methane gap
+ par[0] = 64.8;
+ par[1] = .2;
+ par[2] = 64.8;
+ gMC->Gsvolu("GAP ", "BOX ", idtmed[1008], par, 3);
+
+ // CsI photocathode
+ par[0] = 64.8;
+ par[1] = .25;
+ par[2] = 64.8;
+ gMC->Gsvolu("CSI ", "BOX ", idtmed[1005], par, 3);
+
+ // Anode grid
+ par[0] = 0.;
+ par[1] = .0025;
+ par[2] = 20.;
+ gMC->Gsvolu("GRID", "TUBE", idtmed[1006], par, 3);
+
+ // --- Places the detectors defined with GSVOLU
+ // Place material inside RICH
+ gMC->Gspos("SRIC", 1, "RICH", 0., 0., 0., 0, "ONLY");
+
+ gMC->Gspos("ALUM", 1, "SRIC", 0., -6.075, 0., 0, "ONLY");
+ gMC->Gspos("HONE", 1, "SRIC", 0., -5.862, 0., 0, "ONLY");
+ gMC->Gspos("ALUM", 2, "SRIC", 0., -5.649, 0., 0, "ONLY");
+ gMC->Gspos("OQUA", 1, "SRIC", 0., -5.424, 0., 0, "ONLY");
+
+ AliMatrix(idrotm[1019], 0., 0., 90., 0., 90., 90.);
+
+ for (i = 1; i <= 9; ++i) {
+ zs = (5 - i) * 14.4;
+ gMC->Gspos("SPAC", i, "FREO", 6.7, 0., zs, idrotm[1019], "ONLY");
+ }
+ for (i = 10; i <= 18; ++i) {
+ zs = (14 - i) * 14.4;
+ gMC->Gspos("SPAC", i, "FREO", -6.7, 0., zs, idrotm[1019], "ONLY");
+ }
+
+ gMC->Gspos("FREO", 1, "OQUF", 0., 0., 0., 0, "ONLY");
+ gMC->Gspos("OQUF", 1, "SRIC", 41.3, -4.724, 0., 0, "ONLY");
+ gMC->Gspos("OQUF", 2, "SRIC", 0., -4.724, 0., 0, "ONLY");
+ gMC->Gspos("OQUF", 3, "SRIC", -41.3, -4.724, 0., 0, "ONLY");
+ gMC->Gspos("BARR", 1, "QUAR", 0., 0., -21.65, 0, "ONLY");
+ gMC->Gspos("BARR", 2, "QUAR", 0., 0., 21.65, 0, "ONLY");
+ gMC->Gspos("QUAR", 1, "SRIC", 0., -3.974, 0., 0, "ONLY");
+ gMC->Gspos("GAP ", 1, "META", 0., 4.8, 0., 0, "ONLY");
+ gMC->Gspos("META", 1, "SRIC", 0., 1.276, 0., 0, "ONLY");
+ gMC->Gspos("CSI ", 1, "SRIC", 0., 6.526, 0., 0, "ONLY");
+
+ // Place RICH inside ALICE apparatus
+
+ AliMatrix(idrotm[1000], 90., 0., 70.69, 90., 19.31, -90.);
+ AliMatrix(idrotm[1001], 90., -20., 90., 70., 0., 0.);
+ AliMatrix(idrotm[1002], 90., 0., 90., 90., 0., 0.);
+ AliMatrix(idrotm[1003], 90., 20., 90., 110., 0., 0.);
+ AliMatrix(idrotm[1004], 90., 340., 108.2, 70., 18.2, 70.);
+ AliMatrix(idrotm[1005], 90., 0., 109.31, 90., 19.31, 90.);
+ AliMatrix(idrotm[1006], 90., 20., 108.2, 110., 18.2, 110.);
+
+ gMC->Gspos("RICH", 1, "ALIC", 0., 471.9, 165.26, idrotm[1000], "ONLY");
+ gMC->Gspos("RICH", 2, "ALIC", 171., 470., 0., idrotm[1001], "ONLY");
+ gMC->Gspos("RICH", 3, "ALIC", 0., 500., 0., idrotm[1002], "ONLY");
+ gMC->Gspos("RICH", 4, "ALIC", -171., 470., 0., idrotm[1003], "ONLY");
+ gMC->Gspos("RICH", 5, "ALIC", 161.4, 443.4, -165.3, idrotm[1004], "ONLY");
+ gMC->Gspos("RICH", 6, "ALIC", 0., 471.9, -165.3, idrotm[1005], "ONLY");
+ gMC->Gspos("RICH", 7, "ALIC", -161.4, 443.4, -165.3, idrotm[1006], "ONLY");
+
+}
+
+
+//___________________________________________
+void AliRICHv0::CreateMaterials()
+{
+ //
+ // *** DEFINITION OF AVAILABLE RICH MATERIALS ***
+ // ORIGIN : NICK VAN EIJNDHOVEN
+ // Modified by: N. Colonna (INFN - BARI, Nicola.Colonna@ba.infn.it)
+ // R.A. Fini (INFN - BARI, Rosanna.Fini@ba.infn.it)
+ // R.A. Loconsole (Bari University, loco@riscom.ba.infn.it)
+ //
+ Int_t ISXFLD = gAlice->Field()->Integ();
+ Float_t SXMGMX = gAlice->Field()->Max();
+
+ Float_t ppckov[14] = { 5.63e-9,5.77e-9,5.9e-9,6.05e-9,6.2e-9,6.36e-9,6.52e-9,
+ 6.7e-9,6.88e-9,7.08e-9,7.3e-9,7.51e-9,7.74e-9,8e-9 };
+ Float_t rindex_quarz[14] = { 1.528309,1.533333,
+ 1.538243,1.544223,1.550568,1.55777,
+ 1.565463,1.574765,1.584831,1.597027,
+ 1.611858,1.6277,1.6472,1.6724 };
+ Float_t rindex_quarzo[14] = { 1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1. };
+ Float_t rindex_methane[14] = { 1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1. };
+ Float_t rindex_gri[14] = { 1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1. };
+ Float_t absco_freon[14] = { 179.0987,179.0987,
+ 179.0987,179.0987,179.0987,35.7,12.54,5.92,4.92,3.86,1.42,.336,.134,0. };
+ Float_t absco_quarz[14] = { 20.126,16.27,13.49,11.728,9.224,8.38,7.44,7.17,
+ 6.324,4.483,1.6,.323,.073,0. };
+ Float_t absco_quarzo[14] = { 1e-5,1e-5,1e-5,1e-5,1e-5,1e-5,1e-5,1e-5,1e-5,
+ 1e-5,1e-5,1e-5,1e-5,1e-5 };
+ Float_t absco_csi[14] = { 1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,
+ 1e-4,1e-4,1e-4,1e-4 };
+ Float_t absco_methane[14] = { 1e6,1e6,1e6,1e6,1e6,1e6,1e6,1e6,1e6,1e6,1e6,
+ 1e6,1e6,1e6 };
+ Float_t absco_gri[14] = { 1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,
+ 1e-4,1e-4,1e-4,1e-4 };
+ Float_t effic_all[14] = { 1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1. };
+ Float_t effic_csi[14] = { 4.74e-4,.00438,.009,.0182,.0282,.0653,.1141,.163,
+ .2101,.2554,.293,.376,.3861,.418 };
+ Float_t effic_gri[14] = { 1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1. };
+
+ Float_t afre[2], agri, amet[2], aqua[2], ahon, zfre[2], zgri, zhon,
+ zmet[2], zqua[2];
+ Int_t nlmatfre;
+ Float_t densquao;
+ Int_t nlmatmet, nlmatqua;
+ Float_t wmatquao[2], rindex_freon[14];
+ Int_t i;
+ Float_t aquao[2], epsil, stmin, zquao[2];
+ Int_t nlmatquao;
+ Float_t radlal, densal, tmaxfd, deemax, stemax;
+ Float_t aal, zal, radlgri, densfre, radlhon, densgri, denshon,densqua, densmet, wmatfre[2], wmatmet[2], wmatqua[2];
+
+ Int_t *idtmed = fIdtmed->GetArray()-999;
+
+ TGeant3 *geant3 = (TGeant3*) gMC;
+
+ // --- Photon energy (GeV)
+ // --- Refraction indexes
+ for (i = 0; i < 14; ++i) {
+ rindex_freon[i] = ppckov[i] * .01095 * 1e9 + 1.2177;
+ }
+ // need to be changed
+
+ // --- Absorbtion lenghts (in cm)
+ // DATA ABSCO_QUARZ /
+ // & 5 * 1000000.,
+ // & 29.85, 7.34, 4.134, 1.273, 0.722,
+ // & 0.365, 0.365, 0.365, 0. /
+ // need to be changed
+
+ // --- Detection efficiencies (quantum efficiency for CsI)
+ // --- Define parameters for honeycomb.
+ // Used carbon of equivalent rad. lenght
+
+ ahon = 12.01;
+ zhon = 6.;
+ denshon = 2.265;
+ radlhon = 18.8;
+
+ // --- Parameters to include in GSMIXT, relative to Quarz (SiO2)
+
+ aqua[0] = 28.09;
+ aqua[1] = 16.;
+ zqua[0] = 14.;
+ zqua[1] = 8.;
+ densqua = 2.64;
+ nlmatqua = -2;
+ wmatqua[0] = 1.;
+ wmatqua[1] = 2.;
+
+ // --- Parameters to include in GSMIXT, relative to opaque Quarz (SiO2)
+
+ aquao[0] = 28.09;
+ aquao[1] = 16.;
+ zquao[0] = 14.;
+ zquao[1] = 8.;
+ densquao = 2.64;
+ nlmatquao = -2;
+ wmatquao[0] = 1.;
+ wmatquao[1] = 2.;
+
+ // --- Parameters to include in GSMIXT, relative to Freon (C6F14)
+
+ afre[0] = 12.;
+ afre[1] = 19.;
+ zfre[0] = 6.;
+ zfre[1] = 9.;
+ densfre = 1.7;
+ nlmatfre = -2;
+ wmatfre[0] = 6.;
+ wmatfre[1] = 14.;
+
+ // --- Parameters to include in GSMIXT, relative to methane (CH4)
+
+ amet[0] = 12.01;
+ amet[1] = 1.;
+ zmet[0] = 6.;
+ zmet[1] = 1.;
+ densmet = 7.17e-4;
+ nlmatmet = -2;
+ wmatmet[0] = 1.;
+ wmatmet[1] = 4.;
+
+ // --- Parameters to include in GSMIXT, relative to anode grid (Cu)
+
+ agri = 63.54;
+ zgri = 29.;
+ densgri = 8.96;
+ radlgri = 1.43;
+
+ // --- Parameters to include in GSMATE related to aluminium sheet
+
+ aal = 26.98;
+ zal = 13.;
+ densal = 2.7;
+ radlal = 8.9;
+
+ AliMaterial(1, "Air $", 14.61, 7.3, .001205, 30420., 67500);
+ AliMaterial(6, "HON", ahon, zhon, denshon, radlhon, 0);
+ AliMaterial(16, "CSI", ahon, zhon, denshon, radlhon, 0);
+ AliMixture(20, "QUA", aqua, zqua, densqua, nlmatqua, wmatqua);
+ AliMixture(21, "QUAO", aquao, zquao, densquao, nlmatquao, wmatquao);
+ AliMixture(30, "FRE", afre, zfre, densfre, nlmatfre, wmatfre);
+ AliMixture(40, "MET", amet, zmet, densmet, nlmatmet, wmatmet);
+ AliMixture(41, "METG", amet, zmet, densmet, nlmatmet, wmatmet);
+ AliMaterial(11, "GRI", agri, zgri, densgri, radlgri, 0);
+ AliMaterial(50, "ALUM", aal, zal, densal, radlal, 0);
+
+ tmaxfd = -10.;
+ stemax = -.1;
+ deemax = -.2;
+ epsil = .001;
+ stmin = -.001;
+
+ AliMedium(1, "DEFAULT MEDIUM AIR$", 1, 0, ISXFLD, SXMGMX, tmaxfd, stemax, deemax, epsil, stmin);
+ AliMedium(2, "HONEYCOMB$", 6, 0, ISXFLD, SXMGMX, tmaxfd, stemax, deemax, epsil, stmin);
+ AliMedium(3, "QUARZO$", 20, 1, ISXFLD, SXMGMX, tmaxfd, stemax, deemax, epsil, stmin);
+ AliMedium(4, "FREON$", 30, 1, ISXFLD, SXMGMX, tmaxfd, stemax, deemax, epsil, stmin);
+ AliMedium(5, "METANO$", 40, 1, ISXFLD, SXMGMX, tmaxfd, stemax, deemax, epsil, stmin);
+ AliMedium(6, "CSI$", 16, 1, ISXFLD, SXMGMX,tmaxfd, stemax, deemax, epsil, stmin);
+ AliMedium(7, "GRIGLIA$", 11, 0, ISXFLD, SXMGMX, tmaxfd, stemax, deemax, epsil, stmin);
+ AliMedium(8, "QUARZOO$", 21, 1, ISXFLD, SXMGMX, tmaxfd, stemax, deemax, epsil, stmin);
+ AliMedium(9, "GAP$", 41, 1, ISXFLD, SXMGMX,tmaxfd, .1, -deemax, epsil, -stmin);
+ AliMedium(10, "ALUMINUM$", 50, 1, ISXFLD, SXMGMX, tmaxfd, stemax, deemax, epsil, stmin);
+
+
+ // Switch on delta-ray production in the methane and freon gaps
+
+ gMC->Gstpar(idtmed[1002], "LOSS", 1.);
+ gMC->Gstpar(idtmed[1003], "LOSS", 1.);
+ gMC->Gstpar(idtmed[1004], "LOSS", 1.);
+ gMC->Gstpar(idtmed[1008], "LOSS", 1.);
+ gMC->Gstpar(idtmed[1005], "LOSS", 1.);
+ gMC->Gstpar(idtmed[1002], "HADR", 1.);
+ gMC->Gstpar(idtmed[1003], "HADR", 1.);
+ gMC->Gstpar(idtmed[1004], "HADR", 1.);
+ gMC->Gstpar(idtmed[1008], "HADR", 1.);
+ gMC->Gstpar(idtmed[1005], "HADR", 1.);
+ gMC->Gstpar(idtmed[1002], "DCAY", 1.);
+ gMC->Gstpar(idtmed[1003], "DCAY", 1.);
+ gMC->Gstpar(idtmed[1004], "DCAY", 1.);
+ gMC->Gstpar(idtmed[1008], "DCAY", 1.);
+ gMC->Gstpar(idtmed[1005], "DCAY", 1.);
+ geant3->Gsckov(idtmed[1000], 14, ppckov, absco_methane, effic_all, rindex_methane);
+ geant3->Gsckov(idtmed[1001], 14, ppckov, absco_methane, effic_all, rindex_methane);
+ geant3->Gsckov(idtmed[1002], 14, ppckov, absco_quarz, effic_all,rindex_quarz);
+ geant3->Gsckov(idtmed[1003], 14, ppckov, absco_freon, effic_all,rindex_freon);
+ geant3->Gsckov(idtmed[1004], 14, ppckov, absco_methane, effic_all, rindex_methane);
+ geant3->Gsckov(idtmed[1005], 14, ppckov, absco_csi, effic_csi, rindex_methane);
+ geant3->Gsckov(idtmed[1006], 14, ppckov, absco_gri, effic_gri, rindex_gri);
+ geant3->Gsckov(idtmed[1007], 14, ppckov, absco_quarzo, effic_all, rindex_quarzo);
+ geant3->Gsckov(idtmed[1008], 14, ppckov, absco_methane, effic_all, rindex_methane);
+ geant3->Gsckov(idtmed[1009], 14, ppckov, absco_gri, effic_gri, rindex_gri);
+}
+
+//___________________________________________
+
+void AliRICHv0::Init()
+{
+ printf("\n\n\n Start Init for version 0 - CPC chamber type \n\n\n");
+
+ //
+ // Initialize Tracking Chambers
+ //
+ for (Int_t i=0; i<7; i++) {
+ ( (AliRICHchamber*) (*fChambers)[i])->Init();
+ }
+
+ //
+ // Set the chamber (sensitive region) GEANT identifier
+
+ ((AliRICHchamber*)(*fChambers)[0])->SetGid(1);
+ ((AliRICHchamber*)(*fChambers)[1])->SetGid(2);
+ ((AliRICHchamber*)(*fChambers)[2])->SetGid(3);
+ ((AliRICHchamber*)(*fChambers)[3])->SetGid(4);
+ ((AliRICHchamber*)(*fChambers)[4])->SetGid(5);
+ ((AliRICHchamber*)(*fChambers)[5])->SetGid(6);
+ ((AliRICHchamber*)(*fChambers)[6])->SetGid(7);
+
+ printf("\n\n\n Finished Init for version 0 - CPC chamber type\n\n\n");
+}
+
+//___________________________________________
+void AliRICHv0::StepManager()
+{
+ Int_t copy, id;
+ static Int_t idvol;
+ static Int_t vol[2];
+ Int_t ipart;
+ static Float_t hits[9];
+ TLorentzVector Position;
+ TLorentzVector Momentum;
+ Float_t pos[3];
+ Float_t mom[4];
+ Float_t Localpos[3];
+ Float_t Localmom[4];
+ Float_t Localtheta,Localphi;
+ Float_t theta,phi;
+ Float_t destep, step;
+ static Float_t eloss, xhit, yhit, tlength;
+ const Float_t big=1.e10;
+
+ TClonesArray &lhits = *fHits;
+ TClonesArray &lcerenkovs = *fCerenkovs;
+
+ // Only gas gap inside chamber
+ // Tag chambers and record hits when track enters
+
+ idvol=-1;
+ id=gMC->CurrentVolID(copy);
+ Float_t cherenkov_loss=0.00001;
+
+ // Treat photons produced in Freon and Quartz
+ if (gMC->TrackPid() == 50000050 ) {
+ if (gMC->IsTrackEntering()){
+ if (gMC->VolId("FREO")==gMC->CurrentVolID(copy) || gMC->VolId("QUAR")==gMC->CurrentVolID(copy)){
+ //printf("GOT ONE! Type:%d \n",gMC->TrackPid());
+ }
+ }
+ }
+
+
+ if (gMC->TrackPid() == 50000050 ) {
+ if (gMC->VolId("CSI ")==gMC->CurrentVolID(copy))
+ {
+ if (gMC->Edep() > 0.){
+ gMC->TrackPosition(Position);
+ gMC->TrackMomentum(Momentum);
+ pos[0]=Position(0);
+ pos[1]=Position(1);
+ pos[2]=Position(2);
+ mom[0]=Momentum(0);
+ mom[1]=Momentum(1);
+ mom[2]=Momentum(2);
+ mom[3]=Momentum(3);
+ Double_t tc = mom[0]*mom[0]+mom[1]*mom[1];
+ Double_t rt = TMath::Sqrt(tc);
+ theta = Float_t(TMath::ATan2(rt,Double_t(mom[2])))*kRaddeg;
+ phi = Float_t(TMath::ATan2(Double_t(mom[1]),Double_t(mom[0])))*kRaddeg;
+ gMC->Gmtod(pos,Localpos,1);
+ gMC->Gmtod(mom,Localmom,2);
+
+ gMC->CurrentVolOffID(2,copy);
+ vol[0]=copy;
+ idvol=vol[0]-1;
+
+ ((AliRICHchamber*) (*fChambers)[idvol])
+ ->SigGenInit(Localpos[0], Localpos[2], Localpos[1]);
+ if(idvol<7) {
+ hits[0] = 50000050; // particle type
+ hits[1] = pos[0]; // X-position for hit
+ hits[2] = pos[1]; // Y-position for hit
+ hits[3] = pos[2]; // Z-position for hit
+ hits[4] = theta; // theta angle of incidence
+ hits[5] = phi; // phi angle of incidence
+ hits[8] = (Float_t) fNclusters; // first padhit
+ hits[9] = -1; // last pad hit
+
+ MakePadHits(Localpos[0],Localpos[2],cherenkov_loss,idvol,cerenkov);
+ if (fNclusters > (Int_t)hits[8]) {
+ hits[8]= hits[8]+1;
+ hits[9]= (Float_t) fNclusters;
+ }
+
+ AddHit(gAlice->CurrentTrack(),vol,hits);
+ new(lcerenkovs[fNcerenkovs++]) AliRICHCerenkov(fIshunt,gAlice->CurrentTrack(),vol,hits);
+
+ }
+ }
+ }
+
+
+//
+// treat charged particles
+ } else if (gMC->TrackCharge())
+ {
+
+//If MIP
+ if (gMC->VolId("GAP ")== gMC->CurrentVolID(copy)) {
+// Get current particle id (ipart), track position (pos) and momentum (mom)
+
+ gMC->CurrentVolOffID(3,copy);
+ vol[0]=copy;
+ idvol=vol[0]-1;
+
+ gMC->TrackPosition(Position);
+ gMC->TrackMomentum(Momentum);
+ pos[0]=Position(0);
+ pos[1]=Position(1);
+ pos[2]=Position(2);
+ mom[0]=Momentum(0);
+ mom[1]=Momentum(1);
+ mom[2]=Momentum(2);
+ mom[3]=Momentum(3);
+ gMC->Gmtod(pos,Localpos,1);
+ gMC->Gmtod(mom,Localmom,2);
+
+ ipart = gMC->TrackPid();
+ //
+ // momentum loss and steplength in last step
+ destep = gMC->Edep();
+ step = gMC->TrackStep();
+
+ //
+ // record hits when track enters ...
+ if( gMC->IsTrackEntering()) {
+ gMC->SetMaxStep(fMaxStepGas);
+ Double_t tc = mom[0]*mom[0]+mom[1]*mom[1];
+ Double_t rt = TMath::Sqrt(tc);
+ theta = Float_t(TMath::ATan2(rt,Double_t(mom[2])))*kRaddeg;
+ phi = Float_t(TMath::ATan2(Double_t(mom[1]),Double_t(mom[0])))*kRaddeg;
+
+ Localtheta = Float_t(TMath::ATan2(rt,Double_t(Localmom[2])))*kRaddeg;
+ Localphi = Float_t(TMath::ATan2(Double_t(Localmom[1]),Double_t(Localmom[0])))*kRaddeg;
+
+ hits[0] = Float_t(ipart); // particle type
+ hits[1] = pos[0]; // X-position for hit
+ hits[2] = pos[1]; // Y-position for hit
+ hits[3] = pos[2]; // Z-position for hit
+ hits[4] = theta; // theta angle of incidence
+ hits[5] = phi; // phi angle of incidence
+ hits[8] = (Float_t) fNclusters; // first padhit
+ hits[9] = -1; // last pad hit
+ // phi angle of incidence
+ tlength = 0;
+ eloss = 0;
+
+ Chamber(idvol).LocaltoGlobal(Localpos,hits+1);
+
+ //To make chamber coordinates x-y had to pass LocalPos[0], LocalPos[2]
+ xhit = Localpos[0];
+ yhit = Localpos[2];
+ // Only if not trigger chamber
+ if(idvol<7) {
+ //
+ // Initialize hit position (cursor) in the segmentation model
+ ((AliRICHchamber*) (*fChambers)[idvol])
+ ->SigGenInit(Localpos[0], Localpos[2], Localpos[1]);
+ }
+ }
+
+ //
+ // Calculate the charge induced on a pad (disintegration) in case
+ //
+ // Mip left chamber ...
+ if( gMC->IsTrackExiting() || gMC->IsTrackStop() || gMC->IsTrackDisappeared()){
+ gMC->SetMaxStep(big);
+ eloss += destep;
+ tlength += step;
+
+
+
+ // Only if not trigger chamber
+ if(idvol<7) {
+ if (eloss > 0) MakePadHits(xhit,yhit,eloss,idvol,mip);
+ }
+
+ hits[6]=tlength;
+ hits[7]=eloss;
+ if (fNclusters > (Int_t)hits[8]) {
+ hits[8]= hits[8]+1;
+ hits[9]= (Float_t) fNclusters;
+ }
+
+ new(lhits[fNhits++])
+ AliRICHhit(fIshunt,gAlice->CurrentTrack(),vol,hits);
+ eloss = 0;
+ //
+ // Check additional signal generation conditions
+ // defined by the segmentation
+ // model (boundary crossing conditions)
+ } else if
+ (((AliRICHchamber*) (*fChambers)[idvol])
+ ->SigGenCond(Localpos[0], Localpos[2], Localpos[1]))
+ {
+ ((AliRICHchamber*) (*fChambers)[idvol])
+ ->SigGenInit(Localpos[0], Localpos[2], Localpos[1]);
+ if (eloss > 0) MakePadHits(xhit,yhit,eloss,idvol,mip);
+ xhit = Localpos[0];
+ yhit = Localpos[2];
+ eloss = destep;
+ tlength += step ;
+ //
+ // nothing special happened, add up energy loss
+ } else {
+ eloss += destep;
+ tlength += step ;
+ }
+ }
+ }
+}
+
+
+//___________________________________________
+void AliRICH::MakePadHits(Float_t xhit,Float_t yhit,Float_t eloss, Int_t idvol, Response_t res)
+{
+//
+// Calls the charge disintegration method of the current chamber and adds
+// the simulated cluster to the root treee
+//
+ Int_t clhits[7];
+ Float_t newclust[6][500];
+ Int_t nnew;
+
+//
+// Integrated pulse height on chamber
+
+ clhits[0]=fNhits+1;
+
+ ((AliRICHchamber*) (*fChambers)[idvol])->DisIntegration(eloss, xhit, yhit, nnew, newclust, res);
+ Int_t ic=0;
+
+//
+// Add new clusters
+ for (Int_t i=0; i<nnew; i++) {
+ if (Int_t(newclust[3][i]) > 0) {
+ ic++;
+// Cathode plane
+ clhits[1] = Int_t(newclust[5][i]);
+// Cluster Charge
+ clhits[2] = Int_t(newclust[0][i]);
+// Pad: ix
+ clhits[3] = Int_t(newclust[1][i]);
+// Pad: iy
+ clhits[4] = Int_t(newclust[2][i]);
+// Pad: charge
+ clhits[5] = Int_t(newclust[3][i]);
+// Pad: chamber sector
+ clhits[6] = Int_t(newclust[4][i]);
+
+ AddCluster(clhits);
+ }
+ }
+}
+
+ClassImp(AliRICHchamber)
+
+ AliRICHchamber::AliRICHchamber()
+{
+ fSegmentation = new TObjArray(2);
+ fResponse= new TObjArray(2);
+ fnsec=1;
+ frMin=0.1;
+ frMax=140;
+}
+
+//
+// Get reference to response model
+AliRICHresponse* AliRICHchamber::GetResponseModel(Response_t res)
+{
+ if (res==mip) {
+ return (AliRICHresponse*) (*fResponse)[0];
+ } else if (res==cerenkov) {
+ return (AliRICHresponse*) (*fResponse)[1];
+ }
+ return (AliRICHresponse*) (*fResponse)[0];
+}
+
+// Configure response model
+void AliRICHchamber::ResponseModel(Response_t res, AliRICHresponse* thisResponse)
+{
+
+ if (res==mip) {
+ (*fResponse)[0]=thisResponse;
+ } else if (res==cerenkov) {
+ (*fResponse)[1]=thisResponse;
+ }
+}
+
+void AliRICHchamber::Init()
+{
+
+ ((AliRICHsegmentation *) (*fSegmentation)[0])->Init(this);
+ if (fnsec==2) {
+ ((AliRICHsegmentation *) (*fSegmentation)[1])->Init(this);
+ }
+}
+
+void AliRICHchamber::LocaltoGlobal(Float_t pos[3],Float_t Localpos[3])
+{
+
+ Double_t *fMatrix;
+ fMatrix = fChamberMatrix->GetMatrix();
+ Localpos[0]=pos[0]*fMatrix[0]+pos[1]*fMatrix[3]+pos[2]*fMatrix[6];
+ Localpos[1]=pos[0]*fMatrix[1]+pos[1]*fMatrix[4]+pos[2]*fMatrix[7];
+ Localpos[2]=pos[0]*fMatrix[2]+pos[1]*fMatrix[5]+pos[2]*fMatrix[8];
+ Localpos[0]+=fChamberTrans[0];
+ Localpos[1]+=fChamberTrans[1];
+ Localpos[2]+=fChamberTrans[2];
+}
+
+
+void AliRICHchamber::DisIntegration(Float_t eloss, Float_t xhit, Float_t yhit,
+ Int_t& nnew,Float_t newclust[6][500],Response_t res)
+{
+//
+// Generates pad hits (simulated cluster)
+// using the segmentation and the response model
+
+ Float_t dx, dy;
+ Float_t local[3];
+ Float_t source[3];
+ Int_t Nfp=0;
+
+ //
+ // Width of the integration area
+ //
+ dx=((AliRICHresponse*) (*fResponse)[0])->Nsigma()*((AliRICHresponse*) (*fResponse)[0])->ChwX();
+ dy=((AliRICHresponse*) (*fResponse)[0])->Nsigma()*((AliRICHresponse*) (*fResponse)[0])->ChwY();
+ //
+ // Get pulse height from energy loss
+ Float_t qtot;
+
+ if (res==mip) {
+ qtot = ((AliRICHresponse*) (*fResponse)[0])->IntPH(eloss);
+
+ local[0]=xhit;
+ //Z position of the wires relative to the RICH mother volume
+ local[1]=6.026;
+ local[2]=yhit;
+ //Generate feedback photons
+ Nfp = ((AliRICHresponse*) (*fResponse)[0])->FeedBackPhotons(source,qtot);
+ } else if (res==cerenkov) {
+ qtot = ((AliRICHresponse*) (*fResponse)[1])->IntPH();
+ local[0]=xhit;
+ local[1]=6.026;
+ local[2]=yhit;
+ Nfp = ((AliRICHresponse*) (*fResponse)[1])->FeedBackPhotons(source,qtot);
+ }
+
+ //
+ // Loop Over Pads
+
+ Float_t qcheck=0, qp;
+ nnew=0;
+ for (Int_t i=1; i<=fnsec; i++) {
+ qcheck=0;
+ AliRICHsegmentation * segmentation=(AliRICHsegmentation *) (*fSegmentation)[i-1];
+ for (segmentation->FirstPad(xhit, yhit, dx, dy);
+ segmentation->MorePads();
+ segmentation->NextPad())
+ {
+ if (res==mip) {
+ qp= ((AliRICHresponse*) (*fResponse)[0])->IntXY(segmentation);
+ }
+ if (res==cerenkov) {
+ qp= ((AliRICHresponse*) (*fResponse)[0])->IntXY(segmentation);
+ }
+
+ qp= TMath::Abs(qp);
+
+ if (qp > 1.e-4) {
+ qcheck+=qp;
+ //
+ // --- store signal information
+ newclust[0][nnew]=qtot;
+ newclust[1][nnew]=segmentation->Ix();
+ newclust[2][nnew]=segmentation->Iy();
+ newclust[3][nnew]=qp * qtot;
+ newclust[4][nnew]=segmentation->ISector();
+ newclust[5][nnew]=(Float_t) i;
+ nnew++;
+
+
+ }
+ } // Pad loop
+ } // Cathode plane loop
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
--- /dev/null
+#ifndef RICHv0_H
+#define RICHv0_H
+/////////////////////////////////////////////////////////
+// Manager and hits classes for set:RICH version 0 //
+/////////////////////////////////////////////////////////
+
+#include "AliRICH.h"
+
+class AliRICHv0 : public AliRICH {
+
+ public:
+ AliRICHv0();
+ AliRICHv0(const char *name, const char *title);
+ virtual ~AliRICHv0() {}
+ virtual void CreateGeometry();
+ virtual void CreateMaterials();
+ virtual void Init();
+ virtual Int_t IsVersion() const {return 0;}
+ virtual void StepManager();
+// virtual void Trigger(Float_t (*)[4], Float_t (*)[4], Int_t& iflag);
+ private:
+ ClassDef(AliRICHv0,1) //Hits manager for set:RICH version 0
+
+ };
+
+
+#endif
+
+
+
+
+
+
+
# C++ sources
-SRCS = AliRICH.cxx
+SRCS = AliRICH.cxx AliRICHv0.cxx AliRICHdisplay.cxx AliRICHpoints.cxx AliRICHHitMap.cxx AliRICHSegResV0.cxx AliRICHSegResCkv.cxx
# C++ Headers
#ifdef __CINT__
-
#pragma link off all globals;
#pragma link off all classes;
#pragma link off all functions;
-
+#pragma link C++ enum Response_t;
#pragma link C++ class AliRICH-;
-#pragma link C++ class AliRICHv1;
+#pragma link C++ class AliRICHv0;
#pragma link C++ class AliRICHhit;
-#pragma link C++ class AliRICHmip;
-#pragma link C++ class AliRICHckov;
-#pragma link C++ class AliRICHpadhit;
-
+#pragma link C++ class AliRICHCerenkov;
+#pragma link C++ class AliRICHdigit;
+#pragma link C++ class AliRICHlist;
+#pragma link C++ class AliRICHcluster;
+#pragma link C++ class AliRICHreccluster;
+#pragma link C++ class AliRICHRecCluster;
+#pragma link C++ class AliRICHsegmentation;
+#pragma link C++ class AliRICHresponse;
+#pragma link C++ class AliRICHsegmentationV0;
+#pragma link C++ class AliRICHresponseV0;
+#pragma link C++ class AliRICHresponseCkv;
+#pragma link C++ class AliRICHchamber;
+#pragma link C++ class AliRICHpoints;
+#pragma link C++ class AliRICHdisplay;
+#pragma link C++ class AliRICHHitMap;
+#pragma link C++ class AliRICHHitMapA1;
+#pragma link C++ class Bin-;
+#pragma link C++ class PreCluster-;
#endif
--- /dev/null
+#include "iostream.h"
+
+void RICHdigit (Int_t evNumber1=0,Int_t evNumber2=0)
+{
+/////////////////////////////////////////////////////////////////////////
+// This macro is a small example of a ROOT macro
+// illustrating how to read the output of GALICE
+// and do some analysis.
+//
+/////////////////////////////////////////////////////////////////////////
+
+// Dynamically link some shared libs
+
+ if (gClassTable->GetID("AliRun") < 0) {
+ gROOT->LoadMacro("loadlibs.C");
+ loadlibs();
+ }
+
+
+// Connect the Root Galice file containing Geometry, Kine and Hits
+
+ TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
+ if (file) file->Close();
+ file = new TFile("galice.root","UPDATE");
+ file->ls();
+
+ printf ("I'm after Map \n");
+
+// Get AliRun object from file or create it if not on file
+
+ if (!gAlice) {
+ gAlice = (AliRun*)file->Get("gAlice");
+ if (gAlice) printf("AliRun object found on file\n");
+ if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
+ }
+ printf ("I'm after gAlice \n");
+
+ AliRICH *RICH = (AliRICH*) gAlice->GetDetector("RICH");
+//
+// Event Loop
+//
+ for (int nev=0; nev<= evNumber2; nev++) {
+ Int_t nparticles = gAlice->GetEvent(nev);
+ cout << "nev " <<nev<<endl;
+ cout << "nparticles " <<nparticles<<endl;
+ if (nev < evNumber1) continue;
+ if (nparticles <= 0) return;
+ if (RICH) RICH->Digitise(nev);
+ char hname[30];
+ sprintf(hname,"TreeD%d",nev);
+ gAlice->TreeD()->Write(hname);
+ gAlice->TreeD()->Reset();
+ file->ls();
+ } // event loop
+ file->Close();
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
--- /dev/null
+RICHdisplay (Int_t nevent=0) {
+// Dynamically link some shared libs
+
+ if (gClassTable->GetID("AliRun") < 0) {
+ gROOT->LoadMacro("loadlibs.C");
+ loadlibs();
+ }
+
+
+// Connect the Root Galice file containing Geometry, Kine and Hits
+
+ TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
+ if (file) file->Close();
+ file = new TFile("galice.root","UPDATE");
+ file->ls();
+
+ printf ("I'm after Map \n");
+
+// Get AliRun object from file or create it if not on file
+
+ if (!gAlice) {
+ gAlice = (AliRun*)file->Get("gAlice");
+ if (gAlice) printf("AliRun object found on file\n");
+ if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
+ }
+ printf ("I'm after gAlice \n");
+
+ // Create Event Display object
+ AliRICHdisplay *muondisplay = new AliRICHdisplay(750);
+
+// Display first event
+ gAlice->GetEvent(nevent);
+ muondisplay->ShowNextEvent(0);
+}
--- /dev/null
+void RICHpadtest (Int_t evNumber1=0,Int_t evNumber2=0)
+{
+/////////////////////////////////////////////////////////////////////////
+// This macro is a small example of a ROOT macro
+// illustrating how to read the output of GALICE
+// and do some analysis.
+//
+/////////////////////////////////////////////////////////////////////////
+
+
+ Int_t NpadX = 252; // number of pads on X
+ Int_t NpadY = 374; // number of pads on Y
+
+ Int_t Pad[252][374];
+ for (Int_t i=0;i<NpadX;i++) {
+ for (Int_t j=0;j<NpadY;j++) {
+ Pad[i][j]=0;
+ }
+ }
+
+
+
+// Dynamically link some shared libs
+
+ if (gClassTable->GetID("AliRun") < 0) {
+ gROOT->LoadMacro("loadlibs.C");
+ loadlibs();
+ }
+
+// Connect the Root Galice file containing Geometry, Kine and Hits
+
+ TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
+ if (!file) file = new TFile("galice.root");
+
+// Get AliRun object from file or create it if not on file
+
+ if (!gAlice) {
+ gAlice = (AliRun*)file->Get("gAlice");
+ if (gAlice) printf("AliRun object found on file\n");
+ if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
+ }
+
+// Create some histograms
+
+ Int_t xmin=-NpadX/2;
+ Int_t xmax= NpadX/2;
+ Int_t ymin=-NpadY/2;
+ Int_t ymax= NpadY/2;
+
+ TH2F *hc = new TH2F("hc","Chamber1 signal distribution",NpadX,xmin,xmax,NpadY,ymin,ymax);
+ TH2F *h = new TH2F("h","Chamber1 hit distribution",100,-100,100,100,-100,100);
+ TH1F *charge = new TH1F("charge","Charge distribution",100,0.,1000.);
+
+// Start loop over events
+
+ Int_t Nh=0;
+ Int_t Nh1=0;
+ for (int nev=0; nev<= evNumber2; nev++) {
+ Int_t nparticles = gAlice->GetEvent(nev);
+ //cout<<"nev "<<nev<<endl;
+ printf ("Number of Events : %d\n",nev);
+ //cout<<"nparticles "<<nparticles<<endl;
+ printf ("Number of particles: %d\n",nparticles);
+ if (nev < evNumber1) continue;
+ if (nparticles <= 0) return;
+
+// Get pointers to RICH detector and Hits containers
+
+ AliRICH *RICH = gAlice->GetDetector("RICH");
+ TTree *TH = gAlice->TreeH();
+ Int_t ntracks = TH->GetEntries();
+
+// Start loop on tracks in the hits containers
+
+ // Int_t Nh=0;
+ Int_t Nc=0;
+ for (Int_t track=0; track<ntracks;track++) {
+ // printf("ntracks %d\n",ntracks);
+ gAlice->ResetHits();
+ Int_t nbytes += TH->GetEvent(track);
+ if (RICH) {
+ TClonesArray *Clusters = RICH->Clusters(); // Cluster branch
+ TClonesArray *Hits = RICH->Hits(); // Hits branch
+ //printf("%d %d \n",Clusters,Hits);
+ }
+ //see hits distribution
+ Int_t nhits = Hits->GetEntriesFast();
+ if (nhits) Nh+=nhits;
+ // printf("nhits %d\n",nhits);
+ for (Int_t hit=0;hit<nhits;hit++) {
+ mHit = (AliRICHhit*) Hits->UncheckedAt(hit);
+ Int_t nch = mHit->fChamber; // chamber number
+ Float_t x = mHit->fX; // x-pos of hit
+ Float_t y = mHit->fY; // y-pos
+ // Fill the histograms
+ if( nch==1) {
+ Float_t rhit=TMath::Sqrt(x*x+y*y);
+ if( rhit<= 55 ) Nh1+=nhits;
+ h->Fill(x,y,(float) 1);
+ }
+ }
+ // see signal distribution
+ Int_t nclust = Clusters->GetEntriesFast();
+ // printf("nclust %d\n",nclust);
+ if (nclust) {
+ Nc+=nclust;
+ for (Int_t hit=0;hit<nclust;hit++) {
+ padHit = (AliRICHcluster*) Clusters->UncheckedAt(hit);
+ Int_t nchamber = padHit->fChamber; // chamber number
+ Int_t nhit = padHit->fHitNumber; // hit number
+ Int_t qtot = padHit->fQ; // charge
+ Int_t ipx = padHit->fPadX; // pad number on X
+ Int_t ipy = padHit->fPadY; // pad number on Y
+ Int_t iqpad = padHit->fQpad; // charge per pad
+ Int_t rpad = padHit->fRpad; // R-position of pad
+ if (qtot > 0 && nchamber==1){
+ charge->Fill(qtot,(float) 1);
+ }
+ if(rpad <= 55 && nchamber==1) {
+ // if(nchamber==1) {
+ Pad[ipx+126][ipy+187]+=(iqpad);
+ hc->Fill(ipx,ipy,(float) iqpad);
+ }
+ }
+ }
+ }
+ // cout<<"Nh "<<Nh<<endl;
+ //cout<<"Nc "<<Nc<<endl;
+ printf ("Cluster Number: %d\n",Nc);
+ }
+//Create a canvas, set the view range, show histograms
+ TCanvas *c1 = new TCanvas("c1","Alice RICH pad hits",400,10,600,700);
+ hc->SetXTitle("ix (npads)");
+ hc->Draw();
+ TCanvas *c2 = new TCanvas("c2","Alice RICH hits",400,10,600,700);
+ h->SetFillColor(42);
+ h->SetXTitle("x (cm)");
+ h->Draw();
+ TCanvas *c3 = new TCanvas("c3","Charge distribution",400,10,600,700);
+ charge->SetFillColor(42);
+ charge->SetXTitle("ADC units");
+ charge->Draw();
+
+
+ // calculate the number of pads which give a signal
+
+
+ Int_t Np=0;
+ for (Int_t i=0;i< NpadX;i++) {
+ for (Int_t j=0;j< NpadY;j++) {
+ if (Pad[i][j]>=6){
+ Np+=1;
+ // cout<<"i, j "<<i<<" "<<j<<endl;
+ }
+ }
+ }
+
+ //cout<<"The total number of pads which give a signal "<<Np<<endl;
+ //cout<<"Nh "<<Nh<<endl;
+ //cout<<"Nh1 "<<Nh1<<endl;
+ printf("The total number of pads which give a signal: %d %d\n",Nh,Nh1);
+
+}
--- /dev/null
+void RICHtest (Int_t evNumber1=0,Int_t evNumber2=0)
+{
+/////////////////////////////////////////////////////////////////////////
+// This macro is a small example of a ROOT macro
+// illustrating how to read the output of GALICE
+// and do some analysis.
+//
+/////////////////////////////////////////////////////////////////////////
+
+// Dynamically link some shared libs
+
+ if (gClassTable->GetID("AliRun") < 0) {
+ gROOT->LoadMacro("loadlibs.C");
+ loadlibs();
+ }
+// Connect the Root Galice file containing Geometry, Kine and Hits
+
+ TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
+
+
+ if (!file) {
+ printf("\n Creating galice.root \n");
+ file = new TFile("galice.root");
+ } else {
+ printf("\n galice.root found in file list");
+ }
+ file->ls();
+
+// Get AliRun object from file or create it if not on file
+ if (!gAlice) {
+ gAlice = (AliRun*)(file->Get("gAlice"));
+ if (gAlice) printf("AliRun object found on file\n");
+ if (!gAlice) {
+ printf("\n create new gAlice object");
+ gAlice = new AliRun("gAlice","Alice test program");
+ }
+ }
+
+// Create some histograms
+
+
+ TH1F *ccharge1 = new TH1F("ccharge1","Cluster Charge distribution"
+ ,100,0.,500.);
+ TH1F *pcharge1 = new TH1F("pcharge1","Pad Charge distribution"
+ ,100,0.,200.);
+ TH1F *xresid1 = new TH1F("xresid1","x-Residuals"
+ ,100,-0.5,0.5);
+ TH1F *yresid1 = new TH1F("yresid1","y-Residuals"
+ ,100,-0.2,0.2);
+ TH1F *npads1 = new TH1F("npads1" ,"Pads per Hit"
+ ,20,-0.5,19.5);
+ TH2F *xresys1 = new TH2F("xresys1","x-Residuals systematics"
+ ,50,-2.0,2.0,100,-1.0,1.0);
+ TH2F *yresys1 = new TH2F("yresys1","y-Residuals systematics"
+ ,50,-0.4,0.4,100,-0.1,0.1);
+
+ TH1F *ccharge2 = new TH1F("ccharge2","Cluster Charge distribution"
+ ,100,0.,500.);
+ TH1F *pcharge2 = new TH1F("pcharge2","Pad Charge distribution"
+ ,100,0.,200.);
+ TH1F *xresid2 = new TH1F("xresid2","x-Residuals"
+ ,100,-0.5,0.5);
+ TH1F *yresid2 = new TH1F("yresid2","y-Residuals"
+ ,100,-0.2,0.2);
+ TH1F *npads2 = new TH1F("npads2" ,"Pads per Hit"
+ ,20,-0.5,19.5);
+ TH2F *xresys2 = new TH2F("xresys2","x-Residuals systematics"
+ ,50,-2.0,2.0,100,-1.0,1.0);
+ TH2F *yresys2 = new TH2F("yresys2","y-Residuals systematics"
+ ,50,-0.4,0.4,100,-0.1,0.1);
+
+
+ AliRICHchamber* iChamber;
+//
+// Loop over events
+//
+ Int_t Nh=0;
+ Int_t Nh1=0;
+ for (Int_t nev=0; nev<= evNumber2; nev++) {
+ //cout << "nev " << nev <<endl;
+ printf ("Event Number : %d\n",nev);
+ Int_t nparticles = gAlice->GetEvent(nev);
+ //cout << "nparticles " << nparticles <<endl;
+ printf ("Number of particles: %d\n",nparticles);
+ if (nev < evNumber1) continue;
+ if (nparticles <= 0) return;
+
+ AliRICH *RICH = (AliRICH*) gAlice->GetDetector("RICH");
+ printf("\n track %d %d \n ", nev, RICH);
+
+ TTree *TH = gAlice->TreeH();
+ Int_t ntracks = TH->GetEntries();
+ Int_t Nc=0;
+ Int_t npad[2];
+ Float_t Q[2], xmean[2],ymean[2],xres[2],yres[2], xonpad[2], yonpad[2];
+
+ printf ("Just got it");
+//
+// Loop over events
+//
+ for (Int_t track=0; track<ntracks;track++) {
+
+ gAlice->ResetHits();
+ Int_t nbytes += TH->GetEvent(track);
+ if (RICH) {
+ for(AliRICHhit* mHit=(AliRICHhit*)RICH->FirstHit(-1);
+ mHit;
+ mHit=(AliRICHhit*)RICH->NextHit())
+ {
+ Int_t nch = mHit->fChamber; // chamber number
+ Float_t x = mHit->fX; // x-pos of hit
+ Float_t y = mHit->fY; // y-pos
+//
+//
+ iChamber = & RICH->Chamber(nch-1);
+ response=iChamber->GetResponseModel(mip);
+//
+//
+ if (nch <= 1) {
+ for (Int_t i=0; i<2; i++) {
+ xmean[i]=0;
+ ymean[i]=0;
+ Q[i]=0;
+ npad[i]=0;
+ }
+ for (AliRICHcluster *mPad=(AliRICHcluster*)RICH
+ ->FirstPad(mHit);
+ mPad;
+ mPad=(AliRICHcluster*)RICH->NextPad()
+ )
+ {
+ Int_t nseg = mPad->fCathode; // segmentation
+ Int_t nhit = mPad->fHitNumber; // hit number
+ Int_t qtot = mPad->fQ; // charge
+ Int_t ipx = mPad->fPadX; // pad number on X
+ Int_t ipy = mPad->fPadY; // pad number on Y
+ Int_t iqpad = mPad->fQpad; // charge per pad
+ Int_t secpad = mPad->fRSec; // r-pos of pad
+//
+//
+ segmentation=iChamber->GetSegmentationModel(nseg);
+ Int_t ind=nseg-1;
+// printf("Pad of hit, padx, pady, iqpad, ncha %d %d %d %d %d %d\n",
+// nhit, ipx, ipy, iqpad, nseg, ind);
+
+
+ if (nseg == 1) {
+ pcharge1->Fill(iqpad,(float) 1);
+ ccharge1->Fill(qtot ,(float) 1);
+ } else {
+ pcharge2->Fill(iqpad,(float) 1);
+ ccharge2->Fill(qtot ,(float) 1);
+ }
+// Calculate centre of gravity
+//
+ if (iqpad == 0 && ind==0) {
+ printf("\n attention iqpad 0");
+ }
+
+ if (iqpad > 0 && secpad==1) {
+ Float_t xc,yc;
+ npad[ind]++;
+ segmentation->GetPadCxy(ipx,ipy,xc,yc);
+// printf("\n pad %d %d ", ipx, ipy);
+// printf("\n pos %f %f ", xc, yc);
+ xmean[ind]+=Float_t(iqpad*xc);
+ ymean[ind]+=Float_t(iqpad*yc);
+ Q[ind]+=iqpad;
+ }
+
+ } //pad hit loop
+ for (Int_t i=0; i<2; i++) {
+ segmentation = iChamber->GetSegmentationModel(i+1);
+
+// if (npad[i] ==0) {
+// printf("\n #pads=0 %f %f",x,y);
+// }
+
+ if (Q[i] >0) {
+ xmean[i] = xmean[i]/Q[i];
+ xres[i] = xmean[i]-x;
+ ymean[i] = ymean[i]/Q[i];
+ yres[i] = ymean[i]-y;
+// printf("\n xmean, ymean %f %f",xmean[i],ymean[i]);
+// printf("\n xhit, yhit %f %f",x,y);
+
+// Systematic Error
+//
+ Int_t icx, icy;
+ segmentation->
+ GetPadIxy(x,y,icx,icy);
+ segmentation->
+ GetPadCxy(icx,icy,xonpad[i],yonpad[i]);
+ xonpad[i]-=x;
+ yonpad[i]-=y;
+ }
+ } // plane loop
+ xresid1->Fill(xres[0],(float) 1);
+ yresid1->Fill(yres[0],(float) 1);
+ npads1->Fill(npad[0],(float) 1);
+ if (npad[0] >=2 && Q[0] > 6 ) {
+ xresys1->Fill(xonpad[0],xres[0],(float) 1);
+ yresys1->Fill(yonpad[0],yres[0],(float) 1);
+ }
+
+ xresid2->Fill(xres[1],(float) 1);
+ yresid2->Fill(yres[1],(float) 1);
+ npads2->Fill(npad[1],(float) 1);
+ if (npad[1] >=2 && Q[1] > 6) {
+ xresys2->Fill(xonpad[1],xres[1],(float) 1);
+ yresys2->Fill(yonpad[1],yres[1],(float) 1);
+ }
+ } // chamber 1
+ } // hit loop
+ } // if RICH
+ } // track loop
+ } // event loop
+
+//Create a canvas, set the view range, show histograms
+ TCanvas *c1 = new TCanvas("c1","Charge and Residuals",400,10,600,700);
+ pad11 = new TPad("pad11"," ",0.01,0.51,0.49,0.99);
+ pad12 = new TPad("pad12"," ",0.51,0.51,0.99,0.99);
+ pad13 = new TPad("pad13"," ",0.01,0.01,0.49,0.49);
+ pad14 = new TPad("pad14"," ",0.51,0.01,0.99,0.49);
+ pad11->SetFillColor(11);
+ pad12->SetFillColor(11);
+ pad13->SetFillColor(11);
+ pad14->SetFillColor(11);
+ pad11->Draw();
+ pad12->Draw();
+ pad13->Draw();
+ pad14->Draw();
+
+ pad11->cd();
+ ccharge1->SetFillColor(42);
+ ccharge1->SetXTitle("ADC units");
+ ccharge1->Draw();
+
+ pad12->cd();
+ pcharge1->SetFillColor(42);
+ pcharge1->SetXTitle("ADC units");
+ pcharge1->Draw();
+
+ pad13->cd();
+ xresid1->SetFillColor(42);
+ xresid1->Draw();
+
+ pad14->cd();
+ yresid1->SetFillColor(42);
+ yresid1->Draw();
+
+ TCanvas *c2 = new TCanvas("c2","Cluster Size",400,10,600,700);
+ pad21 = new TPad("pad21"," ",0.01,0.51,0.49,0.99);
+ pad22 = new TPad("pad22"," ",0.51,0.51,0.99,0.99);
+ pad23 = new TPad("pad23"," ",0.01,0.01,0.49,0.49);
+ pad24 = new TPad("pad24"," ",0.51,0.01,0.99,0.49);
+ pad21->SetFillColor(11);
+ pad22->SetFillColor(11);
+ pad23->SetFillColor(11);
+ pad24->SetFillColor(11);
+ pad21->Draw();
+ pad22->Draw();
+ pad23->Draw();
+ pad24->Draw();
+
+ pad21->cd();
+ npads1->SetFillColor(42);
+ npads1->SetXTitle("Cluster Size");
+ npads1->Draw();
+
+ pad23->cd();
+ xresys1->SetXTitle("x on pad");
+ xresys1->SetYTitle("x-xcog");
+ xresys1->Draw();
+
+ pad24->cd();
+ yresys1->SetXTitle("y on pad");
+ yresys1->SetYTitle("y-ycog");
+ yresys1->Draw();
+
+ TCanvas *c3 = new TCanvas("c3","Charge and Residuals",400,10,600,700);
+ pad31 = new TPad("pad31"," ",0.01,0.51,0.49,0.99);
+ pad32 = new TPad("pad32"," ",0.51,0.51,0.99,0.99);
+ pad33 = new TPad("pad33"," ",0.01,0.01,0.49,0.49);
+ pad34 = new TPad("pad34"," ",0.51,0.01,0.99,0.49);
+ pad31->SetFillColor(11);
+ pad32->SetFillColor(11);
+ pad33->SetFillColor(11);
+ pad34->SetFillColor(11);
+ pad31->Draw();
+ pad32->Draw();
+ pad33->Draw();
+ pad34->Draw();
+
+ pad31->cd();
+ ccharge2->SetFillColor(42);
+ ccharge2->SetXTitle("ADC units");
+ ccharge2->Draw();
+
+ pad32->cd();
+ pcharge2->SetFillColor(42);
+ pcharge2->SetXTitle("ADC units");
+ pcharge2->Draw();
+
+ pad33->cd();
+ xresid2->SetFillColor(42);
+ xresid2->Draw();
+
+ pad34->cd();
+ yresid2->SetFillColor(42);
+ yresid2->Draw();
+
+ TCanvas *c4 = new TCanvas("c4","Cluster Size",400,10,600,700);
+ pad41 = new TPad("pad41"," ",0.01,0.51,0.49,0.99);
+ pad42 = new TPad("pad42"," ",0.51,0.51,0.99,0.99);
+ pad43 = new TPad("pad43"," ",0.01,0.01,0.49,0.49);
+ pad44 = new TPad("pad44"," ",0.51,0.01,0.99,0.49);
+ pad41->SetFillColor(11);
+ pad42->SetFillColor(11);
+ pad43->SetFillColor(11);
+ pad44->SetFillColor(11);
+ pad41->Draw();
+ pad42->Draw();
+ pad43->Draw();
+ pad44->Draw();
+
+ pad41->cd();
+ npads2->SetFillColor(42);
+ npads2->SetXTitle("Cluster Size");
+ npads2->Draw();
+
+ pad43->cd();
+ xresys2->SetXTitle("x on pad");
+ xresys2->SetYTitle("x-xcog");
+ xresys2->Draw();
+
+ pad44->cd();
+ yresys2->SetXTitle("y on pad");
+ yresys2->SetYTitle("y-ycog");
+ yresys2->Draw();
+
+}
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff