From ddae0931da828b889cfc1f4041d436a366460aa4 Mon Sep 17 00:00:00 2001 From: fca Date: Tue, 13 Jul 1999 10:08:39 +0000 Subject: [PATCH] New RICH code J.Barbosa, A.Morsch, D.DiBari --- RICH/AliRICH.cxx | 2865 ++++++++++++++----------------------- RICH/AliRICH.h | 721 +++++++--- RICH/AliRICHConst.h | 24 + RICH/AliRICHHitMap.cxx | 99 ++ RICH/AliRICHHitMap.h | 53 + RICH/AliRICHSegResCkv.cxx | 15 + RICH/AliRICHSegResCkv.h | 19 + RICH/AliRICHSegResV0.cxx | 321 +++++ RICH/AliRICHSegResV0.h | 154 ++ RICH/AliRICHdisplay.cxx | 978 +++++++++++++ RICH/AliRICHdisplay.h | 102 ++ RICH/AliRICHpoints.cxx | 272 ++++ RICH/AliRICHpoints.h | 35 + RICH/AliRICHv0.cxx | 830 +++++++++++ RICH/AliRICHv0.h | 34 + RICH/Makefile | 2 +- RICH/RICHLinkDef.h | 27 +- RICH/RICHdigit.C | 70 + RICH/RICHdisplay.C | 34 + RICH/RICHpadtest.C | 163 +++ RICH/RICHtest.C | 342 +++++ 21 files changed, 5189 insertions(+), 1971 deletions(-) create mode 100644 RICH/AliRICHConst.h create mode 100644 RICH/AliRICHHitMap.cxx create mode 100644 RICH/AliRICHHitMap.h create mode 100644 RICH/AliRICHSegResCkv.cxx create mode 100644 RICH/AliRICHSegResCkv.h create mode 100644 RICH/AliRICHSegResV0.cxx create mode 100644 RICH/AliRICHSegResV0.h create mode 100644 RICH/AliRICHdisplay.cxx create mode 100644 RICH/AliRICHdisplay.h create mode 100644 RICH/AliRICHpoints.cxx create mode 100644 RICH/AliRICHpoints.h create mode 100644 RICH/AliRICHv0.cxx create mode 100644 RICH/AliRICHv0.h create mode 100644 RICH/RICHdigit.C create mode 100644 RICH/RICHdisplay.C create mode 100644 RICH/RICHpadtest.C create mode 100644 RICH/RICHtest.C diff --git a/RICH/AliRICH.cxx b/RICH/AliRICH.cxx index 0e632026e00..df0160c6847 100644 --- a/RICH/AliRICH.cxx +++ b/RICH/AliRICH.cxx @@ -1,1924 +1,1199 @@ -/////////////////////////////////////////////////////////////////////////////// -// // -// 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 -/* - -*/ -//End_Html -// // -// // -/////////////////////////////////////////////////////////////////////////////// +//////////////////////////////////////////////// +// Manager and hits classes for set:RICH // +//////////////////////////////////////////////// -#include +#include +#include #include #include +#include +#include +#include -#include #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 +/* + +*/ +//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;iGetGeometry()->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 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)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; digitUncheckedAt(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;digitUncheckedAt(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;IlistAddDigit(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 -/* - -*/ -//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 - /* - - */ - //End_Html - //Begin_Html - /* - - */ - //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;iAt(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 (fCurrentDigitAt(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"<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; trackResetHits(); + 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;trAt(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; trackClear(); + if (fClusters2) fClusters2->Clear(); + + TH1->GetEvent(track); +// +// Loop over hits + AliRICHhit* mHit; + for(int i=0;iGetEntriesFast();++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;trAt(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 "<GetEntriesFast(); + printf(" \n \n nentries %d \n",nentries); + + // start filling the digits + + for (Int_t nent=0;nentAt(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 "< 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;trAt(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"<TreeD()->Fill(); + TTree *TD=gAlice->TreeD(); + Stat_t ndig=TD->GetEntries(); + cout<<"number of digits "<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); } + + + + + + + diff --git a/RICH/AliRICH.h b/RICH/AliRICH.h index 8753ad1e01c..6c28e4a1469 100644 --- a/RICH/AliRICH.h +++ b/RICH/AliRICH.h @@ -3,190 +3,575 @@ //////////////////////////////////////////////// // Manager and hits classes for set:RICH // //////////////////////////////////////////////// - #include "AliDetector.h" #include "AliHit.h" +#include "AliRICHConst.h" +#include "AliDigit.h" +#include +#include +#include +#include +#include -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 + + + + + + + + + + + + diff --git a/RICH/AliRICHConst.h b/RICH/AliRICHConst.h new file mode 100644 index 00000000000..4d83f31cc06 --- /dev/null +++ b/RICH/AliRICHConst.h @@ -0,0 +1,24 @@ +#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 diff --git a/RICH/AliRICHHitMap.cxx b/RICH/AliRICHHitMap.cxx new file mode 100644 index 00000000000..1d2774b1737 --- /dev/null +++ b/RICH/AliRICHHitMap.cxx @@ -0,0 +1,99 @@ +#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; ndigUncheckedAt(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; + } +} + + diff --git a/RICH/AliRICHHitMap.h b/RICH/AliRICHHitMap.h new file mode 100644 index 00000000000..903efad1e74 --- /dev/null +++ b/RICH/AliRICHHitMap.h @@ -0,0 +1,53 @@ +#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 + + diff --git a/RICH/AliRICHSegResCkv.cxx b/RICH/AliRICHSegResCkv.cxx new file mode 100644 index 00000000000..539e7467729 --- /dev/null +++ b/RICH/AliRICHSegResCkv.cxx @@ -0,0 +1,15 @@ +#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; +} + diff --git a/RICH/AliRICHSegResCkv.h b/RICH/AliRICHSegResCkv.h new file mode 100644 index 00000000000..bf61d6714f9 --- /dev/null +++ b/RICH/AliRICHSegResCkv.h @@ -0,0 +1,19 @@ +#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 diff --git a/RICH/AliRICHSegResV0.cxx b/RICH/AliRICHSegResV0.cxx new file mode 100644 index 00000000000..aeba88b0643 --- /dev/null +++ b/RICH/AliRICHSegResV0.cxx @@ -0,0 +1,321 @@ +#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 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); +} diff --git a/RICH/AliRICHSegResV0.h b/RICH/AliRICHSegResV0.h new file mode 100644 index 00000000000..8efae70f9f6 --- /dev/null +++ b/RICH/AliRICHSegResV0.h @@ -0,0 +1,154 @@ +#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 diff --git a/RICH/AliRICHdisplay.cxx b/RICH/AliRICHdisplay.cxx new file mode 100644 index 00000000000..4364c8a5a04 --- /dev/null +++ b/RICH/AliRICHdisplay.cxx @@ -0,0 +1,978 @@ + +////////////////////////////////////////////////////////////////////////// +// // +// AliDisplay // +// // +// Utility class to display ALICE outline, tracks, hits,.. // +// // +////////////////////////////////////////////////////////////////////////// + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#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 +/* + +*/ +//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;trackUncheckedAt(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;digitUncheckedAt(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;trackUncheckedAt(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;digitUncheckedAt(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; trackResetHits(); + 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;hitUncheckedAt(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;pSetMarkerColor(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; trackResetHits(); + 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;hitUncheckedAt(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;pSetMarkerColor(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; + } +} diff --git a/RICH/AliRICHdisplay.h b/RICH/AliRICHdisplay.h new file mode 100644 index 00000000000..592c29acc05 --- /dev/null +++ b/RICH/AliRICHdisplay.h @@ -0,0 +1,102 @@ +#ifndef AliRICHdisplay_H +#define AliRICHdisplay_H + +////////////////////////////////////////////////////////////////////////// +// // +// AliDisplay // +// // +// Utility class to display ALice outline, tracks, hits,.. // +// // +////////////////////////////////////////////////////////////////////////// + +//#ifndef ROOT_TObject +#include +//#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 + diff --git a/RICH/AliRICHpoints.cxx b/RICH/AliRICHpoints.cxx new file mode 100644 index 00000000000..06a0bb96f1a --- /dev/null +++ b/RICH/AliRICHpoints.cxx @@ -0,0 +1,272 @@ +/////////////////////////////////////////////////////////////////////////////// +// // +// This class contains the points for the ALICE event display // +// // +//Begin_Html +/* + +*/ +//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; ndigUncheckedAt(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); + +} + + + diff --git a/RICH/AliRICHpoints.h b/RICH/AliRICHpoints.h new file mode 100644 index 00000000000..ba6ce48bf22 --- /dev/null +++ b/RICH/AliRICHpoints.h @@ -0,0 +1,35 @@ +#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 diff --git a/RICH/AliRICHv0.cxx b/RICH/AliRICHv0.cxx new file mode 100644 index 00000000000..df89ff41f33 --- /dev/null +++ b/RICH/AliRICHv0.cxx @@ -0,0 +1,830 @@ +///////////////////////////////////////////////////////// +// Manager and hits classes for set:RICH version 0 // +///////////////////////////////////////////////////////// + +#include +#include +#include + +#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 + /* + + */ + //End_Html + //Begin_Html + /* + + */ + //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 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 +} + + + + + + + + + + + + + + + + + + + diff --git a/RICH/AliRICHv0.h b/RICH/AliRICHv0.h new file mode 100644 index 00000000000..db448b062bb --- /dev/null +++ b/RICH/AliRICHv0.h @@ -0,0 +1,34 @@ +#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 + + + + + + + diff --git a/RICH/Makefile b/RICH/Makefile index bcfa4920bb5..d291f75f14e 100644 --- a/RICH/Makefile +++ b/RICH/Makefile @@ -9,7 +9,7 @@ PACKAGE = RICH # C++ sources -SRCS = AliRICH.cxx +SRCS = AliRICH.cxx AliRICHv0.cxx AliRICHdisplay.cxx AliRICHpoints.cxx AliRICHHitMap.cxx AliRICHSegResV0.cxx AliRICHSegResCkv.cxx # C++ Headers diff --git a/RICH/RICHLinkDef.h b/RICH/RICHLinkDef.h index ea49fc7efe0..452c95acd6e 100644 --- a/RICH/RICHLinkDef.h +++ b/RICH/RICHLinkDef.h @@ -1,14 +1,27 @@ #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 diff --git a/RICH/RICHdigit.C b/RICH/RICHdigit.C new file mode 100644 index 00000000000..cae6d8fa0c9 --- /dev/null +++ b/RICH/RICHdigit.C @@ -0,0 +1,70 @@ +#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 " <Digitise(nev); + char hname[30]; + sprintf(hname,"TreeD%d",nev); + gAlice->TreeD()->Write(hname); + gAlice->TreeD()->Reset(); + file->ls(); + } // event loop + file->Close(); +} + + + + + + + + + + + + + + diff --git a/RICH/RICHdisplay.C b/RICH/RICHdisplay.C new file mode 100644 index 00000000000..23e765e6114 --- /dev/null +++ b/RICH/RICHdisplay.C @@ -0,0 +1,34 @@ +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); +} diff --git a/RICH/RICHpadtest.C b/RICH/RICHpadtest.C new file mode 100644 index 00000000000..03f3a8f2750 --- /dev/null +++ b/RICH/RICHpadtest.C @@ -0,0 +1,163 @@ +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;iGetID("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 "<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; trackResetHits(); + 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;hitUncheckedAt(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;hitUncheckedAt(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 "<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 "<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 <GetEvent(nev); + //cout << "nparticles " << nparticles <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; trackResetHits(); + 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(); + +} -- 2.39.3