-///////////////////////////////////////////////////////////////////////////////
-// //
-// 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
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * *
+ * Author: The ALICE Off-line Project. *
+ * Contributors are mentioned in the code where appropriate. *
+ * *
+ * Permission to use, copy, modify and distribute this software and its *
+ * documentation strictly for non-commercial purposes is hereby granted *
+ * without fee, provided that the above copyright notice appears in all *
+ * copies and that both the copyright notice and this permission notice *
+ * appear in the supporting documentation. The authors make no claims *
+ * about the suitability of this software for any purpose. It is *
+ * provided "as is" without express or implied warranty. *
+ **************************************************************************/
+
/*
-<img src="picts/AliRICHClass.gif">
+ $Log$
+ Revision 1.18 2000/06/12 15:15:46 jbarbosa
+ Cleaned up version.
+
+ Revision 1.17 2000/06/09 14:58:37 jbarbosa
+ New digitisation per particle type
+
+ Revision 1.16 2000/04/19 12:55:43 morsch
+ Newly structured and updated version (JB, AM)
+
*/
-//End_Html
-// //
-// //
-///////////////////////////////////////////////////////////////////////////////
-#include <TBRIK.h>
+
+////////////////////////////////////////////////
+// Manager and hits classes for set:RICH //
+////////////////////////////////////////////////
+
+#include <TBRIK.h>
+#include <TTUBE.h>
#include <TNode.h>
#include <TRandom.h>
+#include <TObject.h>
+#include <TVector.h>
+#include <TObjArray.h>
+#include <TArrayF.h>
+#include <TFile.h>
+#include <TParticle.h>
+#include <iostream.h>
-#include <TClass.h>
#include "AliRICH.h"
+#include "AliRICHSegmentation.h"
+#include "AliRICHHit.h"
+#include "AliRICHCerenkov.h"
+#include "AliRICHPadHit.h"
+#include "AliRICHDigit.h"
+#include "AliRICHTransientDigit.h"
+#include "AliRICHRawCluster.h"
+#include "AliRICHRecHit.h"
+#include "AliRICHHitMapA1.h"
+#include "AliRICHClusterFinder.h"
#include "AliRun.h"
-#include "TGeant3.h"
+#include "AliMC.h"
+#include "AliPoints.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;
+static TTree *TrH1;
+
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;
+// Default constructor for RICH manager class
+
+ fIshunt = 0;
+ fHits = 0;
+ fPadHits = 0;
+ fNPadHits = 0;
+ fNcerenkovs = 0;
+ fDchambers = 0;
+ fCerenkovs = 0;
+ fNdch = 0;
}
-
-//_____________________________________________________________________________
+
+//___________________________________________
AliRICH::AliRICH(const char *name, const char *title)
- : AliDetector(name,title)
+ : AliDetector(name,title)
{
- //
- // Standard constructor for RICH
- //
- fIshunt = 0;
- fNmips = 0;
- fNckovs = 0;
- fNpadhits = 0;
- //
- // Allocate space for different components
- fHits = new TClonesArray("AliRICHhit", 100);
- fMips = new TClonesArray("AliRICHmip", 100);
- fCkovs = new TClonesArray("AliRICHckov", 100);
- fPadhits = new TClonesArray("AliRICHpadhit", 100);
- //
- // Set parameters to default value
- fChslope = 40;
- fAlphaFeed= 0.04;
- fSxcharge = 0.18;
- fIritri = 0;
- //
- SetMarkerColor(6);
- SetMarkerStyle(20);
- SetMarkerSize(0.5);
+//Begin_Html
+/*
+ <img src="gif/alirich.gif">
+*/
+//End_Html
+
+ fHits = new TClonesArray("AliRICHHit",1000 );
+ gAlice->AddHitList(fHits);
+ fPadHits = new TClonesArray("AliRICHPadHit",100000);
+ fCerenkovs = new TClonesArray("AliRICHCerenkov",1000);
+ gAlice->AddHitList(fCerenkovs);
+ //gAlice->AddHitList(fHits);
+ fNPadHits = 0;
+ fNcerenkovs = 0;
+ fIshunt = 0;
+
+ fNdch = new Int_t[kNCH];
+
+ fDchambers = new TObjArray(kNCH);
+
+ fRecHits = new TObjArray(kNCH);
+
+ Int_t i;
+
+ for (i=0; i<kNCH ;i++) {
+ (*fDchambers)[i] = new TClonesArray("AliRICHDigit",10000);
+ fNdch[i]=0;
+ }
+
+ fNrawch = new Int_t[kNCH];
+
+ fRawClusters = new TObjArray(kNCH);
+ //printf("Created fRwClusters with adress:%p",fRawClusters);
+
+ for (i=0; i<kNCH ;i++) {
+ (*fRawClusters)[i] = new TClonesArray("AliRICHRawCluster",10000);
+ fNrawch[i]=0;
+ }
+
+ fNrechits = new Int_t[kNCH];
+
+ for (i=0; i<kNCH ;i++) {
+ (*fRecHits)[i] = new TClonesArray("AliRICHRecHit",1000);
+ }
+ //printf("Created fRecHits with adress:%p",fRecHits);
+
+
+ SetMarkerColor(kRed);
}
-//_____________________________________________________________________________
+AliRICH::AliRICH(const AliRICH& RICH)
+{
+// Copy Constructor
+}
+
+
+//___________________________________________
AliRICH::~AliRICH()
{
- //
- // Destructor for RICH
- //
- fIshunt = 0;
- delete fHits;
- fMips->Delete(); delete fMips;
- fCkovs->Delete(); delete fCkovs;
- fPadhits->Delete(); delete fPadhits;
+
+// Destructor of RICH manager class
+
+ fIshunt = 0;
+ delete fHits;
+ delete fPadHits;
+ delete fCerenkovs;
}
-//_____________________________________________________________________________
+//___________________________________________
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;
- }
-}
+//
+// Adds a hit to the Hits list
+//
+
+ 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)
+{
+
+//
+// Adds a RICH cerenkov hit to the Cerenkov Hits list
+//
+
+ TClonesArray &lcerenkovs = *fCerenkovs;
+ new(lcerenkovs[fNcerenkovs++]) AliRICHCerenkov(fIshunt,track,vol,cerenkovs);
+ //printf ("Done for Cerenkov %d\n\n\n\n",fNcerenkovs);
}
+//___________________________________________
+void AliRICH::AddPadHit(Int_t *clhits)
+{
+
+//
+// Add a RICH pad hit to the list
+//
+ TClonesArray &lPadHits = *fPadHits;
+ new(lPadHits[fNPadHits++]) AliRICHPadHit(clhits);
+}
//_____________________________________________________________________________
-void AliRICH::AddCkovHit(Int_t track, Int_t *vol, Float_t *hits)
-{
- //
- // Adds a cerenkov hit in the RICH.
+void AliRICH::AddDigits(Int_t id, Int_t *tracks, Int_t *charges, Int_t *digits)
+{
+
//
- TClonesArray &lhits = *fCkovs;
- AliRICHmip *lmip = (AliRICHmip*) fMips->AddrAt(fNmips-1);
+ // Add a RICH digit to the list
//
- // 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);
+
+ TClonesArray &ldigits = *((TClonesArray*)(*fDchambers)[id]);
+ new(ldigits[fNdch[id]++]) AliRICHDigit(tracks,charges,digits);
+}
+
+//_____________________________________________________________________________
+void AliRICH::AddRawCluster(Int_t id, const AliRICHRawCluster& c)
+{
+ //
+ // Add a RICH digit to the list
+ //
+
+ TClonesArray &lrawcl = *((TClonesArray*)(*fRawClusters)[id]);
+ new(lrawcl[fNrawch[id]++]) AliRICHRawCluster(c);
}
//_____________________________________________________________________________
-void AliRICH::AddPadHit(Int_t track, Int_t *vol, Float_t *hits)
-{
+void AliRICH::AddRecHit(Int_t id, Float_t *rechit, Float_t *photons, Int_t *padsx, Int_t* padsy)
+{
+
//
- // Adds pad hits in the RICH
+ // Add a RICH reconstructed hit to the list
//
- 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);
- }
+
+ TClonesArray &lrec = *((TClonesArray*)(*fRecHits)[id]);
+ new(lrec[fNrechits[id]++]) AliRICHRecHit(id,rechit,photons,padsx,padsy);
}
-//_____________________________________________________________________________
+//___________________________________________
void AliRICH::BuildGeometry()
+
{
+
//
// 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);
-}
+ TNode *node, *top;
+
+ const int kColorRICH = kGreen;
+ //
+ top=gAlice->GetGeometry()->GetNode("alice");
+
+
+ new TBRIK("S_RICH","S_RICH","void",71.09999,11.5,73.15);
+
+ top->cd();
+ Float_t pos1[3]={0,471.8999,165.2599};
+ //Chamber(0).SetChamberTransform(pos1[0],pos1[1],pos1[2],
+ new TRotMatrix("rot993","rot993",90,0,70.69,90,19.30999,-90);
+ node = new TNode("RICH1","RICH1","S_RICH",pos1[0],pos1[1],pos1[2],"rot993");
+
-//*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
-//
-// 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;
+ 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");
+
-//_____________________________________________________________________________
-void AliRICH::Init()
-{
- //
- // Initialise RICH after that it has been built
- //
- const Float_t sconv=2*TMath::Sqrt(2*TMath::Log(2));
- const Float_t yK3=1.20;
- Float_t ansp;
- Int_t i;
- //
- sNpx=162;
- sNpy=162;
- sDxp=0.80;
- sDyp=0.80;
- ansp=sDyp/2;
- sDlx=sNpx*sDxp/2;
- sDly=sNpy*sDyp/2;
- sDpad=0.2;
- //
- for(i=0;i<sNpx;i++) {
- sXpad[i]=i*sDxp;
- sYpad[i]=i*sDyp;
- }
- for(i=0;i<2*sNpy+1;i++) sYanode[i]=ansp/2+i*ansp;
- //
-
- sSycharge=4*TMath::ATanH(1/TMath::Sqrt(2*yK3))/TMath::Pi()
- /(1-0.5*TMath::Sqrt(yK3))/sDpad/sconv;
- //
+ 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);
+
}
-
-//___________________________________________________________________________
-void AliRICH::StepManager()
+//___________________________________________
+Int_t AliRICH::DistancetoPrimitive(Int_t , Int_t )
{
- //
- // Called at every step in the RICH
- //
-
- AliMC* pMC = AliMC::GetMC();
- TGeant3 *geant3 = (TGeant3*) pMC;
- const Float_t xshift[3] = { 41.3, 0, -41.3 };
- static Float_t polar[3] = {0, 0, 0};
- 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;
- Int_t ncher;
- 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];
+// Default value
-
-
- // 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
+ return 9999;
+}
- if (geant3->Gckine()->charge) {
+//___________________________________________
+void AliRICH::MakeBranch(Option_t* option)
+{
+ // 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) {
- pMC->Gmtod(geant3->Gckin3()->gpos[i-1], vxloc, 1);
- pMC->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);
- }
- }
+ const Int_t kBufferSize = 4000;
+ char branchname[20];
+
+
+ AliDetector::MakeBranch(option);
+ sprintf(branchname,"%sCerenkov",GetName());
+ if (fCerenkovs && gAlice->TreeH()) {
+ gAlice->TreeH()->Branch(branchname,&fCerenkovs, kBufferSize);
+ printf("Making Branch %s for Cerenkov Hits\n",branchname);
}
- // 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 ?
- pMC->Gmtod(geant3->Gctrak()->vect, sVloc, 1);
- pMC->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 ?
- pMC->Gmtod(geant3->Gctrak()->vect, sVloc, 1);
- pMC->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
- pMC->Gmtod(geant3->Gctrak()->vect, sVloc, 1);
- pMC->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;
+
+ sprintf(branchname,"%sPadHits",GetName());
+ if (fPadHits && gAlice->TreeH()) {
+ gAlice->TreeH()->Branch(branchname,&fPadHits, kBufferSize);
+ printf("Making Branch %s for PadHits\n",branchname);
+ }
+
+// one branch for digits per chamber
+ Int_t i;
+
+ for (i=0; i<kNCH ;i++) {
+ sprintf(branchname,"%sDigits%d",GetName(),i+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;
- }
+ if (fDchambers && gAlice->TreeD()) {
+ gAlice->TreeD()->Branch(branchname,&((*fDchambers)[i]), kBufferSize);
+ printf("Making Branch %s for digits in chamber %d\n",branchname,i+1);
+ }
}
- }
-
- // 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) {
+
+// one branch for raw clusters per chamber
+ for (i=0; i<kNCH ;i++) {
+ sprintf(branchname,"%sRawClusters%d",GetName(),i+1);
- // 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;
+ if (fRawClusters && gAlice->TreeR()) {
+ gAlice->TreeR()->Branch(branchname,&((*fRawClusters)[i]), kBufferSize);
+ printf("Making Branch %s for raw clusters in chamber %d\n",branchname,i+1);
+ }
+ }
+
+ // one branch for rec hits per chamber
+ for (i=0; i<kNCH ;i++) {
+ sprintf(branchname,"%sRecHits%d",GetName(),i+1);
+
+ if (fRecHits && gAlice->TreeR()) {
+ gAlice->TreeR()->Branch(branchname,&((*fRecHits)[i]), kBufferSize);
+ printf("Making Branch %s for rec. hits in chamber %d\n",branchname,i+1);
+ }
+ }
+}
+
+//___________________________________________
+void AliRICH::SetTreeAddress()
+{
+ // Set branch address for the Hits and Digits Tree.
+ char branchname[20];
+ Int_t i;
+
+ AliDetector::SetTreeAddress();
+
+ TBranch *branch;
+ TTree *treeH = gAlice->TreeH();
+ TTree *treeD = gAlice->TreeD();
+ TTree *treeR = gAlice->TreeR();
+
+ if (treeH) {
+ if (fPadHits) {
+ branch = treeH->GetBranch("RICHPadHits");
+ if (branch) branch->SetAddress(&fPadHits);
}
-
- // 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];
+ if (fCerenkovs) {
+ branch = treeH->GetBranch("RICHCerenkov");
+ if (branch) branch->SetAddress(&fCerenkovs);
}
- 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);
- gAlice->SetTrack(0, gAlice->CurrentTrack(), 20000050,
- &geant3->Gckin2()->xphot[i-1][3],geant3->Gckin2()->xphot[i-1],
- polar,geant3->Gctrak()->tofg,"Cherenkov", ncher);
- sMckov[sNphoton - 1] = ncher;
- geant3->Gctrak()->upwght = stwght;
- }
- } else {
- stwght = geant3->Gctrak()->upwght;
- geant3->Gctrak()->upwght = 0.;
- geant3->Gskpho(0);
- geant3->Gctrak()->upwght = stwght;
}
- // 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 ?
- pMC->Gmtod(geant3->Gctrak()->vect, sVloc, 1);
- pMC->Gmtod(&geant3->Gctrak()->vect[3], dir, 2);
- ysx = sVloc[2] + sDly;
- if (geant3->Gctmed()->numed == idtmed[1004-1]) {
- nmodsx = geant3->Gcvolu()->number[geant3->Gcvolu()->nlevel - 4];
- xsx = sVloc[0] + sDlx +
- xshift[geant3->Gcvolu()->number[geant3->Gcvolu()->nlevel - 2] - 1];
- } else if (geant3->Gctmed()->numed == idtmed[1003-1]) {
- nmodsx = geant3->Gcvolu()->number[geant3->Gcvolu()->nlevel - 3];
- xsx = sVloc[0] + sDlx;
- } else {
- nmodsx = 0;
- }
-
- // Momentum and direction of incidence
- psx = geant3->Gctrak()->vect[6];
- // Computing 2nd power
- r1 = dir[0];
- // Computing 2nd power
- r2 = dir[2];
- thetasx = TMath::ATan2(TMath::Sqrt(r1 * r1 + r2 * r2), dir[1]);
- phisx = TMath::ATan2(dir[2], dir[0]);
- if (phisx < 0.) phisx += 2*TMath::Pi();
-
- // Particle type and parent
- //idpartsx = geant3->Gckine()->ipart;
- r1 = geant3->Gctrak()->upwght / 100.;
- idpartgx = Int_t(r1+0.5);
- if (!geant3->Gctrak()->upwght) {
- sPparent = geant3->Gctrak()->vect[6];
- sThincParent = thetasx;
- }
- // new
- for (ll = 0; ll <nrooth; ++ll) rrhh[ll] = 0;
- irivol[0] = geant3->Gcvolu()->number[geant3->Gcvolu()->nlevel - 4];
- irivol[1] = nmodsx;
- ihitrak = gAlice->CurrentTrack();
- rrhh[0] = 1.;
- // Flag to say that this is MIP
- rrhh[1] = (Float_t) geant3->Gckine()->ipart;
- rrhh[2] = xsx;
- rrhh[3] = ysx;
- rrhh[4] = (Float_t) nmodsx;
- // Module Number
- rrhh[5] = thetasx;
- rrhh[6] = geant3->Gctrak()->tofg;
- // in seconds
- rrhh[7] = (Float_t) idpartgx;
- // mips specific
- rrhh[8] = phisx;
- rrhh[9] = psx;
- // charge of current particle in electron charge unit;
- rrhh[10] = geant3->Gckine()->charge;
- rrhh[11] = sVloc[1];
- // Zo of ckov generation
- rrhh[12] = 1.;
- // ckov generation
- AddHit(ihitrak, irivol, rrhh);
- // enew
- }
- }
-
- // Current particle is cherenkov photon
- if (geant3->Gckine()->ipart == 50) {
- pMC->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
- pMC->Gmtod(&geant3->Gctrak()->vect[3], dir, 2);
- cophi = TMath::Cos(TMath::ATan2(dir[0], dir[1]));
- t = (1. - .025 / cophi) * (1. - .05 / cophi);
- pMC->Rndm(ranf, 1);
- if (ranf[0] > t) {
- if (geant3->Gctrak()->upwght && Int_t(geant3->Gctrak()->upwght+0.5)<MAXPH)
- sIloss[Int_t(geant3->Gctrak()->upwght+0.5) - 1] = 15;
- geant3->Gctrak()->istop = 2;
- return;
- }
+ if (treeD) {
+ for (int i=0; i<kNCH; i++) {
+ sprintf(branchname,"%sDigits%d",GetName(),i+1);
+ if (fDchambers) {
+ branch = treeD->GetBranch(branchname);
+ if (branch) branch->SetAddress(&((*fDchambers)[i]));
+ }
+ }
}
-
- // 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;
+ if (treeR) {
+ for (i=0; i<kNCH; i++) {
+ sprintf(branchname,"%sRawClusters%d",GetName(),i+1);
+ if (fRawClusters) {
+ branch = treeR->GetBranch(branchname);
+ if (branch) branch->SetAddress(&((*fRawClusters)[i]));
+ }
}
- // 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 ?
- pMC->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;
- }
- 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;
-
+ for (i=0; i<kNCH; i++) {
+ sprintf(branchname,"%sRecHits%d",GetName(),i+1);
+ if (fRecHits) {
+ branch = treeR->GetBranch(branchname);
+ if (branch) branch->SetAddress(&((*fRecHits)[i]));
}
- }
}
+
}
}
+//___________________________________________
+void AliRICH::ResetHits()
+{
+ // Reset number of clusters and the cluster array for this detector
+ AliDetector::ResetHits();
+ fNPadHits = 0;
+ fNcerenkovs = 0;
+ if (fPadHits) fPadHits->Clear();
+ if (fCerenkovs) fCerenkovs->Clear();
+}
-//_____________________________________________________________________________
-void AliRICH::RichIntegration()
+//____________________________________________
+void AliRICH::ResetDigits()
{
//
- // Integrates over RICH pads
+ // Reset number of digits and the digits array for this detector
//
-
- AliMC* pMC = AliMC::GetMC();
- TGeant3 *geant3 = (TGeant3*) pMC;
-
- 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;
- }
- //
- pMC->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;
- }
- 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;
- }
- }
- // 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];
- }
- // 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
- }
+ for ( int i=0;i<kNCH;i++ ) {
+ if ((*fDchambers)[i]) (*fDchambers)[i]->Clear();
+ if (fNdch) fNdch[i]=0;
}
- }
- // 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;
- pMC->Gdtom(source, source, 1);
- FeedBack(source, qtot);
- return;
}
-//_____________________________________________________________________________
-void AliRICH::AnodicWires(Float_t &y0a)
+//____________________________________________
+void AliRICH::ResetRawClusters()
{
//
- // Position of anodic wires
+ // Reset number of raw clusters and the raw clust array for this detector
//
- 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;
+ for ( int i=0;i<kNCH;i++ ) {
+ if ((*fRawClusters)[i]) ((TClonesArray*)(*fRawClusters)[i])->Clear();
+ if (fNrawch) fNrawch[i]=0;
}
- }
-}
+}
-//_____________________________________________________________________________
-void AliRICH::GetChargeMip(Float_t &qtot)
+//____________________________________________
+void AliRICH::ResetRecHits()
{
//
- // Calculates the charge deposited by a MIP
+ // Reset number of raw clusters and the raw clust array for this detector
//
-
- AliMC* pMC = AliMC::GetMC();
-
- 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) {
- pMC->Rndm(ranf, 1);
- qtot -= fChslope * TMath::Log(ranf[0]);
- }
-}
+ for ( int i=0;i<kNCH;i++ ) {
+ if ((*fRecHits)[i]) ((TClonesArray*)(*fRecHits)[i])->Clear();
+ if (fNrechits) fNrechits[i]=0;
+ }
+}
-//_____________________________________________________________________________
-void AliRICH::GetCharge(Float_t &qtot)
+//___________________________________________
+void AliRICH::SetGeometryModel(Int_t id, AliRICHGeometry *geometry)
{
- //
- // Charge deposited
- //
- AliMC* pMC = AliMC::GetMC();
+//
+// Setter for the RICH geometry model
+//
- Float_t ranf[1];
-
- pMC->Rndm(ranf, 1);
- qtot = -fChslope * TMath::Log(ranf[0]);
-}
-//_____________________________________________________________________________
-void AliRICH::FeedBack(Float_t *source, Float_t qtot)
+ ((AliRICHChamber*) (*fChambers)[id])->GeometryModel(geometry);
+}
+
+//___________________________________________
+void AliRICH::SetSegmentationModel(Int_t id, AliRICHSegmentation *segmentation)
{
- //
- // Generate FeedBack photons
- //
- AliMC* pMC = AliMC::GetMC();
- TGeant3 *geant3 = (TGeant3*) pMC;
+//
+// Setter for the RICH segmentation model
+//
- 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
- pMC->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
- pMC->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);
- pMC->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;
-
- pMC->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;
- pMC->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;
+ ((AliRICHChamber*) (*fChambers)[id])->SegmentationModel(segmentation);
}
-//_____________________________________________________________________________
-Float_t AliRICH::FMathieson(Float_t lambda1, Float_t lambda2)
+//___________________________________________
+void AliRICH::SetResponseModel(Int_t id, AliRICHResponse *response)
{
- //
- // 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;
-
-}
+//
+// Setter for the RICH response model
+//
-//_____________________________________________________________________________
-void AliRICH::UpdateMipHit(Float_t *hits)
-{
- //
- // 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]);
+ ((AliRICHChamber*) (*fChambers)[id])->ResponseModel(response);
}
-//*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
-
-//_____________________________________________________________________________
-void AliRICH::MakeBranch(Option_t *option)
+void AliRICH::SetReconstructionModel(Int_t id, AliRICHClusterFinder *reconst)
{
- //
- // 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");
- }
- }
-}
-//_____________________________________________________________________________
-void AliRICH::SetTreeAddress()
-{
- //
- // 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);
- }
- if (fCkovs) {
- branch = treeH->GetBranch("RICHckov");
- if (branch) branch->SetAddress(&fCkovs);
- }
- if (fPadhits) {
- branch = treeH->GetBranch("RICHpadhit");
- if (branch) branch->SetAddress(&fPadhits);
- }
- }
-}
+//
+// Setter for the RICH reconstruction model (clusters)
+//
-//_____________________________________________________________________________
-void AliRICH::ResetHits()
-{
- //
- // 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();
+ ((AliRICHChamber*) (*fChambers)[id])->ReconstructionModel(reconst);
}
-//_____________________________________________________________________________
-Int_t AliRICH::DistancetoPrimitive(Int_t , Int_t )
-{
- //
- // Distance from mouse RICH on the display
- // Dummy routine
- //
- return 9999;
-}
-
-//_____________________________________________________________________________
-void AliRICH::PreTrack()
+void AliRICH::SetNsec(Int_t id, Int_t nsec)
{
- //
- // Called before transporting a track
- //
- sDetot=0;
- sNsecon=0;
- sNpads=0;
- sNphoton=0;
- sNfeed=0;
- sNfeedd=0;
- sNdir=0;
-}
-//_____________________________________________________________________________
-void AliRICH::Streamer(TBuffer &R__b)
-{
- //
- // 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;
- }
+//
+// Sets the number of padplanes
+//
+
+ ((AliRICHChamber*) (*fChambers)[id])->SetNsec(nsec);
}
-ClassImp(AliRICHv1)
-///////////////////////////////////////////////////////////////////////////////
-// //
-// Ring Imaging Cherenkov //
-// This class contains version 1 of the Ring Imaging Cherenkov //
-// //
-//Begin_Html
-/*
-<img src="picts/AliRICHv1Class.gif">
-*/
-//End_Html
-// //
-///////////////////////////////////////////////////////////////////////////////
+//___________________________________________
-//_____________________________________________________________________________
-AliRICHv1::AliRICHv1() : AliRICH()
+void AliRICH::StepManager()
{
- //
- // Default Constructor RICH for version 1
- //
+
+// Dummy step manager (should never be called)
+
}
-
-//_____________________________________________________________________________
-AliRICHv1::AliRICHv1(const char *name, const char *title)
- : AliRICH(name,title)
+
+void AliRICH::FindClusters(Int_t nev,Int_t lastEntry)
{
- //
- // Standard constructor for RICH version 1
- //
+
+//
+// Loop on chambers and on cathode planes
+//
+ for (Int_t icat=1;icat<2;icat++) {
+ gAlice->ResetDigits();
+ gAlice->TreeD()->GetEvent(1); // spurious +1 ...
+ for (Int_t ich=0;ich<kNCH;ich++) {
+ AliRICHChamber* iChamber=(AliRICHChamber*) (*fChambers)[ich];
+ TClonesArray *pRICHdigits = this->DigitsAddress(ich);
+ if (pRICHdigits == 0)
+ continue;
+ //
+ // Get ready the current chamber stuff
+ //
+ AliRICHResponse* response = iChamber->GetResponseModel();
+ AliRICHSegmentation* seg = iChamber->GetSegmentationModel();
+ AliRICHClusterFinder* rec = iChamber->GetReconstructionModel();
+ if (seg) {
+ rec->SetSegmentation(seg);
+ rec->SetResponse(response);
+ rec->SetDigits(pRICHdigits);
+ rec->SetChamber(ich);
+ if (nev==0) rec->CalibrateCOG();
+ rec->FindRawClusters();
+ }
+ TClonesArray *fRch;
+ fRch=RawClustAddress(ich);
+ fRch->Sort();
+ } // for ich
+
+ gAlice->TreeR()->Fill();
+ TClonesArray *fRch;
+ for (int i=0;i<kNCH;i++) {
+ fRch=RawClustAddress(i);
+ int nraw=fRch->GetEntriesFast();
+ printf ("Chamber %d, raw clusters %d\n",i,nraw);
+ }
+
+ ResetRawClusters();
+
+ } // for icat
+
+ char hname[30];
+ sprintf(hname,"TreeR%d",nev);
+ gAlice->TreeR()->Write(hname,kOverwrite,0);
+ gAlice->TreeR()->Reset();
+
+ //gObjectTable->Print();
}
-//_____________________________________________________________________________
-AliRICHv1::~AliRICHv1()
+
+//______________________________________________________________________________
+void AliRICH::Streamer(TBuffer &R__b)
{
- //
- // Standard destructor for RICH version 1
- //
+ // Stream an object of class AliRICH.
+ AliRICHChamber *iChamber;
+ AliRICHSegmentation *segmentation;
+ AliRICHResponse *response;
+ TClonesArray *digitsaddress;
+ TClonesArray *rawcladdress;
+ TClonesArray *rechitaddress;
+
+ if (R__b.IsReading()) {
+ Version_t R__v = R__b.ReadVersion(); if (R__v) { }
+ AliDetector::Streamer(R__b);
+ R__b >> fNPadHits;
+ R__b >> fPadHits; // diff
+ R__b >> fNcerenkovs;
+ R__b >> fCerenkovs; // diff
+ R__b >> fDchambers;
+ R__b >> fRawClusters;
+ R__b >> fRecHits; //diff
+ R__b >> fDebugLevel; //diff
+ R__b.ReadArray(fNdch);
+ R__b.ReadArray(fNrawch);
+ R__b.ReadArray(fNrechits);
+//
+ R__b >> fChambers;
+// Stream chamber related information
+ for (Int_t i =0; i<kNCH; i++) {
+ iChamber=(AliRICHChamber*) (*fChambers)[i];
+ iChamber->Streamer(R__b);
+ segmentation=iChamber->GetSegmentationModel();
+ segmentation->Streamer(R__b);
+ response=iChamber->GetResponseModel();
+ response->Streamer(R__b);
+ rawcladdress=(TClonesArray*) (*fRawClusters)[i];
+ rawcladdress->Streamer(R__b);
+ rechitaddress=(TClonesArray*) (*fRecHits)[i];
+ rechitaddress->Streamer(R__b);
+ digitsaddress=(TClonesArray*) (*fDchambers)[i];
+ digitsaddress->Streamer(R__b);
+ }
+
+ } else {
+ R__b.WriteVersion(AliRICH::IsA());
+ AliDetector::Streamer(R__b);
+ R__b << fNPadHits;
+ R__b << fPadHits; // diff
+ R__b << fNcerenkovs;
+ R__b << fCerenkovs; // diff
+ R__b << fDchambers;
+ R__b << fRawClusters;
+ R__b << fRecHits; //diff
+ R__b << fDebugLevel; //diff
+ R__b.WriteArray(fNdch, kNCH);
+ R__b.WriteArray(fNrawch, kNCH);
+ R__b.WriteArray(fNrechits, kNCH);
+ //
+ R__b << fChambers;
+// Stream chamber related information
+ for (Int_t i =0; i<kNCH; i++) {
+ iChamber=(AliRICHChamber*) (*fChambers)[i];
+ iChamber->Streamer(R__b);
+ segmentation=iChamber->GetSegmentationModel();
+ segmentation->Streamer(R__b);
+ response=iChamber->GetResponseModel();
+ response->Streamer(R__b);
+ rawcladdress=(TClonesArray*) (*fRawClusters)[i];
+ rawcladdress->Streamer(R__b);
+ rechitaddress=(TClonesArray*) (*fRecHits)[i];
+ rechitaddress->Streamer(R__b);
+ digitsaddress=(TClonesArray*) (*fDchambers)[i];
+ digitsaddress->Streamer(R__b);
+ }
+ }
}
-
-//_____________________________________________________________________________
-void AliRICHv1::CreateGeometry()
+AliRICHPadHit* AliRICH::FirstPad(AliRICHHit* hit,TClonesArray *clusters )
{
- //
- // Create the geometry for RICH version 1
- //
- // Modified by: N. Colonna (INFN - BARI, Nicola.Colonna@ba.infn.it)
- // R.A. Fini (INFN - BARI, Rosanna.Fini@ba.infn.it)
- // R.A. Loconsole (Bari University, loco@riscom.ba.infn.it)
- //
- //Begin_Html
- /*
- <img src="picts/AliRICHv1.gif">
- */
- //End_Html
- //Begin_Html
- /*
- <img src="picts/AliRICHv1Tree.gif">
- */
- //End_Html
-
- AliMC* pMC = AliMC::GetMC();
-
- 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;
- pMC->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;
- pMC->Gsvolu("SRIC", "BOX ", idtmed[1000], par, 3);
-
- // Honeycomb
- par[0] = 63.1;
- par[1] = .188;
- par[2] = 66.55;
- pMC->Gsvolu("HONE", "BOX ", idtmed[1001], par, 3);
-
- // Aluminium sheet
- par[0] = 63.1;
- par[1] = .025;
- par[2] = 66.55;
- pMC->Gsvolu("ALUM", "BOX ", idtmed[1009], par, 3);
-
- // Quartz
- par[0] = 63.1;
- par[1] = .25;
- par[2] = 65.5;
- pMC->Gsvolu("QUAR", "BOX ", idtmed[1002], par, 3);
-
- // Spacers (cylinders)
- par[0] = 0.;
- par[1] = .5;
- par[2] = .5;
- pMC->Gsvolu("SPAC", "TUBE", idtmed[1002], par, 3);
-
- // Opaque quartz
- par[0] = 61.95;
- par[1] = .2;
- par[2] = 66.5;
- pMC->Gsvolu("OQUA", "BOX ", idtmed[1007], par, 3);
-
- // Frame of opaque quartz
- par[0] = 20.65;
- par[1] = .5;
- par[2] = 66.5;
- pMC->Gsvolu("OQUF", "BOX ", idtmed[1007], par, 3);
-
- // Little bar of opaque quartz
- par[0] = 63.1;
- par[1] = .25;
- par[2] = .275;
- pMC->Gsvolu("BARR", "BOX ", idtmed[1007], par, 3);
-
- // Freon
- par[0] = 20.15;
- par[1] = .5;
- par[2] = 65.5;
- pMC->Gsvolu("FREO", "BOX ", idtmed[1003], par, 3);
-
- // Methane
- par[0] = 64.8;
- par[1] = 5.;
- par[2] = 64.8;
- pMC->Gsvolu("META", "BOX ", idtmed[1004], par, 3);
-
- // Methane gap
- par[0] = 64.8;
- par[1] = .2;
- par[2] = 64.8;
- pMC->Gsvolu("GAP ", "BOX ", idtmed[1008], par, 3);
-
- // CsI photocathode
- par[0] = 64.8;
- par[1] = .25;
- par[2] = 64.8;
- pMC->Gsvolu("CSI ", "BOX ", idtmed[1005], par, 3);
-
- // Anode grid
- par[0] = 0.;
- par[1] = .0025;
- par[2] = 20.;
- pMC->Gsvolu("GRID", "TUBE", idtmed[1006], par, 3);
-
- // --- Places the detectors defined with GSVOLU
- // Place material inside RICH
- pMC->Gspos("SRIC", 1, "RICH", 0., 0., 0., 0, "ONLY");
-
- pMC->Gspos("ALUM", 1, "SRIC", 0., -6.075, 0., 0, "ONLY");
- pMC->Gspos("HONE", 1, "SRIC", 0., -5.862, 0., 0, "ONLY");
- pMC->Gspos("ALUM", 2, "SRIC", 0., -5.649, 0., 0, "ONLY");
- pMC->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;
- pMC->Gspos("SPAC", i, "FREO", 6.7, 0., zs, idrotm[1019], "ONLY");
- }
- for (i = 10; i <= 18; ++i) {
- zs = (14 - i) * 14.4;
- pMC->Gspos("SPAC", i, "FREO", -6.7, 0., zs, idrotm[1019], "ONLY");
- }
-
- pMC->Gspos("FREO", 1, "OQUF", 0., 0., 0., 0, "ONLY");
- pMC->Gspos("OQUF", 1, "SRIC", 41.3, -4.724, 0., 0, "ONLY");
- pMC->Gspos("OQUF", 2, "SRIC", 0., -4.724, 0., 0, "ONLY");
- pMC->Gspos("OQUF", 3, "SRIC", -41.3, -4.724, 0., 0, "ONLY");
- pMC->Gspos("BARR", 1, "QUAR", 0., 0., -21.65, 0, "ONLY");
- pMC->Gspos("BARR", 2, "QUAR", 0., 0., 21.65, 0, "ONLY");
- pMC->Gspos("QUAR", 1, "SRIC", 0., -3.974, 0., 0, "ONLY");
- pMC->Gspos("GAP ", 1, "META", 0., 4.8, 0., 0, "ONLY");
- pMC->Gspos("META", 1, "SRIC", 0., 1.276, 0., 0, "ONLY");
- pMC->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.);
-
- pMC->Gspos("RICH", 1, "ALIC", 0., 471.9, 165.26, idrotm[1000], "ONLY");
- pMC->Gspos("RICH", 2, "ALIC", 171., 470., 0., idrotm[1001], "ONLY");
- pMC->Gspos("RICH", 3, "ALIC", 0., 500., 0., idrotm[1002], "ONLY");
- pMC->Gspos("RICH", 4, "ALIC", -171., 470., 0., idrotm[1003], "ONLY");
- pMC->Gspos("RICH", 5, "ALIC", 161.4, 443.4, -165.3, idrotm[1004], "ONLY");
- pMC->Gspos("RICH", 6, "ALIC", 0., 471.9, -165.3, idrotm[1005], "ONLY");
- pMC->Gspos("RICH", 7, "ALIC", -161.4, 443.4, -165.3, idrotm[1006], "ONLY");
-
+//
+ // 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=Int_t(hit->fPHlast);
+ sCurIterPad=Int_t(hit->fPHfirst);
+ return (AliRICHPadHit*) clusters->UncheckedAt(sCurIterPad-1);
+ } else {
+ return 0;
+ }
+
}
-//_____________________________________________________________________________
-void AliRICHv1::DrawModule()
+AliRICHPadHit* AliRICH::NextPad(TClonesArray *clusters)
{
- //
- // Draw a shaded view of the Ring Imaging Cherenkov
- //
- AliMC* pMC = AliMC::GetMC();
- TGeant3 *geant3 = (TGeant3*) pMC;
-
- // Set everything unseen
- pMC->Gsatt("*", "seen", -1);
- //
- // Set ALIC mother transparent
- pMC->Gsatt("ALIC","SEEN",0);
- //
- // Set the volumes visible
-
- pMC->Gsatt("RICH","seen",0);
- pMC->Gsatt("SRIC","seen",0);
- pMC->Gsatt("HONE","seen",1);
- pMC->Gsatt("ALUM","seen",1);
- pMC->Gsatt("QUAR","seen",1);
- pMC->Gsatt("SPAC","seen",1);
- pMC->Gsatt("OQUA","seen",1);
- pMC->Gsatt("OQUF","seen",1);
- pMC->Gsatt("BARR","seen",1);
- pMC->Gsatt("FREO","seen",1);
- pMC->Gsatt("META","seen",1);
- pMC->Gsatt("GAP ","seen",1);
- pMC->Gsatt("CSI ","seen",1);
- pMC->Gsatt("GRID","seen",1);
+ // Iterates over pads
- //
- 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();
- pMC->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");
+ sCurIterPad++;
+ if (sCurIterPad <= sMaxIterPad) {
+ return (AliRICHPadHit*) clusters->UncheckedAt(sCurIterPad-1);
+ } else {
+ return 0;
+ }
}
-//_____________________________________________________________________________
-void AliRICHv1::CreateMaterials()
+
+void AliRICH::Digitise(Int_t nev, Int_t flag, Option_t *option,Text_t *filename)
{
- //
- // *** 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;
-
- AliMC* pMC = AliMC::GetMC();
- TGeant3 *geant3 = (TGeant3*) pMC;
+ // keep galice.root for signal and name differently the file for
+ // background when add! otherwise the track info for signal will be lost !
- // --- 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
-
- pMC->Gstpar(idtmed[1002], "LOSS", 1.);
- pMC->Gstpar(idtmed[1003], "LOSS", 1.);
- pMC->Gstpar(idtmed[1004], "LOSS", 1.);
- pMC->Gstpar(idtmed[1008], "LOSS", 1.);
- pMC->Gstpar(idtmed[1005], "LOSS", 1.);
- pMC->Gstpar(idtmed[1002], "HADR", 1.);
- pMC->Gstpar(idtmed[1003], "HADR", 1.);
- pMC->Gstpar(idtmed[1004], "HADR", 1.);
- pMC->Gstpar(idtmed[1008], "HADR", 1.);
- pMC->Gstpar(idtmed[1005], "HADR", 1.);
- pMC->Gstpar(idtmed[1002], "DCAY", 1.);
- pMC->Gstpar(idtmed[1003], "DCAY", 1.);
- pMC->Gstpar(idtmed[1004], "DCAY", 1.);
- pMC->Gstpar(idtmed[1008], "DCAY", 1.);
- pMC->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);
-}
+ static Bool_t first=kTRUE;
+ static TFile *pFile;
+ char *addBackground = strstr(option,"Add");
-ClassImp(AliRICHhit)
-
-//_____________________________________________________________________________
-AliRICHhit::AliRICHhit(Int_t shunt, Int_t track, Int_t *vol,
- Float_t *hits, Int_t fNpadhits) :
- AliHit(shunt,track)
-{
- //
- // 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];
-}
+ FILE* points; //these will be the digits...
-ClassImp(AliRICHmip)
+ points=fopen("points.dat","w");
-//_____________________________________________________________________________
-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)
-{
- //
- // 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;
- }
+ AliRICHChamber* iChamber;
+ AliRICHSegmentation* segmentation;
+
+ Int_t digitse=0;
+ Int_t trk[50];
+ Int_t chtrk[50];
+ TObjArray *list=new TObjArray;
+ static TClonesArray *pAddress=0;
+ if(!pAddress) pAddress=new TClonesArray("TVector",1000);
+ Int_t digits[5];
+
+ AliRICH *pRICH = (AliRICH *) gAlice->GetDetector("RICH");
+ AliRICHHitMap* pHitMap[10];
+ Int_t i;
+ for (i=0; i<10; i++) {pHitMap[i]=0;}
+ if (addBackground ) {
+ if(first) {
+ fFileName=filename;
+ cout<<"filename"<<fFileName<<endl;
+ pFile=new TFile(fFileName);
+ cout<<"I have opened "<<fFileName<<" file "<<endl;
+ fHits2 = new TClonesArray("AliRICHHit",1000 );
+ fClusters2 = new TClonesArray("AliRICHPadHit",10000);
+ first=kFALSE;
+ }
+ pFile->cd();
+ // Get Hits Tree header from file
+ if(fHits2) fHits2->Clear();
+ if(fClusters2) fClusters2->Clear();
+ if(TrH1) delete TrH1;
+ TrH1=0;
+
+ char treeName[20];
+ sprintf(treeName,"TreeH%d",nev);
+ TrH1 = (TTree*)gDirectory->Get(treeName);
+ if (!TrH1) {
+ printf("ERROR: cannot find Hits Tree for event:%d\n",nev);
+ }
+ // Set branch addresses
+ TBranch *branch;
+ char branchname[20];
+ sprintf(branchname,"%s",GetName());
+ if (TrH1 && fHits2) {
+ branch = TrH1->GetBranch(branchname);
+ if (branch) branch->SetAddress(&fHits2);
+ }
+ if (TrH1 && fClusters2) {
+ branch = TrH1->GetBranch("RICHCluster");
+ if (branch) branch->SetAddress(&fClusters2);
+ }
+ }
+
+ AliRICHHitMap* hm;
+ Int_t countadr=0;
+ Int_t counter=0;
+ for (i =0; i<kNCH; i++) {
+ iChamber=(AliRICHChamber*) (*fChambers)[i];
+ segmentation=iChamber->GetSegmentationModel(1);
+ pHitMap[i] = new AliRICHHitMapA1(segmentation, list);
+ }
+ //
+ // Loop over tracks
+ //
+
+ TTree *treeH = gAlice->TreeH();
+ Int_t ntracks =(Int_t) treeH->GetEntries();
+ for (Int_t track=0; track<ntracks; track++) {
+ gAlice->ResetHits();
+ treeH->GetEvent(track);
+ //
+ // Loop over hits
+ for(AliRICHHit* mHit=(AliRICHHit*)pRICH->FirstHit(-1);
+ mHit;
+ mHit=(AliRICHHit*)pRICH->NextHit())
+ {
+
+ digitse=0;
+
+ Int_t nch = mHit->fChamber-1; // chamber number
+ if (nch >kNCH) continue;
+ iChamber = &(pRICH->Chamber(nch));
+
+ TParticle *current = (TParticle*)(*gAlice->Particles())[track];
+
+ Int_t particle = current->GetPdgCode();
+
+ //printf("Flag:%d\n",flag);
+ //printf("Track:%d\n",track);
+ //printf("Particle:%d\n",particle);
+
+ if (flag == 0)
+ digitse=1;
+
+ if (flag == 1)
+ if(TMath::Abs(particle) == 211 || TMath::Abs(particle) == 111)
+ digitse=1;
+
+ if (flag == 2)
+ if(TMath::Abs(particle)==321 || TMath::Abs(particle)==130 || TMath::Abs(particle)==310
+ || TMath::Abs(particle)==311)
+ digitse=1;
+
+ if (flag == 3 && TMath::Abs(particle)==2212)
+ digitse=1;
+
+ if (flag == 4 && TMath::Abs(particle)==13)
+ digitse=1;
+
+ if (flag == 5 && TMath::Abs(particle)==11)
+ digitse=1;
+
+ if (flag == 6 && TMath::Abs(particle)==2112)
+ digitse=1;
+
+
+ //printf ("Particle: %d, Flag: %d, Digitse: %d\n",particle,flag,digitse);
+
+
+ if (digitse)
+ {
+
+ //
+ // Loop over pad hits
+ for (AliRICHPadHit* mPad=
+ (AliRICHPadHit*)pRICH->FirstPad(mHit,fPadHits);
+ mPad;
+ mPad=(AliRICHPadHit*)pRICH->NextPad(fPadHits))
+ {
+ 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
+ //
+ //
+ //printf("X:%d, Y:%d, Q:%d\n",ipx,ipy,iqpad);
+
+ Float_t thex, they;
+ segmentation=iChamber->GetSegmentationModel(cathode);
+ segmentation->GetPadCxy(ipx,ipy,thex,they);
+ new((*pAddress)[countadr++]) TVector(2);
+ TVector &trinfo=*((TVector*) (*pAddress)[countadr-1]);
+ trinfo(0)=(Float_t)track;
+ trinfo(1)=(Float_t)iqpad;
+
+ digits[0]=ipx;
+ digits[1]=ipy;
+ digits[2]=iqpad;
+
+ AliRICHTransientDigit* pdigit;
+ // build the list of fired pads and update the info
+ if (!pHitMap[nch]->TestHit(ipx, ipy)) {
+ list->AddAtAndExpand(new AliRICHTransientDigit(nch,digits),counter);
+ pHitMap[nch]->SetHit(ipx, ipy, counter);
+ counter++;
+ pdigit=(AliRICHTransientDigit*)list->At(list->GetLast());
+ // list of tracks
+ TObjArray *trlist=(TObjArray*)pdigit->TrackList();
+ trlist->Add(&trinfo);
+ } else {
+ pdigit=(AliRICHTransientDigit*) pHitMap[nch]->GetHit(ipx, ipy);
+ // update charge
+ (*pdigit).fSignal+=iqpad;
+ // update list of tracks
+ TObjArray* trlist=(TObjArray*)pdigit->TrackList();
+ Int_t lastEntry=trlist->GetLast();
+ TVector *ptrkP=(TVector*)trlist->At(lastEntry);
+ TVector &ptrk=*ptrkP;
+ Int_t lastTrack=Int_t(ptrk(0));
+ Int_t lastCharge=Int_t(ptrk(1));
+ if (lastTrack==track) {
+ lastCharge+=iqpad;
+ trlist->RemoveAt(lastEntry);
+ trinfo(0)=lastTrack;
+ trinfo(1)=lastCharge;
+ trlist->AddAt(&trinfo,lastEntry);
+ } else {
+ trlist->Add(&trinfo);
+ }
+ // check the track list
+ Int_t nptracks=trlist->GetEntriesFast();
+ if (nptracks > 2) {
+ printf("Attention - tracks: %d (>2)\n",nptracks);
+ //printf("cat,nch,ix,iy %d %d %d %d \n",icat+1,nch,ipx,ipy);
+ for (Int_t tr=0;tr<nptracks;tr++) {
+ TVector *pptrkP=(TVector*)trlist->At(tr);
+ TVector &pptrk=*pptrkP;
+ trk[tr]=Int_t(pptrk(0));
+ chtrk[tr]=Int_t(pptrk(1));
+ }
+ } // end if nptracks
+ } // end if pdigit
+ } //end loop over clusters
+ }// track type condition
+ } // hit loop
+ } // track loop
+
+ // open the file with background
+
+ if (addBackground ) {
+ ntracks =(Int_t)TrH1->GetEntries();
+ //printf("background - icat,ntracks1 %d %d\n",icat,ntracks);
+ //printf("background - Start loop over tracks \n");
+ //
+ // Loop over tracks
+ //
+ for (Int_t trak=0; trak<ntracks; trak++) {
+ if (fHits2) fHits2->Clear();
+ if (fClusters2) fClusters2->Clear();
+ TrH1->GetEvent(trak);
+ //
+ // Loop over hits
+ AliRICHHit* mHit;
+ for(int j=0;j<fHits2->GetEntriesFast();++j)
+ {
+ mHit=(AliRICHHit*) (*fHits2)[j];
+ Int_t nch = mHit->fChamber-1; // chamber number
+ if (nch >6) continue;
+ iChamber = &(pRICH->Chamber(nch));
+ Int_t rmin = (Int_t)iChamber->RInner();
+ Int_t rmax = (Int_t)iChamber->ROuter();
+ //
+ // Loop over pad hits
+ for (AliRICHPadHit* mPad=
+ (AliRICHPadHit*)pRICH->FirstPad(mHit,fClusters2);
+ mPad;
+ mPad=(AliRICHPadHit*)pRICH->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
+
+ 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;
+ new((*pAddress)[countadr++]) TVector(2);
+ TVector &trinfo=*((TVector*) (*pAddress)[countadr-1]);
+ trinfo(0)=-1; // tag background
+ trinfo(1)=-1;
+ digits[0]=ipx;
+ digits[1]=ipy;
+ digits[2]=iqpad;
+ if (trak <4 && nch==0)
+ printf("bgr - pHitMap[nch]->TestHit(ipx, ipy),trak %d %d\n",
+ pHitMap[nch]->TestHit(ipx, ipy),trak);
+ AliRICHTransientDigit* pdigit;
+ // build the list of fired pads and update the info
+ if (!pHitMap[nch]->TestHit(ipx, ipy)) {
+ list->AddAtAndExpand(new AliRICHTransientDigit(nch,digits),counter);
+
+ pHitMap[nch]->SetHit(ipx, ipy, counter);
+ counter++;
+ printf("bgr new elem in list - counter %d\n",counter);
+
+ pdigit=(AliRICHTransientDigit*)list->At(list->GetLast());
+ // list of tracks
+ TObjArray *trlist=(TObjArray*)pdigit->TrackList();
+ trlist->Add(&trinfo);
+ } else {
+ pdigit=(AliRICHTransientDigit*) pHitMap[nch]->GetHit(ipx, ipy);
+ // update charge
+ (*pdigit).fSignal+=iqpad;
+ // update list of tracks
+ TObjArray* trlist=(TObjArray*)pdigit->TrackList();
+ Int_t lastEntry=trlist->GetLast();
+ TVector *ptrkP=(TVector*)trlist->At(lastEntry);
+ TVector &ptrk=*ptrkP;
+ Int_t lastTrack=Int_t(ptrk(0));
+ if (lastTrack==-1) {
+ continue;
+ } else {
+ trlist->Add(&trinfo);
+ }
+ // check the track list
+ Int_t nptracks=trlist->GetEntriesFast();
+ if (nptracks > 0) {
+ for (Int_t tr=0;tr<nptracks;tr++) {
+ TVector *pptrkP=(TVector*)trlist->At(tr);
+ TVector &pptrk=*pptrkP;
+ 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
+ TTree *fAli=gAlice->TreeK();
+ if (fAli) pFile =fAli->GetCurrentFile();
+ pFile->cd();
+ } // if Add
+
+ Int_t tracks[10];
+ Int_t charges[10];
+ //cout<<"Start filling digits \n "<<endl;
+ Int_t nentries=list->GetEntriesFast();
+ //printf(" \n \n nentries %d \n",nentries);
+
+ // start filling the digits
+
+ for (Int_t nent=0;nent<nentries;nent++) {
+ AliRICHTransientDigit *address=(AliRICHTransientDigit*)list->At(nent);
+ if (address==0) continue;
+
+ Int_t ich=address->fChamber;
+ Int_t q=address->fSignal;
+ iChamber=(AliRICHChamber*) (*fChambers)[ich];
+ AliRICHResponse * response=iChamber->GetResponseModel();
+ Int_t adcmax= (Int_t) response->MaxAdc();
+
+
+ // add white noise and do zero-suppression and signal truncation (new electronics,old electronics gaus 1.2,0.2)
+ //printf("Treshold: %d\n",iChamber->fTresh->GetHitIndex(address->fPadX,address->fPadY));
+ Int_t pedestal = iChamber->fTresh->GetHitIndex(address->fPadX,address->fPadY);
+
+ //printf("Pedestal:%d\n",pedestal);
+ //Int_t pedestal=0;
+ Float_t treshold = (pedestal + 4*1.7);
+
+ Float_t meanNoise = gRandom->Gaus(1.7, 0.25);
+ Float_t noise = gRandom->Gaus(0, meanNoise);
+ q+=(Int_t)(noise + pedestal);
+ //q+=(Int_t)(noise);
+ // magic number to be parametrised !!!
+ if ( q <= treshold)
+ {
+ q = q - pedestal;
+ continue;
+ }
+ q = q - pedestal;
+ if ( q >= adcmax) q=adcmax;
+ 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 - tracks > 10 "<<nptracks<<endl;
+ nptracks=10;
+ }
+ if (nptracks > 2) {
+ printf("Attention - tracks > 2 %d \n",nptracks);
+ //printf("cat,ich,ix,iy,q %d %d %d %d %d \n",
+ //icat,ich,digits[0],digits[1],q);
+ }
+ for (Int_t tr=0;tr<nptracks;tr++) {
+ TVector *ppP=(TVector*)trlist->At(tr);
+ TVector &pp =*ppP;
+ 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 t=nptracks; t<10; t++) {
+ tracks[t]=0;
+ charges[t]=0;
+ }
+ }
+ //write file
+ if (ich==2)
+ fprintf(points,"%4d, %4d, %4d\n",digits[0],digits[1],digits[2]);
+
+ // fill digits
+ pRICH->AddDigits(ich,tracks,charges,digits);
+ }
+ gAlice->TreeD()->Fill();
+
+ list->Delete();
+ for(Int_t ii=0;ii<kNCH;++ii) {
+ if (pHitMap[ii]) {
+ hm=pHitMap[ii];
+ delete hm;
+ pHitMap[ii]=0;
+ }
+ }
+
+ //TTree *TD=gAlice->TreeD();
+ //Stat_t ndig=TD->GetEntries();
+ //cout<<"number of digits "<<ndig<<endl;
+ TClonesArray *fDch;
+ for (int k=0;k<kNCH;k++) {
+ fDch= pRICH->DigitsAddress(k);
+ int ndigit=fDch->GetEntriesFast();
+ printf ("Chamber %d digits %d \n",k,ndigit);
+ }
+ pRICH->ResetDigits();
+ char hname[30];
+ sprintf(hname,"TreeD%d",nev);
+ gAlice->TreeD()->Write(hname,kOverwrite,0);
+
+ // reset tree
+ // gAlice->TreeD()->Reset();
+ delete list;
+ pAddress->Clear();
+ // gObjectTable->Print();
}
-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)
+AliRICH& AliRICH::operator=(const AliRICH& rhs)
{
- //
- // Standard creator for RICH Cherenkov information
- //
- fEnergy = hits[8];
- fStop = (int) hits[9];
-
- //Parent info
- fParent = fNmips;
+// Assignment operator
+ return *this;
+
}
-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)
+Int_t AliRICH::MakePadHits(Float_t xhit,Float_t yhit,Float_t eloss, Int_t idvol, ResponseType res)
{
- //
- // 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
+//
+// Calls the charge disintegration method of the current chamber and adds
+// the simulated cluster to the root treee
+//
+ Int_t clhits[kNCH];
+ Float_t newclust[6][500];
+ Int_t nnew;
+
+//
+// Integrated pulse height on chamber
+
+ clhits[0]=fNhits+1;
+
+ ((AliRICHChamber*) (*fChambers)[idvol])->DisIntegration(eloss, xhit, yhit, nnew, newclust, res);
+ Int_t ic=0;
+
+//
+// Add new clusters
+ for (Int_t i=0; i<nnew; i++) {
+ if (Int_t(newclust[3][i]) > 0) {
+ ic++;
+// Cathode plane
+ clhits[1] = Int_t(newclust[5][i]);
+// Cluster Charge
+ clhits[2] = Int_t(newclust[0][i]);
+// Pad: ix
+ clhits[3] = Int_t(newclust[1][i]);
+// Pad: iy
+ clhits[4] = Int_t(newclust[2][i]);
+// Pad: charge
+ clhits[5] = Int_t(newclust[3][i]);
+// Pad: chamber sector
+ clhits[6] = Int_t(newclust[4][i]);
+
+ AddPadHit(clhits);
+ }
+ }
+return nnew;
}
-
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff