From d15a28e7baaf9ecf573af1e229c75b9527e05793 Mon Sep 17 00:00:00 2001 From: schutz Date: Fri, 17 Dec 1999 16:27:45 +0000 Subject: [PATCH] ---------------------------------------------------------------------- Modified files AliPHOS.cxx AliPHOS.h AliPHOSv0.cxx AliPHOSv0.h Makefile PHOSLinkDef.h Added Files: AliPHOSClusterizer.cxx AliPHOSClusterizer.h AliPHOSClusterizerv1.cxx AliPHOSClusterizerv1.h AliPHOSDigit.cxx AliPHOSDigit.h AliPHOSEmcRecPoint.cxx AliPHOSEmcRecPoint.h AliPHOSGeometry.cxx AliPHOSGeometry.h AliPHOSLink.cxx AliPHOSLink.h AliPHOSPpsdRecPoint.cxx AliPHOSPpsdRecPoint.h AliPHOSRecPoint.cxx AliPHOSRecPoint.h AliPHOSReconstructioner.cxx AliPHOSReconstructioner.h AliPHOSTrackSegment.cxx AliPHOSTrackSegment.h AliPHOSTrackSegmentMaker.cxx AliPHOSTrackSegmentMaker.h AliPHOSTrackSegmentMakerv1.cxx AliPHOSTrackSegmentMakerv1.h canevas.cxx canevas.h testTrackSegment.C Removed Files: Tag: HEAD AliPHOSv1.cxx AliPHOSv1.h AliPHOSv2.cxx AliPHOSv2.h --------------------------------------------------------------------- This is a major revision of the PHOS package undertaken at SUBATECH. It is not anymore compatible with version V3.02. There is an unique but parametrized geometry. The reconstruction is operational until the track segment reconstruction. Use testTrackSegment.C for tests Y. Schutz --- G. Martinez --- D. Peressounko --- PHOS/AliPHOS.cxx | 1246 ++++--------------------- PHOS/AliPHOS.h | 247 +---- PHOS/AliPHOSClusterizer.cxx | 45 + PHOS/AliPHOSClusterizer.h | 43 + PHOS/AliPHOSClusterizerv1.cxx | 252 +++++ PHOS/AliPHOSClusterizerv1.h | 69 ++ PHOS/AliPHOSDigit.cxx | 86 ++ PHOS/AliPHOSDigit.h | 52 ++ PHOS/AliPHOSEmcRecPoint.cxx | 418 +++++++++ PHOS/AliPHOSEmcRecPoint.h | 73 ++ PHOS/AliPHOSGeometry.cxx | 438 +++++++++ PHOS/AliPHOSGeometry.h | 180 ++++ PHOS/AliPHOSLink.cxx | 51 + PHOS/AliPHOSLink.h | 47 + PHOS/AliPHOSPpsdRecPoint.cxx | 258 +++++ PHOS/AliPHOSPpsdRecPoint.h | 55 ++ PHOS/AliPHOSRecPoint.cxx | 61 ++ PHOS/AliPHOSRecPoint.h | 50 + PHOS/AliPHOSReconstructioner.cxx | 60 ++ PHOS/AliPHOSReconstructioner.h | 50 + PHOS/AliPHOSTrackSegment.cxx | 182 ++++ PHOS/AliPHOSTrackSegment.h | 62 ++ PHOS/AliPHOSTrackSegmentMaker.cxx | 53 ++ PHOS/AliPHOSTrackSegmentMaker.h | 48 + PHOS/AliPHOSTrackSegmentMakerv1.cxx | 520 +++++++++++ PHOS/AliPHOSTrackSegmentMakerv1.h | 73 ++ PHOS/AliPHOSv0.cxx | 1346 ++++++++++++++++++++++----- PHOS/AliPHOSv0.h | 112 ++- PHOS/Makefile | 23 +- PHOS/PHOSLinkDef.h | 30 +- PHOS/canevas.cxx | 44 + PHOS/canevas.h | 31 + PHOS/testTrackSegment.C | 95 ++ 33 files changed, 4857 insertions(+), 1543 deletions(-) create mode 100644 PHOS/AliPHOSClusterizer.cxx create mode 100644 PHOS/AliPHOSClusterizer.h create mode 100644 PHOS/AliPHOSClusterizerv1.cxx create mode 100644 PHOS/AliPHOSClusterizerv1.h create mode 100644 PHOS/AliPHOSDigit.cxx create mode 100644 PHOS/AliPHOSDigit.h create mode 100644 PHOS/AliPHOSEmcRecPoint.cxx create mode 100644 PHOS/AliPHOSEmcRecPoint.h create mode 100644 PHOS/AliPHOSGeometry.cxx create mode 100644 PHOS/AliPHOSGeometry.h create mode 100644 PHOS/AliPHOSLink.cxx create mode 100644 PHOS/AliPHOSLink.h create mode 100644 PHOS/AliPHOSPpsdRecPoint.cxx create mode 100644 PHOS/AliPHOSPpsdRecPoint.h create mode 100644 PHOS/AliPHOSRecPoint.cxx create mode 100644 PHOS/AliPHOSRecPoint.h create mode 100644 PHOS/AliPHOSReconstructioner.cxx create mode 100644 PHOS/AliPHOSReconstructioner.h create mode 100644 PHOS/AliPHOSTrackSegment.cxx create mode 100644 PHOS/AliPHOSTrackSegment.h create mode 100644 PHOS/AliPHOSTrackSegmentMaker.cxx create mode 100644 PHOS/AliPHOSTrackSegmentMaker.h create mode 100644 PHOS/AliPHOSTrackSegmentMakerv1.cxx create mode 100644 PHOS/AliPHOSTrackSegmentMakerv1.h create mode 100644 PHOS/canevas.cxx create mode 100644 PHOS/canevas.h create mode 100644 PHOS/testTrackSegment.C diff --git a/PHOS/AliPHOS.cxx b/PHOS/AliPHOS.cxx index 443e1ddc9e5..da4f679e2cf 100644 --- a/PHOS/AliPHOS.cxx +++ b/PHOS/AliPHOS.cxx @@ -13,1137 +13,277 @@ * provided "as is" without express or implied warranty. * **************************************************************************/ -/* -$Log$ -Revision 1.10 1999/12/09 14:57:51 fca -Removal of obsolete PHOS code +//_________________________________________________________________________ +// Base Class of PHOS. +// Only creates the materials +//*-- Author : Laurent Aphecetche SUBATECH +////////////////////////////////////////////////////////////////////////////// -Revision 1.9 1999/11/08 07:12:31 fca -Minor corrections thanks to I.Hrivnacova - -Revision 1.8 1999/09/29 09:24:23 fca -Introduction of the Copyright and cvs Log - -*/ - -//////////////////////////////////////////////// -// Manager and hits classes for set:PHOS // -//////////////////////////////////////////////// - // --- ROOT system --- -#include "TH1.h" -#include "TRandom.h" -#include "TFile.h" -#include "TTree.h" -#include "TBRIK.h" -#include "TNode.h" -#include "TMath.h" // --- Standard library --- -#include -#include -#include -#include -// --- galice header files --- +// --- AliRoot header files --- + #include "AliPHOS.h" +#include "AliMC.h" #include "AliRun.h" -//______________________________________________________________________________ - - ClassImp(AliPHOS) -//______________________________________________________________________________ - -AliPHOS::~AliPHOS(void) +//____________________________________________________________________________ +AliPHOS::AliPHOS(const char* name, const char* title) + : AliDetector(name,title) { - delete fHits; // 28.12.1998 - delete fTreePHOS; // 28.12.1998 - fCradles->Delete(); - delete fCradles; } -//______________________________________________________________________________ - -AliPHOS::AliPHOS() : - fDebugLevel (0), - fTreePHOS (NULL), - fBranchNameOfCradles ("AliPHOSCradles"), - fTreeName ("PHOS") +//____________________________________________________________________________ +AliPHOS::AliPHOS() : AliDetector() { - fIshunt = 0; - - if( NULL==(fCradles=new TObjArray) ) - { - Error("AliPHOS","Can not create fCradles"); - exit(1); - } - DefPars(); } - -//______________________________________________________________________________ - -AliPHOS::AliPHOS(const char *name, const char *title) - : AliDetector (name,title), - fDebugLevel (0), - fTreePHOS (NULL), - fBranchNameOfCradles ("AliPHOSCradles"), - fTreeName ("PHOS") -{ -//Begin_Html -/* - -*/ -//End_Html - - fHits = new TClonesArray("AliPHOShit", 405); - - fIshunt = 0; - - SetMarkerColor(kGreen); - SetMarkerStyle(2); - SetMarkerSize(0.4); - if( NULL==(fCradles=new TObjArray) ) { - Error("AliPHOS","Can not create fCradles"); - exit(1); - } - DefPars(); -} - -//______________________________________________________________________________ - -void AliPHOS::DefPars() -{ - PHOSflags[0]=0; - PHOSflags[1]=1; - PHOSflags[2]=0; - PHOSflags[3]=0; - PHOSflags[4]=0; - PHOSflags[5]=0; - PHOSflags[6]=0; - PHOSflags[7]=0; - PHOSflags[8]=0; - PHOScell[0]=2.2; - PHOScell[1]=18.; - PHOScell[2]=0.01; - PHOScell[3]=0.01; - PHOScell[4]=1.0; - PHOScell[5]=0.1; - PHOScell[6]=0.; - PHOScell[7]=0.; - PHOScell[8]=0.; - PHOSradius=460.; - PHOSsize[0]=104; - PHOSsize[1]=88; - PHOSsize[2]=4; - PHOScradlesA=0.; - PHOSextra[0]=0.001; - PHOSextra[1]=6.95; - PHOSextra[2]=4.; - PHOSextra[3]=5.; - PHOSextra[4]=2.; - PHOSextra[5]=0.06; - PHOSextra[6]=10.; - PHOSextra[7]=3.; - PHOSextra[8]=1.; - PHOSTXW[0]=209.; - PHOSTXW[1]=71.; - PHOSTXW[2]=250.; - PHOSAIR[0]=206.; - PHOSAIR[1]=66.; - PHOSAIR[2]=244.; - PHOSFTI[0]=214.6; - PHOSFTI[1]=80.; - PHOSFTI[2]=260.; - PHOSFTI[3]=467.; -} -//______________________________________________________________________________ - -void AliPHOS::AddHit(Int_t track, Int_t *vol, Float_t *hits) +//____________________________________________________________________________ +AliPHOS::~AliPHOS() { - TClonesArray &lhits = *fHits; - new(lhits[fNhits++]) AliPHOShit(fIshunt,track,vol,hits); + delete fHits ; + delete fDigits ; } - -//___________________________________________ -void AliPHOS::BuildGeometry() -{ - TNode *Node, *Top; - - const int kColorPHOS = kRed; - // - Top=gAlice->GetGeometry()->GetNode("alice"); - - - // PHOS - Float_t pphi=12.9399462; - new TRotMatrix("rot988","rot988",90,-3*pphi,90,90-3*pphi,0,0); - new TRotMatrix("rot989","rot989",90,- pphi,90,90- pphi,0,0); - new TRotMatrix("rot990","rot990",90, pphi,90,90+ pphi,0,0); - new TRotMatrix("rot991","rot991",90, 3*pphi,90,90+3*pphi,0,0); - new TBRIK("S_PHOS","PHOS box","void",107.3,40,130); - Top->cd(); - Node = new TNode("PHOS1","PHOS1","S_PHOS",-317.824921,-395.014343,0,"rot988"); - Node->SetLineColor(kColorPHOS); - fNodes->Add(Node); - Top->cd(); - Node = new TNode("PHOS2","PHOS2","S_PHOS",-113.532333,-494.124908,0,"rot989"); - fNodes->Add(Node); - Node->SetLineColor(kColorPHOS); - Top->cd(); - Node = new TNode("PHOS3","PHOS3","S_PHOS", 113.532333,-494.124908,0,"rot990"); - Node->SetLineColor(kColorPHOS); - fNodes->Add(Node); - Top->cd(); - Node = new TNode("PHOS4","PHOS4","S_PHOS", 317.824921,-395.014343,0,"rot991"); - Node->SetLineColor(kColorPHOS); - fNodes->Add(Node); -} - -//___________________________________________ +//____________________________________________________________________________ void AliPHOS::CreateMaterials() { -// *** DEFINITION OF AVAILABLE PHOS MATERIALS *** - -// CALLED BY : PHOS_MEDIA -// ORIGIN : NICK VAN EIJNDHOVEN - - - - Int_t ISXFLD = gAlice->Field()->Integ(); - Float_t SXMGMX = gAlice->Field()->Max(); - -// --- The PbWO4 crystals --- - Float_t ax[3] = { 207.19,183.85,16. }; - Float_t zx[3] = { 82.,74.,8. }; - Float_t wx[3] = { 1.,1.,4. }; - Float_t dx = 8.28; -// --- Stainless Steel --- - Float_t as[5] = { 55.847,12.011,51.9961,58.69,28.0855 }; - Float_t zs[5] = { 26.,6.,24.,28.,14. }; - Float_t ws[5] = { .6392,8e-4,.2,.14,.02 }; - Float_t ds = 8.; -// --- The polysterene scintillator (CH) --- - Float_t ap[2] = { 12.011,1.00794 }; - Float_t zp[2] = { 6.,1. }; - Float_t wp[2] = { 1.,1. }; - Float_t dp = 1.032; -// --- Tyvek (CnH2n) - Float_t at[2] = { 12.011,1.00794 }; - Float_t zt[2] = { 6.,1. }; - Float_t wt[2] = { 1.,2. }; - Float_t dt = .331; -// --- Polystyrene foam --- - Float_t af[2] = { 12.011,1.00794 }; - Float_t zf[2] = { 6.,1. }; - Float_t wf[2] = { 1.,1. }; - Float_t df = .12; -//--- Foam thermo insulation (actual chemical composition unknown yet!) --- - Float_t ati[2] = { 12.011,1.00794 }; - Float_t zti[2] = { 6.,1. }; - Float_t wti[2] = { 1.,1. }; - Float_t dti = .1; -// --- Textolit (actual chemical composition unknown yet!) --- - Float_t atx[2] = { 12.011,1.00794 }; - Float_t ztx[2] = { 6.,1. }; - Float_t wtx[2] = { 1.,1. }; - Float_t dtx = 1.83; - - Int_t *idtmed = fIdtmed->GetArray()-699; - - AliMixture( 0, "PbWO4$", ax, zx, dx, -3, wx); - AliMixture( 1, "Polystyrene$", ap, zp, dp, -2, wp); - AliMaterial( 2, "Al$", 26.98, 13., 2.7, 8.9, 999); -// --- Absorption length^ is ignored --- - AliMixture( 3, "Tyvek$", at, zt, dt, -2, wt); - AliMixture( 4, "Foam$", af, zf, df, -2, wf); - AliMixture( 5, "Stainless Steel$", as, zs, ds, 5, ws); - AliMaterial( 6, "Si$", 28.09, 14., 2.33, 9.36, 42.3); - AliMixture( 7, "Thermo Insul.$", ati, zti, dti, -2, wti); - AliMixture( 8, "Textolit$", atx, ztx, dtx, -2, wtx); - AliMaterial(99, "Air$", 14.61, 7.3, .001205, 30420., 67500); - - AliMedium(0, "PHOS Xtal $", 0, 1, ISXFLD, SXMGMX, 10., .1, .1, .1, .1); - AliMedium(2, "Al parts $", 2, 0, ISXFLD, SXMGMX, 10., .1, .1, .001, .001); - AliMedium(3, "Tyvek wrapper$", 3, 0, ISXFLD, SXMGMX, 10., .1, .1, .001, .001); - AliMedium(4, "Polyst. foam $", 4, 0, ISXFLD, SXMGMX, 10., .1, .1, .1, .1); - AliMedium(5, "Steel cover $", 5, 0, ISXFLD, SXMGMX, 10., .1, .1, 1e-4, 1e-4); - AliMedium(6, "Si PIN $", 6, 0, ISXFLD, SXMGMX, 10., .1, .1, .01, .01); - AliMedium(7, "Thermo Insul.$", 7, 0, ISXFLD, SXMGMX, 10., .1, .1, .1, .1); - AliMedium(8, "Textolit $", 8, 0, ISXFLD, SXMGMX, 10., .1, .1, .1, .1); - AliMedium(99, "Air $",99, 0, ISXFLD, SXMGMX, 10., 1., .1, .1, 10); - -// --- Generate explicitly delta rays in the steel cover --- - gMC->Gstpar(idtmed[704], "LOSS", 3.); - gMC->Gstpar(idtmed[704], "DRAY", 1.); -// --- and in aluminium parts --- - gMC->Gstpar(idtmed[701], "LOSS", 3.); - gMC->Gstpar(idtmed[701], "DRAY", 1.); -} - -//______________________________________________________________________________ - -void AliPHOS::AddPHOSCradles() -{ - Int_t i; - for(i=0;iGetEntries(); - fCradles->Add(new AliPHOSCradle( IsVersion(), // geometry. - GetCrystalSideSize (), - GetCrystalLength (), - GetWrapThickness (), - GetAirThickness (), - GetPIN_SideSize (), - GetPIN_Length (), - GetRadius (), - GetNz (), - GetNphi (), - GetCradleAngle (i))); - - if( n+1 != fCradles->GetEntries() || NULL == fCradles->At(n) ) - { - cout << " Can not create or add AliPHOSCradle.\n"; - exit(1); - } - } -} - -//______________________________________________________________________________ - -Int_t AliPHOS::DistancetoPrimitive(Int_t , Int_t ) -{ - return 9999; -} - -//___________________________________________ -void AliPHOS::Init() -{ - Int_t i; - // - printf("\n"); - for(i=0;i<35;i++) printf("*"); - printf(" PHOS_INIT "); - for(i=0;i<35;i++) printf("*"); - printf("\n"); - // - // Here the ABSO initialisation code (if any!) - for(i=0;i<80;i++) printf("*"); - printf("\n"); -} + // DEFINITION OF PHOS MATERIALS -//______________________________________________________________________________ + // --- The PbWO4 crystals --- + Float_t AX[3] = {207.19, 183.85, 16.0} ; + Float_t ZX[3] = {82.0, 74.0, 8.0} ; + Float_t WX[3] = {1.0, 1.0, 4.0} ; + Float_t DX = 8.28 ; -void AliPHOS::MakeBranch(Option_t *) -{ -// ROOT output initialization to ROOT file. -// -// AliDetector::MakeBranch() is always called. -// -// There will be also special tree "PHOS" with one branch "AliPHOSCradles" -// if it was set next flag in the galice card file: -// * PHOSflags: YES: X<>0 NO: X=0 -// * PHOSflags(1) : -----X. Create branch for TObjArray of AliPHOSCradle -// Examples: -// PHOSflags 1. -// PHOSflags 636301. -// In that case special bit CradlesBranch_Bit will be set for AliPHOS - - AliDetector::MakeBranch(); - - int i; - float t = GetPHOS_flag(0)/10; - i = (int) t; - i = (int) ((t-i)*10); - if( !i ) - return; - - SetBit(CradlesBranch_Bit); - - if( NULL==(fTreePHOS=new TTree(fTreeName.Data(),"PHOS events tree")) ) - { - Error("MakeBranch","Can not create TTree"); - exit(1); - } - - if( NULL==fTreePHOS->GetCurrentFile() ) - { - Error("MakeBranch","There is no opened ROOT file"); - exit(1); - } - - // Create a new branch in the current Root Tree. - - if( NULL==fTreePHOS->Branch(fBranchNameOfCradles.Data(),"TObjArray",&fCradles,4000,0) ) - { - Error("MakeBranch","Can not create branch"); - exit(1); - } - - printf("The branch %s has been created\n",fBranchNameOfCradles.Data()); -} + AliMixture(0, "PbWO4$", AX, ZX, DX, -3, WX) ; -//______________________________________________________________________________ -void AliPHOS::SetTreeAddress(void) -{ -// ROOT input initialization. -// -// AliDetector::SetTreeAddress() is always called. -// -// If CradlesBranch_Bit is set (see AliPHOS::MakeBranch) than fTreePHOS is -// initilized. - - AliDetector::SetTreeAddress(); - - if( !TestBit(CradlesBranch_Bit) ) - return; - - if( NULL==(fTreePHOS=(TTree*)gDirectory->Get((char*)(fTreeName.Data())) ) ) - { - Error("SetTreeAddress","Can not find Tree \"%s\"\n",fTreeName.Data()); - exit(1); - } - - TBranch *branch = fTreePHOS->GetBranch(fBranchNameOfCradles.Data()); - if( NULL==branch ) - { - Error("SetTreeAddress","Can not find branch %s in TTree:%s",fBranchNameOfCradles.Data(),fTreeName.Data()); - exit(1); - } - - branch->SetAddress(&fCradles); -} + // --- The polysterene scintillator (CH) --- + Float_t AP[2] = {12.011, 1.00794} ; + Float_t ZP[2] = {6.0, 1.0} ; + Float_t WP[2] = {1.0, 1.0} ; + Float_t DP = 1.032 ; -//______________________________________________________________________________ + AliMixture(1, "Polystyrene$", AP, ZP, DP, -2, WP) ; -AliPHOSCradle *AliPHOS::GetCradleOfTheParticle(const TVector3 &p,const TVector3 &v) const -{ -// For a given direction 'p' and source point 'v' returns pointer to AliPHOSCradle -// in that direction or NULL if AliPHOSCradle was not found. + // --- Aluminium --- + AliMaterial(2, "Al$", 26.98, 13., 2.7, 8.9, 999., 0, 0) ; + // --- Absorption length is ignored ^ - for( int m=0; mGetEntries(); m++ ) - { - AliPHOS *PHOS = (AliPHOS *)this; // Removing 'const'... - AliPHOSCradle *cradle = (AliPHOSCradle *)PHOS->fCradles->operator[](m); + // --- Tyvek (CnH2n) --- + Float_t AT[2] = {12.011, 1.00794} ; + Float_t ZT[2] = {6.0, 1.0} ; + Float_t WT[2] = {1.0, 2.0} ; + Float_t DT = 0.331 ; - float x,y,l; - const float d = cradle->GetRadius(); - cradle->GetXY(p,v,d,x,y,l); + AliMixture(3, "Tyvek$", AT, ZT, DT, -2, WT) ; - if( l>0 && TMath::Abs(x)GetNz ()*cradle->GetCellSideSize()/2 - && TMath::Abs(y)GetNphi()*cradle->GetCellSideSize()/2 ) - return cradle; - } + // --- Polystyrene foam --- + Float_t AF[2] = {12.011, 1.00794} ; + Float_t ZF[2] = {6.0, 1.0} ; + Float_t WF[2] = {1.0, 1.0} ; + Float_t DF = 0.12 ; - return NULL; -} + AliMixture(4, "Foam$", AF, ZF, DF, -2, WF) ; -//______________________________________________________________________________ + // --- Titanium --- + Float_t ATIT[3] = {47.88, 26.98, 54.94} ; + Float_t ZTIT[3] = {22.0, 13.0, 25.0} ; + Float_t WTIT[3] = {69.0, 6.0, 1.0} ; + Float_t DTIT = 4.5 ; -void AliPHOS::Reconstruction(Float_t signal_step, UInt_t min_signal_reject) -{ -// Call AliPHOSCradle::Reconstruction(Float_t signal_step, UInt_t min_signal_reject) -// for all AliPHOSCradles. + AliMixture(5, "Titanium$", ATIT, ZTIT, DTIT, -3, WTIT); - for( int i=0; iGetEntries(); i++ ) - GetCradle(i).Reconstruction(signal_step,min_signal_reject); -} + // --- Silicon --- + AliMaterial(6, "Si$", 28.0855, 14., 2.33, 9.36, 42.3, 0, 0) ; -//______________________________________________________________________________ -void AliPHOS::ResetDigits(void) -{ - AliDetector::ResetDigits(); - for( int i=0; iGetEntries(); i++ ) - ((AliPHOSCradle*)(*fCradles)[i]) -> Clear(); -} - -//______________________________________________________________________________ + // --- Foam thermo insulation --- + Float_t ATI[2] = {12.011, 1.00794} ; + Float_t ZTI[2] = {6.0, 1.0} ; + Float_t WTI[2] = {1.0, 1.0} ; + Float_t DTI = 0.1 ; -void AliPHOS::FinishEvent(void) -{ -// Called at the end of each 'galice' event. + AliMixture(7, "Thermo Insul.$", ATI, ZTI, DTI, -2, WTI) ; - if( NULL!=fTreePHOS ) - fTreePHOS->Fill(); -} - -//______________________________________________________________________________ - -void AliPHOS::FinishRun(void) -{ -} + // --- Textolitn --- + Float_t ATX[4] = {16.0, 28.09, 12.011, 1.00794} ; + Float_t ZTX[4] = {8.0, 14.0, 6.0, 1.0} ; + Float_t WTX[4] = {292.0, 68.0, 462.0, 736.0} ; + Float_t DTX = 1.75 ; -//______________________________________________________________________________ + AliMixture(8, "Textolit$", ATX, ZTX, DTX, -4, WTX) ; -void AliPHOS::Print(Option_t *opt) -{ -// Print PHOS information. -// For each AliPHOSCradle the function AliPHOSCradle::Print(opt) is called. + //--- FR4 --- + Float_t AFR[3] = {28.0855, 15.9994, 17.749} ; + Float_t ZFR[3] = {14., 8., 8.875} ; + Float_t WFR[3] = {.28, .32, .4} ; + Float_t DFR = 1.8 ; - AliPHOS &PHOS = *(AliPHOS *)this; // Removing 'const'... + AliMixture(9, "FR4$", AFR, ZFR, DFR, -3, WFR) ; - for( int i=0; iGetEntries(); i++ ) - { - printf("PHOS cradle %d from %d\n",i+1, fCradles->GetEntries()); - PHOS.GetCradle(i).Print(opt); - printf( "---------------------------------------------------\n"); - } -} + // --- The Composite Material for micromegas (so far polyetylene) --- + Float_t ACM[2] = {12.01, 1.} ; + Float_t ZCM[2] = {6., 1.} ; + Float_t WCM[2] = {1., 2.} ; + Float_t DCM = 0.935 ; -//______________________________________________________________________________ -void AliPHOS::SetFlags(Float_t p1,Float_t p2,Float_t p3,Float_t p4, - Float_t p5,Float_t p6,Float_t p7,Float_t p8,Float_t p9) -{ - PHOSflags[0]=p1; - PHOSflags[1]=p2; - PHOSflags[2]=p3; - PHOSflags[3]=p4; - PHOSflags[4]=p5; - PHOSflags[5]=p6; - PHOSflags[6]=p7; - PHOSflags[7]=p8; - PHOSflags[8]=p9; -} + AliMixture(10, "Compo Mat$", ACM, ZCM, DCM, -2, WCM) ; -//______________________________________________________________________________ -void AliPHOS::SetCell(Float_t p1,Float_t p2,Float_t p3,Float_t p4, - Float_t p5,Float_t p6,Float_t p7,Float_t p8,Float_t p9) -{ - PHOScell[0]=p1; - PHOScell[1]=p2; - PHOScell[2]=p3; - PHOScell[3]=p4; - PHOScell[4]=p5; - PHOScell[5]=p6; - PHOScell[6]=p7; - PHOScell[7]=p8; - PHOScell[8]=p9; -} - -//______________________________________________________________________________ -void AliPHOS::SetRadius(Float_t radius) -{ - PHOSradius=radius; -} - -//______________________________________________________________________________ -void AliPHOS::SetCradleSize(Int_t nz, Int_t nphi, Int_t ncradles) -{ - PHOSsize[0]=nz; - PHOSsize[1]=nphi; - PHOSsize[2]=ncradles; -} - -//______________________________________________________________________________ -void AliPHOS::SetCradleA(Float_t angle) -{ - PHOScradlesA=angle; -} - -//______________________________________________________________________________ -void AliPHOS::SetExtra(Float_t p1,Float_t p2,Float_t p3,Float_t p4, - Float_t p5,Float_t p6,Float_t p7,Float_t p8,Float_t p9) -{ - PHOSextra[0] = p1; - PHOSextra[1] = p2; - PHOSextra[2] = p3; - PHOSextra[3] = p4; - PHOSextra[4] = p5; - PHOSextra[5] = p6; - PHOSextra[6] = p7; - PHOSextra[7] = p8; - PHOSextra[8] = p9; -} - -//______________________________________________________________________________ -void AliPHOS::SetTextolitWall(Float_t dx, Float_t dy, Float_t dz) -{ - PHOSTXW[0] = dx; - PHOSTXW[1] = dy; - PHOSTXW[2] = dz; -} + // --- Copper --- + AliMaterial(11, "Cu$", 63.546, 29, 8.96, 1.43, 14.8, 0, 0) ; + + // --- G10 : Printed Circuit material --- + Float_t AG10[4] = { 12., 1., 16., 28.} ; + Float_t ZG10[4] = { 6., 1., 8., 14.} ; + Float_t WG10[4] = { .259, .288, .248, .205} ; + Float_t DG10 = 1.7 ; + + AliMixture(12, "G10$", AG10, ZG10, DG10, -4, WG10); -//______________________________________________________________________________ -void AliPHOS::SetInnerAir(Float_t dx, Float_t dy, Float_t dz) -{ - PHOSAIR[0] = dx; - PHOSAIR[1] = dy; - PHOSAIR[2] = dz; -} + // --- Lead --- + AliMaterial(13, "Pb$", 207.2, 82, 11.35, 0.56, 0., 0, 0) ; -//______________________________________________________________________________ -void AliPHOS::SetFoam(Float_t dx, Float_t dy, Float_t dz, Float_t dr) -{ - PHOSFTI[0] = dx; - PHOSFTI[1] = dy; - PHOSFTI[2] = dz; - PHOSFTI[3] = dr; -} + // --- The gas mixture --- + // Co2 + Float_t ACO[2] = {12.0, 16.0} ; + Float_t ZCO[2] = {6.0, 8.0} ; + Float_t WCO[2] = {1.0, 2.0} ; + Float_t DCO = 0.001977 ; -ClassImp(AliPHOSCradle) - -//______________________________________________________________________________ - -AliPHOSCradle::AliPHOSCradle(void) {} - -//______________________________________________________________________________ - -AliPHOSCradle::AliPHOSCradle( int Geometry , - float CrystalSideSize , - float CrystalLength , - float WrapThickness , - float AirThickness , - float PIN_SideSize , - float PIN_Length , - float Radius , - int Nz , - int Nphi , - float Angle ) : - fGeometry (Geometry), -// fCellEnergy (), -// fChargedTracksInPIN (), - fCrystalSideSize (CrystalSideSize), - fCrystalLength (CrystalLength), - fWrapThickness (WrapThickness), - fAirThickness (AirThickness), - fPIN_SideSize (PIN_SideSize), - fPIN_Length (PIN_Length), - fRadius (Radius), - fNz (Nz), - fNphi (Nphi), - fPhi (Angle) -{ - fCellEnergy = TH2F("CellE","Energy deposition in a cells",fNz,0,fNz,fNphi,0,fNphi); - fCellEnergy .SetDirectory(0); - fChargedTracksInPIN = TH2S("PINCtracks","Amount of charged tracks in PIN",fNz,0,fNz,fNphi,0,fNphi); - fChargedTracksInPIN .SetDirectory(0); -} + AliMixture(14, "CO2$", ACO, ZCO, DCO, -2, WCO); -//______________________________________________________________________________ + // Ar + Float_t DAr = 0.001782 ; + AliMaterial(15, "Ar$", 39.948, 18.0, DAr, 14.0, 0., 0, 0) ; + + // ArCo2 + Char_t namate[21]; + Float_t AGM[2] ; + Float_t ZGM[2] ; + Float_t WGM[2] ; + Float_t DGM ; -AliPHOSCradle::~AliPHOSCradle(void) // 28.12.1998 -{ - fGammasReconstructed.Delete(); - fParticles .Delete(); -} + Float_t AbsL, RadL, Density ; + Float_t buf[1] ; + Int_t nbuf ; -//______________________________________________________________________________ + gMC->Gfmate((*fIdmate)[15], namate, AGM[0], ZGM[0], Density, RadL, AbsL, buf, nbuf) ; // Get properties of Ar + gMC->Gfmate((*fIdmate)[14], namate, AGM[1], ZGM[1], Density, RadL, AbsL, buf, nbuf) ; // Get properties of CO2 -void AliPHOSCradle::Clear(Option_t *) -{ -// Clear digit. information. - fCellEnergy .Reset(); - fChargedTracksInPIN .Reset(); - GetParticles() .Delete(); - GetParticles() .Compress(); - GetGammasReconstructed() .Delete(); - GetGammasReconstructed() .Compress(); + // Create gas mixture -} + Float_t ArContent = 0.80 ; // Ar-content of the Ar/CO2-mixture (80% / 20%) + + WGM[0] = ArContent; + WGM[1] = 1. - ArContent ; + DGM = WGM[0] * DAr + WGM[1] * DCO; -//______________________________________________________________________________ + AliMixture(16, "ArCO2$", AGM, ZGM, DGM, 2, WGM) ; -void AliPHOSCradle::GetXY(const TVector3 &p,const TVector3 &v,float R,float &x,float &y,float &l) const -{ -// This function calculates hit position (x,y) in the CRADLE cells plain from particle in -// the direction given by 'p' (not required to be normalized) and start point -// given by 3-vector 'v'. So the particle trajectory is t(l) = v + p*l -// were 'l' is a number (distance from 'v' to CRADLE cells plain) and 't' is resulting -// three-vector of trajectory point. -// -// After the call to this function user should test that l>=0 (the particle HITED the -// plain) and (x,y) are in the region of CRADLE: -// -// Example: -// AliPHOSCradle cradle(......); -// TVector3 p(....), v(....); -// Float_t x,y,l; -// cradle.GetXY(p,v,x,y,l); -// if( l<0 || TMath::Abs(x)>cradle.GetNz() *cradle.GetCellSideSize()/2 -// || TMath::Abs(y)>cradle.GetNphi()*cradle.GetCellSideSize()/2 ) -// cout << "Outside the CRADLE.\n"; - - // We have to create three vectors: - // s - central point on the PHOS surface - // n1 - first vector in CRADLE plain - // n2 - second vector in CRADLE plain - // This three vectors are orthonormalized. - - double phi = fPhi/180*TMath::Pi(); - TVector3 n1( 0.0 , 0.0 , 1.0 ), // Z direction (X) - n2( -sin(phi) , cos(phi) , 0 ), // around beam (Y) - s ( R*cos(phi) , R*sin(phi) , 0 ); // central point - - const double l1_min = 1e-2; - double l1, - p_n1 = p*n1, // * - scalar product. - p_n2 = p*n2, - v_n1 = v*n1, - v_n2 = v*n2, - s_n1 = s*n1, // 0 - s_n2 = s*n2; // 0 + + // --- Air --- + AliMaterial(99, "Air$", 14.61, 7.3, 0.001205, 30420., 67500., 0, 0) ; - if ( TMath::Abs(l1=p.X()-n1.X()*p_n1-n2.X()*p_n2)>l1_min ) - { l = (-v.X()+s.X()+n1.X()*(v_n1-s_n1)+n2.X()*(v_n2-s_n2))/l1; } - else if ( TMath::Abs(l1=p.Y()-n1.Y()*p_n1-n2.Y()*p_n2)>l1_min ) - { l = (-v.Y()+s.Y()+n1.Y()*(v_n1-s_n1)+n2.Y()*(v_n2-s_n2))/l1; } - else if ( TMath::Abs(l1=p.Z()-n1.Z()*p_n1-n2.Z()*p_n2)>l1_min ) - { l = (-v.Z()+s.Z()+n1.Z()*(v_n1-s_n1)+n2.Z()*(v_n2-s_n2))/l1; } - -// double lx = (-v.X()+s.X()+n1.X()*(v.dot(n1)-s.dot(n1))+n2.X()*(v.dot(n2)-s.dot(n2)))/ -// (p.X()-n1.X()*p.dot(n1)-n2.X()*p.dot(n2)), -// ly = (-v.Y()+s.Y()+n1.Y()*(v.dot(n1)-s.dot(n1))+n2.Y()*(v.dot(n2)-s.dot(n2)))/ -// (p.Y()-n1.Y()*p.dot(n1)-n2.Y()*p.dot(n2)), -// lz = (-v.Z()+s.Z()+n1.Z()*(v.dot(n1)-s.dot(n1))+n2.Z()*(v.dot(n2)-s.dot(n2)))/ -// (p.Z()-n1.Z()*p.dot(n1)-n2.Z()*p.dot(n2)); -// cout.form("x: %g %g %g %g\n",lx,-v.X()+s.X()+n1.X()*(v.dot(n1)-s.dot(n1))+n2.X()*(v.dot(n2)-s.dot(n2)),p.X()-n1.X()*p.dot(n1)-n2.X()*p.dot(n2)); -// cout.form("y: %g %g %g %g\n",lx,-v.Y()+s.Y()+n1.Y()*(v.dot(n1)-s.dot(n1))+n2.Y()*(v.dot(n2)-s.dot(n2)),p.Y()-n1.Y()*p.dot(n1)-n2.Y()*p.dot(n2)); -// cout.form("z: %g %g %g %g\n",lx,-v.Z()+s.Z()+n1.Z()*(v.dot(n1)-s.dot(n1))+n2.Z()*(v.dot(n2)-s.dot(n2)),p.Z()-n1.Z()*p.dot(n1)-n2.Z()*p.dot(n2)); -// cout.form("lx,ly,lz = %g,%g,%g\n",lx,ly,lz); - - x = p_n1*l + v_n1 - s_n1; - y = p_n2*l + v_n2 - s_n2; -} - -//______________________________________________________________________________ + + // DEFINITION OF THE TRACKING MEDIA -void AliPHOSCradle::Print(Option_t *opt) -{ -// Print AliPHOSCradle information. -// -// options: 'd' - print energy deposition for EVERY cell -// 'p' - print particles list that hit the cradle -// 'r' - print list of reconstructed particles - - AliPHOSCradle *cr = (AliPHOSCradle *)this; // Removing 'const'... - - printf("AliPHOSCradle: Nz=%d Nphi=%d, fPhi=%f, E=%g\n",fNz,fNphi,fPhi, - cr->fCellEnergy.GetSumOfWeights()); - - if( NULL!=strchr(opt,'d') ) - { - printf("\n\nCells Energy (in MeV):\n\n |"); - for( int x=0; x=0; y-- ) - { - printf("%3d|",y+1); - for( int x=0; xfCellEnergy.GetBinContent(cr->fCellEnergy.GetBin(x,y))*1000)); - printf("\n"); - } - printf("\n"); - } - - if( NULL!=strchr(opt,'p') ) - { - printf("This cradle was hit by %d particles\n", - ((AliPHOSCradle*)this)->GetParticles().GetEntries()); - TObjArray &p=((AliPHOSCradle*)this)->GetParticles(); - for( int i=0; iPrint(); - } - - if( NULL!=strchr(opt,'p') ) - { - printf("Amount of reconstructed gammas is %d\n", - ((AliPHOSCradle*)this)->GetGammasReconstructed().GetEntries()); - - TObjArray &p=((AliPHOSCradle*)this)->GetGammasReconstructed(); - for( int i=0; iPrint(); - } -} + // for PHOS: idtmed[699->798] equivalent to fIdtmed[0->100] + Int_t * idtmed = fIdtmed->GetArray() - 699 ; + Int_t ISXFLD = gAlice->Field()->Integ() ; + Float_t SXMGMX = gAlice->Field()->Max() ; -//______________________________________________________________________________ + // The scintillator of the calorimeter made of PBW04 -> idtmed[699] + AliMedium(0, "PHOS Xtal $", 0, 1, + ISXFLD, SXMGMX, 10.0, 0.1, 0.1, 0.1, 0.1, 0, 0) ; -void AliPHOSCradle::Distortion(const TH2F *Noise, const TH2F *Stochastic, const TH2F *Calibration) -{ -// This function changes histogram of cell energies fCellEnergy on the base of input -// histograms Noise, Stochastic, Calibration. The histograms must have -// size Nz x Nphi. + // The scintillator of the CPV made of Polystyrene scintillator -> idtmed[700] + AliMedium(1, "CPV scint. $", 1, 1, + ISXFLD, SXMGMX, 10.0, 0.1, 0.1, 0.1, 0.1, 0, 0) ; - ////////////////////////////////// - // Testing the histograms size. // - ////////////////////////////////// - - if( fNz!=fCellEnergy.GetNbinsX() || fNphi!=fCellEnergy.GetNbinsY() ) - { - printf ("Bad size of CellEnergy! Must be: Nz x Nphi = %d x %d\n" - "but size of CellEnergy is: %d x %d\n", - fNz,fNphi,fCellEnergy.GetNbinsX(),fCellEnergy.GetNbinsY()); - exit(1); - } - - if( fNz!=fChargedTracksInPIN.GetNbinsX() || fNphi!=fChargedTracksInPIN.GetNbinsY() ) - { - printf ("Bad size of ChargedTracksInPIN! Must be: Nz x Nphi = %d x %d\n" - "but size of ChargedTracksInPIN is: %d x %d\n", - fNz,fNphi,fChargedTracksInPIN.GetNbinsX(),fChargedTracksInPIN.GetNbinsY()); - exit(1); - } - - if( NULL!=Noise && (fNz!=Noise->GetNbinsX() || fNphi!=Noise->GetNbinsX()) ) - { - printf ("Bad size of Noise! Must be: Nz x Nphi = %d x %d\n" - "but size of Noise is: %d x %d\n", - fNz,fNphi,fChargedTracksInPIN.GetNbinsX(),fChargedTracksInPIN.GetNbinsY()); - exit(1); - } - - if( NULL!=Stochastic && (fNz!=Stochastic->GetNbinsX() || fNphi!=Stochastic->GetNbinsX()) ) - { - printf ("Bad size of Stochastic! Must be: Nz x Nphi = %d x %d\n" - "but size of Stochastic is: %d x %d\n", - fNz,fNphi,fChargedTracksInPIN.GetNbinsX(),fChargedTracksInPIN.GetNbinsY()); - exit(1); - } - - if( NULL!=Calibration && (fNz!=Calibration->GetNbinsX() || fNphi!=Calibration->GetNbinsX()) ) - { - printf ("Bad size of Calibration! Must be: Nz x Nphi = %d x %d\n" - "but size of Calibration is: %d x %d\n", - fNz,fNphi,fChargedTracksInPIN.GetNbinsX(),fChargedTracksInPIN.GetNbinsY()); - exit(1); - } - - //////////////////// - // Do distortion! // - //////////////////// - - for( int y=0; y idtmed[701] + AliMedium(2, "Al parts $", 2, 0, + ISXFLD, SXMGMX, 10.0, 0.1, 0.1, 0.001, 0.001, 0, 0) ; -//////////////////////////////////////////////////////////////////////////////// + // The Tywek which wraps the calorimeter crystals -> idtmed[702] + AliMedium(3, "Tyvek wrapper$", 3, 0, + ISXFLD, SXMGMX, 10.0, 0.1, 0.1, 0.001, 0.001, 0, 0) ; -TH2F* AliPHOSCradle::CreateHistForDistortion(const char *name, const char *title, - Int_t Nx, Int_t Ny, - Float_t MU_mu, Float_t MU_sigma, - Float_t SIGMA_mu, Float_t SIGMA_sigma) -{ -// Create (new TH2F(...)) histogram with information (for every bin) that will -// be used for VALUE creation. -// Two values will be created for each bin: -// MU = TRandom::Gaus(MU_mu,MU_sigma) -// and -// SIGMA = TRandom::Gaus(SIGMA_mu,SIGMA_sigma) -// The VALUE in a particluar bin will be equal -// VALUE = TRandom::Gaus(MU,SIGMA) -// -// Do not forget to delete the histogram at the end of the work. - - TH2F *h = new TH2F( name,title, Nx,1,Nx, Ny,1,Ny ); - if( h==NULL ) - { - Error("CreateHistForDistortion","Can not create the histogram"); - exit(1); - } - h->SetDirectory(0); - - for( int y=0; yGetBin(x,y); - h->SetBinContent(n,r.Gaus( MU_mu, MU_sigma)); - h->SetBinError (n,r.Gaus(SIGMA_mu,SIGMA_sigma)); - } - - return h; -} + // The Polystyrene foam around the calorimeter module -> idtmed[703] + AliMedium(4, "Polyst. foam $", 4, 0, + ISXFLD, SXMGMX, 10.0, 0.1, 0.1, 0.1, 0.1, 0, 0) ; -//////////////////////////////////////////////////////////////////////////////// + // The Titanium around the calorimeter crystal -> idtmed[704] + AliMedium(5, "Titan. cover $", 5, 0, + ISXFLD, SXMGMX, 10.0, 0.1, 0.1, 0.0001, 0.0001, 0, 0) ; -Float_t AliPHOSCradle::GetDistortedValue(const TH2F *h, UInt_t n) -{ - return r.Gaus(((TH2F*)h)->GetBinContent(n),n); -} + // The Silicon of the pin diode to read out the calorimeter crystal -> idtmed[705] + AliMedium(6, "Si PIN $", 6, 0, + ISXFLD, SXMGMX, 10.0, 0.1, 0.1, 0.01, 0.01, 0, 0) ; -//////////////////////////////////////////////////////////////////////////////// -//______________________________________________________________________________ + // The thermo insulating material of the box which contains the calorimeter module -> idtmed[706] + AliMedium(7, "Thermo Insul.$", 7, 0, + ISXFLD, SXMGMX, 10.0, 0.1, 0.1, 0.1, 0.1, 0, 0) ; -#ifdef WIN32 - #define common_for_event_storing COMMON_FOR_EVENT_STORING -#else - #define common_for_event_storing common_for_event_storing_ -#endif + // The Textolit which makes up the box which contains the calorimeter module -> idtmed[707] + AliMedium(8, "Textolit $", 8, 0, + ISXFLD, SXMGMX, 10.0, 0.1, 0.1, 0.1, 0.1, 0, 0) ; -/* extern "C" */ struct -{ - enum { crystals_matrix_amount_max=4, crystals_in_matrix_amount_max=40000 }; - - // Event-independent information - UShort_t crystals_matrix_amount_PHOS, - crystal_matrix_type, - amount_of_crystals_on_Z, - amount_of_crystals_on_PHI; - Float_t radius, - crystal_size, - crystal_length, - matrix_coordinate_Z [crystals_matrix_amount_max], - matrix_coordinate_PHI [crystals_matrix_amount_max]; - UInt_t event_number; - UShort_t crystals_amount_with_amplitudes [crystals_matrix_amount_max], - crystals_amplitudes_Iad [crystals_matrix_amount_max] - [crystals_in_matrix_amount_max][2]; -} common_for_event_storing; - -// integer*4 crystals_amount_max,crystals_in_matrix_amount_max, -// + crystals_matrix_amount_max -// parameter (crystals_matrix_amount_max=4) -// parameter (crystals_in_matrix_amount_max=40000) -// parameter (crystals_amount_max =crystals_matrix_amount_max* -// + crystals_in_matrix_amount_max) -// -// * All units are in GeV, cm, radian -// real crystal_amplitudes_unit, radius_unit, -// + crystal_size_unit, crystal_length_unit, -// + matrix_coordinate_Z_unit, matrix_coordinate_PHI_unit -// integer crystal_amplitudes_in_units_min -// parameter (crystal_amplitudes_in_units_min = 1) -// parameter (crystal_amplitudes_unit = 0.001 ) ! 1.0 MeV -// parameter (radius_unit = 0.1 ) ! 0.1 cm -// parameter (crystal_size_unit = 0.01 ) ! 0.01 cm -// parameter (crystal_length_unit = 0.01 ) ! 0.01 cm -// parameter (matrix_coordinate_Z_unit = 0.1 ) ! 0.1 cm -// parameter (matrix_coordinate_PHI_unit = 1e-4 ) ! 1e-4 radian -// -// integer*2 crystals_matrix_amount_PHOS, crystal_matrix_type, -// + amount_of_crystals_on_Z, amount_of_crystals_on_PHI, -// + crystals_amount_with_amplitudes, crystals_amplitudes_Iad -// integer*4 event_number -// -// real radius, crystal_size, crystal_length, -// + matrix_coordinate_Z, matrix_coordinate_PHI -// -// real crystals_amplitudes, crystals_energy_total -// integer event_file_unit_number -// -// common /common_for_event_storing/ -// + ! Event-independent information -// + crystals_matrix_amount_PHOS, -// + crystal_matrix_type, -// + amount_of_crystals_on_Z, -// + amount_of_crystals_on_PHI, -// + radius, -// + crystal_size, -// + crystal_length, -// + matrix_coordinate_Z (crystals_matrix_amount_max), -// + matrix_coordinate_PHI (crystals_matrix_amount_max), -// + -// + ! Event-dependent information -// + event_number, -// + crystals_amount_with_amplitudes -// + (crystals_matrix_amount_max), -// + crystals_amplitudes_Iad (2,crystals_in_matrix_amount_max, -// + crystals_matrix_amount_max), -// + -// + ! These information don't store in data file -// + crystals_amplitudes (crystals_amount_max), -// + crystals_energy_total, -// + event_file_unit_number - - -// parameter (NGp=1000,nsps=10,nvertmax=1000) -// COMMON /GAMMA/KG,MW(ngp),ID(ngp),JD(ngp),E(ngp),E4(ngp), -// , XW(ngp),YW(ngp),ES(nsps,ngp),ET(nsps,ngp),ISsd(ngp), -// , IGDEV(ngp),ZGDEV(ngp),sigexy(3,ngp),Emimx(2,nsps,ngp), -// , kgfix,igfix(ngp),cgfix(3,ngp),sgfix(3,ngp),hiw(ngp), -// , wsw(nsps,ngp),h1w(ngp),h0w(ngp),raxay(5,ngp), -// , sigmaes0(nsps,ngp),dispeces(nsps,ngp), -// , igamvert(ngp) - - -#ifdef WIN32 -#define rcgamma RCGAMMA -#else -#define rcgamma rcgamma_ -#endif - -/* extern "C" */ struct -{ - enum {NGP=1000, nsps=10, nvertmax=1000}; - int recons_gammas_amount, mw[NGP],ID[NGP],JD[NGP]; - float E[NGP], E4[NGP], XW[NGP], YW[NGP], ES[NGP][nsps],ET[NGP][nsps],ISsd[NGP], - igdev[NGP],Zgdev[NGP]; -// sigexy(3,ngp),Emimx(2,nsps,ngp), -// , kgfix,igfix(ngp),cgfix(3,ngp),sgfix(3,ngp),hiw(ngp), -// , wsw(nsps,ngp),h1w(ngp),h0w(ngp),raxay(5,ngp), -// , sigmaes0(nsps,ngp),dispeces(nsps,ngp), -// , igamvert(ngp) -} rcgamma; - -/* -#ifdef WIN32 -#define reconsfirst RECONSFIRST -#define type_of_call _stdcall -#else -#define reconsfirst reconsfirst_ -#define type_of_call -#endif - -extern "C" void type_of_call reconsfirst(const float &,const float &); -*/ - -void AliPHOSCradle::Reconstruction(Float_t signal_step, UInt_t min_signal_reject) -{ -// Call of PHOS reconstruction program. -// signal_step=0.001 GeV (1MeV) -// min_signal_reject = 15 or 30 MeV - - - common_for_event_storing.event_number = 0; // We do not know event number? - common_for_event_storing.crystals_matrix_amount_PHOS = 1; - common_for_event_storing.crystal_matrix_type = 1; // 1 - rectangular - common_for_event_storing.amount_of_crystals_on_Z = fNz; - common_for_event_storing.amount_of_crystals_on_PHI = fNphi; - - common_for_event_storing.radius = fRadius; - common_for_event_storing.crystal_size = GetCellSideSize(); - common_for_event_storing.crystal_length = fCrystalLength; - - common_for_event_storing.matrix_coordinate_Z [0] = 0; - common_for_event_storing.matrix_coordinate_PHI [0] = fPhi; - - #define k common_for_event_storing.crystals_amount_with_amplitudes[0] - k=0; - - for( int y=0; y=min_signal_reject ) - { - common_for_event_storing.crystals_amplitudes_Iad[0][k][0] = signal; - common_for_event_storing.crystals_amplitudes_Iad[0][k][1] = x + y*fNz; - k++; - } - } - #undef k - - GetGammasReconstructed().Delete(); - GetGammasReconstructed().Compress(); - - const float stochastic_term = 0.03, // per cents over sqrt(E); E in GeV - electronic_noise = 0.01; // GeV -// reconsfirst(stochastic_term,electronic_noise); // Call of reconstruction program. - - for( int i=0; i idtmed[708] + AliMedium(9, "FR4 $", 9, 0, + ISXFLD, SXMGMX, 10.0, 0.1, 0.1, 0.1, 0.0001, 0, 0) ; -//______________________________________________________________________________ -//______________________________________________________________________________ -//______________________________________________________________________________ -//______________________________________________________________________________ -//______________________________________________________________________________ -ClassImp(AliPHOSgamma) + // The Composite Material for micromegas -> idtmed[709] + AliMedium(10, "CompoMat $", 10, 0, + ISXFLD, SXMGMX, 10.0, 0.1, 0.1, 0.1, 0.1, 0, 0) ; -//______________________________________________________________________________ + // Copper -> idtmed[710] + AliMedium(11, "Copper $", 11, 0, + ISXFLD, SXMGMX, 10.0, 0.1, 0.1, 0.1, 0.0001, 0, 0) ; -void AliPHOSgamma::Print(Option_t *) -{ - float mass = fE*fE - fPx*fPx - fPy*fPy - fPz*fPz; + // G10: Printed Circuit material -> idtmed[711] + + AliMedium(12, "G10 $", 12, 0, + ISXFLD, SXMGMX, 10.0, 0.1, 0.1, 0.1, 0.01, 0, 0) ; - if( mass>=0 ) - mass = sqrt( mass); - else - mass = -sqrt(-mass); + // The Lead -> idtmed[712] + + AliMedium(13, "Lead $", 13, 0, + ISXFLD, SXMGMX, 10.0, 0.1, 0.1, 0.1, 0.1, 0, 0) ; - printf("XY=(%+7.2f,%+7.2f) (%+7.2f,%+7.2f,%+7.2f;%7.2f) mass=%8.4f Ipart=%2d\n", - fX,fY,fPx,fPy,fPz,fE,mass,fIpart); -} + // The gas mixture: ArCo2 -> idtmed[715] + + AliMedium(16, "ArCo2 $", 16, 1, + ISXFLD, SXMGMX, 10.0, 0.1, 0.1, 0.1, 0.01, 0, 0) ; + + // Air -> idtmed[798] + AliMedium(99, "Air $", 99, 0, + ISXFLD, SXMGMX, 10.0, 1.0, 0.1, 0.1, 10.0, 0, 0) ; -//______________________________________________________________________________ + // --- Set decent energy thresholds for gamma and electron tracking -AliPHOSgamma &AliPHOSgamma::operator=(const AliPHOSgamma &g) -{ - fX = g.fX; - fY = g.fY; - fE = g.fE; - fPx = g.fPx; - fPy = g.fPy; - fPz = g.fPz; - fIpart = g.fIpart; - - return *this; -} + // Tracking threshold for photons and electrons in the scintillator crystal + gMC->Gstpar(idtmed[699], "CUTGAM",0.5E-4) ; + gMC->Gstpar(idtmed[699], "CUTELE",1.0E-4) ; -//______________________________________________________________________________ -//______________________________________________________________________________ -//______________________________________________________________________________ -//______________________________________________________________________________ -//______________________________________________________________________________ + // Tracking threshold for photons and electrons in the gas + gMC->Gstpar(idtmed[715], "CUTGAM",0.5E-4) ; + gMC->Gstpar(idtmed[715], "CUTELE",1.0E-4) ; -ClassImp(AliPHOShit) + // --- Generate explicitly delta rays in the titan cover --- + gMC->Gstpar(idtmed[704], "LOSS",3.) ; + gMC->Gstpar(idtmed[704], "DRAY",1.) ; -//______________________________________________________________________________ + // --- and in aluminium parts --- + gMC->Gstpar(idtmed[701], "LOSS",3.) ; + gMC->Gstpar(idtmed[701], "DRAY",1.) ; -AliPHOShit::AliPHOShit(Int_t shunt, Int_t track, Int_t *vol, Float_t *hits): -AliHit(shunt, track) -{ - Int_t i; - for (i=0;i<5;i++) fVolume[i] = vol[i]; - fX = hits[0]; - fY = hits[1]; - fZ = hits[2]; - fELOS = hits[3]; } - -//______________________________________________________________________________ diff --git a/PHOS/AliPHOS.h b/PHOS/AliPHOS.h index 7ebf9c28930..b5f3b4ddbc8 100644 --- a/PHOS/AliPHOS.h +++ b/PHOS/AliPHOS.h @@ -1,242 +1,39 @@ -#ifndef PHOS_H -#define PHOS_H +#ifndef ALIPHOS_H +#define ALIPHOS_H /* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * * See cxx source for full Copyright notice */ -/* $Id$ */ - //////////////////////////////////////////////// -// Manager and hits classes for set:PHOS // +// Abstract Base Class for PHOS // +// Version SUBATECH // +// Author Laurent Aphecetche SUBATECH // +// The only provided method here is // +// CreateMaterials, which defines the // +// materials common to all PHOS versions. // +// // //////////////////////////////////////////////// - -// --- ROOT system --- -#include -#include -#include -#include - -// --- galice header files --- -#include "AliDetector.h" -#include "AliHit.h" -#include "AliRun.h" - -class AliPHOSgamma : public TObject { - - public: - virtual ~AliPHOSgamma(void) {} - AliPHOSgamma(void) {} - AliPHOSgamma(const AliPHOSgamma &g) { *this=g; } - AliPHOSgamma(Float_t X, Float_t Y, Float_t E, - Float_t Px, Float_t Py, Float_t Pz, - Int_t Ipart) : - fX(X), fY(Y), fE(E), - fPx(Px), fPy(Py), fPz(Pz), - fIpart(Ipart) - {} - - Float_t fX; // cm. x-coordinate (in beam direction) - Float_t fY; // cm. y-coordinate (around beam) - - Float_t fE; // GeV. energy - - Float_t fPx; // GeV. Gamma momentum Px - Float_t fPy; // GeV. Gamma momentum Py - Float_t fPz; // GeV. Gamma momentum Pz - - Int_t fIpart; // Current particle number (GEANT particle code) - - void Print(Option_t *options=NULL); - AliPHOSgamma &operator=(const AliPHOSgamma &g); - - private: - - ClassDef(AliPHOSgamma,1) // Gamma particle in PHOS cradle -}; - -//______________________________________________________________________________ - -class AliPHOShit : public AliHit { - -public: - Int_t fVolume[5]; //array of volumes - Float_t fELOS; //ELOS - -public: - AliPHOShit() {} - AliPHOShit(Int_t shunt, Int_t track, Int_t *vol, Float_t *hits); - virtual ~AliPHOShit() {} - - ClassDef(AliPHOShit,1) //Hits object for set:PHOS -}; - -//______________________________________________________________________________ - - -class AliPHOSCradle : public TObject { - - public: - - virtual ~AliPHOSCradle(void); - AliPHOSCradle(void); - AliPHOSCradle(int Geometry , - float CrystalSideSize , - float CrystalLength , - float WrapThickness , - float AirThickness , - float PIN_SideSize , - float PIN_Length , - float Radius , - int Nz , - int Nphi , - float Angle ); - Float_t GetCrystalSideSize (void) const {return fCrystalSideSize;} - Float_t GetCellSideSize (void) const {return fCrystalSideSize+2*fWrapThickness+2*fAirThickness;} - Float_t GetCrystalLength (void) const {return fCrystalLength;} - Float_t GetWrapThickness (void) const {return fWrapThickness;} - Float_t GetAirThickness (void) const {return fAirThickness;} - Float_t GetPIN_SideSize (void) const {return fPIN_SideSize;} - Float_t GetPIN_Length (void) const {return fPIN_Length;} - Float_t GetRadius (void) const {return fRadius;} - Int_t GetNz (void) const {return fNz;} - Int_t GetNphi (void) const {return fNphi;} - Float_t GetPhi (void) const {return fPhi;} - - void Clear(Option_t *opt=""); // Clear all data. - void Print(Option_t *opt=""); - void Distortion(const TH2F *Noise=NULL, const TH2F *Stochastic=NULL, const TH2F *Calibration=NULL); - TH2F *CreateHistForDistortion(const char *name, const char *title, Int_t Nx, Int_t Ny, - Float_t MU_mu, Float_t MU_sigma, Float_t SIGMA_mu, Float_t SIGMA_sigma); - Float_t GetDistortedValue(const TH2F *h, UInt_t n); - - void Reconstruction(Float_t signal_step, UInt_t min_signal_reject); - void GetXY(const TVector3 &p,const TVector3 &v,float R,float &x,float &y,float &l) const; - - TObjArray &GetGammasReconstructed (void) {return fGammasReconstructed;} - TObjArray &GetParticles (void) {return fParticles;} - - TH2F fCellEnergy; // GeV. Energy in cells - TH2S fChargedTracksInPIN; // amount. hits in PIN - - - private: - - Int_t fGeometry; // Geometry type: 1 or 2 - Float_t fCrystalSideSize; // cm. - Float_t fCrystalLength; // cm. - Float_t fWrapThickness; // cm. - Float_t fAirThickness; // cm. - Float_t fPIN_SideSize; // cm. - Float_t fPIN_Length; // cm. - - Float_t fRadius; // cm. Distance to PHOS - - - Int_t fNz; // Cells amount in beam direction - Int_t fNphi; // Cells amount around beam - - Float_t fPhi; // degree. Position of CRADLE center - - TObjArray fGammasReconstructed; // List of reconstructed gammas - TObjArray fParticles; // List of particles in the direction of this cradle +// --- ROOT system --- - TRandom r; //! Random number class, do not stream +// --- AliRoot header files --- -// friend class AliPHOS; +#include "AliDetector.h" +#include "AliPHOSGeometry.h" - ClassDef(AliPHOSCradle,1) // PHOS cradle -}; class AliPHOS : public AliDetector { public: - enum {CradlesBranch_Bit=1}; - - AliPHOS(); - AliPHOS(const char *name, const char *title); - virtual ~AliPHOS(); - virtual void AddHit(Int_t, Int_t*, Float_t*); - virtual void BuildGeometry(); - virtual void CreateGeometry() {} - virtual void CreateMaterials(); - Int_t DistancetoPrimitive(Int_t px, Int_t py); - void FinishEvent(void); - - virtual void Init(); - virtual Int_t IsVersion() const =0; - void MakeBranch(Option_t *option); - void SetTreeAddress(void); - void FinishRun(void); - void ResetDigits(void); - void Print(Option_t *opt=""); - AliPHOSCradle *GetCradleOfTheParticle(const TVector3 &p,const TVector3 &v) const; - AliPHOSCradle &GetCradle(int n) {return *(AliPHOSCradle*)fCradles->operator[](n);} - // AliPHOSCradle &GetCradle(int n) {return *((AliPHOSCradle*) (*fCradles)[n]) ;} - void Reconstruction(Float_t signal_step, UInt_t min_signal_reject); - virtual void SetFlags(Float_t p1,Float_t p2=0,Float_t p3=0,Float_t p4=0, - Float_t p5=0,Float_t p6=0,Float_t p7=0,Float_t p8=0,Float_t p9=0); - virtual void SetCell(Float_t p1,Float_t p2=0,Float_t p3=0,Float_t p4=0, - Float_t p5=0,Float_t p6=0,Float_t p7=0,Float_t p8=0,Float_t p9=0); - virtual void SetRadius(Float_t radius); - virtual void SetCradleSize(Int_t nz, Int_t nphi, Int_t ncradles); - virtual void SetCradleA(Float_t angle); - virtual void SetExtra(Float_t p1,Float_t p2=0,Float_t p3=0,Float_t p4=0, - Float_t p5=0,Float_t p6=0,Float_t p7=0,Float_t p8=0,Float_t p9=0); - virtual void SetTextolitWall(Float_t dx, Float_t dy, Float_t dz); - virtual void SetInnerAir(Float_t dx, Float_t dy, Float_t dz); - virtual void SetFoam(Float_t dx, Float_t dy, Float_t dz, Float_t dr); - virtual void StepManager()=0; - virtual void DefPars(); - virtual void AddPHOSCradles(); - - - virtual Int_t GetPHOS_IDTMED_PbWO4 (void){return (*fIdtmed)[0];} - virtual Int_t GetPHOS_IDTMED_Al (void){return (*fIdtmed)[2];} - virtual Int_t GetPHOS_IDTMED_Tyvek (void){return (*fIdtmed)[3];} - virtual Int_t GetPHOS_IDTMED_PIN (void){return (*fIdtmed)[4];} - virtual Int_t GetPHOS_IDTMED_AIR (void){return (*fIdtmed)[99];} - - virtual Int_t &GetPHOS_Ndiv_magic (void) {return PHOS_Ndiv_magic;} - virtual Float_t GetCrystalSideSize (void) const {return PHOScell[0]; } - virtual Float_t GetCrystalLength (void) const {return PHOScell[1]; } - virtual Float_t GetWrapThickness (void) const {return PHOScell[2]; } - virtual Float_t GetAirThickness (void) const {return PHOScell[3]; } - virtual Float_t GetPIN_SideSize (void) const {return PHOScell[4]; } - virtual Float_t GetPIN_Length (void) const {return PHOScell[5]; } - virtual Float_t GetRadius (void) const {return PHOSradius; } - virtual Int_t GetNz (void) const {return PHOSsize[0]; } - virtual Int_t GetNphi (void) const {return PHOSsize[1]; } - virtual Int_t GetCradlesAmount (void) const {return PHOSsize[2]; } - virtual Float_t GetAngleBetweenCradles(void) const {return PHOScradlesA;} - virtual Float_t GetPHOS_flag (Int_t n) const {return PHOSflags[n];} - virtual Float_t GetPHOSextra (Int_t n) const {return PHOSextra[n];} - virtual Float_t GetPHOSFoam (Int_t n) const {return PHOSFTI[n];} - virtual Float_t GetPHOStxwall (Int_t n) const {return PHOSTXW[n];} - virtual Float_t GetPHOSAir (Int_t n) const {return PHOSAIR[n];} - virtual Float_t &GetCradleAngle (Int_t n) {return PHOSangle[n];} - - - TObjArray *fCradles; //! Cradles in PHOS - Int_t fDebugLevel; - - TTree *fTreePHOS; //! Pointer to PHOS tree. - -private: - - TString fBranchNameOfCradles; // - TString fTreeName; // Name of PHOS tree: "PHOS" - -#define MAXCRAD 100 + AliPHOS(const char* name, const char* title) ; + AliPHOS() ; + virtual ~AliPHOS() ; + + virtual void CreateMaterials() ; + virtual AliPHOSGeometry * GetGeometry() = 0 ; - Float_t PHOSflags[9], PHOScell[9], PHOSradius; - Int_t PHOSsize[3]; - Float_t PHOScradlesA,PHOSTXW[3],PHOSAIR[3],PHOSFTI[4],PHOSextra[9], - PHOSangle[MAXCRAD]; - Int_t PHOS_Ndiv_magic; + ClassDef(AliPHOS,2) // Photon Spectrometer Detector - ClassDef(AliPHOS,1) //Hits manager for set:PHOS -}; - -#endif +} ; +#endif // ALIPHOS_H diff --git a/PHOS/AliPHOSClusterizer.cxx b/PHOS/AliPHOSClusterizer.cxx new file mode 100644 index 00000000000..5a8ed337977 --- /dev/null +++ b/PHOS/AliPHOSClusterizer.cxx @@ -0,0 +1,45 @@ +/************************************************************************** + * 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. * + **************************************************************************/ + +//_________________________________________________________________________ +// A brief description of the class +//*-- Author : Yves Schutz SUBATECH +////////////////////////////////////////////////////////////////////////////// + +// --- ROOT system --- + + + +// --- Standard library --- + + + +// --- AliRoot header files --- + +#include "AliPHOSClusterizer.h" + +ClassImp(AliPHOSClusterizer) + +//____________________________________________________________________________ +AliPHOSClusterizer::AliPHOSClusterizer() +{ + // ctor +} + +//____________________________________________________________________________ +AliPHOSClusterizer::~AliPHOSClusterizer() +{ + // dtor +} diff --git a/PHOS/AliPHOSClusterizer.h b/PHOS/AliPHOSClusterizer.h new file mode 100644 index 00000000000..938bf19eb80 --- /dev/null +++ b/PHOS/AliPHOSClusterizer.h @@ -0,0 +1,43 @@ +#ifndef ALIPHOSCLUSTERIZER_H +#define ALIPHOSCLUSTERIZER_H +/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * See cxx source for full Copyright notice */ + +//////////////////////////////////////////////// +// Algorithme class for the clusterization // +// interface class // +// Version SUBATECH // +// Author Yves Schutz SUBATECH // +// // +// pABC // +//////////////////////////////////////////////// + +// --- ROOT system --- + +#include "TObject.h" +#include "TClonesArray.h" + +// --- Standard library --- + +// --- AliRoot header files --- + + +typedef TClonesArray RecPointsList ; // a cluster has a variable size (see ROOT FAQ) +typedef TClonesArray DigitsList ; //for digits saved on disk + +class AliPHOSClusterizer : public TObject { + +public: + + AliPHOSClusterizer() ; // ctor + virtual ~AliPHOSClusterizer() ; // dtor + + virtual Float_t Calibrate(Int_t Amp) = 0 ; + virtual void GetNumberOfClustersFound(Int_t * numb) = 0 ; + virtual void MakeClusters(const DigitsList * dl, RecPointsList * emccl, RecPointsList * ppsdl) = 0 ; + + ClassDef(AliPHOSClusterizer,1) // clusterization interface, version 1 + +} ; + +#endif // AliPHOSCLUSTERIZER_H diff --git a/PHOS/AliPHOSClusterizerv1.cxx b/PHOS/AliPHOSClusterizerv1.cxx new file mode 100644 index 00000000000..0ef7b04ee61 --- /dev/null +++ b/PHOS/AliPHOSClusterizerv1.cxx @@ -0,0 +1,252 @@ +/************************************************************************** + * 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. * + **************************************************************************/ + +//_________________________________________________________________________ +// A brief description of the class +//*-- Author : Yves Schutz SUBATECH +////////////////////////////////////////////////////////////////////////////// + +// --- ROOT system --- + +#include "TMath.h" + +// --- Standard library --- + +#include "iostream.h" + +// --- AliRoot header files --- + +#include "AliPHOSClusterizerv1.h" +#include "AliPHOSDigit.h" +#include "AliPHOSEmcRecPoint.h" +#include "AliPHOSPpsdRecPoint.h" +#include "AliPHOSv0.h" +#include "AliRun.h" + +ClassImp(AliPHOSClusterizerv1) + +//____________________________________________________________________________ +AliPHOSClusterizerv1::AliPHOSClusterizerv1() +{ + fA = 0.; + fB = 0.01 ; + fNumberOfEmcClusters = 0 ; + fNumberOfPpsdClusters = 0 ; + fEmcClusteringThreshold = 0.1; + fEmcEnergyThreshold = 0.01; + fPpsdClusteringThreshold = 0.00000015; + fPpsdEnergyThreshold = 0.0000001; + fW0 = 4.5 ; + fLocMaxCut = 0.06 ; +} + +//____________________________________________________________________________ +Int_t AliPHOSClusterizerv1::AreNeighbours(AliPHOSDigit * d1, AliPHOSDigit * d2) +{ + // neigbours are defined as digits having at least common virtix + // The order of A and B in AreNeighbours(A,B) is important: first (A) should be digit + // in cluster, which compared with digits, which not clasterized yet + Int_t rv = 0 ; + + AliPHOSGeometry * geom = AliPHOSGeometry::GetInstance() ; + + Int_t relid1[4] ; + geom->AbsToRelNumbering(d1->GetId(), relid1) ; + + Int_t relid2[4] ; + geom->AbsToRelNumbering(d2->GetId(), relid2) ; + + if ( (relid1[0] == relid2[0]) && (relid1[1]==relid2[1]) ) { // inside the same PHOS module and the same PPSD Module + Int_t RowDiff = TMath::Abs( relid1[2] - relid2[2] ) ; + Int_t ColDiff = TMath::Abs( relid1[3] - relid2[3] ) ; + + if (( ColDiff<=1 ) && ( RowDiff <= 1 )){ + rv = 1 ; + } + else { + if((relid2[2] > relid1[2]) && (relid2[3] > relid1[3]+1)) + rv = 2; // Difference in row numbers is too large to look further + } + + } + else { + + if( (relid1[0] < relid2[0]) || (relid1[1] < relid2[1]) ) + rv=2 ; + + } + + return rv ; +} + +//____________________________________________________________________________ +void AliPHOSClusterizerv1::FillandSort(const DigitsList * dl, TObjArray * tl) +{ + // copies the digits with energy above thershold and sorts the list + // according to increasing Id number + + + AliPHOSGeometry * geom = AliPHOSGeometry::GetInstance() ; + Int_t relid[4] ; + + TIter next(dl) ; + AliPHOSDigit * digit ; + + while ( (digit = (AliPHOSDigit *)next()) ) { + + Int_t id = digit->GetId() ; + Float_t ene = Calibrate(digit->GetAmp()) ; + geom->AbsToRelNumbering(id, relid) ; + + if(relid[1]==0){ // EMC + if ( ene > fEmcEnergyThreshold ) + tl->Add(digit) ; + } + + else { //Ppsd + if ( ene > fPpsdEnergyThreshold ) + tl->Add(digit) ; + } + + } + tl->Sort() ; +} + +//____________________________________________________________________________ +void AliPHOSClusterizerv1:: GetNumberOfClustersFound(Int_t * numb) +{ + numb[0] = fNumberOfEmcClusters ; + numb[1] = fNumberOfPpsdClusters ; +} + +//____________________________________________________________________________ +Bool_t AliPHOSClusterizerv1::IsInEmc(AliPHOSDigit * digit) +{ + Bool_t rv = kFALSE ; + + AliPHOSGeometry * geom = AliPHOSGeometry::GetInstance() ; + + Int_t relid[4] ; + geom->AbsToRelNumbering(digit->GetId(), relid) ; + + if ( relid[1] == 0 ) + rv = kTRUE; + + return rv ; +} + +//____________________________________________________________________________ +void AliPHOSClusterizerv1::MakeClusters(const DigitsList * dl, RecPointsList * emcl, RecPointsList * ppsdl) +{ + + // Fill and sort the working digits list + TObjArray TempoDigitsList( dl->GetEntries() ) ; + this->FillandSort(dl, &TempoDigitsList) ; + + + + // Clusterization starts + TIter nextdigit(&TempoDigitsList) ; + AliPHOSDigit * digit ; + Bool_t NotRemoved = kTRUE ; + + while ( (digit = (AliPHOSDigit *)nextdigit()) ) { // scan over the list of digits + AliPHOSRecPoint * clu ; + + int * ClusterDigitsList[dl->GetEntries()] ; + Int_t index ; + if (( (this->IsInEmc(digit)) && (Calibrate(digit->GetAmp()) > fEmcClusteringThreshold ) ) || + ( (!this->IsInEmc(digit)) && (Calibrate(digit->GetAmp()) > fPpsdClusteringThreshold ) ) ) { + + Int_t iDigitInCluster = 0 ; + + if (this->IsInEmc(digit) ) { + new ((*emcl)[fNumberOfEmcClusters]) AliPHOSEmcRecPoint(fW0, fLocMaxCut) ;// start a new EMC RecPoint + + clu = (AliPHOSEmcRecPoint *) (*emcl)[fNumberOfEmcClusters++] ; + + clu->AddDigit(*digit,Calibrate(digit->GetAmp())) ; + + ClusterDigitsList[iDigitInCluster++] = (int* ) digit ; + TempoDigitsList.Remove(digit) ; + } + + else { + new ((*ppsdl)[fNumberOfPpsdClusters]) AliPHOSPpsdRecPoint() ;// start a new PPSD cluster + clu = (AliPHOSPpsdRecPoint *) ppsdl->At(fNumberOfPpsdClusters++) ; + clu->AddDigit(*digit,0.) ; + ClusterDigitsList[iDigitInCluster++] = (int* ) digit ; + TempoDigitsList.Remove(digit) ; + nextdigit.Reset() ; + + //Here we remove resting EMC digits, which can not make cluster + if(NotRemoved){ + + while( (digit = (AliPHOSDigit *)nextdigit()) ){ + + if(IsInEmc(digit)) TempoDigitsList.Remove(digit) ; + else + break ; + + }// while digit + + } // if NotRemoved + + } // else + + nextdigit.Reset() ; + + AliPHOSDigit * digitN ; + index = 0 ; + while (index < iDigitInCluster){ // scan over digits already in claster + digit = (AliPHOSDigit *) ClusterDigitsList[index++] ; + + while ( (digitN = (AliPHOSDigit *)nextdigit()) ) { // scan over the reduced list of digits + Int_t ineb = AreNeighbours(digit, digitN); // call (digit,digitN) in THAT oder !!!!! + switch (ineb ) { + case 0 : // not a neibors + break ; + case 1 : // Are neibors + clu->AddDigit(*digitN,Calibrate(digitN->GetAmp())) ; + ClusterDigitsList[iDigitInCluster++] =(int*) digitN ; + TempoDigitsList.Remove(digitN) ; + break ; + case 2 : // to far from each other + goto endofloop; + } // switch + + } // while digitN + + endofloop: ; + nextdigit.Reset() ; + + } // loop over cluster + + } //below energy theshold + + } // while digit + + TempoDigitsList.Clear() ; +} + +//____________________________________________________________________________ +void AliPHOSClusterizerv1::PrintParameters() +{ + cout << "PHOS Clusterizer version 1 :" << endl + << " EMC Clustering threshold = " << fEmcClusteringThreshold << endl + << " EMC Energy threshold = " << fEmcEnergyThreshold << endl + << " PPSD Clustering threshold = " << fPpsdClusteringThreshold << endl + << " PPSD Energy threshold = " << fPpsdEnergyThreshold << endl ; +} diff --git a/PHOS/AliPHOSClusterizerv1.h b/PHOS/AliPHOSClusterizerv1.h new file mode 100644 index 00000000000..cbe1d712bdf --- /dev/null +++ b/PHOS/AliPHOSClusterizerv1.h @@ -0,0 +1,69 @@ +#ifndef ALIPHOSCLUSTERIZERV1_H +#define ALIPHOSCLUSTERIZERV1_H +/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * See cxx source for full Copyright notice */ + +//////////////////////////////////////////////// +// Clusterizer implementation version 1 // +// algorithme class // +// // +// Author Yves Schutz SUBATECH // +// // +// // +//////////////////////////////////////////////// + +// --- ROOT system --- + +// --- Standard library --- + +// --- AliRoot header files --- + +#include "AliPHOSClusterizer.h" +#include "AliPHOSDigit.h" + + + +class AliPHOSClusterizerv1 : public AliPHOSClusterizer { + +public: + + AliPHOSClusterizerv1() ; // ctor + virtual ~AliPHOSClusterizerv1(){} ; // dtor + + Int_t AreNeighbours(AliPHOSDigit * d1, AliPHOSDigit * d2) ; // checks if digits are in neighbour cells + Float_t Calibrate(Int_t Amp){ return fA + fB*Amp ;} //Tranforms Amp to energy + void FillandSort(const DigitsList * dl, TObjArray * tl) ; // sorts the list according to increasing id + Float_t GetLogWeightCut(void){return fW0 ; } + Float_t GetLocalMaxCut(void) {return fLocMaxCut ; } + virtual void GetNumberOfClustersFound(Int_t * numb) ; + Bool_t IsInEmc(AliPHOSDigit * digit) ; // tells id digit is in EMCA + virtual void MakeClusters(const DigitsList * dl, RecPointsList * emcl, RecPointsList * ppsdl) ; // does the job + void PrintParameters() ; + void SetCalibrationParameters(Float_t A,Float_t B){ fA = A ; fB = B;} + void SetEmcClusteringThreshold(Float_t cluth) { fEmcClusteringThreshold = cluth ; } + void SetEmcEnergyThreshold(Float_t enth) { fEmcEnergyThreshold = enth ; } + void SetLocalMaxCut(Float_t cut) { fLocMaxCut = cut ; } + void SetLogWeightCut(Float_t w) { fW0 = w ; } + void SetPpsdClusteringThreshold(Float_t cluth) { fPpsdClusteringThreshold = cluth ; } + void SetPpsdEnergyThreshold(Float_t enth) { fPpsdEnergyThreshold = enth ; } + +private: + + Float_t fA ; + Float_t fB ; + Float_t fEmcClusteringThreshold ; + Float_t fEmcEnergyThreshold ; + Float_t fLocMaxCut ; + Int_t fNumberOfEmcClusters ; + Int_t fNumberOfPpsdClusters ; + Float_t fPpsdClusteringThreshold ; + Float_t fPpsdEnergyThreshold ; + Float_t fW0 ; + +public: + + ClassDef(AliPHOSClusterizerv1,1) // Clusterizer implementation , version 1 + +}; + +#endif // AliPHOSCLUSTERIZERV1_H diff --git a/PHOS/AliPHOSDigit.cxx b/PHOS/AliPHOSDigit.cxx new file mode 100644 index 00000000000..ca0ddd8cc67 --- /dev/null +++ b/PHOS/AliPHOSDigit.cxx @@ -0,0 +1,86 @@ +/************************************************************************** + * 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. * + **************************************************************************/ + +//_________________________________________________________________________ +// Digit class for PHOS that contains an absolute ID and an energy +//*-- Author : Laurent Aphecetche SUBATECH +////////////////////////////////////////////////////////////////////////////// + +// --- ROOT system --- + +// --- Standard library --- + +#include + +// --- AliRoot header files --- + +#include "AliPHOSDigit.h" + + +ClassImp(AliPHOSDigit) + +//____________________________________________________________________________ +AliPHOSDigit::AliPHOSDigit(Int_t id, Int_t DigEnergy) : AliDigitNew( ), fId(id),fAmp(DigEnergy) +{ + // This part should be a true Digitization, but is not for the moment. + // fAmp = energy ; + fId = id; + fAmp = DigEnergy ; +} + +//____________________________________________________________________________ +Bool_t AliPHOSDigit::operator==(AliPHOSDigit const &Digit) const +{ + if ( fId == Digit.fId ) + return kTRUE ; + else + return kFALSE ; +} + +//____________________________________________________________________________ +AliPHOSDigit& AliPHOSDigit::operator+(AliPHOSDigit const &Digit) +{ + fAmp += Digit.fAmp ; + + return *this ; +} + +//____________________________________________________________________________ +ostream& operator << ( ostream& out , const AliPHOSDigit& Digit) +{ + out << "ID " << Digit.fId << " Energy = " << Digit.fAmp ; + + return out ; +} + +//____________________________________________________________________________ +Int_t AliPHOSDigit::Compare(TObject * obj) +{ + Int_t rv ; + + AliPHOSDigit * digit = (AliPHOSDigit *)obj ; + + Int_t iddiff = fId - digit->GetId() ; + + if ( iddiff > 0 ) + rv = 1 ; + else if ( iddiff < 0 ) + rv = -1 ; + else + rv = 0 ; + + return rv ; + +} diff --git a/PHOS/AliPHOSDigit.h b/PHOS/AliPHOSDigit.h new file mode 100644 index 00000000000..d6bada6d210 --- /dev/null +++ b/PHOS/AliPHOSDigit.h @@ -0,0 +1,52 @@ +#ifndef ALIPHOSDIGIT_H +#define ALIPHOSDIGIT_H +/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * See cxx source for full Copyright notice */ + +//////////////////////////////////////////////// +// The digit class: a list of abs Id, energy// +// Version SUBATECH // +// Author Laurent Aphecetche SUBATECH // +// comment: added sortable YS // +// // +//////////////////////////////////////////////// + +// --- ROOT system --- + +#include "TObject.h" + +// --- Standard library --- + +// --- AliRoot header files --- + +#include "AliDigitNew.h" + +class AliPHOSDigit : public AliDigitNew { + +public: + + AliPHOSDigit() {} + AliPHOSDigit(Int_t id, Int_t DigEnergy) ; + virtual ~AliPHOSDigit() {} + + Bool_t operator==(AliPHOSDigit const &rValue) const; + AliPHOSDigit& operator+(AliPHOSDigit const &rValue) ; + + friend ostream& operator << ( ostream& , const AliPHOSDigit&) ; + + Int_t Compare(TObject * obj) ; + Int_t GetId() { return fId ; } + Int_t GetAmp() { return fAmp ; } + Bool_t IsSortable() const{ return kTRUE ; } + +private: + Int_t fId ; // absolute id + Int_t fAmp ; // digitalized energy + +public: + + ClassDef(AliPHOSDigit,1) // Digit in PHOS, version 1 + +} ; + +#endif // ALIPHOSDIGIT_H diff --git a/PHOS/AliPHOSEmcRecPoint.cxx b/PHOS/AliPHOSEmcRecPoint.cxx new file mode 100644 index 00000000000..2c33578d9c7 --- /dev/null +++ b/PHOS/AliPHOSEmcRecPoint.cxx @@ -0,0 +1,418 @@ +/************************************************************************** + * 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. * + **************************************************************************/ + +//_________________________________________________________________________ +// Rec Point in the PHOS EM calorimeter +//*-- Author : Dmitri Peressounko RRC KI +////////////////////////////////////////////////////////////////////////////// + +// --- ROOT system --- + +#include "TMath.h" + +// --- Standard library --- + +#include "iostream.h" + +// --- AliRoot header files --- + +#include "AliPHOSGeometry.h" +#include "AliPHOSEmcRecPoint.h" +#include "AliRun.h" + +ClassImp(AliPHOSEmcRecPoint) + +//____________________________________________________________________________ +AliPHOSEmcRecPoint::AliPHOSEmcRecPoint(Float_t W0, Float_t LocMaxCut) + : AliPHOSRecPoint() +{ + // ctor + + fMulDigit = 0 ; + fAmp = 0. ; + fEnergyList = new Float_t[fMaxDigit]; + AliPHOSGeometry * PHOSGeom = (AliPHOSGeometry *) fGeom ; + fDelta = PHOSGeom->GetCrystalSize(0) ; + fW0 = W0 ; + fLocMaxCut = LocMaxCut ; + fLocPos.SetX(1000000.) ; //Local position should be evaluated +} + +// //____________________________________________________________________________ +// AliPHOSEmcRecPoint::~AliPHOSEmcRecPoint() +// { +// // dtor +// } + +//____________________________________________________________________________ +void AliPHOSEmcRecPoint::AddDigit(AliDigitNew & digit, Float_t Energy) +{ + // adds a digit to the digits list + // and accumulates the total amplitude and the multiplicity + + if ( fMulDigit >= fMaxDigit ) { // increase the size of the lists + int * tempo = new ( int[fMaxDigit*=2] ) ; + Float_t * tempoE = new ( Float_t[fMaxDigit*=2] ) ; + Int_t index ; + + for ( index = 0 ; index < fMulDigit ; index++ ){ + tempo[index] = fDigitsList[index] ; + tempoE[index] = fEnergyList[index] ; + } + + delete fDigitsList ; + delete fEnergyList ; + fDigitsList = tempo ; + fEnergyList = tempoE ; + } + + fDigitsList[fMulDigit] = (int) &digit ; + fEnergyList[fMulDigit++]= Energy ; + fAmp += Energy ; +} + +//____________________________________________________________________________ +Bool_t AliPHOSEmcRecPoint::AreNeighbours(AliPHOSDigit * digit1, AliPHOSDigit * digit2 ) +{ + + Bool_t aren = kFALSE ; + + AliPHOSGeometry * PHOSGeom = (AliPHOSGeometry *) fGeom ; + Int_t relid1[4] ; + PHOSGeom->AbsToRelNumbering(digit1->GetId(), relid1) ; + + Int_t relid2[4] ; + PHOSGeom->AbsToRelNumbering(digit2->GetId(), relid2) ; + + Int_t RowDiff = TMath::Abs( relid1[2] - relid2[2] ) ; + Int_t ColDiff = TMath::Abs( relid1[3] - relid2[3] ) ; + + if (( ColDiff<=1 ) && ( RowDiff <= 1 ) && (ColDiff+RowDiff > 0)) + aren = kTRUE ; + + return aren ; +} + +//____________________________________________________________________________ +Int_t AliPHOSEmcRecPoint::Compare(TObject * obj) +{ + Int_t rv ; + + AliPHOSEmcRecPoint * clu = (AliPHOSEmcRecPoint *)obj ; + + + Int_t PHOSMod1 = this->GetPHOSMod() ; + Int_t PHOSMod2 = clu->GetPHOSMod() ; + + TVector3 LocPos1; + this->GetLocalPosition(LocPos1) ; + TVector3 LocPos2; + clu->GetLocalPosition(LocPos2) ; + + if(PHOSMod1 == PHOSMod2 ) { + Int_t rowdif = (Int_t)TMath::Ceil(LocPos1.X()/fDelta)-(Int_t)TMath::Ceil(LocPos2.X()/fDelta) ; + if (rowdif> 0) + rv = -1 ; + else if(rowdif < 0) + rv = 1 ; + else if(LocPos1.Z()>LocPos2.Z()) + rv = -1 ; + else + rv = 1 ; + } + + else { + if(PHOSMod1 < PHOSMod2 ) + rv = -1 ; + else + rv = 1 ; + } + + return rv ; +} + +//____________________________________________________________________________ +Float_t AliPHOSEmcRecPoint::GetDispersion() +{ + Float_t D = 0 ; + Float_t wtot = 0 ; + + TVector3 LocPos; + GetLocalPosition(LocPos); + Float_t x = LocPos.X() ; + Float_t z = LocPos.Z() ; + // Int_t i = GetPHOSMod() ; + + AliPHOSDigit * digit ; + AliPHOSGeometry * PHOSGeom = (AliPHOSGeometry *) fGeom ; + + Int_t iDigit; + for(iDigit=0; iDigitAbsToRelNumbering(digit->GetId(), relid) ; + PHOSGeom->RelPosInModule(relid, xi, zi); + Float_t w = TMath::Max(0.,fW0+TMath::Log(fEnergyList[iDigit]/fAmp ) ) ; + D += w*((xi-x)*(xi-x) + (zi-z)*(zi-z) ) ; + wtot+=w ; + } + + D /= wtot ; + + return TMath::Sqrt(D) ; +} + +//____________________________________________________________________________ +void AliPHOSEmcRecPoint::GetElipsAxis(Float_t * lambda) +{ + Float_t wtot = 0. ; + Float_t x = 0.; + Float_t z = 0.; + Float_t Dxx = 0.; + Float_t Dzz = 0.; + Float_t Dxz = 0.; + + AliPHOSDigit * digit ; + AliPHOSGeometry * PHOSGeom = (AliPHOSGeometry *) fGeom ; + Int_t iDigit; + + for(iDigit=0; iDigitAbsToRelNumbering(digit->GetId(), relid) ; + PHOSGeom->RelPosInModule(relid, xi, zi); + Float_t w = TMath::Max(0.,fW0+TMath::Log(fEnergyList[iDigit]/fAmp ) ) ; + Dxx += w * xi * xi ; + x += w * xi ; + Dzz += w * zi * zi ; + z += w * zi ; + Dxz += w * xi * zi ; + wtot += w ; + } + + Dxx /= wtot ; + x /= wtot ; + Dxx -= x * x ; + Dzz /= wtot ; + z /= wtot ; + Dzz -= z * z ; + Dxz /= wtot ; + Dxz -= x * z ; + + lambda[0] = TMath::Sqrt( 0.5 * (Dxx + Dzz) + TMath::Sqrt( 0.25 * (Dxx - Dzz) * (Dxx - Dzz) + Dxz * Dxz ) ) ; + lambda[1] = TMath::Sqrt( 0.5 * (Dxx + Dzz) - TMath::Sqrt( 0.25 * (Dxx - Dzz) * (Dxx - Dzz) + Dxz * Dxz ) ) ; +} + +//____________________________________________________________________________ +Float_t AliPHOSEmcRecPoint::GetMaximalEnergy(void) +{ + Float_t menergy = 0. ; + + Int_t iDigit; + + for(iDigit=0; iDigit menergy) + menergy = fEnergyList[iDigit] ; + } + return menergy ; +} + +//____________________________________________________________________________ +Int_t AliPHOSEmcRecPoint::GetMultiplicityAtLevel(Float_t H) +{ + Int_t multipl = 0 ; + Int_t iDigit ; + for(iDigit=0; iDigit H * fAmp) + multipl++ ; + } + return multipl ; +} + +//____________________________________________________________________________ +Int_t AliPHOSEmcRecPoint::GetNumberOfLocalMax(int * maxAt, Float_t * maxAtEnergy) +{ + AliPHOSDigit * digit ; + AliPHOSDigit * digitN ; + + + Int_t iDigitN ; + Int_t iDigit ; + + for(iDigit=0; iDigit fEnergyList[iDigitN] ) { + maxAt[iDigitN] = -1 ; + //but may be digit is not local max too ? + if(fEnergyList[iDigit]fEnergyList[iDigitN]-fLocMaxCut) + maxAt[iDigitN] = -1 ; + } + } // if Areneighbours + } // while digitN + } // slot not empty + } // while digit + + iDigitN = 0 ; + for(iDigit=0; iDigitAbsToRelNumbering(digit->GetId(), relid) ; + PHOSGeom->RelPosInModule(relid, xi, zi); + Float_t w = TMath::Max( 0., fW0 + TMath::Log( fEnergyList[iDigit] / fAmp ) ) ; + x += xi * w ; + z += zi * w ; + wtot += w ; + + } + + x /= wtot ; + z /= wtot ; + fLocPos.SetX(x) ; + fLocPos.SetY(0.) ; + fLocPos.SetZ(z) ; + + LPos = fLocPos ; +} + +// //____________________________________________________________________________ +// AliPHOSEmcRecPoint& AliPHOSEmcRecPoint::operator = (AliPHOSEmcRecPoint Clu) +// { +// int * DL = Clu.GetDigitsList() ; + +// if(fDigitsList) +// delete fDigitsList ; + +// AliPHOSDigit * digit ; + +// Int_t iDigit; + +// for(iDigit=0; iDigit +#include "assert.h" + +// --- AliRoot header files --- + +#include "AliPHOSGeometry.h" +#include "AliPHOSPpsdRecPoint.h" +#include "AliConst.h" + +ClassImp(AliPHOSGeometry) + + AliPHOSGeometry * AliPHOSGeometry::fGeom = 0 ; + +//____________________________________________________________________________ +AliPHOSGeometry::~AliPHOSGeometry(void) +{ + fRotMatrixArray->Delete() ; + delete fRotMatrixArray ; +} + +//____________________________________________________________________________ +Bool_t AliPHOSGeometry::AbsToRelNumbering(const Int_t AbsId, Int_t * RelId) +{ + // RelId[0] = PHOS Module number 1:fNModules + // RelId[1] = 0 if PbW04 + // = PPSD Module number 1:fNumberOfModulesPhi*fNumberOfModulesZ*2 (2->up and bottom level) + // RelId[2] = Row number inside a PHOS or PPSD module + // RelId[3] = Column number inside a PHOS or PPSD module + + Bool_t rv = kTRUE ; + Float_t Id = AbsId ; + + Int_t PHOSModuleNumber = (Int_t)TMath:: Ceil( Id / ( GetNPhi() * GetNZ() ) ) ; + + if ( PHOSModuleNumber > GetNModules() ) { // its a PPSD pad + + Id -= GetNPhi() * GetNZ() * GetNModules() ; + Float_t tempo = 2 * GetNumberOfModulesPhi() * GetNumberOfModulesZ() * GetNumberOfPadsPhi() * GetNumberOfPadsZ() ; + RelId[0] = (Int_t)TMath::Ceil( Id / tempo ) ; + Id -= ( RelId[0] - 1 ) * tempo ; + RelId[1] = (Int_t)TMath::Ceil( Id / ( GetNumberOfPadsPhi() * GetNumberOfPadsZ() ) ) ; + Id -= ( RelId[1] - 1 ) * GetNumberOfPadsPhi() * GetNumberOfPadsZ() ; + RelId[2] = (Int_t)TMath::Ceil( Id / GetNumberOfPadsPhi() ) ; + RelId[3] = (Int_t) ( Id - ( RelId[2] - 1 ) * GetNumberOfPadsPhi() ) ; + } + else { // its a PW04 crystal + + RelId[0] = PHOSModuleNumber ; + RelId[1] = 0 ; + Id -= ( PHOSModuleNumber - 1 ) * GetNPhi() * GetNZ() ; + RelId[2] = (Int_t)TMath::Ceil( Id / GetNPhi() ) ; + RelId[3] = (Int_t)( Id - ( RelId[2] - 1 ) * GetNPhi() ) ; + } + return rv ; +} + +//____________________________________________________________________________ +void AliPHOSGeometry::GetGlobal(const AliRecPoint* RecPoint, TVector3 & gpos, TMatrix & gmat) +{ + + AliPHOSRecPoint * tmpPHOS = (AliPHOSRecPoint *) RecPoint ; + TVector3 LocalPosition ; + + tmpPHOS->GetLocalPosition(gpos) ; + + + if ( tmpPHOS->IsEmc() ) // it is a EMC crystal + { gpos.SetY( -(GetIPtoOuterCoverDistance() + GetUpperPlateThickness() + + GetSecondUpperPlateThickness() + GetUpperCoolingPlateThickness()) ) ; + + } + else + { // it is a PPSD pad + AliPHOSPpsdRecPoint * tmpPpsd = (AliPHOSPpsdRecPoint *) RecPoint ; + if (tmpPpsd->GetUp() ) // it is an upper module + { + gpos.SetY(-( GetIPtoOuterCoverDistance() - GetMicromegas2Thickness() - + GetLeadToMicro2Gap() - GetLeadConverterThickness() - + GetMicro1ToLeadGap() - GetMicromegas1Thickness() / 2.0 ) ) ; + } + else // it is a lower module + gpos.SetY(-( GetIPtoOuterCoverDistance() - GetMicromegas2Thickness() / 2.0) ) ; + } + + Float_t Phi = GetPHOSAngle( tmpPHOS->GetPHOSMod()) ; + Double_t const RADDEG = 180.0 / kPI ; + Float_t rPhi = Phi / RADDEG ; + + TRotation Rot ; + Rot.RotateZ(-rPhi) ; // a rotation around Z by angle + + TRotation dummy = Rot.Invert() ; // to transform from original frame to rotate frame + gpos.Transform(Rot) ; // rotate the baby +} + +//____________________________________________________________________________ +void AliPHOSGeometry::GetGlobal(const AliRecPoint* RecPoint, TVector3 & gpos) +{ + AliPHOSRecPoint * tmpPHOS = (AliPHOSRecPoint *) RecPoint ; + TVector3 LocalPosition ; + tmpPHOS->GetLocalPosition(gpos) ; + + + if ( tmpPHOS->IsEmc() ) // it is a EMC crystal + { gpos.SetY( -(GetIPtoOuterCoverDistance() + GetUpperPlateThickness() + + GetSecondUpperPlateThickness() + GetUpperCoolingPlateThickness()) ) ; + } + else + { // it is a PPSD pad + AliPHOSPpsdRecPoint * tmpPpsd = (AliPHOSPpsdRecPoint *) RecPoint ; + if (tmpPpsd->GetUp() ) // it is an upper module + { + gpos.SetY(-( GetIPtoOuterCoverDistance() - GetMicromegas2Thickness() - + GetLeadToMicro2Gap() - GetLeadConverterThickness() - + GetMicro1ToLeadGap() - GetMicromegas1Thickness() / 2.0 ) ) ; + } + else // it is a lower module + gpos.SetY(-( GetIPtoOuterCoverDistance() - GetMicromegas2Thickness() / 2.0) ) ; + } + + Float_t Phi = GetPHOSAngle( tmpPHOS->GetPHOSMod()) ; + Double_t const RADDEG = 180.0 / kPI ; + Float_t rPhi = Phi / RADDEG ; + + TRotation Rot ; + Rot.RotateZ(-rPhi) ; // a rotation around Z by angle + + TRotation dummy = Rot.Invert() ; // to transform from original frame to rotate frame + gpos.Transform(Rot) ; // rotate the baby +} + +//____________________________________________________________________________ +void AliPHOSGeometry::Init(void) +{ + fRotMatrixArray = new TObjArray(fNModules) ; + + cout << "PHOS geometry setup: parameters for option " << fName << " " << fTitle << endl ; + if ( ((strcmp( fName, "default" )) == 0) || ((strcmp( fName, "GPS2" )) == 0) ) { + fInit = kTRUE ; + this->InitPHOS() ; + this->InitPPSD() ; + this->SetPHOSAngles() ; + } + else { + fInit = kFALSE ; + cout << "PHOS Geometry setup: option not defined " << fName << endl ; + } +} + +//____________________________________________________________________________ +void AliPHOSGeometry::InitPHOS(void) +{ + // PHOS + + fNPhi = 64 ; + fNZ = 64 ; + fNModules = 5 ; + + fPHOSAngle[0] = 0.0 ; // Module position angles are set in CreateGeometry() + fPHOSAngle[1] = 0.0 ; + fPHOSAngle[2] = 0.0 ; + fPHOSAngle[3] = 0.0 ; + + fXtlSize[0] = 2.2 ; + fXtlSize[1] = 18.0 ; + fXtlSize[2] = 2.2 ; + + // all these numbers coming next are subject to changes + + fOuterBoxThickness[0] = 2.8 ; + fOuterBoxThickness[1] = 5.0 ; + fOuterBoxThickness[2] = 5.0 ; + + fUpperPlateThickness = 4.0 ; + + fSecondUpperPlateThickness = 5.0 ; + + fCrystalSupportHeight = 6.95 ; + fCrystalWrapThickness = 0.01 ; + fCrystalHolderThickness = 0.005 ; + fModuleBoxThickness = 2.0 ; + fIPtoOuterCoverDistance = 447.0 ; + fIPtoCrystalSurface = 460.0 ; + + fPinDiodeSize[0] = 1.0 ; + fPinDiodeSize[1] = 0.1 ; + fPinDiodeSize[2] = 1.0 ; + + fUpperCoolingPlateThickness = 0.06 ; + fSupportPlateThickness = 10.0 ; + fLowerThermoPlateThickness = 3.0 ; + fLowerTextolitPlateThickness = 1.0 ; + fGapBetweenCrystals = 0.03 ; + + fTextolitBoxThickness[0] = 1.5 ; + fTextolitBoxThickness[1] = 0.0 ; + fTextolitBoxThickness[2] = 3.0 ; + + fAirThickness[0] = 1.56 ; + fAirThickness[1] = 20.5175 ; + fAirThickness[2] = 2.48 ; + + Float_t XtalModulePhiSize = fNPhi * ( fXtlSize[0] + 2 * fGapBetweenCrystals ) ; + Float_t XtalModuleZSize = fNZ * ( fXtlSize[2] + 2 * fGapBetweenCrystals ) ; + + // The next dimensions are calculated from the above parameters + + fOuterBoxSize[0] = XtalModulePhiSize + 2 * ( fAirThickness[0] + fModuleBoxThickness + + fTextolitBoxThickness[0] + fOuterBoxThickness[0] ) ; + fOuterBoxSize[1] = ( fXtlSize[1] + fCrystalSupportHeight + fCrystalWrapThickness + fCrystalHolderThickness ) + + 2 * (fAirThickness[1] + fModuleBoxThickness + fTextolitBoxThickness[1] + fOuterBoxThickness[1] ) ; + fOuterBoxSize[2] = XtalModuleZSize + 2 * ( fAirThickness[2] + fModuleBoxThickness + + fTextolitBoxThickness[2] + fOuterBoxThickness[2] ) ; + + fTextolitBoxSize[0] = fOuterBoxSize[0] - 2 * fOuterBoxThickness[0] ; + fTextolitBoxSize[1] = fOuterBoxSize[1] - fOuterBoxThickness[1] - fUpperPlateThickness ; + fTextolitBoxSize[2] = fOuterBoxSize[2] - 2 * fOuterBoxThickness[2] ; + + fAirFilledBoxSize[0] = fTextolitBoxSize[0] - 2 * fTextolitBoxThickness[0] ; + fAirFilledBoxSize[1] = fTextolitBoxSize[1] - fSecondUpperPlateThickness ; + fAirFilledBoxSize[2] = fTextolitBoxSize[2] - 2 * fTextolitBoxThickness[2] ; + +} + +//____________________________________________________________________________ +void AliPHOSGeometry::InitPPSD(void) +{ + // PPSD + + fAnodeThickness = 0.0009 ; + fAvalancheGap = 0.01 ; + fCathodeThickness = 0.0009 ; + fCompositeThickness = 0.3 ; + fConversionGap = 0.3 ; + fLeadConverterThickness = 0.56 ; + fLeadToMicro2Gap = 0.1 ; + fLidThickness = 0.2 ; + fMicro1ToLeadGap = 0.1 ; + fMicromegasWallThickness = 0.6 ; + fNumberOfModulesPhi = 4 ; + fNumberOfModulesZ = 4 ; + fNumberOfPadsPhi = 24 ; + fNumberOfPadsZ = 24 ; + fPCThickness = 0.1 ; + fPhiDisplacement = 0.8 ; + fZDisplacement = 0.8 ; + + fMicromegas1Thickness = fLidThickness + 2 * fCompositeThickness + fCathodeThickness + fPCThickness + + fAnodeThickness + fConversionGap + fAvalancheGap ; + fMicromegas2Thickness = fMicromegas1Thickness ; + + + fPPSDModuleSize[0] = 38.0 ; + fPPSDModuleSize[1] = fMicromegas1Thickness ; + fPPSDModuleSize[2] = 38.0 ; + + fPPSDBoxSize[0] = fNumberOfModulesPhi * fPPSDModuleSize[0] + 2 * fPhiDisplacement ; + fPPSDBoxSize[1] = fMicromegas2Thickness + fMicromegas2Thickness + fLeadConverterThickness + fMicro1ToLeadGap + fLeadToMicro2Gap ; + fPPSDBoxSize[2] = fNumberOfModulesZ * fPPSDModuleSize[2] + 2 * fZDisplacement ; + + fIPtoTopLidDistance = fIPtoOuterCoverDistance - fPPSDBoxSize[1] - 1. ; + +} + +//____________________________________________________________________________ +AliPHOSGeometry * AliPHOSGeometry::GetInstance() +{ + assert(fGeom!=0) ; + return (AliPHOSGeometry *) fGeom ; +} + +//____________________________________________________________________________ +AliPHOSGeometry * AliPHOSGeometry::GetInstance(const Text_t* name, const Text_t* title) +{ + AliPHOSGeometry * rv = 0 ; + if ( fGeom == 0 ) { + fGeom = new AliPHOSGeometry(name, title) ; + rv = (AliPHOSGeometry * ) fGeom ; + } + else { + if ( strcmp(fGeom->GetName(), name) != 0 ) { + cout << "AliPHOSGeometry : current geometry is " << fGeom->GetName() << endl + << " you cannot call " << name << endl ; + } + else + rv = (AliPHOSGeometry *) fGeom ; + } + return rv ; +} + +//____________________________________________________________________________ +Bool_t AliPHOSGeometry::RelToAbsNumbering(const Int_t * RelId, Int_t & AbsId) +{ + + // AbsId = 1:fNModules * fNPhi * fNZ -> PbWO4 + // AbsId = 1:fNModules * 2 * (fNumberOfModulesPhi * fNumberOfModulesZ) * fNumberOfPadsPhi * fNumberOfPadsZ -> PPSD + + Bool_t rv = kTRUE ; + + if ( RelId[1] > 0 ) { // its a PPSD pad + + AbsId = GetNPhi() * GetNZ() * GetNModules() // the offset to separate emcal crystals from PPSD pads + + ( RelId[0] - 1 ) * GetNumberOfModulesPhi() * GetNumberOfModulesZ() // the pads offset of PHOS modules + * GetNumberOfPadsPhi() * GetNumberOfPadsZ() * 2 + + ( RelId[1] - 1 ) * GetNumberOfPadsPhi() * GetNumberOfPadsZ() // the pads offset of PPSD modules + + ( RelId[2] - 1 ) * GetNumberOfPadsPhi() // the pads offset of a PPSD row + + RelId[3] ; // the column number + } + else { + if ( RelId[1] == 0 ) { // its a Phos crystal + AbsId = ( RelId[0] - 1 ) * GetNPhi() * GetNZ() // the offset of PHOS modules + + ( RelId[2] - 1 ) * GetNPhi() // the offset of a xtal row + + RelId[3] ; // the column number + } + } + + return rv ; +} + +//____________________________________________________________________________ + +void AliPHOSGeometry::RelPosInAlice(const Int_t Id, TVector3 & pos ) +{ + if (Id > 0) { + + Int_t RelId[4] ; + + AbsToRelNumbering(Id , RelId) ; + + Int_t PHOSModule = RelId[0] ; + + + if ( RelId[1] == 0 ) // it is a PbW04 crystal + { pos.SetY( -(GetIPtoOuterCoverDistance() + GetUpperPlateThickness() + + GetSecondUpperPlateThickness() + GetUpperCoolingPlateThickness()) ) ; + } + if ( RelId[1] > 0 ) { // its a PPSD pad + if ( RelId[1] > GetNumberOfModulesPhi() * GetNumberOfModulesZ() ) // its an bottom module + { + pos.SetY(-( GetIPtoOuterCoverDistance() - GetMicromegas2Thickness() / 2.0) ) ; + } + else // its an upper module + pos.SetY(-( GetIPtoOuterCoverDistance() - GetMicromegas2Thickness() - GetLeadToMicro2Gap() + - GetLeadConverterThickness() - GetMicro1ToLeadGap() - GetMicromegas1Thickness() / 2.0) ) ; + } + + Float_t x, z ; + RelPosInModule(RelId, x, z) ; + + pos.SetX(x); + pos.SetZ(z); + + + Float_t Phi = GetPHOSAngle( PHOSModule) ; + Double_t const RADDEG = 180.0 / kPI ; + Float_t rPhi = Phi / RADDEG ; + + TRotation Rot ; + Rot.RotateZ(-rPhi) ; // a rotation around Z by angle + + TRotation dummy = Rot.Invert() ; // to transform from original frame to rotate frame + + pos.Transform(Rot) ; // rotate the baby + } + else { + pos.SetX(0.); + pos.SetY(0.); + pos.SetZ(0.); + } +} + +//____________________________________________________________________________ +void AliPHOSGeometry::RelPosInModule(const Int_t * RelId, Float_t & x, Float_t & z) +{ + Int_t PPSDModule ; + Int_t Row = RelId[2] ; //offset along z axiz + Int_t Column = RelId[3] ; //offset along x axiz + + Float_t PadSizeZ = GetPPSDModuleSize(2)/ GetNumberOfPadsZ(); + Float_t PadSizeX = GetPPSDModuleSize(0)/ GetNumberOfPadsPhi(); + + if ( RelId[1] == 0 ) { // its a PbW04 crystal + x = -( GetNPhi()/2. - Row + 0.5 ) * GetCrystalSize(0) ; // position ox Xtal with respect + z = -( GetNZ() /2. - Column + 0.5 ) * GetCrystalSize(2) ; // of center of PHOS module + } + else { + if ( RelId[1] > GetNumberOfModulesPhi() * GetNumberOfModulesZ() ) + PPSDModule = RelId[1]-GetNumberOfModulesPhi() * GetNumberOfModulesZ(); + else PPSDModule = RelId[1] ; + Int_t ModRow = 1+(Int_t)TMath::Ceil( (Float_t)PPSDModule / GetNumberOfModulesPhi()-1. ) ; + Int_t ModCol = PPSDModule - ( ModRow-1 ) * GetNumberOfModulesPhi() ; + Float_t x0 = ( GetNumberOfModulesPhi() / 2. - ModRow + 0.5 ) * GetPPSDModuleSize(0) ; + Float_t z0 = ( GetNumberOfModulesZ() / 2. - ModCol + 0.5 ) * GetPPSDModuleSize(2) ; + x = - ( GetNumberOfPadsPhi()/2. - Row - 0.5 ) * PadSizeX + x0 ; // position of pad with respect + z = - ( GetNumberOfPadsZ()/2. - Column - 0.5 ) * PadSizeZ + z0 ; // of center of PHOS module + } +} + +//____________________________________________________________________________ +void AliPHOSGeometry:: SetPHOSAngles() +{ + Double_t const RADDEG = 180.0 / kPI ; + Float_t PPHI = TMath::ATan( fOuterBoxSize[0] / ( 2.0 * fIPtoOuterCoverDistance ) ) ; + PPHI *= RADDEG ; + + for( Int_t i = 1; i <= fNModules ; i++ ) { + Float_t angle = PPHI * 2 * ( i - fNModules / 2.0 - 0.5 ) ; + fPHOSAngle[i-1] = - angle ; + } +} + diff --git a/PHOS/AliPHOSGeometry.h b/PHOS/AliPHOSGeometry.h new file mode 100644 index 00000000000..835d0413818 --- /dev/null +++ b/PHOS/AliPHOSGeometry.h @@ -0,0 +1,180 @@ +#ifndef ALIPHOSGEOMETRY_H +#define ALIPHOSGEOMETRY_H +/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * See cxx source for full Copyright notice */ + +//////////////////////////////////////////////// +// Geometry class for PHOS : singleton // +// Version SUBATECH // +// Author Y. Schutz SUBATECH // +// geometry parametrized for any // +// shape of modules // +//////////////////////////////////////////////// + +// --- ROOT system --- + +#include "TNamed.h" +#include "TString.h" +#include "TObjArray.h" +#include "TVector3.h" + +// --- AliRoot header files --- + +#include "AliGeometry.h" +#include "AliPHOSRecPoint.h" + +class AliPHOSGeometry : public AliGeometry { + +public: + + AliPHOSGeometry() {} ; // must be kept public for root persistency purposes + virtual ~AliPHOSGeometry(void) ; + static AliPHOSGeometry * GetInstance(const Text_t* name, const Text_t* title) ; + static AliPHOSGeometry * GetInstance() ; + virtual void GetGlobal(const AliRecPoint* RecPoint, TVector3 & gpos, TMatrix & gmat) ; + virtual void GetGlobal(const AliRecPoint* RecPoint, TVector3 & gpos) ; + +protected: + + AliPHOSGeometry(const Text_t* name, const Text_t* title) : AliGeometry(name, title) { Init() ; } + void Init(void) ; // steering method for PHOS and CPV + void InitPHOS(void) ; // defines the various PHOS geometry parameters + void InitPPSD(void) ; // defines the various PPSD geometry parameters + +public: + + // General + + Bool_t AbsToRelNumbering(const Int_t AbsId, Int_t * RelId) ; // converts the absolute PHOS numbering to a relative + void RelPosInModule(const Int_t * RelId, Float_t & y, Float_t & z) ; // gets the position of element (pad or Xtal) relative to + // center of PHOS module + void RelPosInAlice(const Int_t AbsId, TVector3 & pos) ; // gets the position of element (pad or Xtal) relative to + // Alice + Bool_t RelToAbsNumbering(const Int_t * RelId, Int_t & AbsId) ; // converts the absolute PHOS numbering to a relative + // inlines + + ///////////// PHOS related parameters + + Bool_t IsInitialized(void) const { return fInit ; } + Float_t GetAirFilledBoxSize(Int_t index) const { return fAirFilledBoxSize[index] ;} + Float_t GetCrystalHolderThickness(void) const { return fCrystalHolderThickness ; } + Float_t GetCrystalSize(Int_t index) const { return fXtlSize[index] ; } + Float_t GetCrystalSupportHeight(void) const { return fCrystalSupportHeight ; } + Float_t GetCrystalWrapThickness(void) const { return fCrystalWrapThickness;} + Float_t GetGapBetweenCrystals(void) const { return fGapBetweenCrystals ; } + Float_t GetIPtoCrystalSurface(void) const { return fIPtoCrystalSurface ; } + Float_t GetIPtoOuterCoverDistance(void) const { return fIPtoOuterCoverDistance ; } + Float_t GetIPtoTopLidDistance(void) const { return fIPtoTopLidDistance ; } + Float_t GetLowerThermoPlateThickness(void) const { return fLowerThermoPlateThickness ; } + Float_t GetLowerTextolitPlateThickness(void) const { return fLowerTextolitPlateThickness ; } + Float_t GetModuleBoxThickness(void) const { return fModuleBoxThickness ; } + Int_t GetNPhi(void) const { return fNPhi ; } + Int_t GetNZ(void) const { return fNZ ; } + Int_t GetNModules(void) const { return fNModules ; } + Float_t GetOuterBoxSize(Int_t index) const { return fOuterBoxSize[index] ; } + Float_t GetOuterBoxThickness(Int_t index) const { return fOuterBoxThickness[index] ; } + Float_t GetPHOSAngle(Int_t index) const { return fPHOSAngle[index-1] ; } + Float_t GetPinDiodeSize(Int_t index) const { return fPinDiodeSize[index] ; } + Float_t GetSecondUpperPlateThickness(void) const { return fSecondUpperPlateThickness ; } + Float_t GetSupportPlateThickness(void) const { return fSupportPlateThickness ; } + Float_t GetTextolitBoxSize(Int_t index) const { return fTextolitBoxSize[index] ; } + Float_t GetTextolitBoxThickness(Int_t index) const { return fTextolitBoxThickness[index]; } + Float_t GetUpperPlateThickness(void) const { return fUpperPlateThickness ; } + Float_t GetUpperCoolingPlateThickness(void) const { return fUpperCoolingPlateThickness ; } + +private: + + void SetPHOSAngles() ; // calculates the PHOS modules PHI angle + +public: + + ///////////// PPSD (PHOS PRE SHOWER DETECTOR) related parameters + + + Float_t GetAnodeThickness(void) const { return fAnodeThickness ; } + Float_t GetAvalancheGap(void) const { return fAvalancheGap ; } + Float_t GetCathodeThickness(void) const { return fCathodeThickness ; } + Float_t GetCompositeThickness(void) const { return fCompositeThickness ; } + Float_t GetConversionGap(void) const { return fConversionGap ; } + Float_t GetLeadConverterThickness(void) const { return fLeadConverterThickness ; } + Float_t GetLeadToMicro2Gap(void) const { return fLeadToMicro2Gap ; } + Float_t GetLidThickness(void) const { return fLidThickness ; } + Float_t GetMicromegas1Thickness(void) const { return fMicromegas1Thickness ; } + Float_t GetMicromegas2Thickness(void) const { return fMicromegas2Thickness ; } + Float_t GetMicromegasWallThickness(void) const { return fMicromegasWallThickness ; } + Float_t GetMicro1ToLeadGap(void) const { return fMicro1ToLeadGap ; } + Int_t GetNumberOfPadsPhi(void) const { return fNumberOfPadsPhi ; } + Int_t GetNumberOfPadsZ(void) const { return fNumberOfPadsZ ; } + Int_t GetNumberOfModulesPhi(void) const { return fNumberOfModulesPhi ; } + Int_t GetNumberOfModulesZ(void) const { return fNumberOfModulesZ ; } + Float_t GetPCThickness(void) const { return fPCThickness ; } + Float_t GetPhiDisplacement(void) const { return fPhiDisplacement ; } + Float_t GetPPSDBoxSize(Int_t index) const { return fPPSDBoxSize[index] ; } + Float_t GetPPSDModuleSize(Int_t index) const { return fPPSDModuleSize[index] ; } + Float_t GetZDisplacement(void) const { return fZDisplacement ; } + +private: + + ///////////// PHOS related parameters + + Float_t fAirFilledBoxSize[3] ; // Air filled box containing one module + Float_t fAirThickness[3] ; // Space filled with air between the module box and the Textolit box + Float_t fCrystalSupportHeight ; // Height of the support of the crystal + Float_t fCrystalWrapThickness ; // Thickness of Tyvek wrapping the crystal + Float_t fCrystalHolderThickness ; // Titanium holder of the crystal + Float_t fGapBetweenCrystals ; // Total Gap between two adjacent crystals + Bool_t fInit ; // Tells if geometry has been succesfully set up + Float_t fIPtoOuterCoverDistance ; // Distances from interaction point to outer cover + Float_t fIPtoCrystalSurface ; // Distances from interaction point to Xtal surface + Float_t fModuleBoxThickness ; // Thickness of the thermo insulating box containing one crystals module + Float_t fLowerTextolitPlateThickness ; // Thickness of lower textolit plate + Float_t fLowerThermoPlateThickness ; // Thickness of lower thermo insulating plate + Int_t fNModules ; // Number of modules constituing PHOS + Int_t fNPhi ; // Number of crystal units in X (phi) direction + Int_t fNZ ; // Number of crystal units in Z direction + Float_t fOuterBoxSize[3] ; // Size of the outer thermo insulating foam box + Float_t fOuterBoxThickness[3] ; // Thickness of the outer thermo insulating foam box + Float_t fPHOSAngle[4] ; // Position angles of modules + Float_t fPinDiodeSize[3] ; // Size of the PIN Diode + TObjArray * fRotMatrixArray ; // Liste of rotation matrices (one per phos module) + Float_t fSecondUpperPlateThickness ; // Thickness of upper polystyrene foam plate + Float_t fSupportPlateThickness ; // Thickness of the Aluminium support plate + Float_t fUpperCoolingPlateThickness ; // Thickness of the upper cooling plate + Float_t fUpperPlateThickness ; // Thickness of the uper thermo insulating foam plate + Float_t fTextolitBoxSize[3] ; // Size of the Textolit box inside the insulating foam box + Float_t fTextolitBoxThickness[3] ; // Thicknesses of th Textolit box + Float_t fXtlSize[3] ; // PWO4 crystal dimensions + + + ///////////// PPSD (PHOS PRE SHOWER DETECTOR) related parameters + + Float_t fAnodeThickness ; // Thickness of the copper layer which makes the anode + Float_t fAvalancheGap ; // Thickness of the gas in the avalanche stage + Float_t fCathodeThickness ; // Thickeness of composite material ensuring rigidity of cathode + Float_t fCompositeThickness ; // Thickeness of composite material ensuring rigidity of anode + Float_t fConversionGap ; // Thickness of the gas in the conversion stage + Float_t fIPtoTopLidDistance ; // Distance from interaction point to top lid of PPSD + Float_t fLeadConverterThickness ; // Thickness of the Lead converter + Float_t fLeadToMicro2Gap ; // Thickness of the air gap between the Lead and Micromegas 2 + Float_t fLidThickness ; // Thickness of top lid + Float_t fMicromegas1Thickness ; // Thickness of the first downstream Micromegas + Float_t fMicromegas2Thickness ; // Thickness of the second downstream Micromegas + Float_t fMicromegasWallThickness ; // Thickness of the Micromegas leak tight box + Float_t fMicro1ToLeadGap ; // Thickness of the air gap between Micromegas 1 and the Lead + Int_t fNumberOfPadsPhi ; // Number of pads on a micromegas module ; + Int_t fNumberOfPadsZ ; // Number of pads on a micromegas module ; + Int_t fNumberOfModulesPhi ; // Number of micromegas modules in phi + Int_t fNumberOfModulesZ ; // Number of micromegas modules in z + Float_t fPCThickness ; // Thickness of the printed circuit board of the anode + Float_t fPhiDisplacement ; // Phi displacement of micromegas1 with respect to micromegas2 + Float_t fPPSDBoxSize[3] ; // Size of large box which contains PPSD; matches PHOS module size + Float_t fPPSDModuleSize[3] ; // Size of an individual micromegas module + Float_t fZDisplacement ; // Z displacement of micromegas1 with respect to micromegas2 + + static AliPHOSGeometry * fGeom ; // pointer to the unique instance of the singleton + + ClassDef(AliPHOSGeometry,1) // PHOS geometry class , version subatech + +} ; + +#endif // AliPHOSGEOMETRY_H diff --git a/PHOS/AliPHOSLink.cxx b/PHOS/AliPHOSLink.cxx new file mode 100644 index 00000000000..1a78b242721 --- /dev/null +++ b/PHOS/AliPHOSLink.cxx @@ -0,0 +1,51 @@ +/************************************************************************** + * 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. * + **************************************************************************/ + +//_________________________________________________________________________ +// Auxiliary class, used ONLY by AliPHOSTrackSegmentMaker +//*-- Author : Dmitri Peressounko SUBATECH +////////////////////////////////////////////////////////////////////////////// + +// --- ROOT system --- + +// --- Standard library --- + +// --- AliRoot header files --- + +#include "AliPHOSLink.h" + +ClassImp(AliPHOSLink) +//____________________________________________________________________________ + AliPHOSLink::AliPHOSLink(Float_t r, Int_t Emc, Int_t Ppsd) +{ + fR = r ; + fEmcN = Emc ; + fPpsdN = Ppsd ; +} + +//____________________________________________________________________________ +Int_t AliPHOSLink::Compare(TObject * obj) +{ + Int_t rv ; + + AliPHOSLink * link = (AliPHOSLink *) obj ; + + if(this->fR < link->GetR() ) + rv = -1 ; + else + rv = 1 ; + + return rv ; +} diff --git a/PHOS/AliPHOSLink.h b/PHOS/AliPHOSLink.h new file mode 100644 index 00000000000..8bac87f9714 --- /dev/null +++ b/PHOS/AliPHOSLink.h @@ -0,0 +1,47 @@ +#ifndef ALIPHOSLINK_H +#define ALIPHOSLINK_H +/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * See cxx source for full Copyright notice */ + +//////////////////////////////////////////////// +// Short description // +// Version SUBATECH // +// Author Dmitri Peressounko SUBATECH // +// comment: auxiliary class used ONLY // +// by AliPHOSTrackSegmentMaker // +//////////////////////////////////////////////// + +// --- ROOT system --- + +#include "TObject.h" + +// --- Standard library --- + +// --- AliRoot header files --- + +class AliPHOSLink : public TObject{ + +public: + + AliPHOSLink( Float_t r, Int_t EMC, Int_t PPSD) ; // ctor + virtual ~AliPHOSLink(){} // dtor + + Int_t Compare(TObject * obj) ; + Int_t GetEmc(void) { return fEmcN; } + Int_t GetPpsd(void) { return fPpsdN ; } + Float_t GetR(void) { return fR ; } + Bool_t IsSortable() const{ return kTRUE ; } + +private: + + Int_t fEmcN ; // Emc index + Int_t fPpsdN ; // Ppsd index + Float_t fR ; // Distance + +public: + + ClassDef(AliPHOSLink,1) // description , version 1 + +}; + +#endif // AliPHOSLINK_H diff --git a/PHOS/AliPHOSPpsdRecPoint.cxx b/PHOS/AliPHOSPpsdRecPoint.cxx new file mode 100644 index 00000000000..3f75e65131a --- /dev/null +++ b/PHOS/AliPHOSPpsdRecPoint.cxx @@ -0,0 +1,258 @@ +/************************************************************************** + * 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. * + **************************************************************************/ + +//_________________________________________________________________________ +// RecPoint in the PHOS PPSD: a list of AliPHOSDigit's +//*-- Author : Yves Schutz SUBATECH +////////////////////////////////////////////////////////////////////////////// + +// --- ROOT system --- + +// --- Standard library --- + +#include + +// --- AliRoot header files --- + +#include "AliPHOSGeometry.h" +#include "AliPHOSPpsdRecPoint.h" +#include "AliRun.h" + +ClassImp(AliPHOSPpsdRecPoint) + +//____________________________________________________________________________ +AliPHOSPpsdRecPoint::AliPHOSPpsdRecPoint(void) +{ + fMulDigit = 0 ; + AliPHOSGeometry * geom = AliPHOSGeometry::GetInstance() ; + fDelta = geom->GetCrystalSize(0) ; + fLocPos.SetX(1000000.) ; //Local position should be evaluated +} + +//____________________________________________________________________________ +AliPHOSPpsdRecPoint::~AliPHOSPpsdRecPoint(void) // dtor +{ + delete fDigitsList ; +} + +//____________________________________________________________________________ +void AliPHOSPpsdRecPoint::AddDigit(AliDigitNew & digit, Float_t Energy) +{ + // adds a digit to the digits list + // and accumulates the total amplitude and the multiplicity + + + if ( fMulDigit >= fMaxDigit ) { // increase the size of the lists + int * tempo = new ( int[fMaxDigit*=2] ) ; + Int_t index ; + + for ( index = 0 ; index < fMulDigit ; index++ ) + tempo[index] = fDigitsList[index] ; + + delete fDigitsList ; + fDigitsList = tempo ; + } + + fDigitsList[fMulDigit++] = (int) &digit ; + fAmp += Energy ; +} + + + + +//____________________________________________________________________________ +Int_t AliPHOSPpsdRecPoint::Compare(TObject * obj) +{ + Int_t rv ; + + AliPHOSPpsdRecPoint * clu = (AliPHOSPpsdRecPoint *)obj ; + + + Float_t x1 , z1 ; + Float_t x2 , z2 ; + + Int_t PhosMod1 ; + Int_t PhosMod2 ; + + Int_t up1 ; + Int_t up2 ; + + if(GetUp()) // upper layer + up1 = 0 ; + else // lower layer + up1 = 1 ; + + if(clu->GetUp()) // upper layer + up2 = 0 ; + else // lower layer + up2 = 1 ; + + TVector3 PosLoc ; + this->GetLocalPosition(PosLoc) ; + x1 = PosLoc.X() ; + z1 = PosLoc.Z() ; + PhosMod1 = this->GetPHOSMod(); + clu->GetLocalPosition(PosLoc) ; + x2 = PosLoc.X() ; + z2 = PosLoc.Z() ; + PhosMod2 = clu->GetPHOSMod(); + + if(PhosMod1 == PhosMod2 ) { + + if(up1 == up2 ){ + Int_t rowdif = (Int_t)TMath::Ceil(x1/fDelta) - (Int_t) TMath::Ceil(x2/fDelta) ; + + if (rowdif> 0) + rv = -1 ; + else if(rowdif < 0) + rv = 1 ; + else if(z1>z2) + rv = -1 ; + else + rv = 1 ; + } + + else { + + if(up1 < up2 ) // Upper level first (up = True or False, True > False) + rv = 1 ; + else + rv = - 1 ; + } + + } // if PhosMod1 == PhosMod2 + + else { + + if(PhosMod1 < PhosMod2 ) + rv = -1 ; + else + rv = 1 ; + +} + + return rv ; +} + +//____________________________________________________________________________ +void AliPHOSPpsdRecPoint::GetLocalPosition(TVector3 &LPos){ + + if( fLocPos.X() < 1000000.) { //allready evaluated + LPos = fLocPos ; + return ; + } + + Int_t relid[4] ; + + Float_t x = 0. ; + Float_t z = 0. ; + + AliPHOSGeometry * PHOSGeom = (AliPHOSGeometry *) fGeom ; + + AliPHOSDigit * digit ; + Int_t iDigit; + + for(iDigit=0; iDigitAbsToRelNumbering(digit->GetId(), relid) ; + PHOSGeom->RelPosInModule(relid, xi, zi); + x += xi ; + z += zi ; + } + + x /= fMulDigit ; + z /= fMulDigit ; + + fLocPos.SetX(x) ; + fLocPos.SetY(0.) ; + fLocPos.SetZ(z) ; + + LPos = fLocPos ; +} + +//____________________________________________________________________________ +Bool_t AliPHOSPpsdRecPoint::GetUp() +{ + Int_t relid[4] ; + + AliPHOSGeometry * PHOSGeom = (AliPHOSGeometry *) fGeom ; + + AliPHOSDigit *digit = (AliPHOSDigit *)fDigitsList[0] ; + + PHOSGeom->AbsToRelNumbering(digit->GetId(),relid); + Bool_t up ; + + if((Int_t)TMath::Ceil((Float_t)relid[1]/ + (PHOSGeom->GetNumberOfModulesPhi()*PHOSGeom->GetNumberOfModulesZ())-0.0001 ) > 1) + up = kFALSE ; + else + up = kTRUE ; + + return up ; +} + +//____________________________________________________________________________ +void AliPHOSPpsdRecPoint::Print(Option_t * option) +{ + cout << "AliPHOSPpsdRecPoint: " << endl ; + + AliPHOSDigit * digit ; + + Int_t iDigit; + + for(iDigit=0; iDigit 0) return fPHOSMod ; + Int_t relid[4] ; + + AliPHOSDigit * digit ; + digit = (AliPHOSDigit *) fDigitsList[0] ; + AliPHOSGeometry * PHOSGeom = (AliPHOSGeometry *) fGeom ; + + PHOSGeom->AbsToRelNumbering(digit->GetId(), relid) ; + fPHOSMod = relid[0]; + return fPHOSMod ; +} diff --git a/PHOS/AliPHOSRecPoint.h b/PHOS/AliPHOSRecPoint.h new file mode 100644 index 00000000000..38f08cc3a95 --- /dev/null +++ b/PHOS/AliPHOSRecPoint.h @@ -0,0 +1,50 @@ +#ifndef ALIPHOSRECPOINT_H +#define ALIPHOSRECPOINT_H +/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * See cxx source for full Copyright notice */ + +//////////////////////////////////////////////// +// Base class for Reconstructed Points // +// in PHOS and PPSD // +// // +// Author Gines MARTINEZ SUBATECH // +// // +// // +//////////////////////////////////////////////// + +// --- ROOT system --- + +// --- Standard library --- + +#include "assert.h" + +// --- AliRoot header files --- + +#include "AliRecPoint.h" + + +class AliPHOSRecPoint : public AliRecPoint { + +public: + + AliPHOSRecPoint() ; // ctor + virtual ~AliPHOSRecPoint() ; // dtor + virtual void AddDigit(AliDigitNew & digit, Float_t Energy) = 0 ; + virtual Int_t GetPHOSMod(void) ; + virtual Bool_t IsEmc(void){return kTRUE ;} + virtual void Print(Option_t * opt = "void") {} + + virtual Int_t Compare(TObject * obj) { assert(0==1) ; } + virtual Bool_t IsSortable() const { return kTRUE ; } + +protected: + + Int_t fPHOSMod; + +public: + + ClassDef(AliPHOSRecPoint,1) + +}; + +#endif // AliPHOSRECPOINT_H diff --git a/PHOS/AliPHOSReconstructioner.cxx b/PHOS/AliPHOSReconstructioner.cxx new file mode 100644 index 00000000000..98d5b553b5e --- /dev/null +++ b/PHOS/AliPHOSReconstructioner.cxx @@ -0,0 +1,60 @@ +/************************************************************************** + * 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. * + **************************************************************************/ + +//_________________________________________________________________________ +// A brief description of the class +//*-- Author : Gines MARTINEZ SUBATECH +////////////////////////////////////////////////////////////////////////////// + +// --- ROOT system --- + +#include "TClonesArray.h" + +// --- Standard library --- + +// --- AliRoot header files --- + +#include "AliPHOSReconstructioner.h" +#include "AliPHOSClusterizer.h" + +ClassImp(AliPHOSReconstructioner) + + +//____________________________________________________________________________ +AliPHOSReconstructioner::AliPHOSReconstructioner() +{ + // ctor +} + +//____________________________________________________________________________ +AliPHOSReconstructioner::AliPHOSReconstructioner(AliPHOSClusterizer & Clusterizer, AliPHOSTrackSegmentMaker & Tracker) +{ + fClusterizer = &Clusterizer ; + fTrackSegmentMaker = &Tracker ; +} + +//____________________________________________________________________________ +AliPHOSReconstructioner::~AliPHOSReconstructioner() +{ + // dtor +} + +//____________________________________________________________________________ + void AliPHOSReconstructioner:: Make(TClonesArray * dl, RecPointsList * emccl, RecPointsList * ppsdl, TrackSegmentsList * trsl) +{ + fClusterizer->MakeClusters(dl, emccl, ppsdl); + + fTrackSegmentMaker->MakeTrackSegments(dl,emccl,ppsdl,trsl) ; +} diff --git a/PHOS/AliPHOSReconstructioner.h b/PHOS/AliPHOSReconstructioner.h new file mode 100644 index 00000000000..f3fb4ec5c31 --- /dev/null +++ b/PHOS/AliPHOSReconstructioner.h @@ -0,0 +1,50 @@ +#ifndef ALIPHOSRECONSTRUCTIONER_H +#define ALIPHOSRECONSTRUCTIONER_H +/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * See cxx source for full Copyright notice */ + +//////////////////////////////////////////////// +// Algorithme class for the reconstruction // +// // +// Author Gines MARTINEZ SUBATECH // +// // +// // +//////////////////////////////////////////////// + +// --- ROOT system --- + +#include "TObject.h" +#include "AliPHOSClusterizer.h" +#include "AliPHOSTrackSegmentMaker.h" +#include "TClonesArray.h" + +// --- Standard library --- + +// --- AliRoot header files --- + +class AliPHOSReconstructioner : public TObject { + +public: + + AliPHOSReconstructioner(); //ctor + AliPHOSReconstructioner(AliPHOSClusterizer& Clusterizer, AliPHOSTrackSegmentMaker& Tracker); //ctor + ~AliPHOSReconstructioner(); // dtor + + AliPHOSClusterizer * GetClusterizer() { return fClusterizer ; } + void Make(TClonesArray * DL, RecPointsList * emccl, RecPointsList * ppsdl, TrackSegmentsList * trsl) ; // does the job + + +private: + + AliPHOSClusterizer * fClusterizer ; // Method of clusterization + + AliPHOSTrackSegmentMaker * fTrackSegmentMaker ; // + + +public: + + ClassDef(AliPHOSReconstructioner,1) // Reconstruction interface , version 1 + +}; + +#endif // ALIPHOSRECONSTRUCTIONER_H diff --git a/PHOS/AliPHOSTrackSegment.cxx b/PHOS/AliPHOSTrackSegment.cxx new file mode 100644 index 00000000000..f1e871bcb26 --- /dev/null +++ b/PHOS/AliPHOSTrackSegment.cxx @@ -0,0 +1,182 @@ +/************************************************************************** + * 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. * + **************************************************************************/ + +//_________________________________________________________________________ +// class of PHOS Sub Track +//*-- Author : Dmitri Peressounko RRC KI +////////////////////////////////////////////////////////////////////////////// + +// --- ROOT system --- + +#include "TVector3.h" + +// --- Standard library --- + +#include "iostream.h" + +// --- AliRoot header files --- + +#include "AliPHOSTrackSegment.h" +#include "AliPHOSv0.h" + +ClassImp(AliPHOSTrackSegment) + +//____________________________________________________________________________ +AliPHOSTrackSegment::AliPHOSTrackSegment( AliPHOSEmcRecPoint * emc , AliPHOSPpsdRecPoint * ppsdRP1, + AliPHOSPpsdRecPoint * ppsdRP2 ) +{ + if( emc ) + fEmcRecPoint = emc ; + + if( ppsdRP1 ) + fPpsdUp = ppsdRP1 ; + + if( ppsdRP2 ) + fPpsdLow = ppsdRP2 ; + + fCutOnDispersion = 1.5 ; +} + +//____________________________________________________________________________ +AliPHOSTrackSegment::~AliPHOSTrackSegment() // dtor +{ +// fEmcRecPoint.Delete() ; Not Owners !!! +// fPpsdUp.Delete() ; +// fPpsdLow.Delete() ; +} + +//____________________________________________________________________________ +Float_t AliPHOSTrackSegment::GetDistanceInPHOSPlane() +{ + + TVector3 vecEmc ; + fEmcRecPoint->GetLocalPosition(vecEmc) ; + + TVector3 vecPpsd ; + if( fPpsdLow->GetMultiplicity() ) + fPpsdLow->GetLocalPosition(vecPpsd) ; + else { + vecPpsd.SetX(10000.) ; + } + vecEmc -= vecPpsd ; + + Float_t R = vecEmc.Mag();; + + return R ; +} + +//____________________________________________________________________________ +Bool_t AliPHOSTrackSegment::GetMomentumDirection( TVector3 & dir ) +{ + // True if determined + Bool_t ifdeterm = kTRUE ; + + if( fPpsdLow ){ + TMatrix mdummy ; + if( fPpsdUp->GetMultiplicity() ) { // draw line trough 2 points + TVector3 posEmc ; + fEmcRecPoint->GetGlobalPosition(posEmc,mdummy) ; + TVector3 posPpsd ; + fPpsdLow->GetGlobalPosition(posPpsd,mdummy) ; + dir = posEmc - posPpsd ; + dir.SetMag(1.) ; + } + + else { // draw line through 3 pionts + TVector3 posEmc ; + fEmcRecPoint->GetGlobalPosition(posEmc,mdummy) ; + TVector3 posPpsdl ; + fPpsdLow->GetGlobalPosition(posPpsdl,mdummy) ; + TVector3 posPpsdup ; + fPpsdUp->GetGlobalPosition(posPpsdup,mdummy) ; + posPpsdl = 0.5*( posPpsdup+posPpsdl ) ; + dir = posEmc - posPpsdl ; + dir.SetMag(1.) ; + } + + } + else + ifdeterm = kFALSE ; + + return ifdeterm ; +} + +//____________________________________________________________________________ +Int_t AliPHOSTrackSegment::GetPartType() +{ + // Returns 0 - gamma + // 1 - e+, e- + // 2 - neutral hadron + // 3 - charged hadron + + Int_t PartType =0; + if( fPpsdUp ){ // Neutral + + if( fPpsdLow ) // Neutral hadron + PartType = 2 ; + else // Gamma + PartType = 0 ; + + } + + else { // Charged + + if( fEmcRecPoint->GetDispersion() > fCutOnDispersion) + PartType = 3 ; + else + PartType = 1 ; + + } + + return PartType ; + +} + +//____________________________________________________________________________ +void AliPHOSTrackSegment::GetPosition( TVector3 & pos ) +{ + // Returns positions of hits + TMatrix Dummy ; + fEmcRecPoint->GetGlobalPosition(pos, Dummy) ; +} + +//____________________________________________________________________________ +void AliPHOSTrackSegment::Print() +{ + cout << "--------AliPHOSTrackSegment-------- "<GetGlobalPosition( pos, Dummy ) ; + + cout << " position " << pos.X() << " " << pos.Y() << " " << pos.Z() << " Energy " << fEmcRecPoint->GetTotalEnergy() << endl ; + cout << "PPSD Low Reconstructed Point: " << endl; + + if(fPpsdLow){ + fPpsdLow->GetGlobalPosition( pos , Dummy ) ; + cout << " position " << pos.X() << " " << pos.Y() << " " << pos.Z() << endl ; + } + + cout << "PPSD Up Reconstructed Point: " << endl; + + if(fPpsdUp ){ + fPpsdUp->GetGlobalPosition( pos, Dummy ) ; + cout << " position " << pos.X() << " " << pos.Y() << " " << pos.Z() << endl ; + } + +} + diff --git a/PHOS/AliPHOSTrackSegment.h b/PHOS/AliPHOSTrackSegment.h new file mode 100644 index 00000000000..1ae2e5594c3 --- /dev/null +++ b/PHOS/AliPHOSTrackSegment.h @@ -0,0 +1,62 @@ +#ifndef ALIPHOSSUBTRACK_H +#define ALIPHOSSUBTRACK_H +/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * See cxx source for full Copyright notice */ + +///////////////////////////////////////////////// +// Short description // +// Version SUBATECH // +// Author Dmitri Peressounko RRC KI // +// comment: contains pairs (triplets) of // +// EMC+PPSD(+PPSD) clusters, and // +// evaluates particle type, // +// energy, etc // +///////////////////////////////////////////////// + +// --- ROOT system --- + +#include "TObject.h" +#include "TVector3.h" + +// --- Standard library --- + +// --- AliRoot header files --- + +#include "AliPHOSEmcRecPoint.h" +#include "AliPHOSPpsdRecPoint.h" + + + +class AliPHOSTrackSegment : public TObject { + +public: + + AliPHOSTrackSegment() {} ; // ctor + AliPHOSTrackSegment(AliPHOSEmcRecPoint * EmcRecPoint , AliPHOSPpsdRecPoint * PpsdUp, + AliPHOSPpsdRecPoint * PpsdLow ) ; + virtual ~AliPHOSTrackSegment() ; // dtor + + Int_t GetPartType() ; // Returns 0 - gamma, 1 - e+, e- ; 2 - neutral hadron ; 3 - charged hadron + Float_t GetEnergy(){ return fEmcRecPoint->GetTotalEnergy() ;} // Returs energy in EMC + Float_t GetDistanceInPHOSPlane(void) ; // computes in PHOS plane the relative position between EMC and PPSD clusters + Bool_t GetMomentumDirection( TVector3 & dir ) ; // True if determined + void GetPosition( TVector3 & pos ) ; // Returns positions of hits + void Print() ; + void SetDispersionCutOff(Float_t Dcut) {fCutOnDispersion = Dcut ; } + + +private: + + AliPHOSEmcRecPoint * fEmcRecPoint ; + AliPHOSPpsdRecPoint * fPpsdLow ; + AliPHOSPpsdRecPoint * fPpsdUp ; + + Float_t fCutOnDispersion ; + +public: + + ClassDef(AliPHOSTrackSegment,1) // description , version 1 + +}; + +#endif // AliPHOSSUBTRACK_H diff --git a/PHOS/AliPHOSTrackSegmentMaker.cxx b/PHOS/AliPHOSTrackSegmentMaker.cxx new file mode 100644 index 00000000000..bcb4b3819bb --- /dev/null +++ b/PHOS/AliPHOSTrackSegmentMaker.cxx @@ -0,0 +1,53 @@ +/************************************************************************** + * 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. * + **************************************************************************/ + +//_________________________________________________________________________ +// A brief description of the class +//*-- Author : Yves Schutz SUBATECH +////////////////////////////////////////////////////////////////////////////// + +// --- ROOT system --- + +#include "TObjArray.h" +#include "TClonesArray.h" + +// --- Standard library --- + +#include "iostream.h" + +// --- AliRoot header files --- + +#include "AliPHOSTrackSegmentMaker.h" +#include "AliPHOSTrackSegment.h" +#include "AliPHOSLink.h" +#include "AliPHOSv0.h" +#include "AliRun.h" + +ClassImp( AliPHOSTrackSegmentMaker) + + +//____________________________________________________________________________ + AliPHOSTrackSegmentMaker:: AliPHOSTrackSegmentMaker() // ctor +{ + fR0 = 4. ; +} + +//____________________________________________________________________________ +void AliPHOSTrackSegmentMaker::MakeTrackSegments(DigitsList * DL, RecPointsList * emcl, RecPointsList * ppsdl, + TrackSegmentsList * trsl) +{ + +} + diff --git a/PHOS/AliPHOSTrackSegmentMaker.h b/PHOS/AliPHOSTrackSegmentMaker.h new file mode 100644 index 00000000000..8422e73932c --- /dev/null +++ b/PHOS/AliPHOSTrackSegmentMaker.h @@ -0,0 +1,48 @@ +#ifndef ALIPHOSSUBTRACKER_H +#define ALIPHOSSUBTRACKER_H +/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * See cxx source for full Copyright notice */ + +/////////////////////////////////////////////////// +// Subtrackin class for PHOS // +// Version SUBATECH // +// Author Dmitri Peressounko RRC Ki // +// comment: finds pairs of clusters EMC+PPSD // +// performs unfolding. // +/////////////////////////////////////////////////// + +// --- ROOT system --- + +// --- Standard library --- + +// --- AliRoot header files --- + +#include "TObjArray.h" +#include "AliPHOSClusterizer.h" +#include "AliPHOSEmcRecPoint.h" +#include "AliPHOSPpsdRecPoint.h" + +typedef TObjArray TrackSegmentsList ; + +class AliPHOSTrackSegmentMaker : public TObject { + +public: + + AliPHOSTrackSegmentMaker() ; + + virtual ~ AliPHOSTrackSegmentMaker(){} // dtor + + virtual void MakeTrackSegments(DigitsList * DL, RecPointsList * emcl, RecPointsList * ppsdl, TrackSegmentsList * trsl ) ; + // does the job + virtual void SetMaxEmcPpsdDistance(Float_t r){ fR0 = r ;} + +private: + Float_t fR0 ; + +public: + +ClassDef( AliPHOSTrackSegmentMaker,1) // subtracking implementation , version 1 + +}; + +#endif // AliPHOSSUBTRACKER_H diff --git a/PHOS/AliPHOSTrackSegmentMakerv1.cxx b/PHOS/AliPHOSTrackSegmentMakerv1.cxx new file mode 100644 index 00000000000..f5e1f5204c5 --- /dev/null +++ b/PHOS/AliPHOSTrackSegmentMakerv1.cxx @@ -0,0 +1,520 @@ +/************************************************************************** + * 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. * + **************************************************************************/ + +//_________________________________________________________________________ +// A brief description of the class +//*-- Author : Yves Schutz SUBATECH +////////////////////////////////////////////////////////////////////////////// + +// --- ROOT system --- + +#include "TObjArray.h" +#include "TClonesArray.h" +#include "TMinuit.h" + +// --- Standard library --- + +#include "iostream.h" + +// --- AliRoot header files --- + +#include "AliPHOSTrackSegmentMakerv1.h" +#include "AliPHOSTrackSegment.h" +#include "AliPHOSLink.h" +#include "AliPHOSv0.h" +#include "AliRun.h" + +ClassImp( AliPHOSTrackSegmentMakerv1) + + +//____________________________________________________________________________ + AliPHOSTrackSegmentMakerv1:: AliPHOSTrackSegmentMakerv1() // ctor +{ + fR0 = 4. ; + AliPHOSGeometry * geom = AliPHOSGeometry::GetInstance() ; + //clusters are sorted in "rows" and "columns" of width geom->GetCrystalSize(0), + fDelta = fR0 + geom->GetCrystalSize(0) ; +} + +//____________________________________________________________________________ +Bool_t AliPHOSTrackSegmentMakerv1::FindFit(AliPHOSEmcRecPoint * emcRP, int * maxAt, Float_t * maxAtEnergy, + Int_t NPar, Float_t * FitParametres) +{ //Calls TMinuit for fitting cluster with several maxima + + AliPHOSGeometry * geom = AliPHOSGeometry::GetInstance() ; + TMinuit *gMinuit = new TMinuit(NPar); //initialize TMinuit with a maximum of 5 params + gMinuit->SetPrintLevel(-1) ; //No PRIntout + gMinuit->SetFCN(UnfoldingChiSquare ); //To set the address of the minimization function + gMinuit->SetObjectFit(emcRP) ; //To tranfer pointer to UnfoldingChiSquare + + //filling initial values for fit parameters + AliPHOSDigit * digit ; + Int_t ierflg = 0; + Int_t index = 0 ; + Int_t NDigits = (Int_t) NPar / 3 ; + Int_t iDigit ; + for(iDigit = 0 ; iDigit < NDigits ; iDigit ++){ + digit = (AliPHOSDigit *) maxAt[iDigit]; + + Int_t RelId[4] ; + Float_t x ; + Float_t z ; + geom->AbsToRelNumbering(digit->GetId(),RelId) ; + geom->RelPosInModule(RelId,x,z) ; + + Float_t Energy = maxAtEnergy[iDigit] ; + + gMinuit->mnparm(index++, " ", x, 0.1, 0, 0, ierflg) ; + if(ierflg != 0){ + cout << "PHOS Unfolding> Unable to set initial value for fit procedure : x = " << x << endl ; + return kFALSE; + } + gMinuit->mnparm(index++, " ", z, 0.1, 0, 0, ierflg) ; + if(ierflg != 0){ + cout << "PHOS Unfolding> Unable to set initial value for fit procedure : z = " << z << endl ; + return kFALSE; + } + gMinuit->mnparm(index++, " ", Energy , 0.05*Energy, 0., 4.*Energy, ierflg) ; + if(ierflg != 0){ + cout << "PHOS Unfolding> Unable to set initial value for fit procedure : Energy = " << Energy << endl ; + return kFALSE; +} + } + + Double_t p0=0.1; //"Tolerance" Evaluation stops when EDM = 0.0001*p0 ; The number of function call slightly + //depends on it. + Double_t p1 = 1.; + + gMinuit->mnexcm("SET STR", 0, 0, ierflg) ; //force TMinuit to reduce function calls + gMinuit->mnexcm("SET GRA", &p1, 1, ierflg) ; //force TMinuit to use my gradient + gMinuit->SetMaxIterations(5); + gMinuit->mnexcm("SET NOW", 0 , 0, ierflg) ; //No Warnings + gMinuit->mnexcm("MIGRAD", &p0, 0, ierflg); // minimize + if(ierflg == 4){ //Minimum not found + cout << "PHOS Unfolding> Fit not converged, cluster abondoned "<< endl ; + return kFALSE ; + } + for(index = 0; index < NPar; index++){ + Double_t err ; + Double_t val ; + gMinuit->GetParameter(index, val, err) ; //Returns value and error of parameter index + FitParametres[index] = val ; + } + gMinuit->Delete() ; + return kTRUE; +} +//____________________________________________________________________________ +void AliPHOSTrackSegmentMakerv1::FillOneModule(DigitsList * Dl, RecPointsList * emcIn, TObjArray * emcOut, + RecPointsList * ppsdIn, TObjArray * ppsdOutUp, + TObjArray * ppsdOutLow, Int_t &PHOSMod, Int_t & emcStopedAt, + Int_t & ppsdStopedAt) +{// Unfold clusters and fill xxxOut arrais with clusters from ome PHOS modeule + AliPHOSEmcRecPoint * emcRecPoint ; + AliPHOSPpsdRecPoint * ppsdRecPoint ; + Int_t index ; + + Int_t NemcUnfolded = emcIn->GetEntries() ; + for(index = emcStopedAt; index < NemcUnfolded; index++){ + + emcRecPoint = (AliPHOSEmcRecPoint *) (*emcIn)[index] ; + + if(emcRecPoint->GetPHOSMod() != PHOSMod ) + break ; + + Int_t NMultipl = emcRecPoint->GetMultiplicity() ; + int maxAt[NMultipl] ; + Float_t maxAtEnergy[NMultipl] ; + Int_t Nmax = emcRecPoint->GetNumberOfLocalMax(maxAt,maxAtEnergy) ; + + if(Nmax<=1) // if cluster is very flat, so that no prononsed maximum, then Nmax = 0 + emcOut->Add(emcRecPoint) ; + else { + UnfoldClusters(Dl, emcIn, emcRecPoint, Nmax, maxAt, maxAtEnergy, emcOut) ; + emcIn->Remove(emcRecPoint); + emcIn->Compress() ; + NemcUnfolded-- ; + index-- ; + } + } + emcStopedAt = index ; + + for(index = ppsdStopedAt; index < ppsdIn->GetEntries(); index++){ + ppsdRecPoint = (AliPHOSPpsdRecPoint *) (*ppsdIn)[index] ; + if(ppsdRecPoint->GetPHOSMod() != PHOSMod ) break ; + if(ppsdRecPoint->GetUp() ) + ppsdOutUp->Add(ppsdRecPoint) ; + else + ppsdOutLow->Add(ppsdRecPoint) ; + } + ppsdStopedAt = index ; + + PHOSMod ++ ; + + emcOut->Sort() ; + ppsdOutUp->Sort() ; + ppsdOutLow->Sort() ; + +} +//____________________________________________________________________________ +Float_t AliPHOSTrackSegmentMakerv1::GetDistanceInPHOSPlane(AliPHOSEmcRecPoint * EmcClu,AliPHOSPpsdRecPoint * PpsdClu, Bool_t &TooFar) +{ + Float_t R = fR0 ; + + TVector3 vecEmc ; + TVector3 vecPpsd ; + + EmcClu->GetLocalPosition(vecEmc) ; + PpsdClu->GetLocalPosition(vecPpsd) ; + if(EmcClu->GetPHOSMod() == PpsdClu->GetPHOSMod()){ + if(vecPpsd.X() >= vecEmc.X() - fDelta ){ + if(vecPpsd.Z() >= vecEmc.Z() - fDelta ){ + AliPHOSGeometry * geom = AliPHOSGeometry::GetInstance() ; + //Correct to difference in CPV and EMC position due to different distance to center. + //we assume, that particle moves from center + Float_t DCPV = geom->GetIPtoOuterCoverDistance(); + Float_t DEMC = geom->GetIPtoCrystalSurface() ; + DEMC = DEMC / DCPV ; + vecPpsd = DEMC * vecPpsd - vecEmc ; + R = vecPpsd.Mag() ; + } // if zPpsd >= zEmc - fDelta + TooFar = kFALSE ; + } // if xPpsd >= xEmc - fDelta + else + TooFar = kTRUE ; + } + else + TooFar = kTRUE ; + + return R ; +} + +//____________________________________________________________________________ +void AliPHOSTrackSegmentMakerv1::MakeLinks(TObjArray * EmcRecPoints, TObjArray * PpsdRecPointsUp, + TObjArray * PpsdRecPointsLow, TClonesArray * LinkLowArray, + TClonesArray *LinkUpArray) +{ + //Finds distanses (links) between all EMC and PPSD clusters, which are not further from each other than fR0 + + TIter nextEmc(EmcRecPoints) ; + Int_t iEmcClu = 0 ; + + AliPHOSPpsdRecPoint * PpsdLow ; + AliPHOSPpsdRecPoint * PpsdUp ; + AliPHOSEmcRecPoint * EmcClu ; + + Int_t iLinkLow = 0 ; + Int_t iLinkUp = 0 ; + + while( (EmcClu = (AliPHOSEmcRecPoint*)nextEmc() ) ) { + Bool_t TooFar ; + TIter nextPpsdLow(PpsdRecPointsLow ) ; + Int_t iPpsdLow = 0 ; + + while( (PpsdLow = (AliPHOSPpsdRecPoint*)nextPpsdLow() ) ) { + Float_t R = GetDistanceInPHOSPlane(EmcClu, PpsdLow, TooFar) ; + + if(TooFar) + break ; + if(R < fR0) + new( (*LinkLowArray)[iLinkLow++]) AliPHOSLink(R, iEmcClu, iPpsdLow) ; + + iPpsdLow++ ; + + } + + TIter nextPpsdUp(PpsdRecPointsUp ) ; + Int_t iPpsdUp = 0 ; + + while( (PpsdUp = (AliPHOSPpsdRecPoint*)nextPpsdUp() ) ) { + Float_t R = GetDistanceInPHOSPlane(EmcClu, PpsdUp, TooFar) ; + + if(TooFar) + break ; + if(R < fR0) + new( (*LinkUpArray)[iLinkUp++]) AliPHOSLink(R, iEmcClu, iPpsdUp) ; + + iPpsdUp++ ; + + } + + iEmcClu++ ; + + } // while nextEmC + + LinkLowArray->Sort() ; //first links with smallest distances + LinkUpArray->Sort() ; +} + +//____________________________________________________________________________ +void AliPHOSTrackSegmentMakerv1::MakePairs(TObjArray * EmcRecPoints, TObjArray * PpsdRecPointsUp, + TObjArray * PpsdRecPointsLow, TClonesArray * LinkLowArray, + TClonesArray * LinkUpArray, TrackSegmentsList * trsl) +{ // Finds the smallest links and makes pairs of PPSD and EMC clusters with smallest distance + TIter nextLow(LinkLowArray) ; + TIter nextUp(LinkUpArray) ; + + AliPHOSLink * linkLow ; + AliPHOSLink * linkUp ; + + AliPHOSEmcRecPoint * emc ; + AliPHOSPpsdRecPoint * ppsdLow ; + AliPHOSPpsdRecPoint * ppsdUp ; + + while ( (linkLow = (AliPHOSLink *)nextLow() ) ){ + emc = (AliPHOSEmcRecPoint *) EmcRecPoints->At(linkLow->GetEmc()) ; + ppsdLow = (AliPHOSPpsdRecPoint *) PpsdRecPointsLow->At(linkLow->GetPpsd()) ; + if((emc)&&(ppsdLow)){ // RecPoints not removed yet + ppsdUp = 0 ; + + while ( (linkUp = (AliPHOSLink *)nextUp() ) ){ + if(linkLow->GetEmc() == linkUp->GetEmc() ){ + ppsdUp = (AliPHOSPpsdRecPoint *) PpsdRecPointsUp->At(linkUp->GetPpsd()) ; + break ; + } + + } // while nextUp + + nextUp.Reset(); + AliPHOSTrackSegment * subtr = new AliPHOSTrackSegment(emc, ppsdUp, ppsdLow ) ; + trsl->Add(subtr) ; + EmcRecPoints->RemoveAt(linkLow->GetEmc()) ; + PpsdRecPointsLow->RemoveAt(linkLow->GetPpsd()) ; + + if(ppsdUp) + PpsdRecPointsUp->RemoveAt(linkUp->GetPpsd()) ; + + } // if NLocMax + } + + TIter nextEmc(EmcRecPoints) ; + nextEmc.Reset() ; + + while( (emc = (AliPHOSEmcRecPoint*)nextEmc()) ){ //to create pairs if no PpsdLow + ppsdLow = NULL ; + ppsdUp = NULL ; + + while ( (linkUp = (AliPHOSLink *)nextUp() ) ){ + + if(EmcRecPoints->IndexOf(emc) == linkUp->GetEmc() ){ + ppsdUp = (AliPHOSPpsdRecPoint *) PpsdRecPointsUp->At(linkUp->GetPpsd()) ; + break ; + } + + } + AliPHOSTrackSegment * subtr = new AliPHOSTrackSegment(emc, ppsdUp, ppsdLow ) ; + trsl->Add(subtr) ; + + if(ppsdUp) + PpsdRecPointsUp->RemoveAt(linkUp->GetPpsd()) ; + + } + +} + +//____________________________________________________________________________ +void AliPHOSTrackSegmentMakerv1::MakeTrackSegments(DigitsList * DL, RecPointsList * emcl, + RecPointsList * ppsdl, TrackSegmentsList * trsl) +{ + //main function, does the job + + Int_t PHOSMod = 1 ; + Int_t emcStopedAt = 0 ; + Int_t ppsdStopedAt = 0 ; + + TObjArray * EmcRecPoints = new TObjArray(100) ; //these arrays keeps pointers + TObjArray * PpsdRecPointsUp = new TObjArray(100) ; //on RecPoints, which are + TObjArray * PpsdRecPointsLow = new TObjArray(100) ; //kept in TClonesArray's emcl and ppsdl + + + TClonesArray * LinkLowArray = new TClonesArray("AliPHOSLink", 100); + TClonesArray * LinkUpArray = new TClonesArray("AliPHOSLink", 100); + + AliPHOSGeometry * geom = AliPHOSGeometry::GetInstance() ; + + while(PHOSMod <= geom->GetNModules() ){ + + FillOneModule(DL, emcl, EmcRecPoints, ppsdl, PpsdRecPointsUp, PpsdRecPointsLow, PHOSMod , emcStopedAt, ppsdStopedAt) ; + + MakeLinks(EmcRecPoints,PpsdRecPointsUp, PpsdRecPointsLow, LinkLowArray, LinkUpArray) ; + + MakePairs(EmcRecPoints,PpsdRecPointsUp, PpsdRecPointsLow, LinkLowArray, LinkUpArray, trsl) ; + + EmcRecPoints->Clear() ; + PpsdRecPointsUp->Clear() ; + PpsdRecPointsLow->Clear() ; + LinkUpArray->Delete(); + LinkLowArray->Delete(); + } +} + +//____________________________________________________________________________ +Double_t AliPHOSTrackSegmentMakerv1::ShowerShape(Double_t r){ +// If you change this function, change also gradiend evaluation in ChiSquare() + Double_t r4 = r*r*r*r ; + Double_t r295 = TMath::Power(r, 2.95) ; + Double_t shape = TMath::Exp( -r4 * (1. / (2.32 + 0.26 * r4) + 0.0316 / (1 + 0.0652 * r295) ) ) ; + return shape ; +} + +//____________________________________________________________________________ +void AliPHOSTrackSegmentMakerv1::UnfoldClusters(DigitsList * DL, RecPointsList * emcIn, AliPHOSEmcRecPoint * iniEmc, + Int_t Nmax, int * maxAt, Float_t * maxAtEnergy, TObjArray * emcList) +{ + //fits cluster with Nmax overlapping showers + + Int_t NPar = 3 * Nmax ; + Float_t FitParameters[NPar] ; + AliPHOSGeometry * geom = AliPHOSGeometry::GetInstance() ; + + if( !FindFit(iniEmc, maxAt, maxAtEnergy, NPar, FitParameters) ) //Fit failed, return and remove cluster + return ; + + Float_t xDigit ; + Float_t zDigit ; + Int_t RelId[4] ; + + Int_t Ndigits = iniEmc->GetMultiplicity() ; + Float_t xpar ; + Float_t zpar ; + Float_t Epar ; + Float_t Distance ; + Float_t Ratio ; + Float_t Efit[Ndigits] ; + Int_t iparam ; + Int_t iDigit ; + + AliPHOSDigit * digit ; + AliPHOSEmcRecPoint * emcRP ; + int * emcDigits = iniEmc->GetDigitsList() ; + Float_t * emcEnergies = iniEmc->GetEnergiesList() ; + + Int_t iRecPoint = emcIn->GetEntries() ; + + for(iDigit = 0 ; iDigit < Ndigits ; iDigit ++){ + digit = (AliPHOSDigit *) emcDigits[iDigit]; + geom->AbsToRelNumbering(digit->GetId(), RelId) ; + geom->RelPosInModule(RelId, xDigit, zDigit) ; + Efit[iDigit] = 0; + iparam = 0 ; + + while(iparam < NPar ){ + xpar = FitParameters[iparam] ; + zpar = FitParameters[iparam+1] ; + Epar = FitParameters[iparam+2] ; + iparam += 3 ; + Distance = (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar) ; + Distance = TMath::Sqrt(Distance) ; + Efit[iDigit] += Epar * ShowerShape(Distance) ; + } + + } + + iparam = 0 ; + Float_t eDigit ; + + while(iparam < NPar ){ + xpar = FitParameters[iparam] ; + zpar = FitParameters[iparam+1] ; + Epar = FitParameters[iparam+2] ; + iparam += 3 ; + new ((*emcIn)[iRecPoint]) AliPHOSEmcRecPoint( iniEmc->GetLogWeightCut(), iniEmc->GetLocMaxCut() ) ; + emcRP = (AliPHOSEmcRecPoint *) emcIn->At(iRecPoint++); + + for(iDigit = 0 ; iDigit < Ndigits ; iDigit ++){ + digit = (AliPHOSDigit *) emcDigits[iDigit]; + geom->AbsToRelNumbering(digit->GetId(), RelId) ; + geom->RelPosInModule(RelId, xDigit, zDigit) ; + Distance = (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar) ; + Distance = TMath::Sqrt(Distance) ; + Ratio = Epar * ShowerShape(Distance) / Efit[iDigit] ; + eDigit = emcEnergies[iDigit] * Ratio ; + emcRP->AddDigit( *digit,eDigit ) ; + } + + emcList->Add(emcRP) ; + + } +} +//______________________________________________________________________________ +void AliPHOSTrackSegmentMakerv1::UnfoldingChiSquare(Int_t &NPar, Double_t *Grad, Double_t & fret, Double_t *x, Int_t iflag) +{ +// NUmber of paramters, Gradient , Chi squared, parameters, what to do + AliPHOSGeometry * geom = AliPHOSGeometry::GetInstance() ; + AliPHOSTrackSegmentMakerv1 TrS ; + + AliPHOSEmcRecPoint * emcRP = (AliPHOSEmcRecPoint *) gMinuit->GetObjectFit() ; //EmcRecPoint to fit + int * emcDigits = emcRP->GetDigitsList() ; + Float_t * emcEnergies = emcRP->GetEnergiesList() ; + fret = 0. ; + Int_t iparam ; + + if(iflag==2) + for(iparam = 0 ; iparam < NPar ; iparam ++) + Grad[iparam] = 0 ; //Will evaluate gradient + + Double_t Efit ; + + AliPHOSDigit * digit ; + Int_t iDigit = 0 ; + while ( (digit = (AliPHOSDigit *)emcDigits[iDigit] )){ + Int_t RelId[4] ; + Float_t xDigit ; + Float_t zDigit ; + geom->AbsToRelNumbering(digit->GetId(), RelId) ; + geom->RelPosInModule(RelId, xDigit, zDigit) ; + + if(iflag == 2){ //calculate gradient + Int_t iParam = 0 ; + Efit = 0 ; + while(iParam < NPar ){ + Double_t Distance = TMath::Sqrt( (xDigit - x[iParam]) * (xDigit - x[iParam]) + + (zDigit - x[++iParam]) * (zDigit - x[iParam]) ) ; + Efit += x[++iParam] * TrS.ShowerShape(Distance) ; + iParam++ ; + } + Double_t sum = 2. * (Efit - emcEnergies[iDigit]) / emcEnergies[iDigit] ; //Here we assume, that sigma = sqrt(E) + iParam = 0 ; + while(iParam < NPar ){ + Double_t xpar = x[iParam] ; + Double_t zpar = x[iParam+1] ; + Double_t Epar = x[iParam+2] ; + Double_t dr = TMath::Sqrt( (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar) ); + Double_t shape = sum * TrS.ShowerShape(dr) ; + Double_t r4 = dr*dr*dr*dr ; + Double_t r295 = TMath::Power(dr,2.95) ; + Double_t deriv =-4. * dr*dr * ( 2.32 / ( (2.32 + 0.26 * r4) * (2.32 + 0.26 * r4) ) + + 0.0316 * (1. + 0.0171 * r295) / ( ( 1. + 0.0652 * r295) * (1. + 0.0652 * r295) ) ) ; + + Grad[iParam++] += Epar * shape * deriv * (xpar - xDigit) ; // Derivative over x + Grad[iParam++] += Epar * shape * deriv * (zpar - zDigit) ; // Derivative over z + Grad[iParam++] += shape ; // Derivative over energy + } + } + Efit = 0; + iparam = 0 ; + while(iparam < NPar ){ + Double_t xpar = x[iparam] ; + Double_t zpar = x[iparam+1] ; + Double_t Epar = x[iparam+2] ; + iparam += 3 ; + Double_t Distance = (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar) ; + Distance = TMath::Sqrt(Distance) ; + Efit += Epar * TrS.ShowerShape(Distance) ; + } + fret += (Efit-emcEnergies[iDigit])*(Efit-emcEnergies[iDigit])/emcEnergies[iDigit] ; + //Here we assume, that sigma = sqrt(E) + iDigit++ ; + } +} diff --git a/PHOS/AliPHOSTrackSegmentMakerv1.h b/PHOS/AliPHOSTrackSegmentMakerv1.h new file mode 100644 index 00000000000..8826c35afca --- /dev/null +++ b/PHOS/AliPHOSTrackSegmentMakerv1.h @@ -0,0 +1,73 @@ +#ifndef ALIPHOSSUBTRACKERV1_H +#define ALIPHOSSUBTRACKERV1_H +/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * See cxx source for full Copyright notice */ + +/////////////////////////////////////////////////// +// Track Segment Maker class for PHOS // +// Version SUBATECH // +// Author Dmitri Peressounko RRC Ki // +// comment: finds pairs of clusters EMC+PPSD // +// performs unfolding. // +/////////////////////////////////////////////////// + +// --- ROOT system --- + +// --- Standard library --- + +// --- AliRoot header files --- + +#include "TObjArray.h" +#include "AliPHOSClusterizer.h" +#include "AliPHOSEmcRecPoint.h" +#include "AliPHOSPpsdRecPoint.h" +#include "AliPHOSTrackSegmentMaker.h" + + +class AliPHOSTrackSegmentMakerv1 : public AliPHOSTrackSegmentMaker { + +public: + + AliPHOSTrackSegmentMakerv1() ; + virtual ~ AliPHOSTrackSegmentMakerv1(){} // dtor + + Bool_t FindFit(AliPHOSEmcRecPoint * emcRP, int * MaxAt, Float_t * maxAtEnergy, + Int_t NPar, Float_t * FitParametres) ; //Used in UnfoldClusters, calls TMinuit + + void FillOneModule(DigitsList * Dl, RecPointsList * emcIn, TObjArray * emcOut, RecPointsList * ppsdIn, + TObjArray * ppsdOutUp, TObjArray * ppsdOutLow, Int_t &PHOSModule, Int_t & emcStopedAt, + Int_t & ppsdStopedAt) ; // Unfolds clusters and fills temporary arrais + + Float_t GetDistanceInPHOSPlane(AliPHOSEmcRecPoint * EmcClu , AliPHOSPpsdRecPoint * Ppsd , Bool_t & TooFar ) ; // see R0 + + void MakeLinks(TObjArray * EmcRecPoints, TObjArray * PpsdRecPointsUp, TObjArray * PpsdRecPointsLow, + TClonesArray * LinkLowArray, TClonesArray *LinkUpArray) ; //Evaluates distances(links) between EMC and PPSD + + void MakePairs(TObjArray * EmcRecPoints, TObjArray * PpsdRecPointsUp, TObjArray * PpsdRecPointsLow, + TClonesArray * LinkLowArray, TClonesArray * LinkUpArray, TrackSegmentsList * trsl) ; + //Finds pairs(triplets) with smallest link + + void MakeTrackSegments(DigitsList * DL, RecPointsList * emcl, RecPointsList * ppsdl, TrackSegmentsList * trsl ) ; // does the job + + void SetMaxEmcPpsdDistance(Float_t r){ fR0 = r ;} //Radius within which we look for ppsd cluster + + Double_t ShowerShape(Double_t r) ; //Shape of shower used in unfolding + + void UnfoldClusters(DigitsList * DL, RecPointsList * emcIn, AliPHOSEmcRecPoint * iniEmc, Int_t Nmax, + int * maxAt, Float_t * maxAtEnergy, TObjArray * emclist) ; //Unfolds overlaping clusters using TMinuit packadge + + void static UnfoldingChiSquare(Int_t &NPar, Double_t *Grad, Double_t & fret, Double_t *x, Int_t iflag); //used in TMinuit + + +private: + + Float_t fDelta ; // parameter used for sorting + Float_t fR0 ; // Maximal distance between EMC and PPSD clusters of one Track Segment in module plane + +public: + +ClassDef( AliPHOSTrackSegmentMakerv1,1) // track segment maker implementation , version 1 + +}; + +#endif // AliPHOSSUBTRACKERV1_H diff --git a/PHOS/AliPHOSv0.cxx b/PHOS/AliPHOSv0.cxx index 0c2e8627942..f77a423bb7f 100644 --- a/PHOS/AliPHOSv0.cxx +++ b/PHOS/AliPHOSv0.cxx @@ -13,253 +13,1169 @@ * provided "as is" without express or implied warranty. * **************************************************************************/ -/* -$Log$ -Revision 1.9 1999/11/08 07:12:31 fca -Minor corrections thanks to I.Hrivnacova +//_________________________________________________________________________ +// Manager class for PHOS version SUBATECH +//*-- Author : Y. Schutz SUBATECH +////////////////////////////////////////////////////////////////////////////// -Revision 1.8 1999/09/29 09:24:23 fca -Introduction of the Copyright and cvs Log - -*/ - -///////////////////////////////////////////////////////// -// Manager and hits classes for set:PHOS version 0 // -///////////////////////////////////////////////////////// - // --- ROOT system --- -#include "TH1.h" -#include "TRandom.h" -#include "TFile.h" -#include "TTree.h" + #include "TBRIK.h" #include "TNode.h" -// --- galice header files --- +// --- Standard library --- + +#include +#include +#include +#include +#include + +// --- AliRoot header files --- + #include "AliPHOSv0.h" +#include "AliPHOSHit.h" +#include "AliPHOSDigit.h" +#include "AliPHOSReconstructioner.h" #include "AliRun.h" -#include "AliMC.h" +#include "AliConst.h" ClassImp(AliPHOSv0) -//______________________________________________________________________________ - - +//____________________________________________________________________________ AliPHOSv0::AliPHOSv0() { - fIdSens=0; + fNTmpHits = 0 ; + fTmpHits = 0 ; } + +//____________________________________________________________________________ +AliPHOSv0::AliPHOSv0(const char *name, const char *title): + AliPHOS(name,title) +{ + + // We use 2 arrays of hits : + // + // - fHits (the "normal" one), which retains the hits associated with + // the current primary particle being tracked + // (this array is reset after each primary has been tracked). + // + // - fTmpHits, which retains all the hits of the current event. It + // is used for the digitization part. + + fHits = new TClonesArray("AliPHOSHit",100) ; + fDigits = new TClonesArray("AliPHOSDigit",100) ; + fTmpHits= new TClonesArray("AliPHOSHit",100) ; + + assert ( fHits != 0 ) ; + assert ( fDigits != 0 ) ; + assert ( fTmpHits != 0 ) ; + + fNTmpHits = fNhits = 0 ; + + fIshunt = 1 ; // All hits are associated with primary particles -//______________________________________________________________________________ + // gets an instance of the geometry parameters class + fGeom = AliPHOSGeometry::GetInstance(title, "") ; + + if (fGeom->IsInitialized() ) + cout << "AliPHOSv0 : PHOS geometry intialized for " << fGeom->GetName() << endl ; + else + cout << "AliPHOSv0 : PHOS geometry initialization failed !" << endl ; +} +//____________________________________________________________________________ +AliPHOSv0::AliPHOSv0(AliPHOSReconstructioner& Reconstructioner, const char *name, const char *title): + AliPHOS(name,title) +{ + + // We use 2 arrays of hits : + // + // - fHits (the "normal" one), which retains the hits associated with + // the current primary particle being tracked + // (this array is reset after each primary has been tracked). + // + // - fTmpHits, which retains all the hits of the current event. It + // is used for the digitization part. + + fHits = new TClonesArray("AliPHOSHit",100) ; + fDigits = new TClonesArray("AliPHOSDigit",100) ; + fTmpHits= new TClonesArray("AliPHOSHit",100) ; + + assert ( fHits != 0 ) ; + assert ( fDigits != 0 ) ; + assert ( fTmpHits != 0 ) ; + + fNTmpHits = fNhits = 0 ; + + fIshunt = 1 ; // All hits are associated with primary particles + + // gets an instance of the geometry parameters class + fGeom = AliPHOSGeometry::GetInstance(title, "") ; + + if (fGeom->IsInitialized() ) + cout << "AliPHOSv0 : PHOS geometry intialized for " << fGeom->GetName() << endl ; + else + cout << "AliPHOSv0 : PHOS geometry initialization failed !" << endl ; + + // Defining the PHOS Reconstructioner + + fReconstructioner = &Reconstructioner; +} + +//____________________________________________________________________________ +AliPHOSv0::~AliPHOSv0() +{ + delete fHits ; + delete fTmpHits ; + delete fDigits ; +} + +//____________________________________________________________________________ +void AliPHOSv0::AddHit(Int_t track, Int_t Id, Float_t * hits) +{ + Int_t hitCounter ; + TClonesArray <mphits = *fTmpHits; + AliPHOSHit *newHit ; + AliPHOSHit *curHit; + bool already = false ; + + // In any case, fills the fTmpHit TClonesArray (with "accumulated hits") + + newHit = new AliPHOSHit(fIshunt, track, Id, hits) ; + + for ( hitCounter = 0 ; hitCounter < fNTmpHits && !already ; hitCounter++ ) { + curHit = (AliPHOSHit*) ltmphits[hitCounter] ; + if( *curHit == *newHit ) { + *curHit = *curHit + *newHit ; + already = true ; + } + } + + if ( !already ) { + new(ltmphits[fNTmpHits]) AliPHOSHit(*newHit) ; + fNTmpHits++ ; + } + + // Please note that the fTmpHits array must survive up to the + // end of the events, so it does not appear e.g. in ResetHits() ( + // which is called at the end of each primary). -AliPHOSv0::AliPHOSv0(const char *name, const char *title) - : AliPHOS(name, title) + // if (IsTreeSelected('H')) { + // And, if we really want raw hits tree, have the fHits array filled also + // TClonesArray &lhits = *fHits; + // new(lhits[fNhits]) AliPHOSHit(*newHit) ; + // fNhits++ ; + // } + + delete newHit; + +} + + +//____________________________________________________________________________ +void AliPHOSv0::BuildGeometry() { - fIdSens=0; + + this->BuildGeometryforPHOS() ; + if ( ( strcmp(fGeom->GetName(), "GPS2" ) == 0 ) ) + this->BuildGeometryforPPSD() ; + else + cout << "AliPHOSv0::BuildGeometry : no charged particle identification system installed" << endl; + } + +//____________________________________________________________________________ +void AliPHOSv0:: BuildGeometryforPHOS(void) +{ + // Build the PHOS geometry for the ROOT display + + const Int_t kColorPHOS = kRed ; + const Int_t kColorXTAL = kBlue ; + + Double_t const RADDEG = 180.0 / kPI ; + + new TBRIK( "OuterBox", "PHOS box", "void", fGeom->GetOuterBoxSize(0)/2, + fGeom->GetOuterBoxSize(1)/2, + fGeom->GetOuterBoxSize(2)/2 ); + + // Textolit Wall box, position inside PHOS + + new TBRIK( "TextolitBox", "PHOS Textolit box ", "void", fGeom->GetTextolitBoxSize(0)/2, + fGeom->GetTextolitBoxSize(1)/2, + fGeom->GetTextolitBoxSize(2)/2); + + // Polystyrene Foam Plate + + new TBRIK( "UpperFoamPlate", "PHOS Upper foam plate", "void", fGeom->GetTextolitBoxSize(0)/2, + fGeom->GetSecondUpperPlateThickness()/2, + fGeom->GetTextolitBoxSize(2)/2 ) ; + + // Air Filled Box -//___________________________________________ -void AliPHOSv0::Init() + new TBRIK( "AirFilledBox", "PHOS air filled box", "void", fGeom->GetAirFilledBoxSize(0)/2, + fGeom->GetAirFilledBoxSize(1)/2, + fGeom->GetAirFilledBoxSize(2)/2 ); + + // Crystals Box + + Float_t XTL_X = fGeom->GetCrystalSize(0) ; + Float_t XTL_Y = fGeom->GetCrystalSize(1) ; + Float_t XTL_Z = fGeom->GetCrystalSize(2) ; + + Float_t XL = fGeom->GetNPhi() * ( XTL_X + 2 * fGeom->GetGapBetweenCrystals() ) / 2.0 + fGeom->GetModuleBoxThickness() ; + Float_t YL = ( XTL_Y + fGeom->GetCrystalSupportHeight() + fGeom->GetCrystalWrapThickness() + fGeom->GetCrystalHolderThickness() ) / 2.0 + + fGeom->GetModuleBoxThickness() / 2.0 ; + Float_t ZL = fGeom->GetNZ() * ( XTL_Z + 2 * fGeom->GetGapBetweenCrystals() ) / 2.0 + fGeom->GetModuleBoxThickness() ; + + new TBRIK( "CrystalsBox", "PHOS crystals box", "void", XL, YL, ZL ) ; + +// position PHOS into ALICE + + Float_t R = fGeom->GetIPtoOuterCoverDistance() + fGeom->GetOuterBoxSize(1) / 2.0 ; + Int_t number = 988 ; + Float_t pphi = TMath::ATan( fGeom->GetOuterBoxSize(0) / ( 2.0 * fGeom->GetIPtoOuterCoverDistance() ) ) ; + pphi *= RADDEG ; + TNode * Top = gAlice->GetGeometry()->GetNode("alice") ; + + char * nodename = new char[20] ; + char * rotname = new char[20] ; + + for( Int_t i = 1; i <= fGeom->GetNModules(); i++ ) { + Float_t angle = pphi * 2 * ( i - fGeom->GetNModules() / 2.0 - 0.5 ) ; + sprintf(rotname, "%s%d", "rot", number++) ; + new TRotMatrix(rotname, rotname, 90, angle, 90, 90 + angle, 0, 0); + Top->cd(); + sprintf(nodename,"%s%d", "Module", i) ; + Float_t X = R * TMath::Sin( angle / RADDEG ) ; + Float_t Y = -R * TMath::Cos( angle / RADDEG ) ; + TNode * OuterBoxNode = new TNode(nodename, nodename, "OuterBox", X, Y, 0, rotname ) ; + OuterBoxNode->SetLineColor(kColorPHOS) ; + fNodes->Add(OuterBoxNode) ; + OuterBoxNode->cd() ; + // now inside the outer box the textolit box + Y = ( fGeom->GetOuterBoxThickness(1) - fGeom->GetUpperPlateThickness() ) / 2. ; + sprintf(nodename,"%s%d", "TexBox", i) ; + TNode * TextolitBoxNode = new TNode(nodename, nodename, "TextolitBox", 0, Y, 0) ; + TextolitBoxNode->SetLineColor(kColorPHOS) ; + fNodes->Add(TextolitBoxNode) ; + // upper foam plate inside outre box + OuterBoxNode->cd() ; + sprintf(nodename, "%s%d", "UFPlate", i) ; + Y = ( fGeom->GetTextolitBoxSize(1) - fGeom->GetSecondUpperPlateThickness() ) / 2.0 ; + TNode * UpperFoamPlateNode = new TNode(nodename, nodename, "UpperFoamPlate", 0, Y, 0) ; + UpperFoamPlateNode->SetLineColor(kColorPHOS) ; + fNodes->Add(UpperFoamPlateNode) ; + // air filled box inside textolit box (not drawn) + TextolitBoxNode->cd(); + Y = ( fGeom->GetTextolitBoxSize(1) - fGeom->GetAirFilledBoxSize(1) ) / 2.0 - fGeom->GetSecondUpperPlateThickness() ; + sprintf(nodename, "%s%d", "AFBox", i) ; + TNode * AirFilledBoxNode = new TNode(nodename, nodename, "AirFilledBox", 0, Y, 0) ; + fNodes->Add(AirFilledBoxNode) ; + // crystals box inside air filled box + AirFilledBoxNode->cd() ; + Y = fGeom->GetAirFilledBoxSize(1) / 2.0 - YL + - ( fGeom->GetIPtoCrystalSurface() - fGeom->GetIPtoOuterCoverDistance() - fGeom->GetModuleBoxThickness() + - fGeom->GetUpperPlateThickness() - fGeom->GetSecondUpperPlateThickness() ) ; + sprintf(nodename, "%s%d", "XTBox", i) ; + TNode * CrystalsBoxNode = new TNode(nodename, nodename, "CrystalsBox", 0, Y, 0) ; + CrystalsBoxNode->SetLineColor(kColorXTAL) ; + fNodes->Add(CrystalsBoxNode) ; + } +} + +//____________________________________________________________________________ +void AliPHOSv0:: BuildGeometryforPPSD(void) { - fIdSens=gMC->VolId("PXTL"); + // Build the PPSD geometry for the ROOT display + + Double_t const RADDEG = 180.0 / kPI ; + + const Int_t kColorPHOS = kRed ; + const Int_t kColorPPSD = kGreen ; + const Int_t kColorGas = kBlue ; + const Int_t kColorAir = kYellow ; + + // Box for a full PHOS module + + new TBRIK( "PPSDBox", "PPSD box", "void", fGeom->GetPPSDBoxSize(0)/2, + fGeom->GetPPSDBoxSize(1)/2, + fGeom->GetPPSDBoxSize(2)/2 ); + + // Box containing one micromegas module + + new TBRIK( "PPSDModule", "PPSD module", "void", fGeom->GetPPSDModuleSize(0)/2, + fGeom->GetPPSDModuleSize(1)/2, + fGeom->GetPPSDModuleSize(2)/2 ); + // top lid + + new TBRIK ( "TopLid", "Micromegas top lid", "void", fGeom->GetPPSDModuleSize(0)/2, + fGeom->GetLidThickness()/2, + fGeom->GetPPSDModuleSize(2)/2 ) ; + // composite panel (top and bottom) + + new TBRIK ( "TopPanel", "Composite top panel", "void", ( fGeom->GetPPSDModuleSize(0) - fGeom->GetMicromegasWallThickness() )/2, + fGeom->GetCompositeThickness()/2, + ( fGeom->GetPPSDModuleSize(2) - fGeom->GetMicromegasWallThickness() )/2 ) ; + + new TBRIK ( "BottomPanel", "Composite bottom panel", "void", ( fGeom->GetPPSDModuleSize(0) - fGeom->GetMicromegasWallThickness() )/2, + fGeom->GetCompositeThickness()/2, + ( fGeom->GetPPSDModuleSize(2) - fGeom->GetMicromegasWallThickness() )/2 ) ; + // gas gap (conversion and avalanche) + + new TBRIK ( "GasGap", "gas gap", "void", ( fGeom->GetPPSDModuleSize(0) - fGeom->GetMicromegasWallThickness() )/2, + ( fGeom->GetConversionGap() + fGeom->GetAvalancheGap() )/2, + ( fGeom->GetPPSDModuleSize(2) - fGeom->GetMicromegasWallThickness() )/2 ) ; + + // anode and cathode + + new TBRIK ( "Anode", "Anode", "void", ( fGeom->GetPPSDModuleSize(0) - fGeom->GetMicromegasWallThickness() )/2, + fGeom->GetAnodeThickness()/2, + ( fGeom->GetPPSDModuleSize(2) - fGeom->GetMicromegasWallThickness() )/2 ) ; + + new TBRIK ( "Cathode", "Cathode", "void", ( fGeom->GetPPSDModuleSize(0) - fGeom->GetMicromegasWallThickness() )/2, + fGeom->GetCathodeThickness()/2, + ( fGeom->GetPPSDModuleSize(2) - fGeom->GetMicromegasWallThickness() )/2 ) ; + // PC + + new TBRIK ( "PCBoard", "Printed Circuit", "void", ( fGeom->GetPPSDModuleSize(0) - fGeom->GetMicromegasWallThickness() )/2, + fGeom->GetPCThickness()/2, + ( fGeom->GetPPSDModuleSize(2) - fGeom->GetMicromegasWallThickness() )/2 ) ; + // Gap between Lead and top micromegas + + new TBRIK ( "LeadToM", "Air Gap top", "void", fGeom->GetPPSDBoxSize(0)/2, + fGeom->GetMicro1ToLeadGap()/2, + fGeom->GetPPSDBoxSize(2)/2 ) ; + +// Gap between Lead and bottom micromegas + + new TBRIK ( "MToLead", "Air Gap bottom", "void", fGeom->GetPPSDBoxSize(0)/2, + fGeom->GetLeadToMicro2Gap()/2, + fGeom->GetPPSDBoxSize(2)/2 ) ; + // Lead converter + + new TBRIK ( "Lead", "Lead converter", "void", fGeom->GetPPSDBoxSize(0)/2, + fGeom->GetLeadConverterThickness()/2, + fGeom->GetPPSDBoxSize(2)/2 ) ; + + // position PPSD into ALICE + + char * nodename = new char[20] ; + char * rotname = new char[20] ; + + Float_t R = fGeom->GetIPtoTopLidDistance() + fGeom->GetPPSDBoxSize(1) / 2.0 ; + Int_t number = 988 ; + TNode * Top = gAlice->GetGeometry()->GetNode("alice") ; + + for( Int_t i = 1; i <= fGeom->GetNModules(); i++ ) { // the number of PHOS modules + Float_t angle = fGeom->GetPHOSAngle(i) ; + sprintf(rotname, "%s%d", "rotg", number++) ; + new TRotMatrix(rotname, rotname, 90, angle, 90, 90 + angle, 0, 0); + Top->cd(); + sprintf(nodename, "%s%d", "Moduleg", i) ; + Float_t X = R * TMath::Sin( angle / RADDEG ) ; + Float_t Y = -R * TMath::Cos( angle / RADDEG ) ; + TNode * PPSDBoxNode = new TNode(nodename , nodename ,"PPSDBox", X, Y, 0, rotname ) ; + PPSDBoxNode->SetLineColor(kColorPPSD) ; + fNodes->Add(PPSDBoxNode) ; + PPSDBoxNode->cd() ; + // inside the PPSD box: + // 1. fNumberOfModulesPhi x fNumberOfModulesZ top micromegas + X = ( fGeom->GetPPSDBoxSize(0) - fGeom->GetPPSDModuleSize(0) ) / 2. ; + for ( Int_t iphi = 1; iphi <= fGeom->GetNumberOfModulesPhi(); iphi++ ) { // the number of micromegas modules in phi per PHOS module + Float_t Z = ( fGeom->GetPPSDBoxSize(2) - fGeom->GetPPSDModuleSize(2) ) / 2. ; + TNode * Micro1Node ; + for ( Int_t iz = 1; iz <= fGeom->GetNumberOfModulesZ(); iz++ ) { // the number of micromegas modules in z per PHOS module + Y = ( fGeom->GetPPSDBoxSize(1) - fGeom->GetMicromegas1Thickness() ) / 2. ; + sprintf(nodename, "%s%d%d%d", "Mic1", i, iphi, iz) ; + Micro1Node = new TNode(nodename, nodename, "PPSDModule", X, Y, Z) ; + Micro1Node->SetLineColor(kColorPPSD) ; + fNodes->Add(Micro1Node) ; + // inside top micromegas + Micro1Node->cd() ; + // a. top lid + Y = ( fGeom->GetMicromegas1Thickness() - fGeom->GetLidThickness() ) / 2. ; + sprintf(nodename, "%s%d%d%d", "Lid", i, iphi, iz) ; + TNode * TopLidNode = new TNode(nodename, nodename, "TopLid", 0, Y, 0) ; + TopLidNode->SetLineColor(kColorPPSD) ; + fNodes->Add(TopLidNode) ; + // b. composite panel + Y = Y - fGeom->GetLidThickness() / 2. - fGeom->GetCompositeThickness() / 2. ; + sprintf(nodename, "%s%d%d%d", "CompU", i, iphi, iz) ; + TNode * CompUpNode = new TNode(nodename, nodename, "TopPanel", 0, Y, 0) ; + CompUpNode->SetLineColor(kColorPPSD) ; + fNodes->Add(CompUpNode) ; + // c. anode + Y = Y - fGeom->GetCompositeThickness() / 2. - fGeom->GetAnodeThickness() / 2. ; + sprintf(nodename, "%s%d%d%d", "Ano", i, iphi, iz) ; + TNode * AnodeNode = new TNode(nodename, nodename, "Anode", 0, Y, 0) ; + AnodeNode->SetLineColor(kColorPHOS) ; + fNodes->Add(AnodeNode) ; + // d. gas + Y = Y - fGeom->GetAnodeThickness() / 2. - ( fGeom->GetConversionGap() + fGeom->GetAvalancheGap() ) / 2. ; + sprintf(nodename, "%s%d%d%d", "GGap", i, iphi, iz) ; + TNode * GGapNode = new TNode(nodename, nodename, "GasGap", 0, Y, 0) ; + GGapNode->SetLineColor(kColorGas) ; + fNodes->Add(GGapNode) ; + // f. cathode + Y = Y - ( fGeom->GetConversionGap() + fGeom->GetAvalancheGap() ) / 2. - fGeom->GetCathodeThickness() / 2. ; + sprintf(nodename, "%s%d%d%d", "Cathode", i, iphi, iz) ; + TNode * CathodeNode = new TNode(nodename, nodename, "Cathode", 0, Y, 0) ; + CathodeNode->SetLineColor(kColorPHOS) ; + fNodes->Add(CathodeNode) ; + // g. printed circuit + Y = Y - fGeom->GetCathodeThickness() / 2. - fGeom->GetPCThickness() / 2. ; + sprintf(nodename, "%s%d%d%d", "PC", i, iphi, iz) ; + TNode * PCNode = new TNode(nodename, nodename, "PCBoard", 0, Y, 0) ; + PCNode->SetLineColor(kColorPPSD) ; + fNodes->Add(PCNode) ; + // h. composite panel + Y = Y - fGeom->GetPCThickness() / 2. - fGeom->GetCompositeThickness() / 2. ; + sprintf(nodename, "%s%d%d%d", "CompDown", i, iphi, iz) ; + TNode * CompDownNode = new TNode(nodename, nodename, "BottomPanel", 0, Y, 0) ; + CompDownNode->SetLineColor(kColorPPSD) ; + fNodes->Add(CompDownNode) ; + Z = Z - fGeom->GetPPSDModuleSize(2) ; + PPSDBoxNode->cd() ; + } // end of Z module loop + X = X - fGeom->GetPPSDModuleSize(0) ; + PPSDBoxNode->cd() ; + } // end of phi module loop + // 2. air gap + PPSDBoxNode->cd() ; + Y = ( fGeom->GetPPSDBoxSize(1) - 2 * fGeom->GetMicromegas1Thickness() - fGeom->GetMicro1ToLeadGap() ) / 2. ; + sprintf(nodename, "%s%d", "GapUp", i) ; + TNode * GapUpNode = new TNode(nodename, nodename, "LeadToM", 0, Y, 0) ; + GapUpNode->SetLineColor(kColorAir) ; + fNodes->Add(GapUpNode) ; + // 3. lead converter + Y = Y - fGeom->GetMicro1ToLeadGap() / 2. - fGeom->GetLeadConverterThickness() / 2. ; + sprintf(nodename, "%s%d", "LeadC", i) ; + TNode * LeadCNode = new TNode(nodename, nodename, "Lead", 0, Y, 0) ; + LeadCNode->SetLineColor(kColorPPSD) ; + fNodes->Add(LeadCNode) ; + // 4. air gap + Y = Y - fGeom->GetLeadConverterThickness() / 2. - fGeom->GetLeadToMicro2Gap() / 2. ; + sprintf(nodename, "%s%d", "GapDown", i) ; + TNode * GapDownNode = new TNode(nodename, nodename, "MToLead", 0, Y, 0) ; + GapDownNode->SetLineColor(kColorAir) ; + fNodes->Add(GapDownNode) ; + // 5. fNumberOfModulesPhi x fNumberOfModulesZ bottom micromegas + X = ( fGeom->GetPPSDBoxSize(0) - fGeom->GetPPSDModuleSize(0) ) / 2. - fGeom->GetPhiDisplacement() ; + for ( Int_t iphi = 1; iphi <= fGeom->GetNumberOfModulesPhi(); iphi++ ) { + Float_t Z = ( fGeom->GetPPSDBoxSize(2) - fGeom->GetPPSDModuleSize(2) ) / 2. - fGeom->GetZDisplacement() ;; + TNode * Micro2Node ; + for ( Int_t iz = 1; iz <= fGeom->GetNumberOfModulesZ(); iz++ ) { + Y = - ( fGeom->GetPPSDBoxSize(1) - fGeom->GetMicromegas2Thickness() ) / 2. ; + sprintf(nodename, "%s%d%d%d", "Mic2", i, iphi, iz) ; + Micro2Node = new TNode(nodename, nodename, "PPSDModule", X, Y, Z) ; + Micro2Node->SetLineColor(kColorPPSD) ; + fNodes->Add(Micro2Node) ; + // inside bottom micromegas + Micro2Node->cd() ; + // a. top lid + Y = ( fGeom->GetMicromegas2Thickness() - fGeom->GetLidThickness() ) / 2. ; + sprintf(nodename, "%s%d", "Lidb", i) ; + TNode * TopLidbNode = new TNode(nodename, nodename, "TopLid", 0, Y, 0) ; + TopLidbNode->SetLineColor(kColorPPSD) ; + fNodes->Add(TopLidbNode) ; + // b. composite panel + Y = Y - fGeom->GetLidThickness() / 2. - fGeom->GetCompositeThickness() / 2. ; + sprintf(nodename, "%s%d", "CompUb", i) ; + TNode * CompUpbNode = new TNode(nodename, nodename, "TopPanel", 0, Y, 0) ; + CompUpbNode->SetLineColor(kColorPPSD) ; + fNodes->Add(CompUpbNode) ; + // c. anode + Y = Y - fGeom->GetCompositeThickness() / 2. - fGeom->GetAnodeThickness() / 2. ; + sprintf(nodename, "%s%d", "Anob", i) ; + TNode * AnodebNode = new TNode(nodename, nodename, "Anode", 0, Y, 0) ; + AnodebNode->SetLineColor(kColorPPSD) ; + fNodes->Add(AnodebNode) ; + // d. conversion gas + Y = Y - fGeom->GetAnodeThickness() / 2. - ( fGeom->GetConversionGap() + fGeom->GetAvalancheGap() ) / 2. ; + sprintf(nodename, "%s%d", "GGapb", i) ; + TNode * GGapbNode = new TNode(nodename, nodename, "GasGap", 0, Y, 0) ; + GGapbNode->SetLineColor(kColorGas) ; + fNodes->Add(GGapbNode) ; + // f. cathode + Y = Y - ( fGeom->GetConversionGap() + fGeom->GetAvalancheGap() ) / 2. - fGeom->GetCathodeThickness() / 2. ; + sprintf(nodename, "%s%d", "Cathodeb", i) ; + TNode * CathodebNode = new TNode(nodename, nodename, "Cathode", 0, Y, 0) ; + CathodebNode->SetLineColor(kColorPPSD) ; + fNodes->Add(CathodebNode) ; + // g. printed circuit + Y = Y - fGeom->GetCathodeThickness() / 2. - fGeom->GetPCThickness() / 2. ; + sprintf(nodename, "%s%d", "PCb", i) ; + TNode * PCbNode = new TNode(nodename, nodename, "PCBoard", 0, Y, 0) ; + PCbNode->SetLineColor(kColorPPSD) ; + fNodes->Add(PCbNode) ; + // h. composite pane + Y = Y - fGeom->GetPCThickness() / 2. - fGeom->GetCompositeThickness() / 2. ; + sprintf(nodename, "%s%d", "CompDownb", i) ; + TNode * CompDownbNode = new TNode(nodename, nodename, "BottomPanel", 0, Y, 0) ; + CompDownbNode->SetLineColor(kColorPPSD) ; + fNodes->Add(CompDownbNode) ; + Z = Z - fGeom->GetPPSDModuleSize(2) ; + PPSDBoxNode->cd() ; + } // end of Z module loop + X = X - fGeom->GetPPSDModuleSize(0) ; + PPSDBoxNode->cd() ; + } // end of phi module loop + } // PHOS modules + delete rotname ; + delete nodename ; } -//___________________________________________ +//____________________________________________________________________________ void AliPHOSv0::CreateGeometry() { -// *** DEFINITION OF THE -0.25GetModule("PHOS") ; + + if ( PHOS_tmp == NULL ) { + + fprintf(stderr, "PHOS detector not found!\n") ; + return; - Int_t *idtmed = fIdtmed->GetArray()-699; - -// --- Dimensions of PbWO4 crystal --- - const Float_t XTL_X=2.2; - const Float_t XTL_Y=18.; - const Float_t XTL_Z=2.2; -// --- Tyvek wrapper thickness - const Float_t PAP_THICK=0.01; -// --- Polystyrene Foam Outer Cover dimensions --- - const Float_t FOC_X=214.6; - const Float_t FOC_Y=80.; - const Float_t FOC_Z=260.; -// --- Inner AIR volume dimensions --- - const Float_t AIR_X=206.; - const Float_t AIR_Y=66.; - const Float_t AIR_Z=244.; -// --- Tyvek Crystal Block dimensions --- - const Float_t TCB_X=198.; - const Float_t TCB_Y=25.; - const Float_t TCB_Z=234.; -// --- Upper Cooling Plate thickness --- - const Float_t UCP_Y=0.06; -// --- Al Support Plate thickness --- - const Float_t ASP_Y=10.; -//--- Distance from IP to Foam Outer Cover top plate (needs to be 447.) --- - const Float_t FOC_R=467.; -//--- Distance from IP to Crystal Block top Surface (needs to be 460.) --- - const Float_t CBS_R=480.; - -// --- Dimensions of volumes --- - - -// --- Define PHOS box volume, fill with Polystyrene foam --- - dphos[0] = FOC_X/2.; - dphos[1] = FOC_Y/2.; - dphos[2] = FOC_Z/2.; - gMC->Gsvolu("PHOS", "BOX ", idtmed[703], dphos, 3); - -// --- Define air-filled box, place inside PHOS --- - dpair[0] = AIR_X/2.; - dpair[1] = AIR_Y/2.; - dpair[2] = AIR_Z/2.; - gMC->Gsvolu("PAIR", "BOX ", idtmed[798], dpair, 3); - gMC->Gspos("PAIR", 1, "PHOS", 0., 0., 0., 0, "ONLY"); - -// --- Define Upper Cooling Panel --- -// --- place it right behind upper foam plate --- - dpucp[0] = TCB_X/2.; - dpucp[1] = UCP_Y/2.; - dpucp[2] = TCB_Z/2.; - gMC->Gsvolu("PUCP", "BOX ", idtmed[701], dpucp, 3); - yo = (AIR_Y-UCP_Y)/2.; - gMC->Gspos("PUCP", 1, "PAIR", 0., yo, 0., 0, "ONLY"); - -// --- Define Crystal Block, fill with Tyvek, position inside PAIR --- - dptcb[0] = TCB_X/2.; - dptcb[1] = TCB_Y/2.; - dptcb[2] = TCB_Z/2.; - gMC->Gsvolu("PTCB", "BOX ", idtmed[702], dptcb, 3); -// --- Divide PTCB in X and Z directions -- - gMC->Gsdvn("PSEC", "PTCB", 11, 1); - gMC->Gsdvn("PMOD", "PSEC", 13, 3); - gMC->Gsdvn("PSTR", "PMOD", 8, 1); - gMC->Gsdvn("PCEL", "PSTR", 8, 3); - yo = (FOC_Y-TCB_Y)/2. -(CBS_R-FOC_R); - gMC->Gspos("PTCB", 1, "PAIR", 0., yo, 0., 0, "ONLY"); - -// --- Define PbWO4 crystal volume, place inside PCEL --- - dpxtl[0] = XTL_X/2.; - dpxtl[1] = XTL_Y/2.; - dpxtl[2] = XTL_Z/2.; - gMC->Gsvolu("PXTL", "BOX ", idtmed[699], dpxtl, 3); - yo = (TCB_Y-XTL_Y)/2. - PAP_THICK; - gMC->Gspos("PXTL", 1, "PCEL", 0., yo, 0., 0, "ONLY"); - -// --- Define Al Support Plate, position it inside PAIR --- -// --- right beneath PTCB --- - dpasp[0] = AIR_X/2.; - dpasp[1] = ASP_Y/2.; - dpasp[2] = AIR_Z/2.; - gMC->Gsvolu("PASP", "BOX ", idtmed[701], dpasp, 3); - yo = (FOC_Y-ASP_Y)/2. - (CBS_R-FOC_R+TCB_Y); - gMC->Gspos("PASP", 1, "PAIR", 0., yo, 0., 0, "ONLY"); - -// --- Divide in X and Z direction (same way as PTCB) --- - // gMC->Gsdvn("PCMO", "PCSE", 13, 3); - // gMC->Gsdvn("PCST", "PCMO", 8, 1); - // gMC->Gsdvn("PCCE", "PCST", 8, 3); - -// --- Position various PHOS units in ALICE setup --- -// --- PHOS itself first --- - r = FOC_R+FOC_Y/2.; - pphi = TMath::ATan(FOC_X/(2.*FOC_R)); - xp1 = -r * TMath::Sin(pphi * 3.); - yp1 = -r * TMath::Cos(pphi * 3.); - xp2 = -r * TMath::Sin(pphi); - yp2 = -r * TMath::Cos(pphi); - pphi *= 180/kPI; - AliMatrix(idrotm[0], 90.,-3*pphi, 90., 90-3*pphi, 0., 0.); - AliMatrix(idrotm[1], 90., -pphi, 90., 90-pphi, 0., 0.); - AliMatrix(idrotm[2], 90., pphi, 90., 90+pphi, 0., 0.); - AliMatrix(idrotm[3], 90., 3*pphi, 90., 90+3*pphi, 0., 0.); - gMC->Gspos("PHOS", 1, "ALIC", xp1, yp1, 0., idrotm[0], "ONLY"); - gMC->Gspos("PHOS", 2, "ALIC", xp2, yp2, 0., idrotm[1], "ONLY"); - gMC->Gspos("PHOS", 3, "ALIC",-xp2, yp2, 0., idrotm[2], "ONLY"); - gMC->Gspos("PHOS", 4, "ALIC",-xp1, yp1, 0., idrotm[3], "ONLY"); - -// --- Set modules seen without tree for drawings --- - gMC->Gsatt("PMOD", "SEEN", -2); - gMC->Gsatt("PCMO", "SEEN", -2); + } + + // Get pointer to the array containing media indeces + Int_t *IDTMED = fIdtmed->GetArray() - 699 ; + + Float_t BigBox[3] ; + BigBox[0] = fGeom->GetOuterBoxSize(0) / 2.0 ; + BigBox[1] = ( fGeom->GetOuterBoxSize(1) + fGeom->GetPPSDBoxSize(1) ) / 2.0 ; + BigBox[2] = fGeom->GetOuterBoxSize(2) / 2.0 ; + + gMC->Gsvolu("PHOS", "BOX ", IDTMED[798], BigBox, 3) ; + + this->CreateGeometryforPHOS() ; + if ( strcmp( fGeom->GetName(), "GPS2") == 0 ) + this->CreateGeometryforPPSD() ; + else + cout << "AliPHOSv0::CreateGeometry : no charged particle identification system installed" << endl; + + // --- Position PHOS mdules in ALICE setup --- + + Int_t IDROTM[99] ; + Double_t const RADDEG = 180.0 / kPI ; + + for( Int_t i = 1; i <= fGeom->GetNModules(); i++ ) { + + Float_t angle = fGeom->GetPHOSAngle(i) ; + AliMatrix(IDROTM[i-1], 90.0, angle, 90.0, 90.0+angle, 0.0, 0.0) ; + + Float_t R = fGeom->GetIPtoOuterCoverDistance() + ( fGeom->GetOuterBoxSize(1) + fGeom->GetPPSDBoxSize(1) ) / 2.0 ; + + Float_t XP1 = R * TMath::Sin( angle / RADDEG ) ; + Float_t YP1 = -R * TMath::Cos( angle / RADDEG ) ; + + gMC->Gspos("PHOS", i, "ALIC", XP1, YP1, 0.0, IDROTM[i-1], "ONLY") ; + + } // for GetNModules + } + +//____________________________________________________________________________ +void AliPHOSv0::CreateGeometryforPHOS() +{ + // Get pointer to the array containing media indeces + Int_t *IDTMED = fIdtmed->GetArray() - 699 ; + + // --- + // --- Define PHOS box volume, fPUFPill with thermo insulating foam --- + // --- Foam Thermo Insulating outer cover dimensions --- + // --- Put it in BigBox = PHOS + + Float_t DPHOS[3] ; + DPHOS[0] = fGeom->GetOuterBoxSize(0) / 2.0 ; + DPHOS[1] = fGeom->GetOuterBoxSize(1) / 2.0 ; + DPHOS[2] = fGeom->GetOuterBoxSize(2) / 2.0 ; + + gMC->Gsvolu("EMCA", "BOX ", IDTMED[706], DPHOS, 3) ; + + Float_t YO = - fGeom->GetPPSDBoxSize(1) / 2.0 ; + + gMC->Gspos("EMCA", 1, "PHOS", 0.0, YO, 0.0, 0, "ONLY") ; + + // --- + // --- Define Textolit Wall box, position inside EMCA --- + // --- Textolit Wall box dimentions --- + + + Float_t DPTXW[3]; + DPTXW[0] = fGeom->GetTextolitBoxSize(0) / 2.0 ; + DPTXW[1] = fGeom->GetTextolitBoxSize(1) / 2.0 ; + DPTXW[2] = fGeom->GetTextolitBoxSize(2) / 2.0 ; + + gMC->Gsvolu("PTXW", "BOX ", IDTMED[707], DPTXW, 3); + + YO = ( fGeom->GetOuterBoxThickness(1) - fGeom->GetUpperPlateThickness() ) / 2. ; + + gMC->Gspos("PTXW", 1, "EMCA", 0.0, YO, 0.0, 0, "ONLY") ; + + // --- + // --- Define Upper Polystyrene Foam Plate, place inside PTXW --- + // --- immediately below Foam Thermo Insulation Upper plate --- + + // --- Upper Polystyrene Foam plate thickness --- + + Float_t DPUFP[3] ; + DPUFP[0] = fGeom->GetTextolitBoxSize(0) / 2.0 ; + DPUFP[1] = fGeom->GetSecondUpperPlateThickness() / 2. ; + DPUFP[2] = fGeom->GetTextolitBoxSize(2) /2.0 ; + + gMC->Gsvolu("PUFP", "BOX ", IDTMED[703], DPUFP, 3) ; + + YO = ( fGeom->GetTextolitBoxSize(1) - fGeom->GetSecondUpperPlateThickness() ) / 2.0 ; + + gMC->Gspos("PUFP", 1, "PTXW", 0.0, YO, 0.0, 0, "ONLY") ; + + // --- + // --- Define air-filled box, place inside PTXW --- + // --- Inner AIR volume dimensions --- -//___________________________________________ -void AliPHOSv0::CreateMaterials() + + Float_t DPAIR[3] ; + DPAIR[0] = fGeom->GetAirFilledBoxSize(0) / 2.0 ; + DPAIR[1] = fGeom->GetAirFilledBoxSize(1) / 2.0 ; + DPAIR[2] = fGeom->GetAirFilledBoxSize(2) / 2.0 ; + + gMC->Gsvolu("PAIR", "BOX ", IDTMED[798], DPAIR, 3) ; + + YO = ( fGeom->GetTextolitBoxSize(1) - fGeom->GetAirFilledBoxSize(1) ) / 2.0 - fGeom->GetSecondUpperPlateThickness() ; + + gMC->Gspos("PAIR", 1, "PTXW", 0.0, YO, 0.0, 0, "ONLY") ; + +// --- Dimensions of PbWO4 crystal --- + + Float_t XTL_X = fGeom->GetCrystalSize(0) ; + Float_t XTL_Y = fGeom->GetCrystalSize(1) ; + Float_t XTL_Z = fGeom->GetCrystalSize(2) ; + + Float_t DPTCB[3] ; + DPTCB[0] = fGeom->GetNPhi() * ( XTL_X + 2 * fGeom->GetGapBetweenCrystals() ) / 2.0 + fGeom->GetModuleBoxThickness() ; + DPTCB[1] = ( XTL_Y + fGeom->GetCrystalSupportHeight() + fGeom->GetCrystalWrapThickness() + fGeom->GetCrystalHolderThickness() ) / 2.0 + + fGeom->GetModuleBoxThickness() / 2.0 ; + DPTCB[2] = fGeom->GetNZ() * ( XTL_Z + 2 * fGeom->GetGapBetweenCrystals() ) / 2.0 + fGeom->GetModuleBoxThickness() ; + + gMC->Gsvolu("PTCB", "BOX ", IDTMED[706], DPTCB, 3) ; + + YO = fGeom->GetAirFilledBoxSize(1) / 2.0 - DPTCB[1] + - ( fGeom->GetIPtoCrystalSurface() - fGeom->GetIPtoOuterCoverDistance() - fGeom->GetModuleBoxThickness() + - fGeom->GetUpperPlateThickness() - fGeom->GetSecondUpperPlateThickness() ) ; + + gMC->Gspos("PTCB", 1, "PAIR", 0.0, YO, 0.0, 0, "ONLY") ; + + // --- + // --- Define Crystal BLock filled with air, position it inside PTCB --- + Float_t DPCBL[3] ; + + DPCBL[0] = fGeom->GetNPhi() * ( XTL_X + 2 * fGeom->GetGapBetweenCrystals() ) / 2.0 ; + DPCBL[1] = ( XTL_Y + fGeom->GetCrystalSupportHeight() + fGeom->GetCrystalWrapThickness() + fGeom->GetCrystalHolderThickness() ) / 2.0 ; + DPCBL[2] = fGeom->GetNZ() * ( XTL_Z + 2 * fGeom->GetGapBetweenCrystals() ) / 2.0 ; + + gMC->Gsvolu("PCBL", "BOX ", IDTMED[798], DPCBL, 3) ; + + // --- Divide PCBL in X (phi) and Z directions -- + gMC->Gsdvn("PROW", "PCBL", Int_t (fGeom->GetNPhi()), 1) ; + gMC->Gsdvn("PCEL", "PROW", Int_t (fGeom->GetNZ()), 3) ; + + YO = -fGeom->GetModuleBoxThickness() / 2.0 ; + + gMC->Gspos("PCBL", 1, "PTCB", 0.0, YO, 0.0, 0, "ONLY") ; + + // --- + // --- Define STeel (actually, it's titanium) Cover volume, place inside PCEL + Float_t DPSTC[3] ; + + DPSTC[0] = ( XTL_X + 2 * fGeom->GetCrystalWrapThickness() ) / 2.0 ; + DPSTC[1] = ( XTL_Y + fGeom->GetCrystalSupportHeight() + fGeom->GetCrystalWrapThickness() + fGeom->GetCrystalHolderThickness() ) / 2.0 ; + DPSTC[2] = ( XTL_Z + 2 * fGeom->GetCrystalWrapThickness() + 2 * fGeom->GetCrystalHolderThickness() ) / 2.0 ; + + gMC->Gsvolu("PSTC", "BOX ", IDTMED[704], DPSTC, 3) ; + + gMC->Gspos("PSTC", 1, "PCEL", 0.0, 0.0, 0.0, 0, "ONLY") ; + + // --- + // --- Define Tyvek volume, place inside PSTC --- + Float_t DPPAP[3] ; + + DPPAP[0] = XTL_X / 2.0 + fGeom->GetCrystalWrapThickness() ; + DPPAP[1] = ( XTL_Y + fGeom->GetCrystalSupportHeight() + fGeom->GetCrystalWrapThickness() ) / 2.0 ; + DPPAP[2] = XTL_Z / 2.0 + fGeom->GetCrystalWrapThickness() ; + + gMC->Gsvolu("PPAP", "BOX ", IDTMED[702], DPPAP, 3) ; + + YO = ( XTL_Y + fGeom->GetCrystalSupportHeight() + fGeom->GetCrystalWrapThickness() ) / 2.0 + - ( XTL_Y + fGeom->GetCrystalSupportHeight() + fGeom->GetCrystalWrapThickness() + fGeom->GetCrystalHolderThickness() ) / 2.0 ; + + gMC->Gspos("PPAP", 1, "PSTC", 0.0, YO, 0.0, 0, "ONLY") ; + + // --- + // --- Define PbWO4 crystal volume, place inside PPAP --- + Float_t DPXTL[3] ; + + DPXTL[0] = XTL_X / 2.0 ; + DPXTL[1] = XTL_Y / 2.0 ; + DPXTL[2] = XTL_Z / 2.0 ; + + gMC->Gsvolu("PXTL", "BOX ", IDTMED[699], DPXTL, 3) ; + + YO = ( XTL_Y + fGeom->GetCrystalSupportHeight() + fGeom->GetCrystalWrapThickness() ) / 2.0 - XTL_Y / 2.0 - fGeom->GetCrystalWrapThickness() ; + + gMC->Gspos("PXTL", 1, "PPAP", 0.0, YO, 0.0, 0, "ONLY") ; + + // --- + // --- Define crystal support volume, place inside PPAP --- + Float_t DPSUP[3] ; + + DPSUP[0] = XTL_X / 2.0 + fGeom->GetCrystalWrapThickness() ; + DPSUP[1] = fGeom->GetCrystalSupportHeight() / 2.0 ; + DPSUP[2] = XTL_Z / 2.0 + fGeom->GetCrystalWrapThickness() ; + + gMC->Gsvolu("PSUP", "BOX ", IDTMED[798], DPSUP, 3) ; + + YO = fGeom->GetCrystalSupportHeight() / 2.0 - ( XTL_Y + fGeom->GetCrystalSupportHeight() + fGeom->GetCrystalWrapThickness() ) / 2.0 ; + + gMC->Gspos("PSUP", 1, "PPAP", 0.0, YO, 0.0, 0, "ONLY") ; + + // --- + // --- Define PIN-diode volume and position it inside crystal support --- + // --- right behind PbWO4 crystal + + // --- PIN-diode dimensions --- + + + Float_t DPPIN[3] ; + DPPIN[0] = fGeom->GetPinDiodeSize(0) / 2.0 ; + DPPIN[1] = fGeom->GetPinDiodeSize(1) / 2.0 ; + DPPIN[2] = fGeom->GetPinDiodeSize(2) / 2.0 ; + + gMC->Gsvolu("PPIN", "BOX ", IDTMED[705], DPPIN, 3) ; + + YO = fGeom->GetCrystalSupportHeight() / 2.0 - fGeom->GetPinDiodeSize(1) / 2.0 ; + + gMC->Gspos("PPIN", 1, "PSUP", 0.0, YO, 0.0, 0, "ONLY") ; + + // --- + // --- Define Upper Cooling Panel, place it on top of PTCB --- + Float_t DPUCP[3] ; + // --- Upper Cooling Plate thickness --- + + DPUCP[0] = DPTCB[0] ; + DPUCP[1] = fGeom->GetUpperCoolingPlateThickness() ; + DPUCP[2] = DPTCB[2] ; + + gMC->Gsvolu("PUCP", "BOX ", IDTMED[701], DPUCP,3) ; + + YO = ( fGeom->GetAirFilledBoxSize(1) - fGeom->GetUpperCoolingPlateThickness() ) / 2. + - ( fGeom->GetIPtoCrystalSurface() - fGeom->GetIPtoOuterCoverDistance() - fGeom->GetModuleBoxThickness() + - fGeom->GetUpperPlateThickness() - fGeom->GetSecondUpperPlateThickness() - fGeom->GetUpperCoolingPlateThickness() ) ; + + gMC->Gspos("PUCP", 1, "PAIR", 0.0, YO, 0.0, 0, "ONLY") ; + + // --- + // --- Define Al Support Plate, position it inside PAIR --- + // --- right beneath PTCB --- + // --- Al Support Plate thickness --- + + Float_t DPASP[3] ; + DPASP[0] = fGeom->GetAirFilledBoxSize(0) / 2.0 ; + DPASP[1] = fGeom->GetSupportPlateThickness() / 2.0 ; + DPASP[2] = fGeom->GetAirFilledBoxSize(2) / 2.0 ; + + gMC->Gsvolu("PASP", "BOX ", IDTMED[701], DPASP, 3) ; + + YO = ( fGeom->GetAirFilledBoxSize(1) - fGeom->GetSupportPlateThickness() ) / 2. + - ( fGeom->GetIPtoCrystalSurface() - fGeom->GetIPtoOuterCoverDistance() + - fGeom->GetUpperPlateThickness() - fGeom->GetSecondUpperPlateThickness() + DPCBL[1] * 2 ) ; + + gMC->Gspos("PASP", 1, "PAIR", 0.0, YO, 0.0, 0, "ONLY") ; + + // --- + // --- Define Thermo Insulating Plate, position it inside PAIR --- + // --- right beneath PASP --- + // --- Lower Thermo Insulating Plate thickness --- + + Float_t DPTIP[3] ; + DPTIP[0] = fGeom->GetAirFilledBoxSize(0) / 2.0 ; + DPTIP[1] = fGeom->GetLowerThermoPlateThickness() / 2.0 ; + DPTIP[2] = fGeom->GetAirFilledBoxSize(2) / 2.0 ; + + gMC->Gsvolu("PTIP", "BOX ", IDTMED[706], DPTIP, 3) ; + + YO = ( fGeom->GetAirFilledBoxSize(1) - fGeom->GetLowerThermoPlateThickness() ) / 2. + - ( fGeom->GetIPtoCrystalSurface() - fGeom->GetIPtoOuterCoverDistance() - fGeom->GetUpperPlateThickness() + - fGeom->GetSecondUpperPlateThickness() + DPCBL[1] * 2 + fGeom->GetSupportPlateThickness() ) ; + + gMC->Gspos("PTIP", 1, "PAIR", 0.0, YO, 0.0, 0, "ONLY") ; + + // --- + // --- Define Textolit Plate, position it inside PAIR --- + // --- right beneath PTIP --- + // --- Lower Textolit Plate thickness --- + + Float_t DPTXP[3] ; + DPTXP[0] = fGeom->GetAirFilledBoxSize(0) / 2.0 ; + DPTXP[1] = fGeom->GetLowerTextolitPlateThickness() / 2.0 ; + DPTXP[2] = fGeom->GetAirFilledBoxSize(2) / 2.0 ; + + gMC->Gsvolu("PTXP", "BOX ", IDTMED[707], DPTXP, 3) ; + + YO = ( fGeom->GetAirFilledBoxSize(1) - fGeom->GetLowerTextolitPlateThickness() ) / 2. + - ( fGeom->GetIPtoCrystalSurface() - fGeom->GetIPtoOuterCoverDistance() - fGeom->GetUpperPlateThickness() + - fGeom->GetSecondUpperPlateThickness() + DPCBL[1] * 2 + fGeom->GetSupportPlateThickness() + + fGeom->GetLowerThermoPlateThickness() ) ; + + gMC->Gspos("PTXP", 1, "PAIR", 0.0, YO, 0.0, 0, "ONLY") ; + +} + +//____________________________________________________________________________ +void AliPHOSv0::CreateGeometryforPPSD() { -// *** DEFINITION OF AVAILABLE PHOS MATERIALS *** -// ORIGIN : NICK VAN EIJNDHOVEN + // Get pointer to the array containing media indeces + Int_t *IDTMED = fIdtmed->GetArray() - 699 ; + + // The box containing all PPSD's for one PHOS module filled with air + Float_t PPSD[3] ; + PPSD[0] = fGeom->GetPPSDBoxSize(0) / 2.0 ; + PPSD[1] = fGeom->GetPPSDBoxSize(1) / 2.0 ; + PPSD[2] = fGeom->GetPPSDBoxSize(2) / 2.0 ; - Int_t ISXFLD = gAlice->Field()->Integ(); - Float_t SXMGMX = gAlice->Field()->Max(); - -// --- The PbWO4 crystals --- - Float_t ax[3] = { 207.19,183.85,16. }; - Float_t zx[3] = { 82.,74.,8. }; - Float_t wx[3] = { 1.,1.,4. }; - Float_t dx = 8.28; -// --- The polysterene scintillator (CH) --- - Float_t ap[2] = { 12.011,1.00794 }; - Float_t zp[2] = { 6.,1. }; - Float_t wp[2] = { 1.,1. }; - Float_t dp = 1.032; -// --- Tyvek (CnH2n) - Float_t at[2] = { 12.011,1.00794 }; - Float_t zt[2] = { 6.,1. }; - Float_t wt[2] = { 1.,2. }; - Float_t dt = .331; -// --- Polystyrene foam --- - Float_t af[2] = { 12.011,1.00794 }; - Float_t zf[2] = { 6.,1. }; - Float_t wf[2] = { 1.,1. }; - Float_t df = .3; - - Int_t *idtmed = fIdtmed->GetArray()-699; + gMC->Gsvolu("PPSD", "BOX ", IDTMED[798], PPSD, 3) ; + + Float_t YO = fGeom->GetOuterBoxSize(1) / 2.0 ; + + gMC->Gspos("PPSD", 1, "PHOS", 0.0, YO, 0.0, 0, "ONLY") ; + + // Now we build a micromegas module + // The box containing the whole module filled with epoxy (FR4) + + Float_t MPPSD[3] ; + MPPSD[0] = fGeom->GetPPSDModuleSize(0) / 2.0 ; + MPPSD[1] = fGeom->GetPPSDModuleSize(1) / 2.0 ; + MPPSD[2] = fGeom->GetPPSDModuleSize(2) / 2.0 ; + + gMC->Gsvolu("MPPS", "BOX ", IDTMED[708], MPPSD, 3) ; + + // Inside MPPSD : + // 1. The Top Lid made of epoxy (FR4) + + Float_t TLPPSD[3] ; + TLPPSD[0] = fGeom->GetPPSDModuleSize(0) / 2.0 ; + TLPPSD[1] = fGeom->GetLidThickness() / 2.0 ; + TLPPSD[2] = fGeom->GetPPSDModuleSize(2) / 2.0 ; + + gMC->Gsvolu("TLPS", "BOX ", IDTMED[708], TLPPSD, 3) ; + + Float_t Y0 = ( fGeom->GetMicromegas1Thickness() - fGeom->GetLidThickness() ) / 2. ; + + gMC->Gspos("TLPS", 1, "MPPS", 0.0, Y0, 0.0, 0, "ONLY") ; + + // 2. the upper panel made of composite material + + Float_t UPPPSD[3] ; + UPPPSD[0] = ( fGeom->GetPPSDModuleSize(0) - fGeom->GetMicromegasWallThickness() ) / 2.0 ; + UPPPSD[1] = fGeom->GetCompositeThickness() / 2.0 ; + UPPPSD[2] = ( fGeom->GetPPSDModuleSize(2) - fGeom->GetMicromegasWallThickness() ) / 2.0 ; + + gMC->Gsvolu("UPPS", "BOX ", IDTMED[709], UPPPSD, 3) ; + + Y0 = Y0 - fGeom->GetLidThickness() / 2. - fGeom->GetCompositeThickness() / 2. ; + + gMC->Gspos("UPPS", 1, "MPPS", 0.0, Y0, 0.0, 0, "ONLY") ; + + // 3. the anode made of Copper + + Float_t ANPPSD[3] ; + ANPPSD[0] = ( fGeom->GetPPSDModuleSize(0) - fGeom->GetMicromegasWallThickness() ) / 2.0 ; + ANPPSD[1] = fGeom->GetAnodeThickness() / 2.0 ; + ANPPSD[2] = ( fGeom->GetPPSDModuleSize(2) - fGeom->GetMicromegasWallThickness() ) / 2.0 ; + + gMC->Gsvolu("ANPS", "BOX ", IDTMED[710], ANPPSD, 3) ; + + Y0 = Y0 - fGeom->GetCompositeThickness() / 2. - fGeom->GetAnodeThickness() / 2. ; + + gMC->Gspos("ANPS", 1, "MPPS", 0.0, Y0, 0.0, 0, "ONLY") ; + + // 4. the conversion gap + avalanche gap filled with gas + + Float_t GGPPSD[3] ; + GGPPSD[0] = ( fGeom->GetPPSDModuleSize(0) - fGeom->GetMicromegasWallThickness() ) / 2.0 ; + GGPPSD[1] = ( fGeom->GetConversionGap() + fGeom->GetAvalancheGap() ) / 2.0 ; + GGPPSD[2] = ( fGeom->GetPPSDModuleSize(2) - fGeom->GetMicromegasWallThickness() ) / 2.0 ; + + gMC->Gsvolu("GGPS", "BOX ", IDTMED[715], GGPPSD, 3) ; + + // --- Divide GGPP in X (phi) and Z directions -- + gMC->Gsdvn("GROW", "GGPS", fGeom->GetNumberOfPadsPhi(), 1) ; + gMC->Gsdvn("GCEL", "GROW", fGeom->GetNumberOfPadsZ() , 3) ; + + Y0 = Y0 - fGeom->GetAnodeThickness() / 2. - ( fGeom->GetConversionGap() + fGeom->GetAvalancheGap() ) / 2. ; + + gMC->Gspos("GGPS", 1, "MPPS", 0.0, Y0, 0.0, 0, "ONLY") ; + + + // 6. the cathode made of Copper + + Float_t CAPPSD[3] ; + CAPPSD[0] = ( fGeom->GetPPSDModuleSize(0) - fGeom->GetMicromegasWallThickness() ) / 2.0 ; + CAPPSD[1] = fGeom->GetCathodeThickness() / 2.0 ; + CAPPSD[2] = ( fGeom->GetPPSDModuleSize(2) - fGeom->GetMicromegasWallThickness() ) / 2.0 ; + + gMC->Gsvolu("CAPS", "BOX ", IDTMED[710], CAPPSD, 3) ; + + Y0 = Y0 - ( fGeom->GetAvalancheGap() + fGeom->GetAvalancheGap() ) / 2. - fGeom->GetCathodeThickness() / 2. ; + + gMC->Gspos("CAPS", 1, "MPPS", 0.0, Y0, 0.0, 0, "ONLY") ; + + // 7. the printed circuit made of G10 + + Float_t PCPPSD[3] ; + PCPPSD[0] = ( fGeom->GetPPSDModuleSize(0) - fGeom->GetMicromegasWallThickness() ) / 2,.0 ; + PCPPSD[1] = fGeom->GetPCThickness() / 2.0 ; + PCPPSD[2] = ( fGeom->GetPPSDModuleSize(2) - fGeom->GetMicromegasWallThickness() ) / 2.0 ; + + gMC->Gsvolu("PCPS", "BOX ", IDTMED[711], CAPPSD, 3) ; + + Y0 = Y0 - fGeom->GetCathodeThickness() / 2. - fGeom->GetPCThickness() / 2. ; + + gMC->Gspos("PCPS", 1, "MPPS", 0.0, Y0, 0.0, 0, "ONLY") ; + + // 8. the lower panel made of composite material + + Float_t LPPPSD[3] ; + LPPPSD[0] = ( fGeom->GetPPSDModuleSize(0) - fGeom->GetMicromegasWallThickness() ) / 2.0 ; + LPPPSD[1] = fGeom->GetCompositeThickness() / 2.0 ; + LPPPSD[2] = ( fGeom->GetPPSDModuleSize(2) - fGeom->GetMicromegasWallThickness() ) / 2.0 ; + + gMC->Gsvolu("LPPS", "BOX ", IDTMED[709], LPPPSD, 3) ; + + Y0 = Y0 - fGeom->GetPCThickness() / 2. - fGeom->GetCompositeThickness() / 2. ; + + gMC->Gspos("LPPS", 1, "MPPS", 0.0, Y0, 0.0, 0, "ONLY") ; + + // Position the fNumberOfModulesPhi x fNumberOfModulesZ modules (MPPSD) inside PPSD to cover a PHOS module + // the top and bottom one's (which are assumed identical) : + + Float_t Yt = ( fGeom->GetPPSDBoxSize(1) - fGeom->GetMicromegas1Thickness() ) / 2. ; + Float_t Yb = - ( fGeom->GetPPSDBoxSize(1) - fGeom->GetMicromegas2Thickness() ) / 2. ; + + Int_t CopyNumbertop = 0 ; + Int_t CopyNumberbot = fGeom->GetNumberOfModulesPhi() * fGeom->GetNumberOfModulesZ() ; + + Float_t X = ( fGeom->GetPPSDBoxSize(0) - fGeom->GetPPSDModuleSize(0) ) / 2. ; + + for ( Int_t iphi = 1; iphi <= fGeom->GetNumberOfModulesPhi(); iphi++ ) { // the number of micromegas modules in phi per PHOS module + Float_t Z = ( fGeom->GetPPSDBoxSize(2) - fGeom->GetPPSDModuleSize(2) ) / 2. ; + + for ( Int_t iz = 1; iz <= fGeom->GetNumberOfModulesZ(); iz++ ) { // the number of micromegas modules in z per PHOS module + gMC->Gspos("MPPS", ++CopyNumbertop, "PPSD", X, Yt, Z, 0, "ONLY") ; + gMC->Gspos("MPPS", ++CopyNumberbot, "PPSD", X, Yb, Z, 0, "ONLY") ; + Z = Z - fGeom->GetPPSDModuleSize(2) ; + } // end of Z module loop + X = X - fGeom->GetPPSDModuleSize(0) ; + } // end of phi module loop + + // The Lead converter between two air gaps + // 1. Upper air gap + + Float_t UAPPSD[3] ; + UAPPSD[0] = fGeom->GetPPSDBoxSize(0) / 2.0 ; + UAPPSD[1] = fGeom->GetMicro1ToLeadGap() / 2.0 ; + UAPPSD[2] = fGeom->GetPPSDBoxSize(2) / 2.0 ; + + gMC->Gsvolu("UAPPSD", "BOX ", IDTMED[798], UAPPSD, 3) ; + + Y0 = ( fGeom->GetPPSDBoxSize(1) - 2 * fGeom->GetMicromegas1Thickness() - fGeom->GetMicro1ToLeadGap() ) / 2. ; + + gMC->Gspos("UAPPSD", 1, "PPSD", 0.0, Y0, 0.0, 0, "ONLY") ; + + // 2. Lead converter + + Float_t LCPPSD[3] ; + LCPPSD[0] = fGeom->GetPPSDBoxSize(0) / 2.0 ; + LCPPSD[1] = fGeom->GetLeadConverterThickness() / 2.0 ; + LCPPSD[2] = fGeom->GetPPSDBoxSize(2) / 2.0 ; + + gMC->Gsvolu("LCPPSD", "BOX ", IDTMED[712], LCPPSD, 3) ; + + Y0 = Y0 - fGeom->GetMicro1ToLeadGap() / 2. - fGeom->GetLeadConverterThickness() / 2. ; + + gMC->Gspos("LCPPSD", 1, "PPSD", 0.0, Y0, 0.0, 0, "ONLY") ; + + // 3. Lower air gap + + Float_t LAPPSD[3] ; + LAPPSD[0] = fGeom->GetPPSDBoxSize(0) / 2.0 ; + LAPPSD[1] = fGeom->GetLeadToMicro2Gap() / 2.0 ; + LAPPSD[2] = fGeom->GetPPSDBoxSize(2) / 2.0 ; + + gMC->Gsvolu("LAPPSD", "BOX ", IDTMED[798], LAPPSD, 3) ; - AliMixture( 0, "PbWO4$", ax, zx, dx, -3, wx); - AliMixture( 1, "Polystyrene$", ap, zp, dp, -2, wp); - AliMaterial(2, "Al$", 26.98, 13., 2.7, 8.9, 999.); -// --- Absorption length^ is ignored --- - AliMixture( 3, "Tyvek$", at, zt, dt, -2, wt); - AliMixture( 4, "Foam$", af, zf, df, -2, wf); - AliMaterial(9, "Air$", 14.61, 7.3, .001205, 30420., 67500); - - AliMedium(0, "PHOS Xtal $", 0, 1, ISXFLD, SXMGMX, 10., .1, .1, .1, .1); - AliMedium(2, "Al parts $", 2, 0, ISXFLD, SXMGMX, 10., .1, .1, .001, .001); - AliMedium(3, "Tyvek wrapper$", 3, 0, ISXFLD, SXMGMX, 10., .1, .1, .001, .001); - AliMedium(4, "Polyst. foam $", 4, 0, ISXFLD, SXMGMX, 10., .1, .1, .1, .1); - AliMedium(99, "Air $", 9, 0, ISXFLD, SXMGMX, 10., 1., .1, .1, 10.); - -// --- Generate explicitly delta rays in aluminium parts --- - gMC->Gstpar(idtmed[701], "LOSS", 3.); - gMC->Gstpar(idtmed[701], "DRAY", 1.); + Y0 = Y0 - fGeom->GetLeadConverterThickness() / 2. - fGeom->GetLeadToMicro2Gap() / 2. ; + + gMC->Gspos("LAPPSD", 1, "PPSD", 0.0, Y0, 0.0, 0, "ONLY") ; + } -void AliPHOSv0::StepManager() +//___________________________________________________________________________ +Int_t AliPHOSv0::Digitize(Float_t Energy){ + Float_t fB = 10000000. ; + Float_t fA = 0. ; + Int_t chan = Int_t(fA + Energy*fB ) ; + return chan ; +} +//___________________________________________________________________________ +void AliPHOSv0::FinishEvent() { + cout << "//_____________________________________________________" << endl ; + cout << " AliPHOSv0::FinishEvent() -- Starting digitalization" << endl ; + Int_t i ; + TClonesArray &lDigits = *fDigits ; + AliPHOSHit * Hit ; + AliPHOSDigit * Digit ; - TClonesArray &lhits = *fHits; - TLorentzVector p; - Int_t copy, i; - Int_t vol[5]; - Float_t hits[4]; - if(gMC->CurrentVolID(copy) == fIdSens) { - // - //We are in the sensitive volume - for(i=0;i<4;i++) { - gMC->CurrentVolOffID(i+1,copy); - vol[4-i]=copy; - } - gMC->CurrentVolOffID(7,copy); - vol[0]=copy; - gMC->TrackPosition(p); - for(i=0;i<3;++i) hits[i]=p[i]; - hits[3]=gMC->Edep(); - new(lhits[fNhits++]) AliPHOShit(fIshunt,gAlice->CurrentTrack(),vol,hits); + for ( i = 0 ; i < fNTmpHits ; i++ ) { + Hit = (AliPHOSHit*)fTmpHits->At(i) ; + assert (Hit!=0) ; + Digit = new AliPHOSDigit(Hit->GetId(),Digitize(Hit->GetEnergy())) ; + new(lDigits[fNdigits]) AliPHOSDigit(* Digit) ; + fNdigits++; delete Digit ; } + + // Reset the array of all the "accumulated hits" of this event. + fNTmpHits = 0 ; + fTmpHits->Delete(); } + +//____________________________________________________________________________ +void AliPHOSv0::Init(void) +{ + + Int_t i; + + printf("\n"); + for(i=0;i<35;i++) printf("*"); + printf(" PHOS_INIT "); + for(i=0;i<35;i++) printf("*"); + printf("\n"); + + // Here the PHOS initialisation code (if any!) + + for(i=0;i<80;i++) printf("*"); + printf("\n"); + +} + +//___________________________________________________________________________ +void AliPHOSv0::MakeBranch(Option_t* opt) +{ + // + // Create a new branch in the current Root Tree + // The branch of fHits is automatically split + // + AliDetector::MakeBranch(opt) ; + + char branchname[10]; + sprintf(branchname,"%s",GetName()); + char *D = strstr(opt,"D"); + + if (fDigits && gAlice->TreeD() && D) { + gAlice->TreeD()->Branch(branchname,&fDigits, fBufferSize); + printf("* AliPHOS::MakeBranch * Making Branch %s for digits\n",branchname); + } +} +//_____________________________________________________________________________ + +void AliPHOSv0::Reconstruction(AliPHOSReconstructioner& Reconstructioner) +{ + fReconstructioner = &Reconstructioner; + if (fEmcClusters) + { fEmcClusters->Delete();} + else + { fEmcClusters= new TClonesArray("AliPHOSEmcRecPoint", 100); } ; + + if (fPpsdClusters) + { fPpsdClusters->Delete(); } + else + { fPpsdClusters = new TClonesArray("AliPHOSPpsdRecPoint", 100) ;} + + if (fTrackSegments) + { fTrackSegments->Delete(); } + else + { fTrackSegments = new TObjArray(100) ;} + + fReconstructioner->Make(fDigits, fEmcClusters, fPpsdClusters, fTrackSegments); + +} + +//____________________________________________________________________________ +void AliPHOSv0::StepManager(void) +{ + Int_t RelId[4] ; // (box, layer, row, column) indices + Float_t xyze[4] ; // position wrt MRS and energy deposited + TLorentzVector pos ; + Int_t copy; + + TString name = fGeom->GetName() ; + + if ( name == "GPS2" ) { // the CPV is a PPSD + if( gMC->CurrentVolID(copy) == gMC->VolId("GCEL") ) + // if( strcmp ( gMC->CurrentVolName(), "GCEL" ) == 0 ) // We are inside a gas cell + { + gMC->TrackPosition(pos) ; + xyze[0] = pos[0] ; + xyze[1] = pos[1] ; + xyze[2] = pos[2] ; + xyze[3] = gMC->Edep() ; + + if ( xyze[3] != 0 ) { // there is deposited energy + gMC->CurrentVolOffID(5, RelId[0]) ; // get the PHOS Module number + gMC->CurrentVolOffID(3, RelId[1]) ; // get the Micromegas Module number + // 1-> Geom->GetNumberOfModulesPhi() * fGeom->GetNumberOfModulesZ() upper + // > fGeom->GetNumberOfModulesPhi() * fGeom->GetNumberOfModulesZ() lower + gMC->CurrentVolOffID(1, RelId[2]) ; // get the row number of the cell + gMC->CurrentVolID(RelId[3]) ; // get the column number + + // get the absolute Id number + + Int_t AbsId ; + fGeom->RelToAbsNumbering(RelId,AbsId) ; + + // add current hit to the hit list + AddHit(gAlice->CurrentTrack(), AbsId, xyze); + + } // there is deposited energy + } // We are inside the gas of the CPV + } // GPS2 configuration + + if(gMC->CurrentVolID(copy) == gMC->VolId("PXTL") ) + // if( strcmp ( gMC->CurrentVolName(), "PXTL" ) == 0 ) { // We are inside a PWO crystal + { + gMC->TrackPosition(pos) ; + xyze[0] = pos[0] ; + xyze[1] = pos[1] ; + xyze[2] = pos[2] ; + xyze[3] = gMC->Edep() ; + + if ( xyze[3] != 0 ) { + gMC->CurrentVolOffID(10, RelId[0]) ; // get the PHOS module number ; + RelId[1] = 0 ; // means PW04 + gMC->CurrentVolOffID(4, RelId[2]) ; // get the row number inside the module + gMC->CurrentVolOffID(3, RelId[3]) ; // get the cell number inside the module + + // get the absolute Id number + + Int_t AbsId ; + fGeom->RelToAbsNumbering(RelId,AbsId) ; + + // add current hit to the hit list + + AddHit(gAlice->CurrentTrack(), AbsId, xyze); + + } // there is deposited energy + } // we are inside a PHOS Xtal +} + diff --git a/PHOS/AliPHOSv0.h b/PHOS/AliPHOSv0.h index 67cb8750d5d..74c0c4f0f73 100644 --- a/PHOS/AliPHOSv0.h +++ b/PHOS/AliPHOSv0.h @@ -1,36 +1,100 @@ -#ifndef PHOSv0_H -#define PHOSv0_H +#ifndef ALIPHOSXXX_H +#define ALIPHOSXXX_H /* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * * See cxx source for full Copyright notice */ -/* $Id$ */ +//////////////////////////////////////////////// +// Short description // +// Author SUBATECH // +// comment // +// // +//////////////////////////////////////////////// -///////////////////////////////////////////////////////// -// Manager and hits classes for set:PHOS version 0 // -///////////////////////////////////////////////////////// +// --- ROOT system --- -// --- galice header files --- +// --- Standard library --- + +// --- AliRoot header files --- + +class AliPHOSxxx { + +public: + + virtual ~AliPHOSxxx() ; // dtor + +private: + +ClassDef(AliPHOSxxx,1) // description , version 1 + +}; + +#endif // AliPHOSXXX_H +//-*-C++-*- +#ifndef ALIPHOSV4_H +#define ALIPHOSV4_H +//////////////////////////////////////////////// +// Manager class for PHOS // +// Version SUBATECH // +// Author Y. Schutz SUBATECH // +// geometry parametrized for any // +// shape of modules // +//////////////////////////////////////////////// + +// --- ROOT system --- +#include "TClonesArray.h" + +// --- AliRoot header files --- #include "AliPHOS.h" - -class AliPHOSv0 : public AliPHOS -{ +#include "AliPHOSGeometry.h" +#include "AliPHOSReconstructioner.h" +#include "AliPHOSTrackSegmentMaker.h" + +class AliPHOSv0 : public AliPHOS { + +public: + + AliPHOSv0(void) ; + AliPHOSv0(const char *name, const char *title="") ; + AliPHOSv0(AliPHOSReconstructioner& Reconstructioner, const char *name, const char *title="") ; + virtual ~AliPHOSv0(void) ; + + virtual void AddHit( Int_t track, Int_t id, Float_t *hits ) ; // adds a pre-digitilized hit to the hit tree + virtual void BuildGeometry(void) ; // creates the geometry for the ROOT display + void BuildGeometryforPHOS(void) ; // creates the PHOS geometry for the ROOT display + void BuildGeometryforPPSD(void) ; // creates the PPSD geometry for the ROOT display + virtual void CreateGeometry(void) ; // creates the geometry for GEANT + void CreateGeometryforPHOS(void) ; // creates the PHOS geometry for GEANT + void CreateGeometryforPPSD(void) ; // creates the PPSD geometry for GEANT + Int_t Digitize(Float_t Energy); + RecPointsList* EmcClusters() {return fEmcClusters;} // gets TClonesArray of cluster in the crystals + void FinishEvent(void) ; // makes the digits from the hits + virtual void Init(void) ; // does nothing + void MakeBranch(Option_t* opt) ; + RecPointsList* PpsdClusters() {return fPpsdClusters;} // gets TClonesArray of clusters in the PPSD + void Reconstruction(AliPHOSReconstructioner& Reconstructioner) ; + void ResetClusters(){} ; + void SetReconstructioner(AliPHOSReconstructioner& Reconstructioner) {fReconstructioner = &Reconstructioner;} // + virtual void StepManager(void) ; // does the tracking through PHOS and a preliminary digitalization + TObjArray * TrackSegments(){return fTrackSegments ;} + // inlines + + virtual AliPHOSGeometry * GetGeometry() { return fGeom ; } + Int_t IsVersion(void) const { return 4 ; } - protected: +private: - Int_t fIdSens; //Sensitive volume for phos + AliPHOSGeometry * fGeom ; // geometry definition + RecPointsList * fEmcClusters; //! + Int_t fNTmpHits ; //! used internally for digitalization (!=do not stream) + RecPointsList * fPpsdClusters; //! + TObjArray * fTrackSegments ;//! + TClonesArray * fTmpHits ; //! idem + AliPHOSReconstructioner * fReconstructioner ; // Reconstrutioner of the PHOS event: Clusterization and subtracking procedures + AliPHOSTrackSegmentMaker * fTrackSegmentMaker ; +public: - public: - AliPHOSv0(); - AliPHOSv0(const char *name, const char *title); - virtual ~AliPHOSv0(){} - virtual void CreateGeometry(); - virtual void CreateMaterials(); - virtual void Init(); - virtual Int_t IsVersion() const {return 0;} - virtual void StepManager(); + ClassDef(AliPHOSv0,1) // PHOS main class , version subatech - ClassDef(AliPHOSv0,1) //Hits manager for set:PHOS version 0 }; - -#endif +#endif // AliPHOSV4_H diff --git a/PHOS/Makefile b/PHOS/Makefile index 5c2dc8111eb..69bd9293cdb 100644 --- a/PHOS/Makefile +++ b/PHOS/Makefile @@ -9,7 +9,14 @@ PACKAGE = PHOS # C++ sources -SRCS = AliPHOS.cxx AliPHOSv0.cxx AliPHOSv1.cxx AliPHOSv2.cxx +SRCS = AliPHOS.cxx AliPHOSv0.cxx AliPHOSHit.cxx \ + AliPHOSGeometry.cxx \ + AliPHOSDigit.cxx \ + AliPHOSRecPoint.cxx AliPHOSEmcRecPoint.cxx AliPHOSPpsdRecPoint.cxx \ + AliPHOSClusterizer.cxx AliPHOSClusterizerv1.cxx AliPHOSLink.cxx \ + AliPHOSReconstructioner.cxx AliPHOSTrackSegment.cxx \ + AliPHOSTrackSegmentMaker.cxx AliPHOSTrackSegmentMakerv1.cxx + # C++ Headers @@ -21,12 +28,6 @@ DICT = PHOSCint.cxx DICTH = $(DICT:.cxx=.h) DICTO = $(patsubst %.cxx,tgt_$(ALICE_TARGET)/%.o,$(DICT)) -# FORTRAN sources - -# FORTRAN Objects - -# C Objects - # C++ Objects OBJS = $(patsubst %.cxx,tgt_$(ALICE_TARGET)/%.o,$(SRCS)) $(DICTO) @@ -34,12 +35,10 @@ OBJS = $(patsubst %.cxx,tgt_$(ALICE_TARGET)/%.o,$(SRCS)) $(DICTO) # C++ compilation flags CXXFLAGS = $(CXXOPTS) -I$(ROOTSYS)/include -I. -I$(ALICE_ROOT)/include/ - -# FORTRAN compilation flags -ALSRCS = $(SRCS) +ALSRCS = $(SRCS) $(SHSRCS) $(RCSRCS) $(DUSRCS) dummies.c -ALOBJS = $(OBJS) +ALOBJS = $(SHOBJS) $(RCOBJS) $(DUOBJS) ##### COMMANDS ##### @@ -53,7 +52,7 @@ $(DICT): $(HDRS) depend: $(SRCS) -TOCLEAN = $(OBJS) *Cint.cxx *Cint.h +TOCLEAN = $(OBJS) *Cint.cxx *Cint.h ############################### General Macros ################################ diff --git a/PHOS/PHOSLinkDef.h b/PHOS/PHOSLinkDef.h index 434f9fe09bd..cdd80376a39 100644 --- a/PHOS/PHOSLinkDef.h +++ b/PHOS/PHOSLinkDef.h @@ -1,20 +1,22 @@ #ifdef __CINT__ -/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * - * See cxx source for full Copyright notice */ - -/* $Id$ */ - + #pragma link off all globals; #pragma link off all classes; #pragma link off all functions; - -#pragma link C++ class AliPHOS; -#pragma link C++ class AliPHOSv0; -#pragma link C++ class AliPHOSv1; -#pragma link C++ class AliPHOSv2; -#pragma link C++ class AliPHOShitv2; -#pragma link C++ class AliPHOSCradle; -#pragma link C++ class AliPHOSgamma; -#pragma link C++ class AliPHOShit; +#pragma link C++ class AliPHOS ; +#pragma link C++ class AliPHOSClusterizer ; +#pragma link C++ class AliPHOSClusterizerv1 ; +#pragma link C++ class AliPHOSDigit ; +#pragma link C++ class AliPHOSEmcRecPoint- ; +#pragma link C++ class AliPHOSGeometry ; +#pragma link C++ class AliPHOSHit ; +#pragma link C++ class AliPHOSLink ; +#pragma link C++ class AliPHOSPpsdRecPoint ; +#pragma link C++ class AliPHOSReconstructioner ; +#pragma link C++ class AliPHOSRecPoint ; +#pragma link C++ class AliPHOSv0 ; +#pragma link C++ class AliPHOSTrackSegment ; +#pragma link C++ class AliPHOSTrackSegmentMaker ; +#pragma link C++ class AliPHOSTrackSegmentMakerv1 ; #endif diff --git a/PHOS/canevas.cxx b/PHOS/canevas.cxx new file mode 100644 index 00000000000..43548cd8671 --- /dev/null +++ b/PHOS/canevas.cxx @@ -0,0 +1,44 @@ +/************************************************************************** + * 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. * + **************************************************************************/ + +/* +//_________________________________________________________________________ +// Short description +//*-- Author : SUBATECH +////////////////////////////////////////////////////////////////////////////// + +// --- ROOT system --- + +// --- Standard library --- + +// --- AliRoot header files --- + +#include "AliPHOSxxx.h" + +ClassImp(AliPHOSxxx) + + +//____________________________________________________________________________ +void AliPHOSxxx::AliPHOSxxx() +{ + // ctor +} + +//____________________________________________________________________________ +void ~AliPHOSxxx::AliPHOSxxx() +{ + // dtor +} + diff --git a/PHOS/canevas.h b/PHOS/canevas.h new file mode 100644 index 00000000000..1ddaa6a2029 --- /dev/null +++ b/PHOS/canevas.h @@ -0,0 +1,31 @@ +#ifndef ALIPHOSXXX_H +#define ALIPHOSXXX_H +/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * See cxx source for full Copyright notice */ + +//////////////////////////////////////////////// +// Short description // +// Author SUBATECH // +// comment // +// // +//////////////////////////////////////////////// + +// --- ROOT system --- + +// --- Standard library --- + +// --- AliRoot header files --- + +class AliPHOSxxx { + +public: + + virtual ~AliPHOSxxx() ; // dtor + +private: + +ClassDef(AliPHOSxxx,1) // description , version 1 + +}; + +#endif // AliPHOSXXX_H diff --git a/PHOS/testTrackSegment.C b/PHOS/testTrackSegment.C new file mode 100644 index 00000000000..142cca84dd2 --- /dev/null +++ b/PHOS/testTrackSegment.C @@ -0,0 +1,95 @@ +void testTrackSegment (Int_t evt = 0) +{ + +//========= Dynamically link some shared libs + if (gClassTable->GetID("AliRun") < 0) + { + gROOT->LoadMacro("loadlibs.C"); + loadlibs(); + } + +//========== Opening galice.root file + TFile * file = new TFile("galice.root"); + +//========== Get AliRun object from file + gAlice = (AliRun*) file->Get("gAlice"); + +//=========== Gets the PHOS object and associated geometry from the file + AliPHOSv0 * PHOS = (AliPHOSv0 *)gAlice->GetDetector("PHOS"); + AliPHOSGeometry * Geom = AliPHOSGeometry::GetInstance(PHOS->GetGeometry()->GetName(),PHOS->GetGeometry()->GetTitle()); + +//========== Creates the Clusterizer + AliPHOSClusterizerv1 clusterizer; + clusterizer.SetEmcEnergyThreshold(0.01) ; + clusterizer.SetEmcClusteringThreshold(0.1) ; + clusterizer.SetPpsdEnergyThreshold(0.0000001) ; + clusterizer.SetPpsdClusteringThreshold(0.0000002) ; + clusterizer.SetLocalMaxCut(0.03) ; + clusterizer.SetCalibrationParameters(0.,0.0000001) ; + + +//========== Creates the track segment maker + AliPHOSTrackSegmentMakererv1 tracksegmentmaker ; + +//========== Creates the Reconstructioner + AliPHOSReconstructioner Reconstructioner(clusterizer,tracksegmentmaker); + + //=========== Connects the various Tree's for evt + gAlice->GetEvent(evt); + //=========== Gets the Digit TTree + gAlice->TreeD()->GetEvent(0) ; + //=========== Gets the number of entries in the Digits array + Int_t nId = PHOS->Digits()->GetEntries(); + // printf("AnaPHOSv0.C> Number of entries in the Digit array is %d \n",nId); + + //=========== Do the reconstruction + + AliPHOSDigit * digit ; + TIter next(PHOS->Digits()) ; + Float_t Etot=0 ; + while((digit = (AliPHOSDigit *)next())) Etot+=clusterizer.Calibrate(digit->GetAmp()) ; + cout <<"Found " << nId << " digits with total energy " << Etot << endl ; + + + + PHOS->Reconstruction(Reconstructioner); + + //================Make checks=========================== + AliPHOSDigit * digit ; + TIter next(PHOS->Digits()) ; + Float_t Etot=0 ; + while((digit = (AliPHOSDigit *)next())) Etot+=clusterizer.Calibrate(digit->GetAmp() ); + cout <<"Found " << nId << " digits with total energy " << Etot << endl ; + + TClonesArray * EmcRP = PHOS->EmcClusters() ; + Etot = 0.; + TIter nextemc(EmcRP) ; + AliPHOSEmcRecPoint * emc ; + while((emc = (AliPHOSEmcRecPoint *)nextemc())) { + Etot+=emc->GetTotalEnergy() ; + TVector3 pos ; + emc->GetLocalPosition(pos ) ; + TMatrix Dummy ; + emc->GetGlobalPosition(pos,Dummy) ; + } + + cout << "Found " << EmcRP->GetEntries() << " EMC Clusters with total energy "<PpsdClusters() ; + cout << "Found " << PpsdRP->GetEntries() << " Ppsd Clusters " << endl ; + + TObjArray * trsegl = PHOS->TrackSegments() ; + AliPHOSTrackSegment trseg ; + + Int_t NTrackSegments = trsegl->GetEntries() ; + Int_t index ; + Etot = 0 ; + for(index = 0; index < NTrackSegments ; index++){ + trseg = (AliPHOSTrackSegment * )trsegl->At(index) ; + Etot+= trseg->GetEnergy() ; + trseg->Print() ; + } + cout << "Found " << trsegl->GetEntries() << " Track segments with total energy "<< Etot << endl ; + +} + + -- 2.39.3