EdbShowAlg_CA Class Reference

#include <EdbShowAlg.h>

Inheritance diagram for EdbShowAlg_CA:

Collaboration diagram for EdbShowAlg_CA:

Public Member Functions
	ClassDef (EdbShowAlg_CA, 1)
	EdbShowAlg_CA ()
void	Execute ()
void	Finalize ()
void	Initialize ()
virtual	~EdbShowAlg_CA ()
Public Member Functions inherited from EdbShowAlg
void	AddRecoShowerArray (EdbShowerP *shower)
	ClassDef (EdbShowAlg, 1)
Double_t	DeltaR_NoPropagation (EdbSegP *s, EdbSegP *stest)
Double_t	DeltaR_WithoutPropagation (EdbSegP *s, EdbSegP *stest)
Double_t	DeltaR_WithPropagation (EdbSegP *s, EdbSegP *stest)
Double_t	DeltaTheta (EdbSegP *s1, EdbSegP *s2)
Double_t	DeltaThetaComponentwise (EdbSegP *s1, EdbSegP *s2)
Double_t	DeltaThetaSingleAngles (EdbSegP *s1, EdbSegP *s2)
	EdbShowAlg ()
	EdbShowAlg (TString AlgName, Int_t AlgValue)
TString	GetAlgName () const
Int_t	GetAlgValue () const
Double_t	GetMinimumDist (EdbSegP *seg1, EdbSegP *seg2)
TObjArray *	GetRecoShowerArray () const
Int_t	GetRecoShowerArrayN () const
EdbShowerP *	GetShower (Int_t i) const
Double_t	GetSpatialDist (EdbSegP *s1, EdbSegP *s2)
void	Help ()
Bool_t	IsInConeTube (EdbSegP *sTest, EdbSegP *sStart, Double_t CylinderRadius, Double_t ConeAngle)
void	Print ()
void	PrintAll ()
void	PrintMore ()
void	PrintParameters ()
void	PrintParametersShort ()
void	SetActualAlgParameterset (Int_t ActualAlgParametersetNr)
void	SetEdbPVRec (EdbPVRec *Ali)
void	SetEdbPVRecPIDNumbers (Int_t FirstPlate_eAliPID, Int_t LastPlate_eAliPID, Int_t MiddlePlate_eAliPID, Int_t NumberPlate_eAliPID)
void	SetInBTArray (TObjArray *InBTArray)
void	SetParameter (Int_t parNr, Float_t par)
void	SetParameters (Float_t *par)
void	SetRecoShowerArray (TObjArray *RecoShowerArray)
void	SetRecoShowerArrayN (Int_t RecoShowerArrayN)
void	SetUseAliSub (Bool_t UseAliSub)
void	Transform_eAli (EdbSegP *InitiatorBT, Float_t ExtractSize)
void	UpdateShowerIDs ()
void	UpdateShowerMetaData ()
virtual	~EdbShowAlg ()

Private Member Functions
Bool_t	FindPrecedingBTs (EdbSegP *s, EdbSegP *InBT, EdbPVRec *gAli, EdbShowerP *shower)
Bool_t	GetConeTubeDistanceToInBT (EdbSegP *sa, EdbSegP *s, Double_t CylinderRadius, Double_t ConeAngle)

Additional Inherited Members
Protected Member Functions inherited from EdbShowAlg
void	Set0 ()
Protected Attributes inherited from EdbShowAlg
Int_t	eActualAlgParametersetNr
TString	eAlgName
Int_t	eAlgValue
EdbPVRec *	eAli
EdbPVRec *	eAli_Sub
Int_t	eAli_SubNpat
Int_t	eAliNpat
Int_t	eFirstPlate_eAliPID
TObjArray *	eInBTArray
Int_t	eInBTArrayN
Int_t	eLastPlate_eAliPID
Int_t	eMiddlePlate_eAliPID
Int_t	eNumberPlate_eAliPID
Int_t	eParaN
TString	eParaString [10]
Float_t	eParaValue [10]
EdbShowerP *	eRecoShower
TObjArray *	eRecoShowerArray
Int_t	eRecoShowerArrayN
Int_t	eUseAliSub

Constructor & Destructor Documentation

EdbShowAlg_CA()

EdbShowAlg_CA::EdbShowAlg_CA

(

)

{
// Default Constructor
    cout << "EdbShowAlg_CA::EdbShowAlg_CA()   Default Constructor"<<endl;
 
// Reset all:
    Set0();
 
    eAlgName="CA";
    eAlgValue=2; // see default.par_SHOWREC for labelling
    eParaN =4; // depends on the algorithm, how much parameters it has got.
 
    eParaValue[0]=700;
    eParaString[0]="CylinderRadius";  // micrometer
    eParaValue[1]=0.025;
    eParaString[1]="ConeAngle";       // radiant
    eParaValue[2]=110;
    eParaString[2]="DeltaRToPreceedingBT";  // micrometer
    eParaValue[3]=0.05;
    eParaString[3]="DeltaThetaToPreceedingBT";  // radiant
 
//   700  0.025  110  0.05
}

~EdbShowAlg_CA()

EdbShowAlg_CA::~EdbShowAlg_CA ( )

virtual

{
    cout << "EdbShowAlg_CA::~EdbShowAlg_CA()"<<endl;
// Default Destructor
}

Member Function Documentation

ClassDef()

EdbShowAlg_CA::ClassDef	(	EdbShowAlg_CA	,
		1
	)

Execute()

void EdbShowAlg_CA::Execute ( )

virtual

Reimplemented from EdbShowAlg.

{
    cout << "EdbShowAlg_CA::Execute()" << endl;
    cout << "EdbShowAlg_CA::Execute   DOING MAIN SHOWER RECONSTRUCTION HERE" << endl;
 
    EdbSegP* InBT;
    EdbSegP* Segment;
    EdbShowerP* RecoShower=0;
 
    Bool_t    StillToLoop=kTRUE;
    Int_t     ActualPID;
    Int_t     newActualPID;
    Int_t     STEP=-1;
    Int_t     NLoopedPattern=0;
 
    if (gEDBDEBUGLEVEL>3) {
        cout << "======================================================"<<endl;
        cout << "  eFirstPlate_eAliPID =  " <<  eFirstPlate_eAliPID  << endl;
        cout << "  eLastPlate_eAliPID =  " <<   eLastPlate_eAliPID << endl;
        cout << "  eNumberPlate_eAliPID =  " <<   eNumberPlate_eAliPID << endl;
        cout << "  eAli->GetPattern(eFirstPlate_eAliPID)->Z()   =  " <<  eAli->GetPattern(eFirstPlate_eAliPID)->Z()   << endl;
        cout << "  eAli->GetPattern(eLastPlate_eAliPID)->Z()   =  " <<  eAli->GetPattern(eLastPlate_eAliPID)->Z()   << endl;
        cout << "  eInBTArrayN =  " <<  eInBTArrayN  << endl;
        cout << "======================================================"<<endl;
    }
 
    if (eFirstPlate_eAliPID-eLastPlate_eAliPID<0) STEP=1;
    if (gEDBDEBUGLEVEL>3) cout << "EdbShowAlg_CA::Execute--- STEP for patternloop direction =  " << STEP << endl;
 
 
//--- Loop over InBTs:
 
// Since eInBTArray is filled in ascending ordering by zpositon
// We use the descending loop to begin with BT with lowest z first.
    for (Int_t i=eInBTArrayN-1; i>=0; --i) {
 
// CounterOutPut
        if (gEDBDEBUGLEVEL==2) if ((i%50)==0) cout << eInBTArrayN <<" InBT in total, still to do:"<<Form("%4d",i)<< "\r\r\r\r"<<flush;
        if (gEDBDEBUGLEVEL>=4) cout << "Doing  i=" << i << endl;
//-----------------------------------
// 1) Make local_gAli with cut parameters:
//-----------------------------------
 
// Create new EdbShowerP Object for storage;
// See EdbShowerP why I have to call the Constructor as "unique" ideficable value
        // cout << "Create new EdbShowerP Object for storage."<<endl;
        // cout << i << endl;
        // cout << eAlgValue << endl;
        // cout << eActualAlgParametersetNr << endl;
        RecoShower = new EdbShowerP(i,eAlgValue,eActualAlgParametersetNr);
        // cout << "Create new EdbShowerP Object for storage;...done."<<endl;
 
// Get InitiatorBT from eInBTArray
        InBT=(EdbSegP*)eInBTArray->At(i);
        if (gEDBDEBUGLEVEL>2) InBT->PrintNice();
        Float_t X0=InBT->X();
        Float_t Y0=InBT->Y();
 
// Add InBT to RecoShower:
// This has to be done, since by definition the first BT in the RecoShower is the InBT.
// Otherwise, also the definition of shower axis and transversal profiles is wrong!
        RecoShower -> AddSegment(InBT);
        if (gEDBDEBUGLEVEL>4)  cout << "Segment  (InBT) " << Segment << " was added to RecoShower." <<  endl;
 
 
// Transform (make size smaller, extract only events having same MC) the  eAli  object:
        Transform_eAli(InBT);
        if (gEDBDEBUGLEVEL>3) eAli_Sub->Print();
 
        if (gEDBDEBUGLEVEL>2) cout << "  eAli_Sub->GetPattern(eFirstPlate_eAliPID)->Z()   =  " <<  eAli_Sub->GetPattern(eFirstPlate_eAliPID)->Z()   << endl;
        if (gEDBDEBUGLEVEL>2) cout << "  eAli_Sub->GetPattern(eLastPlate_eAliPID)->Z()   =  " <<  eAli_Sub->GetPattern(eLastPlate_eAliPID)->Z()   << endl;
 
 
//-----------------------------------
// 2) Loop over (whole) eAli, check BT for Cuts
// eAli_Sub
//-----------------------------------
        ActualPID= InBT->PID() ;
        newActualPID= InBT->PID() ;
 
        while (StillToLoop) {
            if (gEDBDEBUGLEVEL>3) cout << "EdbShowAlg_CA::Execute--- --- Doing patterloop id= " << ActualPID << " for pattern with Z position= " << eAli_Sub->GetPattern(ActualPID)->Z() << endl;
 
            for (Int_t btloop_cnt=0; btloop_cnt<eAli_Sub->GetPattern(ActualPID)->GetN(); ++btloop_cnt) {
 
                Segment = (EdbSegP*)eAli_Sub->GetPattern(ActualPID)->GetSegment(btloop_cnt);
                if (gEDBDEBUGLEVEL>4) Segment->PrintNice();
 
 
// Now apply cut conditions: CA ConeAdvanced Alg  --------------------
                if ( Segment->MCEvt() > 0 ) if ( Segment->MCEvt()!=InBT->MCEvt() ) continue; // MCEvtNr (>0) or BgMCNr (-999)
                if ( Abs(Segment->X()-X0) > 7000 ) continue;
                if ( Abs(Segment->Y()-Y0) > 7000 ) continue;
                if ( !GetConeTubeDistanceToInBT(Segment, InBT, eParaValue[0], eParaValue[1]) ) continue;
//  if ( !FindPrecedingBTs(Segment, InBT, eAli, RecoShower)) continue;
                if ( !FindPrecedingBTs(Segment, InBT, eAli_Sub, RecoShower)) continue;
// end of    cut conditions: CA ConeAdvanced Alg  --------------------
 
// If we arrive here, Basetrack  Segment  has passed criteria
// and is then added to the RecoShower:
// Check if its not the InBT which is already added:
                if (Segment->X()==InBT->X()&&Segment->Y()==InBT->Y()) {
                    ;    // is InBT, do nothing;
                }
                else {
                    RecoShower -> AddSegment(Segment);
                }
                if (gEDBDEBUGLEVEL>4)  cout << "Segment  " << Segment << " was added at  &Segment : " << &Segment  <<  endl;
 
 
            } // of btloop_cnt
 
            if (gEDBDEBUGLEVEL>3) cout << "EdbShowAlg_CA::Execute--- --- ActualPID= " << newActualPID << "  done. Reconstructed shower has up to now: " << RecoShower->N()  << " Segments." << endl;
 
//------------
            newActualPID=ActualPID+STEP;
            ++NLoopedPattern;
 
            if (gEDBDEBUGLEVEL>3) cout << "EdbShowAlg_CA::Execute--- --- newActualPID= " << newActualPID << endl;
            if (gEDBDEBUGLEVEL>3) cout << "EdbShowAlg_CA::Execute--- --- NLoopedPattern= " << NLoopedPattern << endl;
            if (gEDBDEBUGLEVEL>3) cout << "EdbShowAlg_CA::Execute--- --- eNumberPlate_eAliPID= " << eNumberPlate_eAliPID << endl;
            if (gEDBDEBUGLEVEL>3) cout << "EdbShowAlg_CA::Execute--- --- StillToLoop= " << StillToLoop << endl;
 
// This if holds in the case of STEP== +1
            if (STEP==1) {
                if (newActualPID>eLastPlate_eAliPID) StillToLoop=kFALSE;
                if (newActualPID>eLastPlate_eAliPID) cout << "EdbShowAlg_CA::Execute--- ---Stop Loop since: newActualPID>eLastPlate_eAliPID"<<endl;
            }
// This if holds in the case of STEP== -1
            if (STEP==-1) {
                if (newActualPID<eLastPlate_eAliPID) StillToLoop=kFALSE;
                if (newActualPID<eLastPlate_eAliPID && gEDBDEBUGLEVEL>3) cout << "EdbShowAlg_CA::Execute--- ---Stop Loop since: newActualPID<eLastPlate_eAliPID"<<endl;
            }
// This if holds  general, since eNumberPlate_eAliPID is not dependent of the structure of the gAli subject:
            if (NLoopedPattern>eNumberPlate_eAliPID) StillToLoop=kFALSE;
            if (NLoopedPattern>eNumberPlate_eAliPID && gEDBDEBUGLEVEL>3) cout << "EdbShowAlg_CA::Execute--- ---Stop Loop since: NLoopedPattern>eNumberPlate_eAliPID"<<endl;
 
            ActualPID=newActualPID;
        } // of // while (StillToLoop)
 
// Obligatory when Shower Reconstruction is finished!
        RecoShower->Update();
 
 
        if (gEDBDEBUGLEVEL>3) RecoShower ->PrintBasics();
        if (gEDBDEBUGLEVEL>3) RecoShower ->PrintNice();
        if (gEDBDEBUGLEVEL>3) RecoShower ->PrintSegments();
 
//     if (gEDBDEBUGLEVEL>3) cout << "EdbShowAlg_CA::Execute--- Before adding to array delete the histograms...by finalize() of shower."<<endl;
//     RecoShower ->Finalize();
//     if (gEDBDEBUGLEVEL>3) cout << "EdbShowAlg_CA::Execute--- Finalize done..."<<endl;
 
 
// Add Shower Object to Shower Reco Array.
// Not, if its empty or containing only one BT:
        if (RecoShower->N()>1) {
      AddRecoShowerArray(RecoShower);
    }
    else {
      // Delete the created shower object. Saves _lots_ of memory and speeds up reconstruction for large #InBTs!
      // TODO Also Include this in the other RecoAlgs ...
      Log(3,"EdbShowAlg_CA::Execute()","InBT # %d RecoShower->N()<=1: Dont add shower to RecoShowerArray. Delete this shower for memory safing.");
      delete RecoShower;
      RecoShower=NULL;
    }
 
//cout << " eRecoShowerArray->GetEntries()  " <<  eRecoShowerArray->GetEntries()   << endl;
 
 
// Set back loop values:
        StillToLoop=kTRUE;
        NLoopedPattern=0;
    } // of  //   for (Int_t i=eInBTArrayN-1; i>=0; --i) {
 
 
// Set new value for  eRecoShowerArrayN  (may now be < eInBTArrayN).
    SetRecoShowerArrayN(eRecoShowerArray->GetEntries());
 
 
    cout << "EdbShowAlg_CA::eRecoShowerArray() Entries: " << eRecoShowerArray->GetEntries() << endl;
    cout << "EdbShowAlg_CA::Execute()...done." << endl;
    return;
}

Finalize()

void EdbShowAlg_CA::Finalize ( )

virtual

Reimplemented from EdbShowAlg.

{
    return;
}

FindPrecedingBTs()

Bool_t EdbShowAlg_CA::FindPrecedingBTs ( EdbSegP *  s, EdbSegP *  InBT, EdbPVRec *  gAli, EdbShowerP *  shower  )

private

{
    EdbSegP* s_TestBT;
//       Int_t nentries=shower->GetEntries();
    Int_t nentries=shower->N();
    Double_t dZ;
 
// In case the shower is empty we do not have to search for a preceeding BT:
    if (nentries==0) return kTRUE;
 
// We do not test BaseTracks which are before the initiator BaseTrack!
    if (s->Z()<InBT->Z()) {
        cout << "--- --- EdbShowAlg_CA::FindPrecedingBTs(): s->Z()<InBT->Z()..: We do not test BaseTracks which are before the initiator BaseTrack!  return kFALSE;" << endl;
        return kFALSE;
    }
 
// For the very first Z position we do not test
// if testBT has Preceeders, only if it it has a BT around (case for e+e- coming from gammma):
// Take 50microns and 80mrad in (dR/dT) around.
// This does not affect the normal results, but helps for
// events which may have a second BT close to InBT (like in e+e-)
    if (s->Z()==InBT->Z()) {
//cout << DeltaTheta(s, InBT) << "  " << DeltaR_WithPropagation(s, InBT) << endl;
        if (DeltaTheta(s, InBT) < 0.08 && DeltaR_WithoutPropagation(s, InBT) < 50.0 ) return kTRUE;
    }
 
    if (gEDBDEBUGLEVEL>3) cout << "--- --- EdbShowAlg_CA::FindPrecedingBTs(): Testing " << nentries << " entries of the shower." << endl;
 
    for (Int_t i=nentries-1; i>=0; --i) {
        s_TestBT = (EdbSegP*)( shower->GetSegment(i) );
 
        if (gEDBDEBUGLEVEL>3) cout << "--- --- EdbShowAlg_CA::FindPrecedingBTs(): Do  s_TestBT->ID()  s->ID()  s_TestBT->MCEvt() s_TestBT->Z()  s->Z() "<< s_TestBT->ID() << " " << s->ID() << " " << s_TestBT->MCEvt() <<"  " << s_TestBT->Z() << "  " <<  s->Z()  <<endl;
        if (gEDBDEBUGLEVEL>3) s_TestBT->PrintNice();
        if (gEDBDEBUGLEVEL>3) s->PrintNice();
 
 
        dZ=TMath::Abs(s_TestBT->Z()-s->Z());
        if (dZ<30) continue;                  // Exclude the case of same Zpositions...
        if (dZ>(3.5*(1300.0+30.0))) continue;     // Exclude the case of more than 3 plates before...  // modified
        // often, experimentally observed, dZ approx 1350 ... so the upper cut really does not
        // look for 3 plates, but only for two ...: seems to work now. (24.Oct.2016)
 
        if (gEDBDEBUGLEVEL>3) cout << "--- --- EdbShowAlg_CA::FindPrecedingBTs(): Checking dT,dR and dZ for i:  " << i << "  " << DeltaTheta(s, s_TestBT)  << "  " << DeltaR_WithPropagation(s, s_TestBT) << "  "<<dZ << endl;
 
        if (DeltaTheta(s, s_TestBT) >  eParaValue[3] ) continue;
        if (DeltaR_WithPropagation(s, s_TestBT) > eParaValue[2]) continue;
 
        if (gEDBDEBUGLEVEL>3) cout << "--- --- EdbShowAlg_CA::FindPrecedingBTs(): Checking dT,dR and dZ for i:  " << i << "  " << DeltaTheta(s, s_TestBT)  << "  " << DeltaR_WithPropagation(s, s_TestBT) << "  "<<dZ << "   ok!"<<endl;
        return kTRUE;
    }
//---------------------------------------------
    return kFALSE;
}

GetConeTubeDistanceToInBT()

Bool_t EdbShowAlg_CA::GetConeTubeDistanceToInBT ( EdbSegP *  sa, EdbSegP *  s, Double_t  CylinderRadius, Double_t  ConeAngle  )

private

Outside if angle greater than ConeAngle (to be fulfilled for Cone and Tube in both cases)

if angle smaller than ConeAngle, then you can differ between Tuberadius and CylinderRadius

{
// This is an updated version since the TMath::ACos(cosangle) does not really refer to a coneangle:
// We rather use from now on the ROOT implementation of TVector3::Angle()
 
    TVector3 x1(InBT->X(),InBT->Y(),InBT->Z());
    TVector3 x2(sa->X(),sa->Y(),sa->Z());
    TVector3 direction_x1(InBT->TX()*1300,InBT->TY()*1300,1300);
    TVector3 direction_x2(sa->TX()*1300,sa->TY()*1300,1300);
    TVector3 u1=x2-x1;
 
    Double_t direction_x1_norm= direction_x1.Mag();
    Double_t cosangle=  (direction_x1*u1)/(u1.Mag()*direction_x1_norm);
    Double_t angle = TMath::ACos(cosangle);  // NO THIS IS NOT THE CONE ANGLE!! (this is rather an old relict)
 
    TVector3 direction_1(InBT->TX()*1300,InBT->TY()*1300,1300);
    TVector3 direction_2(sa->TX()*1300,sa->TY()*1300,1300);
 
    /*
    // DEBUG STUFF ... TO BE REMOVED IF EVERYTHING IS CLEAR
    cout <<"------"<<endl;
    cout << "angle =  " << angle << endl;
    cout << "angle V2 =  " << x1.Angle(x2) << endl;
    cout << "angle V3 =  " << direction_x1.Angle(direction_x2) << endl;
    cout << "angle V4 =  " << direction_1.Angle(direction_2) << "   " <<  direction_1.Angle(direction_2)/3.14*180.0 << endl;
 
    cout << "angle V5   direction_x1.Angle(u1  =  " << direction_x1.Angle(u1) << endl;
    cout << "angle V5   u1.Angle(direction_x1  =  " << u1.Angle(direction_x1) << endl;
    */
    angle=u1.Angle(direction_x1);
 
// For the case where the two basetracks have same z position
// the angle is about 90 degree so it makes no sense to calculate it...
// therefore we set it artificially to zero:
    if (InBT->Z()==sa->Z() ) {
        angle=0.0;
//if (gEDBDEBUGLEVEL>3) //cout << "same z position, set angle artificially to zero" << endl;
    }
 
    if (angle>ConeAngle) {
        return kFALSE;
    }
 
    Double_t TubeDistance = 1.0/direction_x1_norm  *  ( (x2-x1).Cross(direction_x1) ).Mag();
//   cout << "CylinderRadius= " <<  CylinderRadius << "  TubeDistance= " << TubeDistance << endl;
    if (TubeDistance>CylinderRadius) {
        return kFALSE;
    }
 
    return kTRUE;
}

Initialize()

void EdbShowAlg_CA::Initialize ( )

virtual

Reimplemented from EdbShowAlg.

{
    return;
}

---

The documentation for this class was generated from the following files:

• src/libShowRec/EdbShowAlg.h
• src/libShowRec/EdbShowAlg.cxx