X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=EMCAL%2FAliEMCALv1.cxx;h=f800fe0dc50bb4a465400ccba289f8f1d2f06efa;hb=d569f6a40956200044d1c45a2ad35d635c713097;hp=f6d88c3a9d600c5c4e010778744be65909a773b8;hpb=d75bea678edb2ec6d7f562693b8a9bf8e4dbc943;p=u%2Fmrichter%2FAliRoot.git diff --git a/EMCAL/AliEMCALv1.cxx b/EMCAL/AliEMCALv1.cxx index f6d88c3a9d6..f800fe0dc50 100644 --- a/EMCAL/AliEMCALv1.cxx +++ b/EMCAL/AliEMCALv1.cxx @@ -16,146 +16,247 @@ /* $Id$ */ //_________________________________________________________________________ -// Implementation version v1 of EMCAL Manager class -// An object of this class does not produce digits -// It is the one to use if you do want to produce outputs in TREEH -// +//*-- Implementation version v1 of EMCAL Manager class +//*-- An object of this class does not produce digits +//*-- It is the one to use if you do want to produce outputs in TREEH +//*-- //*-- Author: Sahal Yacoob (LBL /UCT) //*-- : Jennifer Klay (LBL) +// This Class not stores information on all particles prior to EMCAL entry - in order to facilitate analysis. +// This is done by setting fIShunt =2, and flagging all parents of particles entering the EMCAL. +// 15/02/2002 .... Yves Schutz +// 1. fSamplingFraction and fLayerToPreshowerRatio have been removed +// 2. Timing signal is collected and added to hit // --- ROOT system --- -#include "TPGON.h" -#include "TTUBS.h" -#include "TNode.h" -#include "TRandom.h" -#include "TTree.h" -#include "TGeometry.h" - +#include +#include +#include // --- Standard library --- -#include -#include -#include -#include -#include -#include // --- AliRoot header files --- - #include "AliEMCALv1.h" #include "AliEMCALHit.h" #include "AliEMCALGeometry.h" -#include "AliConst.h" #include "AliRun.h" #include "AliMC.h" +#include "AliStack.h" ClassImp(AliEMCALv1) //______________________________________________________________________ -AliEMCALv1::AliEMCALv1():AliEMCALv0(){ - // ctor +AliEMCALv1::AliEMCALv1() + : AliEMCALv0(), + fCurPrimary(-1), + fCurParent(-1), + fCurTrack(-1), + fTimeCut(30e-09) +{ + // default ctor } + //______________________________________________________________________ -AliEMCALv1::AliEMCALv1(const char *name, const char *title): - AliEMCALv0(name,title){ +AliEMCALv1::AliEMCALv1(const char *name, const char *title) + : AliEMCALv0(name,title), + fCurPrimary(-1), + fCurParent(-1), + fCurTrack(-1), + fTimeCut(30e-09) +{ // Standard Creator. fHits= new TClonesArray("AliEMCALHit",1000); - gAlice->AddHitList(fHits); + // gAlice->GetMCApp()->AddHitList(fHits); // 20-dec-04 - advice of Andreas fNhits = 0; - fSamplingFraction = 12.9 ; - fLayerToPreshowerRatio = 5.0/6.0; - fIshunt = 1; // All hits are associated with primary particles + fIshunt = 2; // All hits are associated with particles entering the calorimeter } + //______________________________________________________________________ AliEMCALv1::~AliEMCALv1(){ - // dtor - - if ( fHits) { - fHits->Delete(); - delete fHits; - fHits = 0; - } + // dtor + + if ( fHits ) { + fHits->Clear(); + delete fHits; + fHits = 0; + } } + //______________________________________________________________________ void AliEMCALv1::AddHit(Int_t shunt, Int_t primary, Int_t tracknumber, Int_t iparent, Float_t ienergy, Int_t id, Float_t * hits,Float_t * p){ // Add a hit to the hit list. - // An EMCAL hit is the sum of all hits in a single segment - // originating from the same enterring particle - Int_t hitCounter; + // An EMCAL hit is the sum of all hits in a tower section + // originating from the same entering particle + Int_t hitCounter=0; - AliEMCALHit *newHit; - AliEMCALHit *curHit; + AliEMCALHit *newHit=0; + AliEMCALHit *curHit=0; Bool_t deja = kFALSE; newHit = new AliEMCALHit(shunt, primary, tracknumber, iparent, ienergy, id, hits, p); - for ( hitCounter = fNhits-1; hitCounter >= 0 && !deja; hitCounter-- ) { + for ( hitCounter = fNhits-1; hitCounter >= 0 && !deja; hitCounter-- ) { curHit = (AliEMCALHit*) (*fHits)[hitCounter]; // We add hits with the same tracknumber, while GEANT treats // primaries succesively - if(curHit->GetPrimary() != primary) break; + if(curHit->GetPrimary() != primary) + break; if( *curHit == *newHit ) { *curHit = *curHit + *newHit; deja = kTRUE; - } // end if + } // end if } // end for hitCounter - + if ( !deja ) { new((*fHits)[fNhits]) AliEMCALHit(*newHit); fNhits++; } // end if - + delete newHit; } //______________________________________________________________________ void AliEMCALv1::StepManager(void){ - // Accumulates hits as long as the track stays in a single - // crystal or PPSD gas Cell - Int_t id[2]; // (layer, phi, Eta) indices - Int_t absid; - // position wrt MRS and energy deposited - Float_t xyze[4]={0.,0.,0.,0.}; - Float_t pmom[4]={0.,0.,0.,0.}; - TLorentzVector pos; // Lorentz vector of the track current position. - TLorentzVector mom; // Lorentz vector of the track current momentum. - Int_t tracknumber = gAlice->CurrentTrack(); - Int_t primary = 0; - static Int_t iparent = 0; - static Float_t ienergy = 0; - Int_t copy = 0; - - - if(gMC->IsTrackEntering() && (strcmp(gMC->CurrentVolName(),"XALU") == 0)) // This Particle in enterring the Calorimeter - { - iparent = tracknumber; - gMC->TrackMomentum(mom); - ienergy = mom[3]; -} - if(gMC->CurrentVolID(copy) == gMC->VolId("XPHI") ) { // We are in a Scintillator Layer - - gMC->CurrentVolOffID(1, id[0]); // get the POLY copy number; - gMC->CurrentVolID(id[1]); // get the phi number inside the layer - primary = gAlice->GetPrimary(tracknumber); - gMC->TrackPosition(pos); - gMC->TrackMomentum(mom); - xyze[0] = pos[0]; - xyze[1] = pos[1]; - xyze[2] = pos[2]; - xyze[3] = fSamplingFraction*(gMC->Edep()); // Correct for sampling calorimeter - pmom[0] = mom[0]; - pmom[1] = mom[1]; - pmom[2] = mom[2]; - pmom[3] = mom[3]; + // Accumulates hits as long as the track stays in a tower + + Int_t id[2]={0,0}; // (phi, Eta) indices + // position wrt MRS and energy deposited + Float_t xyzte[5]={0.,0.,0.,0.,0.};// position wrt MRS, time and energy deposited + Float_t pmom[4]={0.,0.,0.,0.}; + TLorentzVector pos; // Lorentz vector of the track current position. + TLorentzVector mom; // Lorentz vector of the track current momentum. + Int_t tracknumber = gAlice->GetMCApp()->GetCurrentTrackNumber(); + static Float_t ienergy = 0; + Int_t copy = 0; + + AliEMCALGeometry * geom = GetGeometry() ; + + TParticle *part=0; + Int_t parent=-1; + + static Int_t idXPHI = gMC->VolId("XPHI"); + if(gMC->CurrentVolID(copy) == idXPHI ) { // We are in a Scintillator Layer + Float_t depositedEnergy ; - if(xyze[3] > 0.){// Track is inside the crystal and deposits some energy - absid = (id[0]-1)*(fGeom->GetNPhi()) + id[1]; - if((absid/fGeom->GetNPhi()) < (2*fGeom->GetNZ())) - {xyze[3] = fLayerToPreshowerRatio*xyze[3] ;} // Preshower readout must be scaled - AddHit(fIshunt, primary,tracknumber, iparent, ienergy, absid, xyze, pmom); + if( ((depositedEnergy = gMC->Edep()) > 0.) && (gMC->TrackTime() < fTimeCut)){// Track is inside a scintillator and deposits some energy + if (fCurPrimary==-1) + fCurPrimary=gAlice->GetMCApp()->GetPrimary(tracknumber); + + if (fCurParent==-1 || tracknumber != fCurTrack) { + // Check parentage + //Int_t parent=tracknumber; + parent=tracknumber; + if (fCurParent != -1) { + while (parent != fCurParent && parent != -1) { + //TParticle *part=gAlice->GetMCApp()->Particle(parent); + part=gAlice->GetMCApp()->Particle(parent); + parent=part->GetFirstMother(); + } + } + if (fCurParent==-1 || parent==-1) { + //Int_t parent=tracknumber; + //TParticle *part=gAlice->GetMCApp()->Particle(parent); + parent=tracknumber; + part=gAlice->GetMCApp()->Particle(parent); + while (parent != -1 && geom->IsInEMCAL(part->Vx(),part->Vy(),part->Vz())) { + parent=part->GetFirstMother(); + if (parent!=-1) + part=gAlice->GetMCApp()->Particle(parent); + } + fCurParent=parent; + if (fCurParent==-1) + Error("StepManager","Cannot find parent"); + else { + //TParticle *part=gAlice->GetMCApp()->Particle(fCurParent); + part=gAlice->GetMCApp()->Particle(fCurParent); + ienergy = part->Energy(); + } + while (parent != -1) { + part=gAlice->GetMCApp()->Particle(parent); + part->SetBit(kKeepBit); + parent=part->GetFirstMother(); + } + } + fCurTrack=tracknumber; + } + gMC->TrackPosition(pos); + xyzte[0] = pos[0]; + xyzte[1] = pos[1]; + xyzte[2] = pos[2]; + xyzte[3] = gMC->TrackTime() ; + + gMC->TrackMomentum(mom); + pmom[0] = mom[0]; + pmom[1] = mom[1]; + pmom[2] = mom[2]; + pmom[3] = mom[3]; + + gMC->CurrentVolOffID(1, id[0]); // get the POLY copy number; + gMC->CurrentVolID(id[1]); // get the phi number inside the layer + + Int_t tower = (id[0]-1) % geom->GetNZ() + 1 + (id[1] - 1) * geom->GetNZ() ; + Int_t layer = static_cast((id[0]-1)/(geom->GetNZ())) + 1 ; + Int_t absid = tower; + Int_t nlayers = geom->GetNECLayers(); + if ((layer > nlayers)||(layer<1)) + Fatal("StepManager", "Wrong calculation of layer number: layer = %d > %d\n", layer, nlayers) ; + + Float_t lightYield = depositedEnergy ; + // Apply Birk's law (copied from G3BIRK) + + if (gMC->TrackCharge()!=0) { // Check + Float_t birkC1Mod = 0; + if (fBirkC0==1){ // Apply correction for higher charge states + if (TMath::Abs(gMC->TrackCharge())>=2) + birkC1Mod=fBirkC1*7.2/12.6; + else + birkC1Mod=fBirkC1; + } + Float_t dedxcm=0.; + if (gMC->TrackStep()>0) + dedxcm=1000.*gMC->Edep()/gMC->TrackStep(); + else + dedxcm=0; + lightYield=lightYield/(1.+birkC1Mod*dedxcm+fBirkC2*dedxcm*dedxcm); + } + + // use sampling fraction to get original energy --HG + xyzte[4] = lightYield * geom->GetSampling(); + + if (gDebug == 2) + printf("StepManager: id0 = %d, id1 = %d, absid = %d tower = %d layer = %d energy = %f\n", id[0], id[1], absid, tower, layer, xyzte[4]) ; + + AddHit(fIshunt, fCurPrimary,tracknumber, fCurParent, ienergy, absid, xyzte, pmom); } // there is deposited energy + } } + +//___________________________________________________________ +void AliEMCALv1::RemapTrackHitIDs(Int_t *map) { + // remap track index numbers for primary and parent indices + // (Called by AliStack::PurifyKine) + if (Hits()==0) + return; + TIter hitIter(Hits()); + Int_t iHit=0; + while (AliEMCALHit *hit=dynamic_cast(hitIter()) ) { + if (map[hit->GetIparent()]==-99) + cout << "Remapping, found -99 for parent id " << hit->GetIparent() << ", " << map[hit->GetIparent()] << ", iHit " << iHit << endl; + hit->SetIparent(map[hit->GetIparent()]); + if (map[hit->GetPrimary()]==-99) + cout << "Remapping, found -99 for primary id " << hit->GetPrimary() << ", " << map[hit->GetPrimary()] << ", iHit " << iHit << endl; + hit->SetPrimary(map[hit->GetPrimary()]); + iHit++; + } +} + +//___________________________________________________________ +void AliEMCALv1::FinishPrimary() { + // finish primary + fCurPrimary=-1; + fCurParent=-1; + fCurTrack=-1; }