X-Git-Url: http://git.uio.no/git/?p=u%2Fmrichter%2FAliRoot.git;a=blobdiff_plain;f=EMCAL%2FAliEMCALv1.cxx;h=2bdd6a8d1a4a4e0b5a6449e5a5d0ca5312a0e684;hp=d0fd8e2c4f9e4f961c34acc20528ca7db867ffa8;hb=2942f542d4b5bec82bbddfadd400306661f8dc8d;hpb=fdebddebb521c2f316af46d2034286ce5100c285 diff --git a/EMCAL/AliEMCALv1.cxx b/EMCAL/AliEMCALv1.cxx index d0fd8e2c4f9..2bdd6a8d1a4 100644 --- a/EMCAL/AliEMCALv1.cxx +++ b/EMCAL/AliEMCALv1.cxx @@ -30,8 +30,9 @@ // 2. Timing signal is collected and added to hit // --- ROOT system --- -#include "TParticle.h" -#include "TVirtualMC.h" +#include +#include +#include // --- Standard library --- @@ -41,23 +42,37 @@ #include "AliEMCALGeometry.h" #include "AliRun.h" #include "AliMC.h" +#include "AliStack.h" +using std::cout; +using std::endl; 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, + const Bool_t checkGeoAndRun) + : AliEMCALv0(name,title,checkGeoAndRun), + fCurPrimary(-1), + fCurParent(-1), + fCurTrack(-1), + fTimeCut(30e-09) +{ // Standard Creator. fHits= new TClonesArray("AliEMCALHit",1000); - gAlice->GetMCApp()->AddHitList(fHits); + // gAlice->GetMCApp()->AddHitList(fHits); // 20-dec-04 - advice of Andreas fNhits = 0; fIshunt = 2; // All hits are associated with particles entering the calorimeter @@ -65,24 +80,25 @@ AliEMCALv1::AliEMCALv1(const char *name, const char *title): //______________________________________________________________________ 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 tower section // originating from the same entering particle - Int_t hitCounter; + 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); @@ -102,67 +118,87 @@ void AliEMCALv1::AddHit(Int_t shunt, Int_t primary, Int_t tracknumber, Int_t ipa new((*fHits)[fNhits]) AliEMCALHit(*newHit); fNhits++; } // end if - + delete newHit; } //______________________________________________________________________ void AliEMCALv1::StepManager(void){ // Accumulates hits as long as the track stays in a tower - Int_t id[2]; // (phi, Eta) indices + 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(); - Int_t primary = 0; - static Int_t iparent = 0; static Float_t ienergy = 0; Int_t copy = 0; AliEMCALGeometry * geom = GetGeometry() ; - if(gMC->IsTrackEntering() && (strcmp(gMC->CurrentVolName(),"XALU") == 0)){ // This Particle in enterring the Calorimeter - gMC->TrackPosition(pos) ; - xyzte[0] = pos[0] ; - xyzte[1] = pos[1] ; - xyzte[2] = pos[2] ; - if ( (xyzte[0]*xyzte[0] + xyzte[1]*xyzte[1]) - < (geom->GetEnvelop(0)+geom->GetGap2Active()+1.5 )*(geom->GetEnvelop(0)+geom->GetGap2Active()+1.5 ) ) { - iparent = tracknumber; - gMC->TrackMomentum(mom); - ienergy = mom[3]; - TParticle * part = 0 ; - Int_t parent = iparent ; - while ( parent != -1 ) { // <------------- flags this particle to be kept and - //all the ancestors of this particle - part = gAlice->GetMCApp()->Particle(parent) ; - part->SetBit(kKeepBit); - parent = part->GetFirstMother() ; - } - } - } - if(gMC->CurrentVolID(copy) == gMC->VolId("XPHI") ) { // We are in a Scintillator Layer - + TParticle *part=0; + Int_t parent=-1; + + static Int_t idXPHI = TVirtualMC::GetMC()->VolId("XPHI"); + if(TVirtualMC::GetMC()->CurrentVolID(copy) == idXPHI ) { // We are in a Scintillator Layer Float_t depositedEnergy ; - if( (depositedEnergy = gMC->Edep()) > 0.){// Track is inside a scintillator and deposits some energy - - gMC->TrackPosition(pos); + if( ((depositedEnergy = TVirtualMC::GetMC()->Edep()) > 0.) && (TVirtualMC::GetMC()->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; + } + TVirtualMC::GetMC()->TrackPosition(pos); xyzte[0] = pos[0]; xyzte[1] = pos[1]; xyzte[2] = pos[2]; - xyzte[3] = gMC->TrackTime() ; + xyzte[3] = TVirtualMC::GetMC()->TrackTime() ; - gMC->TrackMomentum(mom); + TVirtualMC::GetMC()->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 + TVirtualMC::GetMC()->CurrentVolOffID(1, id[0]); // get the POLY copy number; + TVirtualMC::GetMC()->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 ; @@ -171,16 +207,59 @@ void AliEMCALv1::StepManager(void){ if ((layer > nlayers)||(layer<1)) Fatal("StepManager", "Wrong calculation of layer number: layer = %d > %d\n", layer, nlayers) ; - Float_t lightYield = depositedEnergy ; - xyzte[4] = lightYield ; - - primary = gAlice->GetMCApp()->GetPrimary(tracknumber); + // Apply Birk's law (copied from G3BIRK) + if (TVirtualMC::GetMC()->TrackCharge()!=0) { // Check + Float_t birkC1Mod = 0; + if (fBirkC0==1){ // Apply correction for higher charge states + if (TMath::Abs(TVirtualMC::GetMC()->TrackCharge())>=2) + birkC1Mod=fBirkC1*7.2/12.6; + else + birkC1Mod=fBirkC1; + } + Float_t dedxcm=0.; + if (TVirtualMC::GetMC()->TrackStep()>0) + dedxcm=1000.*TVirtualMC::GetMC()->Edep()/TVirtualMC::GetMC()->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, primary,tracknumber, iparent, ienergy, absid, xyzte, pmom); + 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; +}