From ab37d09c293539da00fc7924d281e50fdfdd8f0a Mon Sep 17 00:00:00 2001 From: mhorner Date: Wed, 9 Jun 2004 19:54:37 +0000 Subject: [PATCH] Fixed bug in parent assignment and implemented Birk's law in the Step Manager --- EMCAL/AliEMCAL.cxx | 6 ++ EMCAL/AliEMCAL.h | 5 +- EMCAL/AliEMCALGeometry.cxx | 28 +++++++++ EMCAL/AliEMCALGeometry.h | 2 + EMCAL/AliEMCALHit.h | 3 + EMCAL/AliEMCALv1.cxx | 117 ++++++++++++++++++++++++++----------- EMCAL/AliEMCALv1.h | 7 ++- 7 files changed, 132 insertions(+), 36 deletions(-) diff --git a/EMCAL/AliEMCAL.cxx b/EMCAL/AliEMCAL.cxx index de0a62278b9..5ecae675616 100644 --- a/EMCAL/AliEMCAL.cxx +++ b/EMCAL/AliEMCAL.cxx @@ -152,6 +152,12 @@ void AliEMCAL::CreateMaterials() gMC->Gstpar(idtmed[1601],"CUTELE",0.001) ; gMC->Gstpar(idtmed[1601],"BCUTE",0.0001) ; + //set constants for Birk's Law implentation + fBirkC0 = 1; + fBirkC1 = 0.013/dP; + fBirkC2 = 9.6e-6/(dP * dP); + + } //____________________________________________________________________________ diff --git a/EMCAL/AliEMCAL.h b/EMCAL/AliEMCAL.h index ef44e103fd4..3e7a20b0aa0 100644 --- a/EMCAL/AliEMCAL.h +++ b/EMCAL/AliEMCAL.h @@ -53,8 +53,11 @@ class AliEMCAL : public AliDetector { protected: AliEMCALGeometry * fGeom ; // the geometry object + Int_t fBirkC0; // constants for Birk's Law implementation + Double_t fBirkC1; + Double_t fBirkC2; - ClassDef(AliEMCAL,5) // Electromagnetic calorimeter (base class) + ClassDef(AliEMCAL,6) // Electromagnetic calorimeter (base class) } ; diff --git a/EMCAL/AliEMCALGeometry.cxx b/EMCAL/AliEMCALGeometry.cxx index 02463440e26..f7d6713bbb7 100644 --- a/EMCAL/AliEMCALGeometry.cxx +++ b/EMCAL/AliEMCALGeometry.cxx @@ -403,3 +403,31 @@ void AliEMCALGeometry::XYZFromIndex(Int_t absid, TVector3 &v) const { return; } + +Bool_t AliEMCALGeometry::IsInEMCAL(Double_t x, Double_t y, Double_t z) const { + // Checks whether point is inside the EMCal volume + // + // Code uses cylindrical approximation made of inner radius (for speed) + // + // Points behind EMCAl, i.e. R > outer radius, but eta, phi in acceptance + // are considered to inside + + Double_t r=sqrt(x*x+y*y); + + if ( r > fEnvelop[0] ) { + Double_t theta; + theta = TMath::ATan2(r,z); + Double_t eta; + if(theta == 0) + eta = 9999; + else + eta = -TMath::Log(TMath::Tan(theta/2.)); + if (eta < fArm1EtaMin || eta > fArm1EtaMax) + return 0; + + Double_t phi = TMath::ATan2(y,x) * 180./TMath::Pi(); + if (phi > fArm1PhiMin && phi < fArm1PhiMax) + return 1; + } + return 0; +} diff --git a/EMCAL/AliEMCALGeometry.h b/EMCAL/AliEMCALGeometry.h index 81af6074606..3dd0be00586 100644 --- a/EMCAL/AliEMCALGeometry.h +++ b/EMCAL/AliEMCALGeometry.h @@ -48,6 +48,8 @@ public: virtual void GetGlobal(const AliRecPoint *, TVector3 &, TMatrix &) const {} virtual void GetGlobal(const AliRecPoint *, TVector3 &) const {} virtual Bool_t Impact(const TParticle *) const {return kTRUE;} + + Bool_t IsInEMCAL(Double_t x, Double_t y, Double_t z) const; // General Bool_t IsInitialized(void) const { return fgInit ; } // Return EMCA geometrical parameters diff --git a/EMCAL/AliEMCALHit.h b/EMCAL/AliEMCALHit.h index c7b16850556..3c7d8e1c736 100644 --- a/EMCAL/AliEMCALHit.h +++ b/EMCAL/AliEMCALHit.h @@ -44,6 +44,9 @@ public: Float_t GetPz(void) const{return fPz;} Float_t GetPe(void) const{return fPe;} + void SetIparent(Int_t i_parent) {fIparent=i_parent;} + void SetPrimary(Int_t primary) {fPrimary=primary;} + Bool_t operator == (AliEMCALHit const &rValue) const; AliEMCALHit operator + (const AliEMCALHit& rValue); diff --git a/EMCAL/AliEMCALv1.cxx b/EMCAL/AliEMCALv1.cxx index 5a9f82b699d..88d973f34b2 100644 --- a/EMCAL/AliEMCALv1.cxx +++ b/EMCAL/AliEMCALv1.cxx @@ -46,14 +46,14 @@ ClassImp(AliEMCALv1) //______________________________________________________________________ -AliEMCALv1::AliEMCALv1():AliEMCALv0(){ +AliEMCALv1::AliEMCALv1():AliEMCALv0(), fCurPrimary(-1), fCurParent(-1), fCurTrack(-1){ // ctor } //______________________________________________________________________ AliEMCALv1::AliEMCALv1(const char *name, const char *title): - AliEMCALv0(name,title){ + AliEMCALv0(name,title), fCurPrimary(-1), fCurParent(-1), fCurTrack(-1) { // Standard Creator. fHits= new TClonesArray("AliEMCALHit",1000); @@ -73,6 +73,7 @@ AliEMCALv1::~AliEMCALv1(){ 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){ @@ -116,42 +117,50 @@ void AliEMCALv1::StepManager(void){ 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 - Float_t depositedEnergy ; if( (depositedEnergy = gMC->Edep()) > 0.){// Track is inside a scintillator and deposits some energy - - // use sampling fraction to get original energy --HG - depositedEnergy = depositedEnergy * geom->GetSampling(); - + if (fCurPrimary==-1) + fCurPrimary=gAlice->GetMCApp()->GetPrimary(tracknumber); + + if (fCurParent==-1 || tracknumber != fCurTrack) { + // Check parentage + Int_t parent=tracknumber; + if (fCurParent != -1) { + while (parent != fCurParent && parent != -1) { + TParticle *part=gAlice->GetMCApp()->Particle(parent); + parent=part->GetFirstMother(); + } + } + if (fCurParent==-1 || parent==-1) { + Int_t parent=tracknumber; + TParticle *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); + 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]; @@ -174,16 +183,56 @@ 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 (gMC->TrackCharge()!=0) { // Check + Float_t BirkC1_mod = 0; + if (fBirkC0==1){ // Apply correction for higher charge states + if (abs(gMC->TrackCharge())>=2) + BirkC1_mod=fBirkC1*7.2/12.6; + else + BirkC1_mod=fBirkC1; + } + Float_t dedxcm; + if (gMC->TrackStep()>0) + dedxcm=1000.*gMC->Edep()/gMC->TrackStep(); + else + dedxcm=0; + lightYield=lightYield/(1.+BirkC1_mod*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 hit_it(Hits()); + Int_t i_hit=0; + while (AliEMCALHit *hit=dynamic_cast(hit_it()) ) { + if (map[hit->GetIparent()]==-99) + cout << "Remapping, found -99 for parent id " << hit->GetIparent() << ", " << map[hit->GetIparent()] << ", i_hit " << i_hit << endl; + hit->SetIparent(map[hit->GetIparent()]); + if (map[hit->GetPrimary()]==-99) + cout << "Remapping, found -99 for primary id " << hit->GetPrimary() << ", " << map[hit->GetPrimary()] << ", i_hit " << i_hit << endl; + hit->SetPrimary(map[hit->GetPrimary()]); + i_hit++; + } +} + +void AliEMCALv1::FinishPrimary() { + fCurPrimary=-1; + fCurParent=-1; + fCurTrack=-1; +} diff --git a/EMCAL/AliEMCALv1.h b/EMCAL/AliEMCALv1.h index a9ac230d787..bd34611f253 100644 --- a/EMCAL/AliEMCALv1.h +++ b/EMCAL/AliEMCALv1.h @@ -36,6 +36,8 @@ public: // Gives the version number virtual Int_t IsVersion(void) const {return 1;} virtual void StepManager(void) ; + virtual void RemapTrackHitIDs(Int_t *map); + virtual void FinishPrimary(); virtual const TString Version(void)const {return TString("v0");} // assignement operator requested by coding convention but not needed AliEMCALv1 & operator = (const AliEMCALv0 & /*rvalue*/){ @@ -44,8 +46,11 @@ public: private: + Int_t fCurPrimary; + Int_t fCurParent; + Int_t fCurTrack; - ClassDef(AliEMCALv1,6)//Implementation of EMCAL manager class to produce hits in a Central Calorimeter + ClassDef(AliEMCALv1,7)//Implementation of EMCAL manager class to produce hits in a Central Calorimeter }; -- 2.43.0