kaliveda-dev team mailing list archive
-
kaliveda-dev team
-
Mailing list archive
-
Message #01003
[Merge] lp:~frankland/kaliveda/phd_rustine into lp:kaliveda
John Frankland has proposed merging lp:~frankland/kaliveda/phd_rustine into lp:kaliveda.
Requested reviews:
John Frankland (frankland)
For more details, see:
https://code.launchpad.net/~frankland/kaliveda/phd_rustine/+merge/156519
corrected implementation of PHD for silicon detectors
corrects KVDetector::GetCorrectedEnergy: we increase first the mass, then the charge, of ions in case their
measured dE is > the maximum possible for their apparent (Z,A)
in KVSelector analysis, we automatically apply a correction of PHD for transmitted ions with Z>10 which reach the CsI on rings 1-9 (in fact, we re-calibrate them)
--
https://code.launchpad.net/~frankland/kaliveda/phd_rustine/+merge/156519
Your team KaliVeda Development Team is subscribed to branch lp:kaliveda.
=== modified file 'KVIndra/analysis/KVINDRAReconDataAnalyser.cpp'
--- KVIndra/analysis/KVINDRAReconDataAnalyser.cpp 2012-07-17 17:31:04 +0000
+++ KVIndra/analysis/KVINDRAReconDataAnalyser.cpp 2013-04-02 10:59:30 +0000
@@ -312,7 +312,6 @@
void KVINDRAReconDataAnalyser::preAnalysis()
{
// Read and set raw data for the current reconstructed event
-
if(!theRawData) return;
// all recon events are numbered 1, 2, ... : therefore entry number is N-1
Long64_t rawEntry = fSelector->GetEventNumber() - 1;
@@ -322,7 +321,7 @@
theRawData->GetEntry(rawEntry);
for(int i=0; i<NbParFired; i++){
KVACQParam* par = gIndra->GetACQParam((*parList)[ParNum[i]]->GetName());
- if(par) par->SetData(ParVal[i]);
+ if(par) {par->SetData(ParVal[i]);}
}
}
=== modified file 'KVIndra/analysis/KVSelector.cpp'
--- KVIndra/analysis/KVSelector.cpp 2012-06-14 14:10:31 +0000
+++ KVIndra/analysis/KVSelector.cpp 2013-04-02 10:59:30 +0000
@@ -349,8 +349,18 @@
pid.fCpuUser, pid.fMemVirtual/1024., disk.Data());
}
}
+
+ // read event
fChain->GetTree()->GetEntry(fTreeEntry);
- gDataAnalyser->preAnalysis();
+ // read raw data associated to event
+ gDataAnalyser->preAnalysis();
+ if(KVINDRAReconNuc::PHDNeedCorrection){
+ // correct PHD/silicon energies (rustine)
+ KVINDRAReconNuc* nn=0;
+ while( (nn = (KVINDRAReconNuc*)GetEvent()->GetNextParticle()) ){
+ nn->CorrectPHD();
+ }
+ }
//additional selection criteria ?
if(fPartCond){
=== modified file 'KVIndra/detectors/KVCsI.cpp'
--- KVIndra/detectors/KVCsI.cpp 2013-03-08 15:00:51 +0000
+++ KVIndra/detectors/KVCsI.cpp 2013-04-02 10:59:30 +0000
@@ -424,7 +424,7 @@
//______________________________________________________________________________
-Double_t KVCsI::GetCorrectedEnergy(const KVNucleus* nuc, Double_t lum, Bool_t trans)
+Double_t KVCsI::GetCorrectedEnergy(KVNucleus* nuc, Double_t lum, Bool_t trans)
{
// Calculate calibrated energy loss for a nucleus (Z,A) giving total light output "lum".
// If "lum" is not given, the total light of the detector
=== modified file 'KVIndra/detectors/KVCsI.h'
--- KVIndra/detectors/KVCsI.h 2013-03-08 15:00:51 +0000
+++ KVIndra/detectors/KVCsI.h 2013-04-02 10:59:30 +0000
@@ -86,7 +86,7 @@
void SetACQParams();
void SetCalibrators();
- Double_t GetCorrectedEnergy(const KVNucleus*, Double_t lum = -1., Bool_t transmission=kTRUE);
+ Double_t GetCorrectedEnergy(KVNucleus *, Double_t lum = -1., Bool_t transmission=kTRUE);
Double_t GetLightFromEnergy(Int_t Z, Int_t A, Double_t E = -1.);
void SetPinLaser(Int_t n){ if(n>0&&n<255) fPinLaser = (Char_t)n; };
=== modified file 'KVIndra/particles/KVINDRAReconNuc.cpp'
--- KVIndra/particles/KVINDRAReconNuc.cpp 2013-01-17 10:16:09 +0000
+++ KVIndra/particles/KVINDRAReconNuc.cpp 2013-04-02 10:59:30 +0000
@@ -104,6 +104,8 @@
//IN ALL CASES THE RETURNED VALUE OF GetA() IS POSITIVE
//
+Bool_t KVINDRAReconNuc::PHDNeedCorrection = kFALSE;
+
void KVINDRAReconNuc::init()
{
//default initialisations
@@ -113,6 +115,8 @@
fPileup=kFALSE;
fUseFullChIoEnergyForCalib=kTRUE;
fECsI=fESi=fEChIo=0.;
+ fESi_old=fEnergy_old=0.;
+ fCorrectPHD=kFALSE;
}
KVINDRAReconNuc::KVINDRAReconNuc():fCodes()
@@ -807,6 +811,7 @@
else
fECsI = GetCsI()->GetCorrectedEnergy(this, -1., kFALSE);
if(fECsI<=0){
+ //Info("Calib", "ECsI = %f",fECsI);
SetECode(kECode15);// bad - no CsI energy
return;
}
@@ -939,3 +944,65 @@
return idtel->GetIDSubCodeString(code);
}
+
+//______________________________________________________________________________
+
+void KVINDRAReconNuc::Streamer(TBuffer &R__b)
+{
+ // Stream an object of class KVINDRAReconNuc.
+ //
+ // Sets flag for correction of PHD for transmitted ions with Z>10
+ // which reach the CsI on rings 1-9 (see CorrectPHD)
+
+ UInt_t R__s, R__c;
+ if (R__b.IsReading()) {
+ fESi_old=fEnergy_old=0.;
+ fCorrectPHD=kFALSE;
+ Version_t R__v = R__b.ReadVersion(&R__s, &R__c);
+ R__b.ReadClassBuffer(KVINDRAReconNuc::Class(),this,R__v,R__s,R__c);
+ if( IsIdentified() && IsCalibrated()
+ && GetRingNumber()<10 && GetZ()>10 && StoppedInCsI() && R__v < 11 ) {
+ fCorrectPHD=kTRUE;
+ if(!PHDNeedCorrection){
+ Info("Streamer", "APPLYING PHD CORRECTION TO DATA WRITTEN WITH KALIVEDA <v1.8.10");
+ PHDNeedCorrection=kTRUE;
+ }
+ }
+ } else {
+ R__b.WriteClassBuffer(KVINDRAReconNuc::Class(),this);
+ }
+}
+
+void KVINDRAReconNuc::CorrectPHD()
+{
+ // Implements a 'rustine' for correction of PHD in silicon
+ // For KaliVeda < v1.8.10 (KVINDRAReconNuc class version < 11)
+ // the PHD was calculated incorrectly for transmitted ions,
+ // the Moulton formula was applied to the incident energy of the ion,
+ // instead of the energy loss of the ion in the silicon.
+ // Here we perform a 'correction' for all ions with Z>10
+ // which reach the CsI on rings 1-9 (no PHD correction was applied
+ // for ions with Z<=10). The 'correct' PHD correction
+ // is applied to the Si energy, leading to a new deduced total energy
+ // of the ion with a new correction for target energy losses.
+ //
+ // This method is called by KVSelector::Process, so that fragments
+ // used in data analysis have correct energies
+
+ if(!fCorrectPHD) return;
+ KVTarget* t = gMultiDetArray->GetTarget();
+ if(t){
+ t->SetIncoming(kFALSE); t->SetOutgoing(kTRUE);
+ }
+ //Info("CorrectPHD", "PHD CORRECTION RUSTINE");
+ //Print();
+ fESi_old=GetEnergySi();
+ fEnergy_old=GetEnergy();
+ /*Info("CorrectPHD",
+ "OLD : Etot = %f Echio = %f Esi = %f Ecsi = %f Etarg = %f",
+ fEnergy_old, GetEnergyChIo(), fESi_old, GetEnergyCsI(), GetTargetEnergyLoss());*/
+ Calibrate();
+ /*Info("CorrectPHD",
+ "NEW : Etot = %f Echio = %f Esi = %f Ecsi = %f Etarg = %f ecod=%d\n",
+ GetEnergy(), GetEnergyChIo(), GetEnergySi(), GetEnergyCsI(), GetTargetEnergyLoss(),(int)GetCodes().GetVedaECode());*/
+}
=== modified file 'KVIndra/particles/KVINDRAReconNuc.h'
--- KVIndra/particles/KVINDRAReconNuc.h 2013-02-13 10:52:19 +0000
+++ KVIndra/particles/KVINDRAReconNuc.h 2013-04-02 10:59:30 +0000
@@ -41,10 +41,16 @@
Float_t fECsI;//csi contribution to energy
Float_t fESi;//si contribution to energy
Float_t fEChIo;//chio contribution to energy
+ /* for PHD corrections */
+ Bool_t fCorrectPHD;//!set to kTRUE in Streamer if PHD needs correction
+ Float_t fESi_old;//!silicon energy before PHD correction
+ Float_t fEnergy_old;//!total energy before PHD correction, including target losses
void CheckCsIEnergy();
public:
+ static Bool_t PHDNeedCorrection;//static flag set when PHD of analysed events need correcting
+ void CorrectPHD();
Int_t GetIDSubCode(const Char_t * id_tel_type,
KVIDSubCode & code) const;
const Char_t *GetIDSubCodeString(const Char_t * id_tel_type,
@@ -153,7 +159,20 @@
Int_t GetIDSubCode(const Char_t * id_tel_type = "") const;
const Char_t *GetIDSubCodeString(const Char_t * id_tel_type = "") const;
- ClassDef(KVINDRAReconNuc, 10) //Nucleus identified by INDRA array
+ Float_t GetOldSiliconEnergy() const
+ {
+ // If silicon energy is corrected by 'rustine' in Streamer method,
+ // this will return the energy before correction of PHD
+ return fESi_old;
+ }
+ Float_t GetOldTotalEnergy() const
+ {
+ // If silicon energy is corrected by 'rustine' in Streamer method,
+ // this will return the total energy (including target losses)
+ // before correction of PHD
+ return fEnergy_old;
+ }
+ ClassDef(KVINDRAReconNuc, 11) //Nucleus identified by INDRA array
};
//____________________________________________________________________________________________//
=== modified file 'KVIndra/particles/KVIndraparticlesLinkDef.h'
--- KVIndra/particles/KVIndraparticlesLinkDef.h 2011-11-29 09:12:57 +0000
+++ KVIndra/particles/KVIndraparticlesLinkDef.h 2013-04-02 10:59:30 +0000
@@ -2,5 +2,5 @@
#pragma link off all globals;
#pragma link off all classes;
#pragma link off all functions;
-#pragma link C++ class KVINDRAReconNuc+;
+#pragma link C++ class KVINDRAReconNuc-;
#endif
=== modified file 'KVMultiDet/detectors/KVDetector.cpp'
--- KVMultiDet/detectors/KVDetector.cpp 2012-10-17 09:33:13 +0000
+++ KVMultiDet/detectors/KVDetector.cpp 2013-04-02 10:59:30 +0000
@@ -803,7 +803,7 @@
//______________________________________________________________________________//
-Double_t KVDetector::GetCorrectedEnergy(const KVNucleus *nuc, Double_t e, Bool_t transmission)
+Double_t KVDetector::GetCorrectedEnergy( KVNucleus *nuc, Double_t e, Bool_t transmission)
{
// Returns the total energy loss in the detector for a given nucleus
// including inactive absorber layers.
@@ -833,6 +833,32 @@
enum SolType solution = kEmax;
if(!transmission) solution = kEmin;
+ // check that apparent energy loss in detector is compatible with a & z
+ // if it is too big, we will adjust first a, then z, until it is
+ if(e > GetMaxDeltaE(z,a)){
+ bool ok=kFALSE;
+ //Info("GetCorrectedEnergy", "Looking for solution for de=%f Z=%d A=%d",e,z,a);
+ for(int znew=z;znew<z+10;znew++){
+ //Info("GetCorrectedEnergy", "Try Z=%d",z);
+ if(znew>z) nuc->SetZ(znew);
+ KVNumberList arange = nuc->GetKnownARange();
+ arange.Begin();
+ while(!ok && !arange.End()){
+ int anew = arange.Next();
+ //Info("GetCorrectedEnergy", "Try A=%d?",anew);
+ if(z==znew && anew<=a) continue;
+ //Info("GetCorrectedEnergy", "ok try it",anew);
+ nuc->SetA(a=anew);
+ if(e<=GetMaxDeltaE(znew,a)){
+ // check consistency with transmission status
+ if(transmission||(!transmission && (e<GetPunchThroughEnergy(znew,a)))) ok=kTRUE;
+ }
+ }
+ if(ok) {z=znew; /*Info("GetCorrectedEnergy", "1.Found solution: Z=%d A=%d",z,a);*/break;}
+ }
+ //if(ok) Info("GetCorrectedEnergy", "2.Found solution: Z=%d A=%d",z,a);
+
+ }
Double_t EINC, ERES = GetEResAfterDetector();
if(transmission && ERES>0.){
// if residual energy is known we use it to calculate EINC.
@@ -843,44 +869,10 @@
// the corrected dE of this detector would not depend on the
// measured dE !
}
- EINC = GetIncidentEnergy(z, a, e, solution);//will be <0 if deltaE is too big
- Bool_t einc_neg = kFALSE;
- if(EINC<0.){
- einc_neg=kTRUE;
- /*Info("GetCorrectedEnergy",
- "%s : (%d,%d) dE=%f > max theoretical dE=%f",
- GetName(),z,a,e,GetMaxDeltaE(z,a));*/
- // deltaE is bigger than max theoretical dE for (Z,A)
- // increase Z until we find a solution
- KVNucleus tmpnuc;
- tmpnuc.SetMassFormula(nuc->GetMassFormula());
- while(EINC<0. && z<(nuc->GetZ()+10)){
- tmpnuc.SetZ(++z);
- a = tmpnuc.GetA();
- EINC = GetIncidentEnergy(z,a,e,solution);
- }
- /*if(EINC>0){
- Info("GetCorrectedEnergy",
- "%s : energy loss compatible with (%d,%d), Einc=%f",
- GetName(),z,a,EINC);
- }
- else
- {
- Info("GetCorrectedEnergy",
- "%s : still no solution found even with (%d,%d)",
- GetName(),z,a);
- }*/
- }
-
-// ERES = GetERes(nuc->GetZ(),nuc->GetA(),EINC);
+ EINC = GetIncidentEnergy(z, a, e, solution);
ERES = GetERes(z,a,EINC);
SetEResAfterDetector(-1.);
- //incident energy - residual energy = total real energy loss
- /*if(einc_neg && EINC>0)
- Info("GetCorrectedEnergy",
- "%s : (%d,%d) dE=%f --> (%d,%d) corrected dE=%f",
- GetName(),nuc->GetZ(),nuc->GetA(),e,z,a,EINC-ERES);*/
return (EINC - ERES);
}
@@ -1381,10 +1373,19 @@
// Overrides KVMaterial::GetEIncOfMaxDeltaE
// Returns incident energy corresponding to maximum energy loss in the
// active layer of the detector, for a given nucleus.
-
+
return GetELossFunction(Z,A)->GetMaximumX();
}
+Double_t KVDetector::GetMaxDeltaE(Int_t Z, Int_t A)
+{
+ // Overrides KVMaterial::GetMaxDeltaE
+ // Returns maximum energy loss in the
+ // active layer of the detector, for a given nucleus.
+
+ return GetELossFunction(Z,A)->GetMaximum();
+}
+
Double_t KVDetector::GetDeltaE(Int_t Z, Int_t A, Double_t Einc)
{
// Overrides KVMaterial::GetDeltaE
=== modified file 'KVMultiDet/detectors/KVDetector.h'
--- KVMultiDet/detectors/KVDetector.h 2012-10-17 09:33:13 +0000
+++ KVMultiDet/detectors/KVDetector.h 2013-04-02 10:59:30 +0000
@@ -192,7 +192,7 @@
virtual void SetEnergyLoss(Double_t e) {
SetEnergy(e);
};
- virtual Double_t GetCorrectedEnergy(const KVNucleus*, Double_t e =
+ virtual Double_t GetCorrectedEnergy(KVNucleus *, Double_t e =
-1., Bool_t transmission=kTRUE);
virtual Int_t FindZmin(Double_t ELOSS = -1., Char_t mass_formula = -1);
@@ -324,6 +324,7 @@
// detector implementations: this version just returns -1.
return -1;
};
+ virtual Double_t GetMaxDeltaE(Int_t Z, Int_t A);
virtual Double_t GetEIncOfMaxDeltaE(Int_t Z, Int_t A);
virtual Double_t GetDeltaE(Int_t Z, Int_t A, Double_t Einc);
virtual Double_t GetTotalDeltaE(Int_t Z, Int_t A, Double_t Einc);
@@ -399,7 +400,7 @@
virtual void DeduceACQParameters(Int_t /*zz*/=-1,Int_t /*aa*/=-1){};
- ClassDef(KVDetector, 8) //Base class for the description of detectors in multidetector arrays
+ ClassDef(KVDetector, 8) //Base class for the description of detectors in multidetector arrays
};
inline KVCalibrator *KVDetector::GetCalibrator(const Char_t * name,
=== modified file 'KVMultiDet/events/KVWilckeReactionParameters.h'
--- KVMultiDet/events/KVWilckeReactionParameters.h 2013-03-08 15:00:51 +0000
+++ KVMultiDet/events/KVWilckeReactionParameters.h 2013-04-02 10:59:30 +0000
@@ -104,8 +104,8 @@
// instead of (RCT**3+RCP**3)**1/3
// The Wilcke formula gives negative V0 values and an incorrect form of potential.
// The actual values in Wilcke correspond to the correct Bondorf formula
- Double_t v0 = pow((zp+zt),2)/pow(pow(ChargeRadius_Myers(zt,at),3)+pow(ChargeRadius_Myers(zp,ap),3),THIRD);
- v0 -= (pow(zt,2)/ChargeRadius_Myers(zt,at) + pow(zp,2)/ChargeRadius_Myers(zp,ap));
+ Double_t v0 = pow((zp+zt),2.)/pow(pow(ChargeRadius_Myers(zt,at),3.)+pow(ChargeRadius_Myers(zp,ap),3.),THIRD);
+ v0 -= (pow(zt,2.)/ChargeRadius_Myers(zt,at) + pow(zp,2.)/ChargeRadius_Myers(zp,ap));
v0*=3.*e2_Wilcke/5.;
return v0;
}
Follow ups
