/* $Id$ */
+/* History of cvs commits:
+ *
+ * $Log$
+ * Revision 1.111 2007/07/24 09:41:19 morsch
+ * AliStack included for kKeepBit.
+ *
+ * Revision 1.110 2007/03/10 08:58:52 kharlov
+ * Protection for noCPV geometry
+ *
+ * Revision 1.109 2007/03/01 11:37:37 kharlov
+ * Strip units changed from 8x1 to 8x2 (T.Pocheptsov)
+ *
+ * Revision 1.108 2007/02/02 09:40:50 alibrary
+ * Includes required by ROOT head
+ *
+ * Revision 1.107 2007/02/01 10:34:47 hristov
+ * Removing warnings on Solaris x86
+ *
+ * Revision 1.106 2006/11/14 17:11:15 hristov
+ * Removing inheritances from TAttLine, TAttMarker and AliRndm in AliModule. The copy constructor and assignment operators are moved to the private part of the class and not implemented. The corresponding changes are propagated to the detectors
+ *
+ * Revision 1.105 2006/09/13 07:31:01 kharlov
+ * Effective C++ corrections (T.Pocheptsov)
+ *
+ * Revision 1.104 2005/05/28 14:19:05 schutz
+ * Compilation warnings fixed by T.P.
+ *
+ */
+
//_________________________________________________________________________
// Implementation version v1 of PHOS Manager class
//---
-// Layout EMC + PPSD has name GPS2:
-// Produces cumulated hits
//---
// Layout EMC + CPV has name IHEP:
// Produces hits for CPV, cumulated hits
//---
-// Layout EMC + CPV + PPSD has name GPS:
-// Produces hits for CPV, cumulated hits
//---
//*-- Author: Yves Schutz (SUBATECH)
// --- ROOT system ---
-
-#include "TBRIK.h"
-#include "TNode.h"
-#include "TRandom.h"
-#include "TTree.h"
-
+#include <TClonesArray.h>
+#include <TParticle.h>
+#include <TVirtualMC.h>
// --- Standard library ---
-#include <stdio.h>
-#include <string.h>
-#include <stdlib.h>
-#include <strstream.h>
// --- AliRoot header files ---
-
-#include "AliPHOSv1.h"
+#include "AliPHOSCPVDigit.h"
+#include "AliPHOSGeometry.h"
#include "AliPHOSHit.h"
-#include "AliPHOSDigit.h"
-#include "AliPHOSReconstructioner.h"
+#include "AliPHOSv1.h"
#include "AliRun.h"
-#include "AliConst.h"
#include "AliMC.h"
+#include "AliStack.h"
ClassImp(AliPHOSv1)
//____________________________________________________________________________
-AliPHOSv1::AliPHOSv1()
+AliPHOSv1::AliPHOSv1():
+ fLightYieldMean(0.),
+ fIntrinsicPINEfficiency(0.),
+ fLightYieldAttenuation(0.),
+ fRecalibrationFactor(0.),
+ fElectronsPerGeV(0.),
+ fAPDGain(0.),
+ fLightFactor(0.),
+ fAPDFactor(0.)
{
- // ctor
-
- // Create an empty array of AliPHOSCPVModule to satisfy
- // AliPHOSv1::Streamer when reading root file
-
- fReconstructioner = 0;
- fTrackSegmentMaker = 0;
-
- if ( 0==(fEMCModules=new TClonesArray("AliPHOSCPVModule",0)) ) {
- Error("AliPHOSv1","Can not create array of EMC modules");
- exit(1);
- }
-
- if ( 0==(fCPVModules=new TClonesArray("AliPHOSCPVModule",0)) ) {
- Error("AliPHOSv1","Can not create array of CPV modules");
- exit(1);
- }
-
+ //Def ctor.
}
//____________________________________________________________________________
AliPHOSv1::AliPHOSv1(const char *name, const char *title):
-AliPHOSv0(name,title)
+ AliPHOSv0(name,title),
+ fLightYieldMean(0.),
+ fIntrinsicPINEfficiency(0.),
+ fLightYieldAttenuation(0.),
+ fRecalibrationFactor(0.),
+ fElectronsPerGeV(0.),
+ fAPDGain(0.),
+ fLightFactor(0.),
+ fAPDFactor(0.)
{
- // ctor : title is used to identify the layout
- // GPS2 = 5 modules (EMC + PPSD)
- // IHEP = 5 modules (EMC + CPV )
- // MIXT = 4 modules (EMC + CPV ) and 1 module (EMC + PPSD)
//
// We store hits :
// - fHits (the "normal" one), which retains the hits associated with
// (this array is reset after each primary has been tracked).
//
- fPinElectronicNoise = 0.010 ;
- fDigitThreshold = 0.1 ; // 1 GeV
+
// We do not want to save in TreeH the raw hits
// But save the cumulated hits instead (need to create the branch myself)
// It is put in the Digit Tree because the TreeH is filled after each primary
- // and the TreeD at the end of the event (branch is set in FinishEvent() ).
+ // and the TreeD at the end of the event (branch is set in FinishEvent() ).
fHits= new TClonesArray("AliPHOSHit",1000) ;
- gAlice->AddHitList(fHits);
-
- fNhits = 0 ;
-
- fReconstructioner = 0;
- fTrackSegmentMaker = 0;
-
- fIshunt = 1 ; // All hits are associated with primary particles
-
- // Create array of EMC modules of the size of PHOS modules number
-
- if ( 0==(fEMCModules=new TClonesArray("AliPHOSCPVModule",fGeom->GetNModules())) ) {
- Error("AliPHOSv1","Can not create array of EMC modules");
- exit(1);
- }
- TClonesArray &lemcmodule = *fEMCModules;
- for (Int_t i=0; i<fGeom->GetNModules(); i++) new(lemcmodule[i]) AliPHOSCPVModule();
-
- // Create array of CPV modules for the IHEP's version of CPV
-
- if ( strcmp(fGeom->GetName(),"IHEP") == 0 || strcmp(fGeom->GetName(),"MIXT") == 0 ) {
- // Create array of CPV modules of the size of PHOS modules number
-
- if ( 0==(fCPVModules=new TClonesArray("AliPHOSCPVModule",fGeom->GetNCPVModules())) ) {
- Error("AliPHOSv1","Can not create array of CPV modules");
- exit(1);
- }
- TClonesArray &lcpvmodule = *fCPVModules;
- for (Int_t i=0; i<fGeom->GetNCPVModules(); i++) new(lcpvmodule[i]) AliPHOSCPVModule();
- }
- else {
- // Create an empty array of AliPHOSCPVModule to satisfy
- // AliPHOSv1::Streamer when writing root file
+ gAlice->GetMCApp()->AddHitList(fHits) ;
- fCPVModules=new TClonesArray("AliPHOSCPVModule",0);
-
- }
-}
-
-//____________________________________________________________________________
-AliPHOSv1::AliPHOSv1(AliPHOSReconstructioner * Reconstructioner, const char *name, const char *title):
- AliPHOSv0(name,title)
-{
- // ctor : title is used to identify the layout
- // GPS2 = 5 modules (EMC + PPSD)
- // 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.
-
- fPinElectronicNoise = 0.010 ;
-
- // We do not want to save in TreeH the raw hits
- //fHits = new TClonesArray("AliPHOSHit",100) ;
-
- fDigits = 0 ;
- fHits= new TClonesArray("AliPHOSHit",1000) ;
- gAlice->AddHitList(fHits);
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 << "AliPHOS" << Version() << " : PHOS geometry intialized for " << fGeom->GetName() << endl ;
- else
- cout << "AliPHOS" << Version() << " : PHOS geometry initialization failed !" << endl ;
-
- // Defining the PHOS Reconstructioner
-
- fReconstructioner = Reconstructioner ;
-
+ fIshunt = 2 ; // All hits are associated with primary particles
+
+ //Photoelectron statistics:
+ // The light yield is a poissonian distribution of the number of
+ // photons created in the PbWo4 crystal, calculated using following formula
+ // NumberOfPhotons = EnergyLost * LightYieldMean* APDEfficiency *
+ // exp (-LightYieldAttenuation * DistanceToPINdiodeFromTheHit);
+ // LightYieldMean is parameter calculated to be over 47000 photons per GeV
+ // APDEfficiency is 0.02655
+ // k_0 is 0.0045 from Valery Antonenko
+ // The number of electrons created in the APD is
+ // NumberOfElectrons = APDGain * LightYield
+ // The APD Gain is 300
+ fLightYieldMean = 47000;
+ fIntrinsicPINEfficiency = 0.02655 ; //APD= 0.1875/0.1271 * 0.018 (PIN)
+ fLightYieldAttenuation = 0.0045 ;
+ fRecalibrationFactor = 13.418/ fLightYieldMean ;
+ fElectronsPerGeV = 2.77e+8 ;
+ fAPDGain = 300. ;
+ fLightFactor = fLightYieldMean * fIntrinsicPINEfficiency ;
+ fAPDFactor = (fRecalibrationFactor/100.) * fAPDGain ;
}
//____________________________________________________________________________
AliPHOSv1::~AliPHOSv1()
{
// dtor
-
- if ( fHits) {
+ if ( fHits) {
fHits->Delete() ;
delete fHits ;
fHits = 0 ;
- }
-
- if ( fDigits) {
- fDigits->Delete() ;
- delete fDigits ;
- fDigits = 0 ;
- }
-
- if ( fEmcRecPoints ) {
- fEmcRecPoints->Delete() ;
- delete fEmcRecPoints ;
- fEmcRecPoints = 0 ;
- }
-
- if ( fPpsdRecPoints ) {
- fPpsdRecPoints->Delete() ;
- delete fPpsdRecPoints ;
- fPpsdRecPoints = 0 ;
- }
-
- if ( fTrackSegments ) {
- fTrackSegments->Delete() ;
- delete fTrackSegments ;
- fTrackSegments = 0 ;
- }
-
+ }
}
//____________________________________________________________________________
-void AliPHOSv1::AddHit(Int_t shunt, Int_t primary, Int_t tracknumber, Int_t Id, Float_t * hits, Int_t trackpid)
+void AliPHOSv1::AddHit(Int_t shunt, Int_t primary, Int_t Id, Float_t * hits)
{
// Add a hit to the hit list.
- // A PHOS hit is the sum of all hits in a single crystal
- // or in a single PPSD gas cell
+ // A PHOS hit is the sum of all hits in a single crystal from one primary and within some time gate
Int_t hitCounter ;
AliPHOSHit *newHit ;
AliPHOSHit *curHit ;
Bool_t deja = kFALSE ;
+ AliPHOSGeometry * geom = GetGeometry() ;
- newHit = new AliPHOSHit(shunt, primary, tracknumber, Id, hits, trackpid) ;
+ newHit = new AliPHOSHit(shunt, primary, Id, hits) ;
- for ( hitCounter = 0 ; hitCounter < fNhits && !deja ; hitCounter++ ) {
- curHit = (AliPHOSHit*) (*fHits)[hitCounter] ;
+ for ( hitCounter = fNhits-1 ; hitCounter >= 0 && !deja ; hitCounter-- ) {
+ curHit = dynamic_cast<AliPHOSHit*>((*fHits)[hitCounter]) ;
+ if(curHit->GetPrimary() != primary) break ;
+ // We add hits with the same primary, while GEANT treats primaries succesively
if( *curHit == *newHit ) {
- *curHit = *curHit + *newHit ;
+ *curHit + *newHit ;
deja = kTRUE ;
}
}
if ( !deja ) {
new((*fHits)[fNhits]) AliPHOSHit(*newHit) ;
+ // get the block Id number
+ Int_t relid[4] ;
+ geom->AbsToRelNumbering(Id, relid) ;
+
fNhits++ ;
}
-
- delete newHit;
-}
-
-//___________________________________________________________________________
-Int_t AliPHOSv1::Digitize(Float_t Energy)
-{
- // Applies the energy calibration
- Float_t fB = 100000000. ;
- Float_t fA = 0. ;
- Int_t chan = Int_t(fA + Energy*fB ) ;
- return chan ;
+ delete newHit;
}
//____________________________________________________________________________
-void AliPHOSv1::Hit2Digit(Int_t ntracks){
- //Collects all hits in the same active volume into digits
-
- if(fDigits!= 0)
- fDigits->Clear() ;
- else
- fDigits = new TClonesArray("AliPHOSDigit",1000) ;
-
- // Branch address for digit tree
- char branchname[20];
- sprintf(branchname,"%s",GetName());
- gAlice->TreeD()->Branch(branchname,&fDigits,fBufferSize);
-
- gAlice->TreeD()->GetEvent(0);
-
-
- Int_t i ;
- Int_t relid[4];
- Int_t j ;
- AliPHOSHit * hit ;
- AliPHOSDigit * newdigit ;
- AliPHOSDigit * curdigit ;
- Bool_t deja = kFALSE ;
-
- Int_t itrack ;
- for (itrack=0; itrack<ntracks; itrack++){
-
- //=========== Get the Hits Tree for the Primary track itrack
- gAlice->ResetHits();
- gAlice->TreeH()->GetEvent(itrack);
-
- for ( i = 0 ; i < fHits->GetEntries() ; i++ ) {
- hit = (AliPHOSHit*)fHits->At(i) ;
-
- // Assign primary number only if contribution is significant
- if( hit->GetEnergy() > fDigitThreshold)
- newdigit = new AliPHOSDigit( hit->GetPrimary(), hit->GetId(), Digitize( hit->GetEnergy() ) ) ;
- else
- newdigit = new AliPHOSDigit( -1 , hit->GetId(), Digitize( hit->GetEnergy() ) ) ;
- deja =kFALSE ;
-
-
- for ( j = 0 ; j < fNdigits ; j++) {
- curdigit = (AliPHOSDigit*) fDigits->At(j) ;
- if ( *curdigit == *newdigit) {
- *curdigit = *curdigit + *newdigit ;
- deja = kTRUE ;
- }
- }
-
- if ( !deja ) {
- new((*fDigits)[fNdigits]) AliPHOSDigit(* newdigit) ;
- fNdigits++ ;
- }
-
- delete newdigit ;
- }
-
- } // loop over tracks
-
- // Noise induced by the PIN diode of the PbWO crystals
-
- Float_t energyandnoise ;
- for ( i = 0 ; i < fNdigits ; i++ ) {
- newdigit = (AliPHOSDigit * ) fDigits->At(i) ;
-
- fGeom->AbsToRelNumbering(newdigit->GetId(), relid) ;
-
- if (relid[1]==0){ // Digits belong to EMC (PbW0_4 crystals)
- energyandnoise = newdigit->GetAmp() + Digitize(gRandom->Gaus(0., fPinElectronicNoise)) ;
-
- if (energyandnoise < 0 )
- energyandnoise = 0 ;
-
- if ( newdigit->GetAmp() < fDigitThreshold ) // if threshold not surpassed, remove digit from list
- fDigits->RemoveAt(i) ;
- }
- }
-
- fDigits->Compress() ;
-
- fNdigits = fDigits->GetEntries() ;
- fDigits->Expand(fNdigits) ;
-
- for (i = 0 ; i < fNdigits ; i++) {
- newdigit = (AliPHOSDigit *) fDigits->At(i) ;
- newdigit->SetIndexInList(i) ;
- }
-
- gAlice->TreeD()->Fill() ;
-
- gAlice->TreeD()->Write(0,TObject::kOverwrite) ;
-
-}
-//___________________________________________________________________________
-void AliPHOSv1::MakeBranch(Option_t* opt)
-{
- // Create new branche in the current Root Tree in the digit Tree
- AliDetector::MakeBranch(opt) ;
-
- // Create new branches EMC<i> for hits in EMC modules
-
- for( Int_t i=0; i<fGeom->GetNModules(); i++ ) GetEMCModule(i).MakeBranch("EMC",i+1);
-
- // Create new branches CPV<i> for hits in CPV modules for IHEP geometry
-
- if ( strcmp(fGeom->GetName(),"IHEP") == 0 || strcmp(fGeom->GetName(),"MIXT") == 0 ) {
- for( Int_t i=0; i<fGeom->GetNCPVModules(); i++ ) GetCPVModule(i).MakeBranch("CPV",i+1);
- }
-
-}
-
-//_____________________________________________________________________________
-void AliPHOSv1::Reconstruction(AliPHOSReconstructioner * Reconstructioner)
-{
- // 1. Reinitializes the existing RecPoint, TrackSegment, and RecParticles Lists and
- // 2. Creates TreeR with a branch for each list
- // 3. Steers the reconstruction processes
- // 4. Saves the 3 lists in TreeR
- // 5. Write the Tree to File
-
- fReconstructioner = Reconstructioner ;
-
- char branchname[10] ;
-
- // 1.
-
- // gAlice->MakeTree("R") ;
- Int_t splitlevel = 0 ;
-
- fEmcRecPoints->Delete() ;
-
- if ( fEmcRecPoints && gAlice->TreeR() ) {
- sprintf(branchname,"%sEmcRP",GetName()) ;
- gAlice->TreeR()->Branch(branchname, "TObjArray", &fEmcRecPoints, fBufferSize, splitlevel) ;
- }
-
- fPpsdRecPoints->Delete() ;
-
- if ( fPpsdRecPoints && gAlice->TreeR() ) {
- sprintf(branchname,"%sPpsdRP",GetName()) ;
- gAlice->TreeR()->Branch(branchname, "TObjArray", &fPpsdRecPoints, fBufferSize, splitlevel) ;
- }
-
- fTrackSegments->Delete() ;
-
- if ( fTrackSegments && gAlice->TreeR() ) {
- sprintf(branchname,"%sTS",GetName()) ;
- gAlice->TreeR()->Branch(branchname, &fTrackSegments, fBufferSize) ;
- }
-
- fRecParticles->Delete() ;
-
- if (strcmp(fGeom->GetName(),"GPS2") == 0 || strcmp(fGeom->GetName(),"MIXT") == 0) {
- if ( fRecParticles && gAlice->TreeR() ) {
- sprintf(branchname,"%sRP",GetName()) ;
- gAlice->TreeR()->Branch(branchname, &fRecParticles, fBufferSize) ;
- }
- }
-
- // 3.
- if (strcmp(fGeom->GetName(),"GPS2") == 0 || strcmp(fGeom->GetName(),"MIXT") == 0)
- fReconstructioner->MakePPSD(fDigits, fEmcRecPoints, fPpsdRecPoints, fTrackSegments, fRecParticles);
- if (strcmp(fGeom->GetName(),"IHEP") == 0 || strcmp(fGeom->GetName(),"MIXT") == 0)
- fReconstructioner->MakeCPV (fDigits, fEmcRecPoints, fPpsdRecPoints);
-
- // 4. Expand or Shrink the arrays to the proper size
-
- Int_t size ;
-
- size = fEmcRecPoints->GetEntries() ;
- fEmcRecPoints->Expand(size) ;
-
- size = fPpsdRecPoints->GetEntries() ;
- fPpsdRecPoints->Expand(size) ;
-
- size = fTrackSegments->GetEntries() ;
- fTrackSegments->Expand(size) ;
-
- size = fRecParticles->GetEntries() ;
- fRecParticles->Expand(size) ;
-
- gAlice->TreeR()->Fill() ;
- // 5.
+void AliPHOSv1::FinishPrimary()
+{
+ // called at the end of each track (primary) by AliRun
+ // hits are reset for each new track
+ // accumulate the total hit-multiplicity
- gAlice->TreeR()->Write(0,TObject::kOverwrite) ;
-
- // Deleting reconstructed objects
- ResetReconstruction();
-
}
//____________________________________________________________________________
-void AliPHOSv1::ResetHits()
-{
- // Reset hit tree for EMC and CPV
- // Yuri Kharlov, 28 September 2000
-
- AliDetector::ResetHits();
- for (Int_t i=0; i<fGeom->GetNModules(); i++) ((AliPHOSCPVModule*)(*fEMCModules)[i]) -> Clear();
- if ( strcmp(fGeom->GetName(),"IHEP") == 0 || strcmp(fGeom->GetName(),"MIXT") == 0 ) {
- for (Int_t i=0; i<fGeom->GetNCPVModules(); i++) ((AliPHOSCPVModule*)(*fCPVModules)[i]) -> Clear();
- }
-
-}
-//____________________________________________________________________________
-void AliPHOSv1::ResetReconstruction()
-{
- // Deleting reconstructed objects
-
- if ( fEmcRecPoints ) fEmcRecPoints->Delete();
- if ( fPpsdRecPoints ) fPpsdRecPoints->Delete();
- if ( fTrackSegments ) fTrackSegments->Delete();
- if ( fRecParticles ) fRecParticles->Delete();
+void AliPHOSv1::FinishEvent()
+{
+ // called at the end of each event by AliRun
+ // accumulate the hit-multiplicity and total energy per block
+ // if the values have been updated check it
+ AliDetector::FinishEvent();
}
-
-//____________________________________________________________________________
-void AliPHOSv1::SetTreeAddress()
-{
- // TBranch *branch;
- AliPHOS::SetTreeAddress();
-
-// //Branch address for TreeR: RecPpsdRecPoint
-// TTree *treeR = gAlice->TreeR();
-// if ( treeR && fPpsdRecPoints ) {
-// branch = treeR->GetBranch("PHOSPpsdRP");
-// if (branch) branch->SetAddress(&fPpsdRecPoints) ;
-// }
-
- // Set branch address for the Hits Tree for hits in EMC modules
- // Yuri Kharlov, 23 November 2000.
-
- for( Int_t i=0; i<fGeom->GetNModules(); i++ ) GetEMCModule(i).SetTreeAddress("EMC",i+1);
-
- // Set branch address for the Hits Tree for hits in CPV modules for IHEP geometry
- // Yuri Kharlov, 28 September 2000.
-
- if ( strcmp(fGeom->GetName(),"IHEP") == 0 || strcmp(fGeom->GetName(),"MIXT") == 0 ) {
- for( Int_t i=0; i<fGeom->GetNCPVModules(); i++ ) GetCPVModule(i).SetTreeAddress("CPV",i+1);
- }
-
-}
-
//____________________________________________________________________________
-
void AliPHOSv1::StepManager(void)
{
- // Accumulates hits as long as the track stays in a single crystal or PPSD gas Cell
-
-// if (gMC->IsTrackEntering())
-// cout << "Track enters the volume " << gMC->CurrentVolName() << endl;
-// if (gMC->IsTrackExiting())
-// cout << "Track leaves the volume " << gMC->CurrentVolName() << endl;
+ // Accumulates hits as long as the track stays in a single crystal or CPV gas Cell
- Int_t relid[4] ; // (box, layer, row, column) indices
- Int_t absid ; // absolute cell ID number
- Float_t xyze[4]={0,0,0,0} ; // position wrt MRS and energy deposited
- TLorentzVector pos ; // Lorentz vector of the track current position
+ Int_t relid[4] ; // (box, layer, row, column) indices
+ Int_t absid ; // absolute cell ID number
+ Float_t xyzte[5]={-1000.,-1000.,-1000.,0.,0.} ; // position wrt MRS, time and energy deposited
+ TLorentzVector pos ; // Lorentz vector of the track current position
Int_t copy ;
- Int_t i ;
-
- Int_t tracknumber = gAlice->CurrentTrack() ;
- Int_t primary = gAlice->GetPrimary( gAlice->CurrentTrack() );
- TString name = fGeom->GetName() ;
- Int_t trackpid = gMC->TrackPid() ;
- if ( name == "GPS2" || name == "MIXT" ) { // ======> CPV is a GPS' PPSD
- if( gMC->CurrentVolID(copy) == gMC->VolId("GCEL") ) // 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
- if ( name == "MIXT" && strcmp(gMC->CurrentVolOffName(5),"PHO1") == 0 )
- relid[0] += fGeom->GetNModules() - fGeom->GetNPPSDModules();
- gMC->CurrentVolOffID(3, relid[1]) ; // get the Micromegas Module number
- // 1-> fGeom->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
-
- fGeom->RelToAbsNumbering(relid, absid) ;
-
- // add current hit to the hit list
- AddHit(fIshunt, primary, tracknumber, absid, xyze, trackpid);
+ Int_t moduleNumber ;
+
+ static Int_t idPCPQ = -1;
+ if (strstr(fTitle.Data(),"noCPV") == 0)
+ idPCPQ = gMC->VolId("PCPQ");
- } // there is deposited energy
- } // We are inside the gas of the CPV
- } // GPS2 configuration
+ if( gMC->CurrentVolID(copy) == idPCPQ &&
+ (gMC->IsTrackEntering() ) &&
+ gMC->TrackCharge() != 0) {
+
+ gMC -> TrackPosition(pos);
- if ( name == "IHEP" || name == "MIXT" ) { // ======> CPV is a IHEP's one
+ Float_t xyzm[3], xyzd[3] ;
+ Int_t i;
+ for (i=0; i<3; i++) xyzm[i] = pos[i];
+ gMC -> Gmtod (xyzm, xyzd, 1); // transform coordinate from master to daughter system
- // Yuri Kharlov, 28 September 2000
- if( gMC->CurrentVolID(copy) == gMC->VolId("CPVQ") &&
- gMC->IsTrackEntering() &&
- gMC->TrackCharge() != 0) {
+ Float_t xyd[3]={0,0,0} ; //local position of the entering
+ xyd[0] = xyzd[0];
+ xyd[1] =-xyzd[2];
+ xyd[2] =-xyzd[1];
+
+ // Current momentum of the hit's track in the local ref. system
+ TLorentzVector pmom ; //momentum of the particle initiated hit
+ gMC -> TrackMomentum(pmom);
+ Float_t pm[3], pd[3];
+ for (i=0; i<3; i++)
+ pm[i] = pmom[i];
+
+ gMC -> Gmtod (pm, pd, 2); // transform 3-momentum from master to daughter system
+ pmom[0] = pd[0];
+ pmom[1] =-pd[1];
+ pmom[2] =-pd[2];
- // Charged track has just entered to the CPV sensitive plane
-
- AliPHOSv1 &phos = *(AliPHOSv1*)gAlice->GetModule("PHOS");
-
- Int_t moduleNumber;
- gMC->CurrentVolOffID(3,moduleNumber);
- moduleNumber--;
-
- // Current position of the hit in the CPV module ref. system
-
- gMC -> TrackPosition(pos);
- Float_t xyzm[3], xyzd[3], xyd[3]={0,0,0};
- Int_t i;
- for (i=0; i<3; i++) xyzm[i] = pos[i];
- gMC -> Gmtod (xyzm, xyzd, 1); // transform coordinate from master to daughter system
- xyd[0] = xyzd[0];
- xyd[1] =-xyzd[2];
-
- // Current momentum of the hit's track in the CPV module ref. system
+ // Digitize the current CPV hit:
+
+ // 1. find pad response and
+ gMC->CurrentVolOffID(3,moduleNumber);
+ moduleNumber--;
+
+ TClonesArray *cpvDigits = new TClonesArray("AliPHOSCPVDigit",0); // array of digits for current hit
+ CPVDigitize(pmom,xyd,cpvDigits);
- TLorentzVector pmom;
- gMC -> TrackMomentum(pmom);
- Float_t pm[3], pd[3];
- for (i=0; i<3; i++) pm[i] = pmom[i];
- gMC -> Gmtod (pm, pd, 2); // transform 3-momentum from master to daughter system
- pmom[0] = pd[0];
- pmom[1] =-pd[1];
- pmom[2] =-pd[2];
-
- // Current particle type of the hit's track
-
- Int_t ipart = gMC->TrackPid();
-
- // Add the current particle in the list of the CPV hits.
-
- phos.GetCPVModule(moduleNumber).AddHit(fIshunt,primary,pmom,xyd,ipart);
-
- if (fDebug == 1) {
- printf("CPV hit added to module #%2d: p = (% .4f, % .4f, % .4f, % .4f) GeV,\n",
- moduleNumber+1,pmom.Px(),pmom.Py(),pmom.Pz(),pmom.E());
- printf( " xy = (%8.4f, %8.4f) cm, ipart = %d, primary = %d\n",
- xyd[0],xyd[1],ipart,primary);
+ Float_t xmean = 0;
+ Float_t zmean = 0;
+ Float_t qsum = 0;
+ Int_t idigit,ndigits;
+
+ // 2. go through the current digit list and sum digits in pads
+
+ ndigits = cpvDigits->GetEntriesFast();
+ for (idigit=0; idigit<ndigits-1; idigit++) {
+ AliPHOSCPVDigit *cpvDigit1 = dynamic_cast<AliPHOSCPVDigit*>(cpvDigits->UncheckedAt(idigit));
+ Float_t x1 = cpvDigit1->GetXpad() ;
+ Float_t z1 = cpvDigit1->GetYpad() ;
+ for (Int_t jdigit=idigit+1; jdigit<ndigits; jdigit++) {
+ AliPHOSCPVDigit *cpvDigit2 = dynamic_cast<AliPHOSCPVDigit*>(cpvDigits->UncheckedAt(jdigit));
+ Float_t x2 = cpvDigit2->GetXpad() ;
+ Float_t z2 = cpvDigit2->GetYpad() ;
+ if (x1==x2 && z1==z2) {
+ Float_t qsumpad = cpvDigit1->GetQpad() + cpvDigit2->GetQpad() ;
+ cpvDigit2->SetQpad(qsumpad) ;
+ cpvDigits->RemoveAt(idigit) ;
+ }
}
-
- // Digitize the current CPV hit:
-
- // 1. find pad response and
+ }
+ cpvDigits->Compress() ;
+
+ // 3. add digits to temporary hit list fTmpHits
+
+ ndigits = cpvDigits->GetEntriesFast();
+ for (idigit=0; idigit<ndigits; idigit++) {
+ AliPHOSCPVDigit *cpvDigit = dynamic_cast<AliPHOSCPVDigit*>(cpvDigits->UncheckedAt(idigit));
+ relid[0] = moduleNumber + 1 ; // CPV (or PHOS) module number
+ relid[1] =-1 ; // means CPV
+ relid[2] = cpvDigit->GetXpad() ; // column number of a pad
+ relid[3] = cpvDigit->GetYpad() ; // row number of a pad
- TClonesArray *cpvDigits = new TClonesArray("AliPHOSCPVDigit",0); // array of digits for current hit
- CPVDigitize(pmom,xyd,moduleNumber,cpvDigits);
+ // get the absolute Id number
+ GetGeometry()->RelToAbsNumbering(relid, absid) ;
- Float_t xmean = 0;
- Float_t zmean = 0;
- Float_t qsum = 0;
- Int_t idigit,ndigits;
-
- // 2. go through the current digit list and sum digits in pads
-
- ndigits = cpvDigits->GetEntriesFast();
- for (idigit=0; idigit<ndigits-1; idigit++) {
- AliPHOSCPVDigit *cpvDigit1 = (AliPHOSCPVDigit*) cpvDigits->UncheckedAt(idigit);
- Float_t x1 = cpvDigit1->GetXpad() ;
- Float_t z1 = cpvDigit1->GetYpad() ;
- for (Int_t jdigit=idigit+1; jdigit<ndigits; jdigit++) {
- AliPHOSCPVDigit *cpvDigit2 = (AliPHOSCPVDigit*) cpvDigits->UncheckedAt(jdigit);
- Float_t x2 = cpvDigit2->GetXpad() ;
- Float_t z2 = cpvDigit2->GetYpad() ;
- if (x1==x2 && z1==z2) {
- Float_t qsum = cpvDigit1->GetQpad() + cpvDigit2->GetQpad() ;
- cpvDigit2->SetQpad(qsum) ;
- cpvDigits->RemoveAt(idigit) ;
- }
- }
- }
- cpvDigits->Compress() ;
+ // add current digit to the temporary hit list
- // 3. add digits to temporary hit list fTmpHits
+ xyzte[3] = gMC->TrackTime() ;
+ xyzte[4] = cpvDigit->GetQpad() ; // amplitude in a pad
- ndigits = cpvDigits->GetEntriesFast();
- for (idigit=0; idigit<ndigits; idigit++) {
- AliPHOSCPVDigit *cpvDigit = (AliPHOSCPVDigit*) cpvDigits->UncheckedAt(idigit);
- relid[0] = moduleNumber + 1 ; // CPV (or PHOS) module number
- relid[1] =-1 ; // means CPV
- relid[2] = cpvDigit->GetXpad() ; // column number of a pad
- relid[3] = cpvDigit->GetYpad() ; // row number of a pad
-
- // get the absolute Id number
- fGeom->RelToAbsNumbering(relid, absid) ;
-
- // add current digit to the temporary hit list
- xyze[0] = 0. ;
- xyze[1] = 0. ;
- xyze[2] = 0. ;
- xyze[3] = cpvDigit->GetQpad() ; // amplitude in a pad
- primary = -1; // No need in primary for CPV
- AddHit(fIshunt, primary, tracknumber, absid, xyze, trackpid);
-
- if (cpvDigit->GetQpad() > 0.02) {
- xmean += cpvDigit->GetQpad() * (cpvDigit->GetXpad() + 0.5);
- zmean += cpvDigit->GetQpad() * (cpvDigit->GetYpad() + 0.5);
- qsum += cpvDigit->GetQpad();
- }
+ Int_t primary = gAlice->GetMCApp()->GetPrimary( gAlice->GetMCApp()->GetCurrentTrackNumber() );
+ AddHit(fIshunt, primary, absid, xyzte);
+
+ if (cpvDigit->GetQpad() > 0.02) {
+ xmean += cpvDigit->GetQpad() * (cpvDigit->GetXpad() + 0.5);
+ zmean += cpvDigit->GetQpad() * (cpvDigit->GetYpad() + 0.5);
+ qsum += cpvDigit->GetQpad();
}
+ }
+ if (cpvDigits) {
+ cpvDigits->Delete();
delete cpvDigits;
+ cpvDigits=0;
}
- } // end of IHEP configuration
-
- if(gMC->CurrentVolID(copy) == gMC->VolId("PXTL") ) { // We are inside a PBWO crystal
- gMC->TrackPosition(pos) ;
- xyze[0] = pos[0] ;
- xyze[1] = pos[1] ;
- xyze[2] = pos[2] ;
- xyze[3] = gMC->Edep() ;
-
- // Track enters to the crystal from the top edge
-
- if (gMC->IsTrackEntering()) {
- Float_t posloc[3];
- gMC -> Gmtod (xyze, posloc, 1);
- if (posloc[1] > fGeom->GetCrystalSize(1)/2-0.01) {
- Int_t row,cel;
- Float_t xyd[3]={0,0,0};
- AliPHOSv1 &phos = *(AliPHOSv1*)gAlice->GetModule("PHOS");
-
- Int_t moduleNumber;
- gMC->CurrentVolOffID(10,moduleNumber);
- if ( name == "MIXT" && strcmp(gMC->CurrentVolOffName(10),"PHO1") == 0 )
- moduleNumber += fGeom->GetNModules() - fGeom->GetNPPSDModules();
- moduleNumber--;
-
- gMC->CurrentVolOffID(4, row) ; // get the row number inside the module
- gMC->CurrentVolOffID(3, cel) ; // get the cell number inside the module
- xyd[0] = -(posloc[2] + (cel-0.5-fGeom->GetNZ() /2) *
- (fGeom->GetCrystalSize(2) + 2 * fGeom->GetGapBetweenCrystals()));
- xyd[1] = posloc[0] + (row-0.5-fGeom->GetNPhi()/2) *
- (fGeom->GetCrystalSize(0) + 2 * fGeom->GetGapBetweenCrystals());
-
- // Current momentum of the hit's track in the CPV module ref. system
-
- TLorentzVector pmom;
- gMC -> TrackMomentum(pmom);
- Float_t pm[3], pd[3];
- for (i=0; i<3; i++) pm[i] = pmom[i];
- gMC -> Gmtod (pm, pd, 2); // transform 3-momentum from master to daughter system
- pmom[0] = pd[0];
- pmom[1] =-pd[1];
- pmom[2] =-pd[2];
-
- // Current particle type of the hit's track
-
- Int_t ipart = gMC->TrackPid();
+ }
- // Add the current particle in the list of the EMC hits.
+
+ static Int_t idPXTL = gMC->VolId("PXTL");
+ if(gMC->CurrentVolID(copy) == idPXTL ) { // We are inside a PBWO crystal
- phos.GetEMCModule(moduleNumber).AddHit(fIshunt,primary,pmom,xyd,ipart);
+ gMC->TrackPosition(pos) ;
+ xyzte[0] = pos[0] ;
+ xyzte[1] = pos[1] ;
+ xyzte[2] = pos[2] ;
- if (fDebug == 1) {
- printf("EMC hit added to module #%2d: p = (% .4f, % .4f, % .4f, % .4f) GeV,\n",
- moduleNumber+1,pmom.Px(),pmom.Py(),pmom.Pz(),pmom.E());
- printf( " xy = (%8.4f, %8.4f) cm, ipart = %d, primary = %d\n",
- xyd[0],xyd[1],ipart,primary);
+ Float_t lostenergy = gMC->Edep();
+
+ //Put in the TreeK particle entering PHOS and all its parents
+ if ( gMC->IsTrackEntering() ){
+ Float_t xyzd[3] ;
+ gMC -> Gmtod (xyzte, xyzd, 1); // transform coordinate from master to daughter system
+ if (xyzd[1] < -GetGeometry()->GetCrystalSize(1)/2.+0.1){ //Entered close to forward surface
+ Int_t parent = gAlice->GetMCApp()->GetCurrentTrackNumber() ;
+ TParticle * part = gAlice->GetMCApp()->Particle(parent) ;
+ Float_t vert[3],vertd[3] ;
+ vert[0]=part->Vx() ;
+ vert[1]=part->Vy() ;
+ vert[2]=part->Vz() ;
+ gMC -> Gmtod (vert, vertd, 1); // transform coordinate from master to daughter system
+ if(vertd[1]<-GetGeometry()->GetCrystalSize(1)/2.-0.1){ //Particle is created in foront of PHOS
+ //0.1 to get rid of numerical errors
+ part->SetBit(kKeepBit);
+ while ( parent != -1 ) {
+ part = gAlice->GetMCApp()->Particle(parent) ;
+ part->SetBit(kKeepBit);
+ parent = part->GetFirstMother() ;
+ }
}
}
}
+ if ( lostenergy != 0 ) { // Track is inside the crystal and deposits some energy
+ xyzte[3] = gMC->TrackTime() ;
+
+ gMC->CurrentVolOffID(10, moduleNumber) ; // get the PHOS module number ;
+
+ Int_t strip ;
+ gMC->CurrentVolOffID(3, strip);
+ Int_t cell ;
+ gMC->CurrentVolOffID(2, cell);
- // Track is inside the crystal and deposits some energy
+ //Old formula for row is wrong. For example, I have strip 56 (28 for 2 x 8), row must be 1.
+ //But row == 1 + 56 - 56 % 56 == 57 (row == 1 + 28 - 28 % 28 == 29)
+ //Int_t row = 1 + GetGeometry()->GetEMCAGeometry()->GetNStripZ() - strip % (GetGeometry()->GetEMCAGeometry()->GetNStripZ()) ;
+ Int_t row = GetGeometry()->GetEMCAGeometry()->GetNStripZ() - (strip - 1) % (GetGeometry()->GetEMCAGeometry()->GetNStripZ()) ;
+ Int_t col = (Int_t) TMath::Ceil((Double_t) strip/(GetGeometry()->GetEMCAGeometry()->GetNStripZ())) -1 ;
+
+ // Absid for 8x2-strips. Looks nice :)
+ absid = (moduleNumber-1)*GetGeometry()->GetNCristalsInModule() +
+ row * 2 + (col*GetGeometry()->GetEMCAGeometry()->GetNCellsXInStrip() + (cell - 1) / 2)*GetGeometry()->GetNZ() - (cell & 1 ? 1 : 0);
- if ( xyze[3] != 0 ) {
- gMC->CurrentVolOffID(10, relid[0]) ; // get the PHOS module number ;
- if ( name == "MIXT" && strcmp(gMC->CurrentVolOffName(5),"PHO1") == 0 )
- relid[0] += fGeom->GetNModules() - fGeom->GetNPPSDModules();
- relid[1] = 0 ; // means PBW04
- 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
- fGeom->RelToAbsNumbering(relid, absid) ;
+ //Calculates the light yield, the number of photons produced in the
+ //crystal
+ //There is no dependence of reponce on distance from energy deposition to APD
+ Float_t lightYield = gRandom->Poisson(fLightFactor * lostenergy) ;
- // add current hit to the hit list
+ //Calculates de energy deposited in the crystal
+ xyzte[4] = fAPDFactor * lightYield ;
- AddHit(fIshunt, primary,tracknumber, absid, xyze, trackpid);
+ Int_t primary ;
+ if(fIshunt == 2){
+ primary = gAlice->GetMCApp()->GetCurrentTrackNumber() ;
+ TParticle * part = gAlice->GetMCApp()->Particle(primary) ;
+ while ( !part->TestBit(kKeepBit) ) {
+ primary = part->GetFirstMother() ;
+ if(primary == -1){
+ primary = gAlice->GetMCApp()->GetPrimary( gAlice->GetMCApp()->GetCurrentTrackNumber() );
+ break ; //there is a possibility that particle passed e.g. thermal isulator and hits a side
+ //surface of the crystal. In this case it may have no primary at all.
+ //We can not easily separate this case from the case when this is part of the shower,
+ //developed in the neighboring crystal.
+ }
+ part = gAlice->GetMCApp()->Particle(primary) ;
+ }
+ }
+ else
+ primary = gAlice->GetMCApp()->GetPrimary( gAlice->GetMCApp()->GetCurrentTrackNumber() );
+
+
+ // add current hit to the hit list
+ // Info("StepManager","%d %d", primary, tracknumber) ;
+ AddHit(fIshunt, primary, absid, xyzte);
+
} // there is deposited energy
} // we are inside a PHOS Xtal
+
}
//____________________________________________________________________________
-void AliPHOSv1::CPVDigitize (TLorentzVector p, Float_t *zxhit, Int_t moduleNumber, TClonesArray *cpvDigits)
+void AliPHOSv1::CPVDigitize (TLorentzVector p, Float_t *zxhit, TClonesArray *cpvDigits)
{
// ------------------------------------------------------------------------
// Digitize one CPV hit:
// 2 October 2000
// ------------------------------------------------------------------------
- const Float_t kCelWr = fGeom->GetPadSizePhi()/2; // Distance between wires (2 wires above 1 pad)
+ const Float_t kCelWr = GetGeometry()->GetPadSizePhi()/2; // Distance between wires (2 wires above 1 pad)
const Float_t kDetR = 0.1; // Relative energy fluctuation in track for 100 e-
const Float_t kdEdx = 4.0; // Average energy loss in CPV;
const Int_t kNgamz = 5; // Ionization size in Z
Float_t pNorm = p.Py();
Float_t eloss = kdEdx;
-// cout << "CPVDigitize: YVK : "<<hitX<<" "<<hitZ<<" | "<<pX<<" "<<pZ<<" "<<pNorm<<endl;
+//Info("CPVDigitize", "YVK : %f %f | %f %f %d", hitX, hitZ, pX, pZ, pNorm) ;
- Float_t dZY = pZ/pNorm * fGeom->GetCPVGasThickness();
- Float_t dXY = pX/pNorm * fGeom->GetCPVGasThickness();
+ Float_t dZY = pZ/pNorm * GetGeometry()->GetCPVGasThickness();
+ Float_t dXY = pX/pNorm * GetGeometry()->GetCPVGasThickness();
gRandom->Rannor(rnor1,rnor2);
eloss *= (1 + kDetR*rnor1) *
- TMath::Sqrt((1 + ( pow(dZY,2) + pow(dXY,2) ) / pow(fGeom->GetCPVGasThickness(),2)));
- Float_t zhit1 = hitZ + fGeom->GetCPVActiveSize(1)/2 - dZY/2;
- Float_t xhit1 = hitX + fGeom->GetCPVActiveSize(0)/2 - dXY/2;
+ TMath::Sqrt((1 + ( pow(dZY,2) + pow(dXY,2) ) / pow(GetGeometry()->GetCPVGasThickness(),2)));
+ Float_t zhit1 = hitZ + GetGeometry()->GetCPVActiveSize(1)/2 - dZY/2;
+ Float_t xhit1 = hitX + GetGeometry()->GetCPVActiveSize(0)/2 - dXY/2;
Float_t zhit2 = zhit1 + dZY;
Float_t xhit2 = xhit1 + dXY;
// Finite size of ionization region
- Int_t nCellZ = fGeom->GetNumberOfCPVPadsZ();
- Int_t nCellX = fGeom->GetNumberOfCPVPadsPhi();
+ Int_t nCellZ = GetGeometry()->GetNumberOfCPVPadsZ();
+ Int_t nCellX = GetGeometry()->GetNumberOfCPVPadsPhi();
Int_t nz3 = (kNgamz+1)/2;
Int_t nx3 = (kNgamx+1)/2;
cpvDigits->Expand(nIter*kNgamx*kNgamz);
- TClonesArray &ldigits = *(TClonesArray *)cpvDigits;
+ TClonesArray &ldigits = *(static_cast<TClonesArray *>(cpvDigits));
for (Int_t iter=0; iter<nIter; iter++) {
Float_t zhit = zxe[0][iter];
Float_t xhit = zxe[1][iter];
Float_t qhit = zxe[2][iter];
- Float_t zcell = zhit / fGeom->GetPadSizeZ();
- Float_t xcell = xhit / fGeom->GetPadSizePhi();
+ Float_t zcell = zhit / GetGeometry()->GetPadSizeZ();
+ Float_t xcell = xhit / GetGeometry()->GetPadSizePhi();
if ( zcell<=0 || xcell<=0 ||
zcell>=nCellZ || xcell>=nCellX) return;
Int_t izcell = (Int_t) zcell;
// 3 October 2000
// ------------------------------------------------------------------------
- Double_t dz = fGeom->GetPadSizeZ() / 2;
- Double_t dx = fGeom->GetPadSizePhi() / 2;
- Double_t z = zhit * fGeom->GetPadSizeZ();
- Double_t x = xhit * fGeom->GetPadSizePhi();
+ Double_t dz = GetGeometry()->GetPadSizeZ() / 2;
+ Double_t dx = GetGeometry()->GetPadSizePhi() / 2;
+ Double_t z = zhit * GetGeometry()->GetPadSizeZ();
+ Double_t x = xhit * GetGeometry()->GetPadSizePhi();
Double_t amplitude = qhit *
(CPVCumulPadResponse(z+dz,x+dx) - CPVCumulPadResponse(z+dz,x-dx) -
CPVCumulPadResponse(z-dz,x+dx) + CPVCumulPadResponse(z-dz,x-dx));
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff