* 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 <stdio.h>
-#include <string.h>
-#include <stdlib.h>
-#include <iostream.h>
-// --- 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
-/*
-<img src="picts/aliphos.gif">
-*/
-//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;i<GetCradlesAmount();i++) {
-
- int n = fCradles->GetEntries();
- 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; m<fCradles->GetEntries(); 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)<cradle->GetNz ()*cradle->GetCellSideSize()/2
- && TMath::Abs(y)<cradle->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; i<fCradles->GetEntries(); 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; i<fCradles->GetEntries(); 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; i<fCradles->GetEntries(); 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<fNz; x++ )
- printf(" %4d|",x+1);
- printf("\n");
-
- for( int y=fNphi-1; y>=0; y-- )
- {
- printf("%3d|",y+1);
- for( int x=0; x<fNz; x++ )
- printf("%6d",(int)(cr->fCellEnergy.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; i<p.GetEntries(); i++ )
- ((AliPHOSgamma*)(p[i]))->Print();
- }
-
- 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; i<p.GetEntries(); i++ )
- ((AliPHOSgamma*)(p[i]))->Print();
- }
-}
+ // 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<fNphi; y++ )
- for( int x=0; x<fNz; x++ )
- {
- const int n = fCellEnergy.GetBin(x,y); // Bin number
- static TRandom r;
-
- Float_t E_old=fCellEnergy.GetBinContent(n), E_new=E_old;
-
- if( NULL!=Stochastic )
- E_new = r.Gaus(E_old,sqrt(E_old)*GetDistortedValue(Stochastic,n));
-
- if( NULL!=Calibration )
- E_new *= GetDistortedValue(Calibration,n);
-
- if( NULL!=Noise )
- E_new += GetDistortedValue(Noise,n);
-
- fCellEnergy.SetBinContent(n,E_new);
- }
-}
+ // Various Aluminium parts made of Al -> 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; y<Ny; y++ )
- for( int x=0; x<Nx; x++ )
- {
- const int n = h->GetBin(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<fNphi; y++ )
- for( int x=0; x<fNz; x++ )
- {
- UInt_t n = fCellEnergy.GetBin(x,y);
- UInt_t signal = (int) (fCellEnergy.GetBinContent(n)/signal_step);
- if( signal>=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<rcgamma.recons_gammas_amount; i++ )
- {
-// new (GetGammasReconstructed().UncheckedAt(i) ) AliPHOSgamma;
-// AliPHOSgamma &g = *(AliPHOSgamma*)(GetGammasReconstructed().UncheckedAt(i));
-
- AliPHOSgamma *gggg = new AliPHOSgamma;
- if( NULL==gggg )
- {
- Error("Reconstruction","Can not create AliPHOSgamma");
- exit(1);
- }
-
- GetGammasReconstructed().Add(gggg);
- AliPHOSgamma &g=*gggg;
-
- Float_t thetta, alpha, betta, R=fRadius+rcgamma.Zgdev[i]/10;
-
- g.fX = rcgamma.YW[i]/10;
- g.fY = rcgamma.XW[i]/10;
- g.fE = rcgamma.E [i];
-
- thetta = atan(g.fX/R);
-
- alpha = atan(g.fY/R);
- betta = fPhi/180*TMath::Pi() + alpha;
-
- g.fPx = g.fE * cos(thetta) * cos(betta);
- g.fPy = g.fE * cos(thetta) * sin(betta);
- g.fPz = g.fE * sin(thetta);
- }
-}
+ // FR4: The Plastic which makes up the frame of micromegas -> 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];
}
-
-//______________________________________________________________________________
-#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 <TArray.h>
-#include <TRandom.h>
-#include <TH2.h>
-#include <TVector3.h>
-
-// --- 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
--- /dev/null
+/**************************************************************************
+ * 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
+}
--- /dev/null
+#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
--- /dev/null
+/**************************************************************************
+ * 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 ;
+}
--- /dev/null
+#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
--- /dev/null
+/**************************************************************************
+ * 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 <iostream>
+
+// --- 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 ;
+
+}
--- /dev/null
+#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
--- /dev/null
+/**************************************************************************
+ * 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; iDigit<fMulDigit; iDigit++) {
+ digit = (AliPHOSDigit *) fDigitsList[iDigit];
+ Int_t relid[4] ;
+ Float_t xi ;
+ Float_t zi ;
+ PHOSGeom->AbsToRelNumbering(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; iDigit<fMulDigit; iDigit++) {
+ digit = (AliPHOSDigit *) fDigitsList[iDigit];
+ Int_t relid[4] ;
+ Float_t xi ;
+ Float_t zi ;
+ PHOSGeom->AbsToRelNumbering(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<fMulDigit; iDigit++) {
+
+ if(fEnergyList[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<fMulDigit; iDigit++) {
+
+ if(fEnergyList[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<fMulDigit; iDigit++) {
+ maxAt[iDigit] = fDigitsList[iDigit] ;
+ }
+
+ for(iDigit=0 ; iDigit<fMulDigit; iDigit++) {
+ if(maxAt[iDigit]!= -1) {
+ digit = (AliPHOSDigit *) maxAt[iDigit];
+
+ for(iDigitN=0; iDigitN<fMulDigit; iDigitN++) {
+ digitN = (AliPHOSDigit *) fDigitsList[iDigitN];
+
+ if ( AreNeighbours(digit,digitN) ) {
+ if (fEnergyList[iDigit] > fEnergyList[iDigitN] ) {
+ maxAt[iDigitN] = -1 ;
+ //but may be digit is not local max too ?
+ if(fEnergyList[iDigit]<fEnergyList[iDigitN]+fLocMaxCut)
+ maxAt[iDigit] = -1 ;
+ }
+ else {
+ maxAt[iDigit] = -1 ;
+ //but may be digitN 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; iDigit<fMulDigit; iDigit++) {
+ if(maxAt[iDigit] != -1){
+ maxAt[iDigitN] = maxAt[iDigit] ;
+ maxAtEnergy[iDigitN++] = fEnergyList[iDigit] ;
+ }
+ }
+ return iDigitN ;
+}
+
+//____________________________________________________________________________
+void AliPHOSEmcRecPoint::GetLocalPosition(TVector3 &LPos)
+{
+ if( fLocPos.X() < 1000000.) { //allready evaluated
+ LPos = fLocPos ;
+ return ;
+ }
+
+ Float_t wtot = 0. ;
+
+ Int_t relid[4] ;
+
+ Float_t x = 0. ;
+ Float_t z = 0. ;
+
+ AliPHOSDigit * digit ;
+
+ AliPHOSGeometry * PHOSGeom = (AliPHOSGeometry *) fGeom ;
+
+ Int_t iDigit;
+
+
+ for(iDigit=0; iDigit<fMulDigit; iDigit++) {
+ digit = (AliPHOSDigit *) fDigitsList[iDigit];
+
+ Float_t xi ;
+ Float_t zi ;
+ PHOSGeom->AbsToRelNumbering(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<fMulDigit; iDigit++) {
+// digit = (AliPHOSDigit *) DL[iDigit];
+// AddDigit(*digit) ;
+// }
+
+// fAmp = Clu.GetTotalEnergy() ;
+// fGeom = Clu.GetGeom() ;
+// TVector3 LocPos;
+// Clu.GetLocalPosition(LocPos) ;
+// fLocPos = LocPos;
+// fMulDigit = Clu.GetMultiplicity() ;
+// fMaxDigit = Clu.GetMaximumMultiplicity() ;
+// fPHOSMod = Clu.GetPHOSMod() ;
+// fW0 = Clu.GetLogWeightCut() ;
+// fDelta = Clu.GetDelta() ;
+// fLocMaxCut = Clu.GetLocMaxCut() ;
+
+// delete DL ;
+
+// return *this ;
+// }
+
+//____________________________________________________________________________
+void AliPHOSEmcRecPoint::Print(Option_t * option)
+{
+ cout << "AliPHOSEmcRecPoint: " << endl ;
+
+ AliPHOSDigit * digit ;
+ Int_t iDigit;
+
+ for(iDigit=0; iDigit<fMulDigit; iDigit++) {
+ digit = (AliPHOSDigit *) fDigitsList[iDigit];
+ cout << " digit Id = " << digit->GetId()
+ << " digit Energy = " << fEnergyList[iDigit] << endl ;
+ }
+ cout << " Multiplicity = " << fMulDigit << endl ;
+ cout << " Cluster Energy = " << fAmp << endl ;
+
+}
+
+//______________________________________________________________________________
+void AliPHOSEmcRecPoint::Streamer(TBuffer &R__b)
+{
+ // Stream an object of class AliPHOSEmcRecPoint.
+
+ if (R__b.IsReading()) {
+ Version_t R__v = R__b.ReadVersion(); if (R__v) { }
+ AliPHOSRecPoint::Streamer(R__b);
+ R__b >> fDelta;
+ R__b >> fLocMaxCut;
+ R__b.ReadArray(fEnergyList);
+ R__b >> fW0;
+ } else {
+ R__b.WriteVersion(AliPHOSEmcRecPoint::IsA());
+ AliPHOSRecPoint::Streamer(R__b);
+ R__b << fDelta;
+ R__b << fLocMaxCut;
+ R__b.WriteArray(fEnergyList, GetMaximumDigitMultiplicity() );
+ R__b << fW0;
+ }
+}
+
+
--- /dev/null
+#ifndef ALIPHOSEMCRECPOINT_H
+#define ALIPHOSEMCRECPOINT_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice */
+
+//////////////////////////////////////////////////
+// Rec Point in the EM calorimeter of PHOS //
+// //
+// Author Dmitri Peressounko RRC KI //
+// comment: contains list of AliPHOSDigit's * //
+// and evaluates a few average values //
+//////////////////////////////////////////////////
+
+// --- ROOT system ---
+
+#include "TObject.h"
+
+// --- Standard library ---
+
+// --- AliRoot header files ---
+
+#include "AliPHOSDigit.h"
+#include "AliPHOSRecPoint.h"
+
+class AliPHOSEmcRecPoint : public AliPHOSRecPoint {
+
+public:
+
+ AliPHOSEmcRecPoint() ;
+ AliPHOSEmcRecPoint(Float_t W0, Float_t LocMaxCut) ;
+ // virtual ~AliPHOSEmcRecPoint() ;
+ void AddDigit(AliDigitNew & digit, Float_t Energy) ; // add a digit to the digits list
+ Int_t Compare(TObject * obj) ; // method for sorting
+
+ Float_t GetDelta (){ return fDelta ; }
+ Float_t GetDispersion() ; // computes the dispersion of the shower
+ void GetElipsAxis(Float_t * lambda) ; // computes the axis of shower ellipsoide
+ Float_t * GetEnergiesList(){return fEnergyList ;}
+ Float_t GetLocMaxCut () {return fLocMaxCut ; }
+ Float_t GetLogWeightCut (){return fW0 ; }
+ Float_t GetMaximalEnergy(void) ; // get the highest energy in the cluster
+ Int_t GetMaximumMultiplicity() { return fMaxDigit ; }
+ Int_t GetMultiplicity(void) const { return fMulDigit ; }
+ Int_t GetMultiplicityAtLevel(const Float_t level) ; // computes multiplicity of digits with energy above relative level
+ Int_t GetNumberOfLocalMax(int * maxAt, Float_t * maxAtEnergy) ; // searches for the local maxima
+
+ Float_t GetTotalEnergy(void) const { return fAmp ; } // in EMC RecPoint Amp = Energy
+ void GetLocalPosition(TVector3 &Lpos) ; // computes the position in the PHOS module
+ Bool_t IsEmc(void) {return kTRUE ; }
+ Bool_t IsSortable() const { return kTRUE ; }
+ void Print(Option_t * opt = "void") ;
+
+private:
+
+ Bool_t AreNeighbours(AliPHOSDigit * digit1, AliPHOSDigit * digit2 ) ;
+
+public:
+
+ // AliPHOSEmcRecPoint& operator = (AliPHOSEmcRecPoint clu) ;
+
+private:
+ Float_t fDelta ; // parameter used to sort the clusters
+ Float_t fLocMaxCut ; // parameter used for local maximum searc
+ Float_t * fEnergyList ; //energy of digits
+ Float_t fW0 ; // logarithmic weight factor for center of gravity calculation
+
+public:
+
+ClassDef(AliPHOSEmcRecPoint,1) // EMC cluster, version 1
+
+};
+
+#endif // AliPHOSEMCRECPOINT_H
--- /dev/null
+/**************************************************************************
+ * 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. *
+ **************************************************************************/
+
+//_________________________________________________________________________
+// Geometry class for PHOS version SUBATECH
+//*-- Author : Y. Schutz SUBATECH
+//////////////////////////////////////////////////////////////////////////////
+
+// --- ROOT system ---
+
+#include "TVector3.h"
+#include "TRotation.h"
+
+// --- Standard library ---
+
+#include <iostream.h>
+#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 <E> : 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 ;
+ }
+}
+
--- /dev/null
+#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
--- /dev/null
+/**************************************************************************
+ * 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 ;
+}
--- /dev/null
+#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
--- /dev/null
+/**************************************************************************
+ * 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 <iostream.h>
+
+// --- 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; iDigit<fMulDigit; iDigit++) {
+ digit = (AliPHOSDigit *) fDigitsList[iDigit];
+
+ Float_t xi ;
+ Float_t zi ;
+ PHOSGeom->AbsToRelNumbering(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<fMulDigit; iDigit++) {
+ digit = (AliPHOSDigit *) fDigitsList[iDigit];
+ cout << " digit Id = " << digit->GetId()
+ << " digit Energy = " << digit->GetAmp() << endl ;
+ }
+ cout << " Multiplicity = " << fMulDigit << endl ;
+}
+
+
+// //____________________________________________________________________________
+// AliPHOSPpsdRecPoint& AliPHOSPpsdRecPoint::operator = (AliPHOSPpsdRecPoint Clu)
+// {
+// int* DL = Clu.GetDigitsList() ;
+
+// if(fDigitsList)
+// delete fDigitsList ;
+
+// AliPHOSDigit * digit ;
+
+// Int_t iDigit;
+
+// for(iDigit=0; iDigit<fMulDigit; iDigit++) {
+// digit = (AliPHOSDigit *) DL[iDigit];
+// AddDigit(*digit) ;
+// }
+
+// fDelta = Clu.GetDelta() ;
+// delete DL ;
+
+// fAmp = Clu.GetEnergy() ;
+// fGeom = Clu.GetGeom() ;
+// TVector3 LocPos;
+// Clu.GetLocalPosition(LocPos) ;
+// fLocPos = LocPos;
+// fMulDigit = Clu.GetMultiplicity() ;
+// fMaxDigit = Clu.GetMaximumMultiplicity() ;
+// fPHOSMod = Clu.GetPHOSMod() ;
+
+
+// return *this ;
+// }
--- /dev/null
+#ifndef ALIPHOSPPSDRECPOINT_H
+#define ALIPHOSPPSDRECPOINT_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice */
+
+////////////////////////////////////////////////
+// Cluster in the PPSD of PHOS //
+// Version SUBATECH //
+// Author Yves Schutz SUBATECH //
+// comment: its a list of AliPHOSDigit's //
+// //
+////////////////////////////////////////////////
+
+// --- ROOT system ---
+
+// --- Standard library ---
+
+// --- AliRoot header files ---
+
+#include "AliPHOSDigit.h"
+#include "AliPHOSRecPoint.h"
+
+class AliPHOSPpsdRecPoint : public AliPHOSRecPoint {
+
+public:
+
+ AliPHOSPpsdRecPoint() ; // ctor
+ virtual ~AliPHOSPpsdRecPoint() ; // dtor
+ void AddDigit(AliDigitNew & digit, Float_t Energy) ;
+ Int_t Compare(TObject * obj) ; // method to sort clusters
+
+ Float_t GetDelta(void) {return fDelta ;}
+ Int_t GetMultiplicity(void) const { return fMulDigit ; }
+ Int_t GetMaximumMultiplicity() { return fMaxDigit ; }
+ void GetLocalPosition(TVector3 &LPos) ; // computes the position in the module of the cluster center
+ //projection of ALICE axes on PHOS Module, y = 0 .
+ Bool_t GetUp() ; // true if cluster is in upper ppsd
+ Bool_t IsEmc(void) {return kFALSE ; }
+ Bool_t IsSortable() const { return kTRUE ; }
+ void Print(Option_t * opt = "void") ;
+
+ // AliPHOSPpsdRecPoint& operator = (AliPHOSPpsdRecPoint Clu) ;
+
+private:
+
+ Float_t fDelta ; // parameter used for sorting
+
+
+public:
+
+ ClassDef(AliPHOSPpsdRecPoint,1) // PPSD RecPoint, version 1
+
+};
+
+#endif // AliPHOSPPSDRECPOINT_H
--- /dev/null
+/**************************************************************************
+ * 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. *
+ **************************************************************************/
+
+//_________________________________________________________________________
+// PHOSRecPoint base class deriving from AliRecPoint
+//*-- Author : Gines MARTINEZ SUBATECH
+//////////////////////////////////////////////////////////////////////////////
+
+// --- ROOT system ---
+
+// --- Standard library ---
+
+// --- AliRoot header files ---
+
+#include "AliPHOSGeometry.h"
+#include "AliPHOSDigit.h"
+#include "AliPHOSRecPoint.h"
+
+ClassImp(AliPHOSRecPoint)
+
+
+//____________________________________________________________________________
+AliPHOSRecPoint::AliPHOSRecPoint()
+ : AliRecPoint()
+{
+ fGeom = AliPHOSGeometry::GetInstance() ;
+ fPHOSMod = 0;
+}
+
+//____________________________________________________________________________
+AliPHOSRecPoint::~AliPHOSRecPoint()
+{
+ // dtor
+}
+
+//____________________________________________________________________________
+Int_t AliPHOSRecPoint::GetPHOSMod()
+{
+ if(fPHOSMod > 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 ;
+}
--- /dev/null
+#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
--- /dev/null
+/**************************************************************************
+ * 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) ;
+}
--- /dev/null
+#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
--- /dev/null
+/**************************************************************************
+ * 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-------- "<<endl ;
+ cout << "EMC Reconstructed Point: " << fEmcRecPoint << endl;
+
+ TVector3 pos ;
+ TMatrix Dummy ;
+
+ fEmcRecPoint->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 ;
+ }
+
+}
+
--- /dev/null
+#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
--- /dev/null
+/**************************************************************************
+ * 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)
+{
+
+}
+
--- /dev/null
+#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
--- /dev/null
+/**************************************************************************
+ * 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++ ;
+ }
+}
--- /dev/null
+#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
* 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 <stdio.h>
+#include <string.h>
+#include <stdlib.h>
+#include <strstream.h>
+#include <assert.h>
+
+// --- 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.25<y<0.25 TILTED GEOMETRY OF THE PHOS ***
-// ORIGIN : NICK VAN EIJNDHOVEN
-
- Float_t pphi;
- Float_t r, dptcb[3], dpair[3], dphos[3], dpucp[3], dpasp[3];
- Float_t dpxtl[3];
- Float_t yo;
- Int_t idrotm[99];
- Float_t xp1, yp1, xp2, yp2;
+
+ AliPHOSv0 *PHOS_tmp = (AliPHOSv0*)gAlice->GetModule("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 << "<I> 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
+}
+
-#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
# 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
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)
# 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 #####
depend: $(SRCS)
-TOCLEAN = $(OBJS) *Cint.cxx *Cint.h
+TOCLEAN = $(OBJS) *Cint.cxx *Cint.h
############################### General Macros ################################
#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
--- /dev/null
+/**************************************************************************
+ * 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
+}
+
--- /dev/null
+#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
--- /dev/null
+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 "<<Etot << endl ;
+ TClonesArray * PpsdRP = PHOS->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 ;
+
+}
+
+