EdbShowerAlg_GS Class Reference

#include <EdbShowerAlg.h>

Inheritance diagram for EdbShowerAlg_GS:

Collaboration diagram for EdbShowerAlg_GS:

Public Member Functions
void	AddInVtx (EdbVertex *vtx)
Double_t	CalcIP (EdbSegP *s, double x, double y, double z)
Double_t	CalcIP (EdbSegP *s, EdbVertex *v)
TObjArray *	CheckCleanPairs (EdbSegP *InBT, TObjArray *RecoShowerArrayFromFindPairs)
Bool_t	CheckInput ()
Bool_t	CheckPairDuplications (Int_t SegPID, Int_t SegID, Int_t Seg2PID, Int_t Seg2ID, TArrayI *SegmentPIDArray, TArrayI *SegmentIDArray, TArrayI *Segment2PIDArray, TArrayI *Segment2IDArray, Int_t RecoShowerArrayN)
	ClassDef (EdbShowerAlg_GS, 1)
void	Convert_InVtxArray_To_InBTArray ()
void	CreateANNPair ()
void	CreateANNPlots ()
	EdbShowerAlg_GS ()
void	Execute ()
void	FillANNTree (TObjArray *RecoShowerArray, EdbSegP *Inbt)
void	Finalize ()
TObjArray *	FindPairs (EdbSegP *InBT, EdbPVRec *eAli_Sub)
TObjArray *	FindPairsPreselected (EdbPVRec *eAli_Sub)
Bool_t	FindPairsPreselected (EdbSegP *InitiatorBT, EdbPVRec *eAli_Sub)
TObjArray *	GetInVtxArray () const
Int_t	GetInVtxArrayN () const
void	Init ()
void	Initialize ()
Bool_t	IsPossibleFakeDoublet (EdbSegP *s1, EdbSegP *s2)
void	ReloadANNs (Int_t RecoMode)
TObjArray *	SelectHighestPInMCArray (TObjArray *BTArray)
void	Set0 ()
void	SetCleanPairs (Bool_t CleanPairs)
void	SetCutModeFull (Bool_t CutModeFull)
void	SetFindPairsPreselected (Bool_t FindPairsPreselected)
void	SetInVtx (EdbVertex *vtx)
void	SetInVtxArray (TObjArray *InVtxArray)
void	SetRecoMode (Int_t RecoMode)
virtual	~EdbShowerAlg_GS ()
Public Member Functions inherited from EdbShowerAlg
void	AddRecoShowerArray (TObjArray *RecoShowerArray)
EdbVertex *	CalcVertex (TObjArray *segments)
	ClassDef (EdbShowerAlg, 1)
void	Convert_EdbPVRec_To_InBTArray ()
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)
	EdbShowerAlg ()
	EdbShowerAlg (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
Double_t	GetSpatialDist (EdbSegP *s1, EdbSegP *s2)
Bool_t	IsInConeTube (EdbSegP *sTest, EdbSegP *sStart, Double_t CylinderRadius, Double_t ConeAngle)
void	Print ()
void	PrintAll ()
void	PrintMore ()
void	PrintParameters ()
void	PrintParametersShort ()
void	PrintRecoShowerArray ()
void	SetActualAlgParameterset (Int_t ActualAlgParametersetNr)
void	SetDebug ()
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 parvalue)
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)
virtual	~EdbShowerAlg ()

Private Attributes
TMultiLayerPerceptron *	eANNPair
TMultiLayerPerceptron *	eANNPairCaseA
TMultiLayerPerceptron *	eANNPairCaseB
TMultiLayerPerceptron *	eANNPairCaseC
Float_t	eANNPairCut [3]
TTree *	eANNPairTree
Bool_t	eCutModeFull
Bool_t	eFindPairsPreselected
TObjArray *	eInVtxArray
Int_t	eInVtxArrayN
Bool_t	eInVtxArraySet
Int_t	eRecoMode
TArrayI *	eSegment2IDArray
TArrayI *	eSegment2PIDArray
TArrayI *	eSegmentIDArray
TArrayI *	eSegmentPIDArray
Bool_t	eSetCleanPairs
Float_t	eValueGSNN_var00
Float_t	eValueGSNN_var01
Float_t	eValueGSNN_var02
Float_t	eValueGSNN_var03
Float_t	eValueGSNN_var04
Float_t	eValueGSNN_var05
Float_t	eValueGSNN_var06
Float_t	eValueGSNN_varInput
Float_t	eValueGSNN_varOutput

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

Constructor & Destructor Documentation

EdbShowerAlg_GS()

EdbShowerAlg_GS::EdbShowerAlg_GS

(

)

{
    // Default Constructor
    //
    // This Shower Algorithm finds all compatible basetrack pairings that could
    // come from an ePlus-eMinus pair.
    // It needs an EdbPVRec* object.
    // Additionally -but not necessary- an array (or one) vertex, w.r.t. which the
    // basetracks should be looked at.
    //   If vertex array is not given, then the basetracks of the volume will be
    //   used to create interim helper vertices.
    // Usage:
    //   Alg = new EdbShowerAlg_GS();
    //   Alg->SetEdbPVRec(EdbPVRec* pvr);
    //   Alg->SetInVtx(EdbVertex* vtx);  // additionally...
    //   Alg->Execute();
    //   Alg->GetRecoShowerArray();
    //
    // Returns the compatible Pair Segments (stored as EdbTrackP*).
 
    Log(2,"EdbShowerAlg_GS::EdbShowerAlg_GS","EdbShowerAlg_GS:: Default Constructor");
 
    cout << "Howto do:  " << endl;
    cout << "   Alg = new EdbShowerAlg_GS();" << endl;
    cout << "   Alg->SetEdbPVRec(EdbPVRec* pvr);" << endl;
    cout << "   Alg->SetInVtx(EdbVertex* vtx)  // additionally... " << endl;
    cout << "   Alg->Execute();" << endl;
    cout << "   Alg->GetRecoShowerArray();   // .... Returns you the compatible Pair Segments (stored as EdbTrackP*)" << endl;
    cout << "   " << endl;
 
 
    //  Init with values according to GS Alg:
    //  Create essential objects:
    Init();
 
    // Reset all internal variables.
    Set0();
 
    // Create Neural Network for better fake BT rejection
    CreateANNPair();
}

~EdbShowerAlg_GS()

EdbShowerAlg_GS::~EdbShowerAlg_GS ( )

virtual

{
    // Default Destructor
    Log(2,"EdbShowerAlg_GS::~EdbShowerAlg_GS","EdbShowerAlg_GS:: Default Destructor");
}

Member Function Documentation

AddInVtx()

void EdbShowerAlg_GS::AddInVtx

(

EdbVertex *

vtx

)

{
    Log(3,"EdbShowerAlg_GS::AddInVtx","AddInVtx()");
    eInVtxArray->Add(vtx);
    ++eInVtxArrayN;
    eInVtxArraySet=kTRUE;
    if (gEDBDEBUGLEVEL>2) {
        cout << "EdbShowerAlg_GS::AddInVtx(): Added vtx (XYZ,MC): " <<  vtx->X() << " " <<vtx->Y() << " " << vtx->Z() << " " << vtx->MCEvt() << " " << endl;
        cout << "EdbShowerAlg_GS::AddInVtx(): Added One Vertex. Now there are " << eInVtxArrayN << " InVtx stored." << endl;
    }
    Log(3,"EdbShowerAlg_GS::AddInVtx","AddInVtx()...done");
    return;
}

CalcIP() [1/2]

double EdbShowerAlg_GS::CalcIP	(	EdbSegP *	s,
		double	x,
		double	y,
		double	z
	)

                                                                       {
    // Calculate IP between a given segment and a given x,y,z.
    // return the IP value.
    double ax = s->TX();
    double ay = s->TY();
    double bx = s->X()-ax*s->Z();
    double by = s->Y()-ay*s->Z();
    double a;
    double r;
    double xx,yy,zz;
    a = (ax*(x-bx)+ay*(y-by)+1.*(z-0.))/(ax*ax+ay*ay+1.);
    xx = bx +ax*a;
    yy = by +ay*a;
    zz = 0. +1.*a;
    r = sqrt((xx-x)*(xx-x)+(yy-y)*(yy-y)+(zz-z)*(zz-z));
    return r;
}

CalcIP() [2/2]

double EdbShowerAlg_GS::CalcIP	(	EdbSegP *	s,
		EdbVertex *	v
	)

                                                       {
    // calculate IP between a given segment and a given vertex.
    // return the IP value.
    // this is used for IP cut.
 
    // if vertex is not given, use the selected vertex.
    // if(v==NULL) v=gEDA->GetSelectedVertex();
 
    if (v==NULL) return -1.;
    return CalcIP(s, v->X(), v->Y(), v->Z());
}

CheckCleanPairs()

TObjArray * EdbShowerAlg_GS::CheckCleanPairs	(	EdbSegP *	InBT,
		TObjArray *	RecoShowerArrayFromFindPairs
	)

{
    // Check Clean Pairs Function.
    // It looks for basetrack pairs that come in a row directly behind each other.
    // Then it will take the one most upstream in this line.
    // More precise to explain:: WHAT DOES THIS FUNCTION DO ???
 
    Log(3,"EdbShowerAlg_GS::CheckCleanPairs","Starting CheckCleanPairs() now...");
 
    if (NULL==RecoShowerArray) return NULL;
 
    TObjArray* NewRecoShowerArray= new TObjArray();
    NewRecoShowerArray->Clear();
    int NewRecoShowerArrayN=0;
 
    EdbTrackP* TrackPair1=NULL;
    EdbTrackP* TrackPair2=NULL;
 
    int ntrack=RecoShowerArray->GetEntries();
    // cout << "EdbShowerAlg_GS::CheckCleanPairs   ntrack = "  << ntrack << endl;
 
    for (Int_t pat_one_cnt=0; pat_one_cnt<ntrack; ++pat_one_cnt) {
        if (gEDBDEBUGLEVEL>2)  cout << "CheckCleanPairs   Doing pat_one_cnt = " << pat_one_cnt << endl;
        // (if -1) then the last one is not checked
        TrackPair1=(EdbTrackP*)RecoShowerArray->At(pat_one_cnt);
        bool taketrack1=true;
        for (Int_t pat_two_cnt=pat_one_cnt; pat_two_cnt<ntrack; ++pat_two_cnt) {
            if (gEDBDEBUGLEVEL>3) cout << "CheckCleanPairs      Doing pat_two_cnt = " << pat_two_cnt << endl;
            //if only one track at all take it anyway:
            if (ntrack==1) continue;
            TrackPair2=(EdbTrackP*)RecoShowerArray->At(pat_two_cnt);
 
            // Check  if track1 has a track before (smaller Z value)
            // and if so, if dtheta is smaller than 0.1:
            // If both is so, then track 1 is NOT taken for final array:
            //       cout << TrackPair1->Z() << "  " <<  TrackPair2->Z()  << endl;
            EdbSegP* s1=(EdbSegP* )TrackPair1->GetSegment(0);
            EdbSegP* s2=(EdbSegP* )TrackPair2->GetSegment(0);
            if (TrackPair1->Z()>TrackPair2->Z() && DeltaThetaSingleAngles((EdbSegP*)s1,(EdbSegP*)s2)<0.1) taketrack1=false;
            if (!taketrack1) break;
        }
 
        if (!taketrack1) continue;
 
        // Add TrackPair1
        NewRecoShowerArray->Add(TrackPair1);
        ++NewRecoShowerArrayN;
 
        // Print
        if (TrackPair1->N()<2) cout << "CheckCleanPairs   This Track has only ONE entry" << endl;
 
    }
    if (gEDBDEBUGLEVEL>2) cout << "CheckCleanPairs   From " << ntrack << "  originally, there are now after Zposition and overlap cuts: " << NewRecoShowerArrayN << " left." << endl;
 
    Log(3,"EdbShowerAlg_GS::CheckCleanPairs","Starting CheckCleanPairs() now...done.");
    return NewRecoShowerArray;
}

CheckInput()

Bool_t EdbShowerAlg_GS::CheckInput

(

)

{
    // Before execution of the main reco routine, we perform a
    // check, if either an Initiator Vertex (array) was set,
    // or if an Initiator Basetrack (array) was set.
    // In this case, return kTRUE, else return kFALSE.
 
    Log(2,"EdbShowerAlg_GS::CheckInput","CheckInput()");
 
    Bool_t IsInput=kFALSE;
    Bool_t VtxArray_Or_InBTArray=kFALSE;
    Log(2,"EdbShowerAlg_GS::CheckInput","CheckInput()   Check for eInBTArray:");
 
 
    if (eInBTArrayN==0) {
        Log(2,"EdbShowerAlg_GS::CheckInput","CheckInput()   No eInBTArray. Check for eInVtxArray:");
 
        // Check if there is an Convert_InVtxArray_To_InBTArray();
        if (eInVtxArrayN==0) {
            Log(2,"EdbShowerAlg_GS::CheckInput","CheckInput()   No eInVtxArray.");
            Log(2,"EdbShowerAlg_GS::CheckInput","CheckInput()   Set ALL BTs from the volume as InBTs:");
            Log(2,"EdbShowerAlg_GS::CheckInput","CheckInput()   Set BTs from the volume as InBTs:");
            if (eFindPairsPreselected==kTRUE) {
                Log(2,"EdbShowerAlg_GS::CheckInput","CheckInput()   Take Preselected Pairs found from");
                Log(2,"EdbShowerAlg_GS::CheckInput","CheckInput()   the function FindPairsPreselected():");
            }
            VtxArray_Or_InBTArray=kFALSE;
            // Convert_EdbPVRec_To_InBTArray();
            // Taking all BTs takes very long time, we rather skip it.
            //______________________________________________________________________________
            // New we try finding already preselcted InBTs.... might speed up the style...
            if (eFindPairsPreselected==kTRUE) {
                SetInBTArray(this->FindPairsPreselected(eAli));
            }
            else {
                Convert_EdbPVRec_To_InBTArray();
            }
            eInBTArrayN=eInBTArray->GetEntries();
            Log(2,"EdbShowerAlg_GS::CheckInput","CheckInput()   eInBTArray->GetEntries() = %d .",eInBTArray->GetEntries());
            //______________________________________________________________________________
            Log(2,"EdbShowerAlg_GS::CheckInput","CheckInput()   FindPairsPreselected() done.");
            //eRecoMode=1;
//             eRecoMode=2;
            Log(2,"EdbShowerAlg_GS::CheckInput","CheckInput()   Use default eRecoMode=%d",eRecoMode);
            Log(2,"EdbShowerAlg_GS::CheckInput","CheckInput()   Do pair search now, but take care, this is gonna take  _______VERY LONG_______ time.");
        }
        else {
            Log(2,"EdbShowerAlg_GS::CheckInput","CheckInput()   eInVtxArray there. Converting to InBTArray() now...");
            Convert_InVtxArray_To_InBTArray();
            VtxArray_Or_InBTArray=kTRUE;
            eRecoMode=0;
            Log(2,"EdbShowerAlg_GS::CheckInput","CheckInput()   eRecoMode=0");
            Log(2,"EdbShowerAlg_GS::CheckInput","CheckInput()   eInVtxArray there. Converting to InBTArray() now....done.");
        }
    }
 
    cout << "EdbShowerAlg_GS::CheckInput: Finally take eRecoMode = " << eRecoMode << endl;
 
    // Reload ANN weight file:
    // If vertices are given we will take the trained neural net for basetrack pairs
    // with the known vertex.
    // If no vertices are given then we have to assume that we dont know
    // what the vertex position yet is and then we look for basetrack pairings
    // by themselves. Therefore we need the ANN weight file trained on no vertex info.
 
    ReloadANNs(eRecoMode);
 
 
    // For the MC Training, we use only highest P basetracks, if they come from signal
    TObjArray* arr = SelectHighestPInMCArray(eInBTArray);
//      SetInBTArray(arr);
    // For the MC Training, we use only highest P basetracks, if they come from signal
 
 
    Log(2,"EdbShowerAlg_GS::CheckInput","CheckInput()...done.");
    return IsInput;
}

CheckPairDuplications()

Bool_t EdbShowerAlg_GS::CheckPairDuplications	(	Int_t	SegPID,
		Int_t	SegID,
		Int_t	Seg2PID,
		Int_t	Seg2ID,
		TArrayI *	SegmentPIDArray,
		TArrayI *	SegmentIDArray,
		TArrayI *	Segment2PIDArray,
		TArrayI *	Segment2IDArray,
		Int_t	RecoShowerArrayN
	)

{
    Log(3,"EdbShowerAlg_GS::CheckPairDuplications","CheckPairDuplications()");
    if (gEDBDEBUGLEVEL>2) cout << "EdbShowerAlg_GS::CheckPairDuplications for " << SegPID << "," <<  SegID <<","<< Seg2PID <<"," << Seg2ID << " compare with " << RecoShowerArrayN << " pairs:" <<  endl;
    for (Int_t i=0; i<RecoShowerArrayN; i++) {
        // PID and ID of Seg and Seg2 to be exchanged for duplications
        if ( SegPID==eSegment2PIDArray->At(i) && Seg2PID==eSegmentPIDArray->At(i) && SegID==eSegment2IDArray->At(i) && Seg2ID==eSegmentIDArray->At(i)) {
            if (gEDBDEBUGLEVEL>2) cout << "EdbShowerAlg_GS::CheckPairDuplications Found duplication for ... return true"<<endl;
            return kTRUE;
        }
    }
    return kFALSE;
}

ClassDef()

EdbShowerAlg_GS::ClassDef	(	EdbShowerAlg_GS	,
		1
	)

Convert_InVtxArray_To_InBTArray()

void EdbShowerAlg_GS::Convert_InVtxArray_To_InBTArray

(

)

{
    // This function takes vertices from the InVtxArray object,
    // and fills the InBTArray with those: X/Y/Z position as
    // from the Vertices, TX(),TY() equal to 0
    //   (this doesnt matter when having InVtxArray with vertices filled
    //    for this algo, since it does rely only on the point
    //    position here).
    // This function is useful, if one wants to start from
    // already calculated vertices.
    //
    // If vtx-array is not filled, on the other hand, it is recommended
    // to use the full  EdbPVRec object as possible starting points,
    // look in EdbShowerAlg::Convert_EdbPVRec_To_InBTArray().
 
    Log(2,"EdbShowerAlg_GS::EdbShowerAlg_GS","Convert_InVtxArray_To_InBTArray()");
    if (eInBTArrayN!=0) {
        eInBTArray->Clear();
        cout << " eInBTArray->Clear();" << endl;
    }
 
    if (eInVtxArray==NULL || eInVtxArrayN==0 || eInVtxArraySet==kFALSE) cout << "NO  eInVtxArray  " << endl;
    EdbVertex* vtx;
 
    Log(2,"EdbShowerAlg_GS::EdbShowerAlg_GS","Convert_InVtxArray_To_InBTArray()...start loop (N=%d):",eInVtxArrayN);
    for (Int_t i=0; i<eInVtxArrayN; i++) {
        vtx= (EdbVertex*)eInVtxArray->At(i);
        EdbSegP* seg = new EdbSegP(i,vtx->X(),vtx->Y(),0,0);
        seg->SetZ(vtx->Z());
        // vtx can have by default initialization MCEvt==0,
        // this we dont want, in case: we set MCEvt of vtx from 0 to -999
        if (vtx->MCEvt()==0) vtx->SetMC(-999);
        seg->SetMC(vtx->MCEvt(),vtx->MCEvt());
        seg->SetFlag(0);
        eInBTArray->Add(seg);
 
        //cout <<"DEBUG   vtx->MCEvt() " << vtx->MCEvt() << endl;
        // seg->PrintNice();
    }
 
    eInBTArrayN=eInBTArray->GetEntries();
    if (gEDBDEBUGLEVEL>2) cout << "EdbShowerAlg_GS::Convert_InVtxArray_To_InBTArray   Converted " << eInBTArrayN << "InVtx to InBT." << endl;
    Log(2,"EdbShowerAlg_GS::EdbShowerAlg_GS","Convert_InVtxArray_To_InBTArray()...done.");
    return;
}

CreateANNPair()

void EdbShowerAlg_GS::CreateANNPair

(

)

{
    // Create the neural network, which is used as last layer in the
    // GS algo reconstruction.
    // Load its weights from the $FEDRASYS path.
    // Initial testings of the method have been done using the ROOT
    // TMVA framework, espacially the multivariate method "Boosted Decision Trees"
    // and the two given Neural Network Implementations and linear cuts.
    // ROC curves can be found in the manual. The two NN implementation dont
    // differ too much, and we can take -with good  conscience- the standard
    // ROOT TMultiLayerPerceptron one.
 
    Log(2,"EdbShowerAlg_GS::EdbShowerAlg_GS","CreateANNPair()");
 
    // Just to avoid the root errormessages of creating a tree non memory resident.
    TFile* localANNPairTreeFile = new TFile("localANNPairTreeFile.root","RECREATE");
 
    // Create the tree which will store the Inputvariables for the GS ANN discrimination:
    eANNPairTree= new TTree("eANNPairTree","eANNPairTree");
    cout << "EdbShowerAlg_GS::CreateANNPair()   eANNPairTree  Memory creation done." << endl;
    eANNPairTree->Branch("eValueGSNN_varInput",&eValueGSNN_varInput,"eValueGSNN_varInput/F");
    eANNPairTree->Branch("eValueGSNN_varOutput",&eValueGSNN_varOutput,"eValueGSNN_varOutput/F");
    eANNPairTree->Branch("eValueGSNN_var00",&eValueGSNN_var00,"eValueGSNN_var00/F");
    eANNPairTree->Branch("eValueGSNN_var01",&eValueGSNN_var01,"eValueGSNN_var01/F");
    eANNPairTree->Branch("eValueGSNN_var02",&eValueGSNN_var02,"eValueGSNN_var02/F");
    eANNPairTree->Branch("eValueGSNN_var03",&eValueGSNN_var03,"eValueGSNN_var03/F");
    eANNPairTree->Branch("eValueGSNN_var04",&eValueGSNN_var04,"eValueGSNN_var04/F");
    eANNPairTree->Branch("eValueGSNN_var05",&eValueGSNN_var05,"eValueGSNN_var05/F");
    eANNPairTree->Branch("eValueGSNN_var06",&eValueGSNN_var06,"eValueGSNN_var06/F");
    cout << "EdbShowerAlg_GS::CreateANNPair()   eANNPairTree  SetBranchAddress done." << endl << endl;
 
    // :@eValueGSNN_var06 is not used, because Flag() is a MC information only.
    TString layout="";
    TString weights="";
 
    // Attention: the three cases have a different number of input neurons:
    // Case A and Case C have 6 input neurons
    // Case B has only 5 input neurons (because dZ is a constant and we cannot use ANN
    //                                  with constant input neurons).
    //
    // CASE A : Vertex Info out of "RealVertex"
    // CASE B : No Vertex Info, but "FakeVertex"    made out of one BT,  dZ==650
    // CASE C : No Vertex Info, but "InterimVertex" made out of two BTs, dZ variable
 
    // CASE A :
    layout="@eValueGSNN_var00,@eValueGSNN_var01,@eValueGSNN_var02,@eValueGSNN_var03,@eValueGSNN_var04,@eValueGSNN_var05:7:6:eValueGSNN_varInput";
    // Create the Multilayerperceptron
    eANNPairCaseA  = new TMultiLayerPerceptron(layout,eANNPairTree,"","");
    cout << "EdbShowerAlg_GS::CreateANNPair()   eANNPairCaseA  TMultiLayerPerceptron creation done." << endl;
    cout << "EdbShowerAlg_GS::CreateANNPair()   Layout = " << endl;
    cout << layout.Data() << endl;
    weights= TString(gSystem->ExpandPathName("$FEDRA_ROOT"))+TString("/src/libShower/weights/Reco/ShowerAlg_GS/weights_CASE_A.txt");
    eANNPairCaseA->LoadWeights(weights);
    cout << "EdbShowerAlg_GS::CreateANNPair()   eANNPair: weights = " << endl;
    cout << weights.Data() << endl << endl;
 
    // CASE B :
    layout="@eValueGSNN_var01,@eValueGSNN_var02,@eValueGSNN_var04,@eValueGSNN_var05:7:6:eValueGSNN_varInput";
    // Create the Multilayerperceptron
    eANNPairCaseB  = new TMultiLayerPerceptron(layout,eANNPairTree,"","");
    cout << "EdbShowerAlg_GS::CreateANNPair()   eANNPairCaseB  TMultiLayerPerceptron creation done." << endl;
    cout << "EdbShowerAlg_GS::CreateANNPair()   Layout = " << endl;
    cout << layout.Data() << endl;
    weights= TString(gSystem->ExpandPathName("$FEDRA_ROOT"))+TString("/src/libShower/weights/Reco/ShowerAlg_GS/weights_CASE_B.txt");
    eANNPairCaseB->LoadWeights(weights);
    cout << "EdbShowerAlg_GS::CreateANNPair()   eANNPair: weights = " << endl;
    cout << weights.Data() << endl << endl;
 
    // CASE C :
    layout="@eValueGSNN_var00,@eValueGSNN_var01,@eValueGSNN_var02,@eValueGSNN_var03,@eValueGSNN_var04,@eValueGSNN_var05:7:6:eValueGSNN_varInput";
    // Create the Multilayerperceptron
    eANNPairCaseC  = new TMultiLayerPerceptron(layout,eANNPairTree,"","");
    cout << "EdbShowerAlg_GS::CreateANNPair()   eANNPairCaseC  TMultiLayerPerceptron creation done." << endl;
    cout << "EdbShowerAlg_GS::CreateANNPair()   Layout = " << endl;
    cout << layout.Data() << endl;
    weights= TString(gSystem->ExpandPathName("$FEDRA_ROOT"))+TString("/src/libShower/weights/Reco/ShowerAlg_GS/weights_CASE_C.txt");
    eANNPairCaseC->LoadWeights(weights);
    cout << "EdbShowerAlg_GS::CreateANNPair()   eANNPair: weights = " << endl;
    cout << weights.Data() << endl << endl;
 
    // By default we take the weights for the Case C (explained in the manual)
    // since we dont know a priori, if vertices are given.
    eANNPair=eANNPairCaseC;
 
    Log(2,"EdbShowerAlg_GS::EdbShowerAlg_GS","CreateANNPair()...done.");
    return;
}

CreateANNPlots()

void EdbShowerAlg_GS::CreateANNPlots

(

)

{
    // Just  a plot creating function for the GS_ANN Input variables.
 
    Log(3,"EdbShowerAlg_GS::CreateANNPlots","CreateANNPlots...");
 
    gROOT->SetStyle("Bold");
    // Draw Histograms for the last cut on Gamma Search
    TCanvas* canv_input = new TCanvas();
    canv_input->Divide(3,2);
    canv_input->cd(1);
    eANNPairTree->Draw("eValueGSNN_var00");
    canv_input->cd(2);
    eANNPairTree->Draw("eValueGSNN_var01");
    canv_input->cd(3);
    eANNPairTree->Draw("eValueGSNN_var02");
    canv_input->cd(4);
    eANNPairTree->Draw("eValueGSNN_var03");
    canv_input->cd(5);
    eANNPairTree->Draw("eValueGSNN_var04");
    canv_input->cd(6);
    eANNPairTree->Draw("eValueGSNN_var05");
 
    TCanvas* canv_InOutput = new TCanvas();
    canv_InOutput->Divide(3,1);
    canv_InOutput->cd(1);
    eANNPairTree->Draw("eValueGSNN_varInput");
    canv_InOutput->cd(2);
    eANNPairTree->Draw("eValueGSNN_varOutput");
    canv_InOutput->cd(3);
    eANNPairTree->Draw("eValueGSNN_varInput:eValueGSNN_varOutput","","colz");
 
    Log(3,"EdbShowerAlg_GS::CreateANNPlots","CreateANNPlots...done.");
    return;
}

Execute()

void EdbShowerAlg_GS::Execute ( )

virtual

Debug Purpose:

Reimplemented from EdbShowerAlg.

{
    // Main shower reco function.
    // Loops over all filled Initiator basetracks.
    // Fills the TObjArray with reconstructed showers.
    // GS Reco algo is in function FindPairs().
 
    Log(2,"EdbShowerAlg_GS::Execute","Execute()");
    Log(2,"EdbShowerAlg_GS::Execute","Execute()   DOING MAIN SHOWER RECONSTRUCTION HERE");
    Log(2,"EdbShowerAlg_GS::Execute","Execute()   Check for eInBTArray:");
 
    // Check Input which BT filling should be done.
    Bool_t VtxArray_Or_InBTArray = kFALSE;
    VtxArray_Or_InBTArray = CheckInput();
 
 
    //--- Needed interim objects:
    EdbSegP* InBT=NULL;
//   EdbShowerP* RecoShower; //libShowRec
    EdbTrackP* RecoShower;   //libShower
 
    Int_t PermilleCount=eInBTArrayN/1000;
    if (PermilleCount==0) PermilleCount=1; // floating point exception failsafe.
 
    Int_t     STEP=-1;
    if (eFirstPlate_eAliPID-eLastPlate_eAliPID<0) STEP=1;
    if (gEDBDEBUGLEVEL>2) cout << "EdbShowerAlg_GS::Execute--- STEP for patternloop direction =  " << STEP << endl;
    if (gEDBDEBUGLEVEL>2) cout << "EdbShowerAlg_GS::Execute--- eRecoShowerArrayN =  " << eRecoShowerArrayN << endl;
 
    //--- Loop over InBTs:
    if (gEDBDEBUGLEVEL>=2) {
        cout << "EdbShowerAlg_GS::Execute    Loop over InBTs: " << endl;
        cout << "EdbShowerAlg_GS::Execute    (N=" << eInBTArrayN << "):" << endl;
        cout << "EdbShowerAlg_GS::Execute    (each dot means a progress of 1 permille of N):" << endl;
        cout << "EdbShowerAlg_GS::Execute    (each   K means a progress of 1 K InBTs):" << endl;
    }
 
    //-----------------------------------
    // Since eInBTArray is filled in ascending ordering by zpositon
    // We use the descending loop to begin with BT with lowest z first.
    // InBT doest have to be necessary a real BaseTrack, it can also be a vertex (with its positions and angle zero) !!!
    for (Int_t i=eInBTArrayN-1; i>=0; --i) {
 
        if (eDebug) if (i<eInBTArrayN-1) continue;
 
        //-----------------------------------
        // 0) Set counters...
        //-----------------------------------
        Log(3,"EdbShowerAlg_GS::Execute","//-----------------------------------");
        Log(3,"EdbShowerAlg_GS::Execute","// 0) Set counters...:");
        Log(3,"EdbShowerAlg_GS::Execute","//-----------------------------------");
 
        // CounterOutPut
//         if (gEDBDEBUGLEVEL==2) if ((i%100)==0) cout << eInBTArrayN <<" InBT in total, still to do:"<<Form("%4d",i)<< "\r\r\r\r"<<flush;
        if (gEDBDEBUGLEVEL==2)  if ((i%PermilleCount)==0) cout << "." << flush;
        if (gEDBDEBUGLEVEL==2)  if ((i%1000)==0) cout << "K" << flush;
 
        // Get InitiatorBT from eInBTArray  InBT
        InBT=(EdbSegP*)eInBTArray->At(i);
        if (gEDBDEBUGLEVEL>2) {
            Log(3,"EdbShowerAlg_GS::Execute","//---- Print Initiator BT(Vtx)   ------------------");
            InBT->PrintNice();
        }
 
        //-----------------------------------
        // 1) Make local_gAli with cut parameters:
        //-----------------------------------
        Log(3,"EdbShowerAlg_GS::Execute","//-----------------------------------");
        Log(3,"EdbShowerAlg_GS::Execute","// 1) Make local_gAli with cut parameters:");
        Log(3,"EdbShowerAlg_GS::Execute","//-----------------------------------");
        // Transform (make size smaller, extract only events having same MC) the  eAli  object:
        Transform_eAli(InBT,999999);
        if (gEDBDEBUGLEVEL>2)  eAli_Sub->Print();
 
 
        //-----------------------------------
        // 2) FindPairs
        //-----------------------------------
        Log(3,"EdbShowerAlg_GS::Execute","//-----------------------------------");
        Log(3,"EdbShowerAlg_GS::Execute","// 2) FindPairs");
        Log(3,"EdbShowerAlg_GS::Execute","//-----------------------------------");
        TObjArray* Pairs = FindPairs(InBT,eAli_Sub);
        if (gEDBDEBUGLEVEL>2) cout << "TObjArray* Pairs = FindPairs(InBT,eAli_Sub); done. Entries= " << Pairs->GetEntries() << endl;
 
 
 
        //-----------------------------------
        // 3) Clear Found Pairs
        //-----------------------------------
        Log(3,"EdbShowerAlg_GS::Execute","//-----------------------------------");
        Log(3,"EdbShowerAlg_GS::Execute","// 3) Clear Found Pairs");
        Log(3,"EdbShowerAlg_GS::Execute","//-----------------------------------");
        TObjArray* CleanPairs;
        if (eSetCleanPairs==kTRUE) {
            CleanPairs = CheckCleanPairs( InBT, Pairs );
            // After we found pairs (after   clear) we add the internal RecoShowerArray:
            AddRecoShowerArray(CleanPairs);
            FillANNTree(CleanPairs,InBT);
        }
        else {
            // After we found pairs (without clear) we add the internal RecoShowerArray:
            AddRecoShowerArray(Pairs);
            FillANNTree(Pairs,InBT);
        }
 
        //-----------------------------------
        // 4) Finally Apply GS_ANN cut on Pairs
        // ---> Is already done in the FindPairs() routine!!
        //-----------------------------------
 
 
    } // Loop over InBT
    //-----------------------------------
 
    cout << "EdbShowerAlg_GS::Execute  Loop over InBT finished."<< endl;
    cout << "EdbShowerAlg_GS::Execute  eRecoShowerArray=  " <<  eRecoShowerArray << endl;
    cout << "EdbShowerAlg_GS::Execute  eRecoShowerArrayN=  " <<  eRecoShowerArrayN << endl;
 
 
    if (gEDBDEBUGLEVEL>2) {
        cout << "EdbShowerAlg_GS::Execute  Print all BT pairs now:" << endl;
        for (Int_t i=0; i <eRecoShowerArray->GetEntries();  i++) {
//   EdbShowerP* sh = (EdbShowerP* )eRecoShowerArray->At(i);  // libShowRec
            EdbTrackP* sh = (EdbTrackP* )eRecoShowerArray->At(i); // libShower
            sh->PrintNice();
        }
    }
 
    if (eRecoShowerArrayN>0) {
        EdbTrackP* sh = (EdbTrackP* )eRecoShowerArray->At(0); // libShower
        cout << "EdbShowerAlg_GS::Execute  Shower    0: PrintNice: " << endl;
        sh->PrintNice();
        sh = (EdbTrackP* )eRecoShowerArray->At(eRecoShowerArray->GetEntries()-1); // libShower
        cout << "EdbShowerAlg_GS::Execute  Shower last("<< eRecoShowerArrayN <<"): PrintNice: " << endl;
        sh->PrintNice();
    }
 
    eANNPair->Write();
    eANNPairTree->Write();
 
    Log(2,"EdbShowerAlg_GS::Execute","Execute()...done.");
    return;
}

FillANNTree()

void EdbShowerAlg_GS::FillANNTree	(	TObjArray *	RecoShowerArray,
		EdbSegP *	Inbt
	)

only here we distinguish, since the function DeltaR_WithPropagation is sensitive to the Z-ordering which basetrack comes first.

{
    // This Function will be used to fill the ANN tree,
    // when the routine is called _after_  FindPairs.
    // Value calulation is the same as in the FindPairs function.
 
    Log(3,"EdbShowerAlg_GS::FillANNTree","FillANNTree...");
    Log(3,"EdbShowerAlg_GS::FillANNTree","RecoShowerArray->GetEntries() = %d", RecoShowerArray->GetEntries());
 
    if (RecoShowerArray->GetEntries()==0) return;
 
    for (Int_t i=0; i< RecoShowerArray->GetEntries(); ++i) {
 
        EdbTrackP*  shower = (EdbTrackP*)RecoShowerArray->At(i);
        EdbSegP*  Segment = (EdbSegP*)shower->GetSegment(0);
        EdbSegP*  Segment2 = (EdbSegP*)shower->GetSegment(1);
 
        EdbVertex* vtx=new EdbVertex();
        if (eRecoMode==1 || eRecoMode==2 ) {
            vtx->SetXYZ(InBT->X()-650*InBT->TX(),InBT->Y()-650*InBT->TY(),InBT->Z()-650);
        }
        else {
            vtx->SetXYZ(InBT->X(),InBT->Y(),InBT->Z());
        }
 
        TObjArray *segments = new TObjArray(2);
        segments->Add(Segment);
        segments->Add(Segment2);
        EdbVertex* vetex = CalcVertex(segments);
        delete segments;
        Float_t IP_Seg1_To_Vtx = CalcIP(Segment, vtx);
        Float_t IP_Seg2_To_Vtx = CalcIP(Segment2,vtx);
        Float_t IP_InBT_To_Vtx = CalcIP(InBT,vtx);
        Float_t IP_Seg1ToVtxSeg1Seg2=CalcIP(Segment ,vetex);
        Float_t IP_Seg2ToVtxSeg1Seg2=CalcIP(Segment2,vetex);
 
        eValueGSNN_varInput=-1;
        // Purity 1:    Input =  1.0;
        // Purity 0.5:  Input =  0.5;
        // Purity else: Input =  0.0;
        if (Segment2->Flag()+Segment->Flag()==0&&TMath::Abs(Segment2->Flag())==11&&Segment->MCEvt()>0) {
            eValueGSNN_varInput=1;
        }
        else if (Segment2->Flag()+Segment->Flag()!=0&&TMath::Abs(Segment2->Flag())==11&&Segment->MCEvt()>0) {
            eValueGSNN_varInput=0.5;
        }
        else {
            eValueGSNN_varInput=0;
        }
 
        // Set Tree Input Variable Values:
        // Calculation of these values is EXACTLY the same as in the
        // FindPairs() function.
        eValueGSNN_var00=TMath::Min(IP_Seg1_To_Vtx,IP_Seg2_To_Vtx);
        if (eRecoMode==1) eValueGSNN_var00=IP_InBT_To_Vtx;
        if (eRecoMode==2) {
            eValueGSNN_var00=TMath::Min(IP_Seg1ToVtxSeg1Seg2,IP_Seg2ToVtxSeg1Seg2);
        }
        eValueGSNN_var01=GetMinimumDist(Segment,Segment2);
        eValueGSNN_var02=DeltaR_WithPropagation(Segment,Segment2);
        if (Segment2->Z()>Segment->Z()) eValueGSNN_var02=DeltaR_WithPropagation(Segment2,Segment);
        eValueGSNN_var03=InBT->Z()-vtx->Z();
        if (eRecoMode==2) {
            eValueGSNN_var03=TMath::Min(Segment->Z(),Segment2->Z())-vetex->Z();
        }
        if (eRecoMode==0) {
            eValueGSNN_var03=TMath::Min(Segment->Z(),Segment2->Z())-vtx->Z();
        }
        eValueGSNN_var04=DeltaThetaSingleAngles(Segment,Segment2);
        eValueGSNN_var05=TMath::Abs(Segment->PID()-Segment2->PID());
        eValueGSNN_var06=Segment2->Flag()+Segment->Flag();
 
        if (gEDBDEBUGLEVEL>2) {
            cout <<"-------------   FillANNTree Print eValueGSNN   first, pring seg1 and seg2------------- " << endl;
            Segment->PrintNice();
            Segment2->PrintNice();
            cout <<"-------------   FillANNTree Print eValueGSNN_var00= " << eValueGSNN_var00 << "   ---- " << endl;
            cout <<"-------------   FillANNTree Print eValueGSNN_var01= " << eValueGSNN_var01 << "   ---- " << endl;
            cout <<"-------------   FillANNTree Print eValueGSNN_var02= " << eValueGSNN_var02 << "   ---- " << endl;
            cout <<"-------------   FillANNTree Print eValueGSNN_var03= " << eValueGSNN_var03 << "   ---- " << endl;
            cout <<"-------------   FillANNTree Print eValueGSNN_var04= " << eValueGSNN_var04 << "   ---- " << endl;
            cout <<"-------------   FillANNTree Print eValueGSNN_var05= " << eValueGSNN_var05 << "   ---- " << endl;
        }
 
        eANNPairTree->Fill();
 
    } // for (Int_t i=0; i< RecoShowerArray->GetEntries(); ++i) {
 
    Log(3,"EdbShowerAlg_GS::FillANNTree","FillANNTree...done.");
    return;
}

Finalize()

void EdbShowerAlg_GS::Finalize ( )

virtual

Reimplemented from EdbShowerAlg.

{
    // do nothing yet.
}

FindPairs()

TObjArray * EdbShowerAlg_GS::FindPairs	(	EdbSegP *	InBT,
		EdbPVRec *	eAli_Sub
	)

if (eRecoMode==1 && distZ>10) continue; in this reco mode we loop only over segment2 , since InBT and vtx are correlated!!!

{
    // Function to find pairs of an EdbPVRec object.
    // with respect to a Initiator Basetrack (either a real vertex
    // or a helper vertex from the BT itself).
    //
    // ATTENTION: CURRENTLY (23.07.2011) STILL IN MODIFYING PHASE!!!
    // ATTENTION: CURRENTLY (19.09.2011) STILL IN MODIFYING PHASE!!!
    // ATTENTION: CURRENTLY (01.10.2011) STILL IN MODIFYING PHASE!!!
    // ATTENTION: CURRENTLY (05.10.2011) STILL IN MODIFYING PHASE!!!
    // ATTENTION: CURRENTLY (11.10.2011) STILL IN MODIFYING PHASE!!!
 
    Log(3,"EdbShowerAlg_GS::FindPairs","Starting FindPairs(InBT,eAli_Sub) now....");
 
    TObjArray* RecoShowerArray= new TObjArray();
    RecoShowerArray->Clear();
    Int_t RecoShowerArrayN=0;
 
    eSegmentPIDArray->Reset();
    eSegmentIDArray->Reset();
    eSegment2PIDArray->Reset();
    eSegment2IDArray->Reset();
 
    EdbSegP* Segment=NULL;
    EdbSegP* Segment2=NULL;
    Float_t x_av,y_av,z_av,tx_av,ty_av,distZ;
    EdbSegP* Segment_Sum=new EdbSegP(0,0,0,0,0,0);
    Float_t  IP_Pair_To_InBT=0;
    Float_t IP_Pair_To_InBT_Seg=0;
    Float_t IP_Pair_To_InBT_Seg2=0;
    Float_t IP_Pair_To_InBT_SegSum=0;
    Float_t IP_Pair_To_InBT_SegSmaller=0;
    Float_t  IP_InBT_To_Vtx=0;
    Float_t  IP_Seg1_To_Vtx=0;
    Float_t  IP_Seg2_To_Vtx=0;
 
    Int_t NPairTriesTotal=0;
    Double_t paramsAC[6]; // for Case A and C
    Double_t paramsB[4]; // for Case B
 
 
    EdbSegP* InBT=NULL;
    if (NULL==InitiatorBT) {
        InBT= new EdbSegP();
        EdbPattern* pat = eAli_Sub->GetPatternZLowestHighest(1);
        InBT->SetX(pat->X());
        InBT->SetY(pat->Y());
        InBT->SetZ(pat->Z());
        InBT->SetTX(0);
        InBT->SetTY(0);
        InBT->SetMC(-999,-999);
        cout << "WARNING   EdbShowerAlg_GS::FindPairs   InBT==NULL. Created a dummy InBT:" << endl;
        InBT->PrintNice();
    }
    else {
        InBT=InitiatorBT;
    }
 
    // Create Helper Vertex:
    EdbVertex* vtx=new EdbVertex();
    vtx->SetMC(InBT->MCEvt());
 
    // Now here distinguish the Case A and Case B for the gamma reco:
    // Case A(0): vertex is the InBT.
    // Case B(1): vertex is made out of the InBT, propagated back half a plate.
    // Case C(2): like case B, but for the IP cutvariable, the IP to new vertex, made
    //            out of BT 1 & 2 is used.
    if (eRecoMode==1 || eRecoMode==2 ) {
        if (gEDBDEBUGLEVEL>2) cout << "EdbShowerAlg_GS::FindPairs   eRecoMode==1: Propagate InBT half a plate back for the vertex." << endl;
        vtx->SetXYZ(InBT->X()-650*InBT->TX(),InBT->Y()-650*InBT->TY(),InBT->Z()-650);
    }
    else {
        if (gEDBDEBUGLEVEL>2) cout << "EdbShowerAlg_GS::FindPairs   eRecoMode==0: Take InBT position for the vertex." << endl;
        //cout << "DEBUG     EdbShowerAlg_GS::FindPairs   eRecoMode==0: Take InBT position for the vertex." << endl;
        vtx->SetXYZ(InBT->X(),InBT->Y(),InBT->Z());
    }
 
    IP_InBT_To_Vtx= CalcIP(InBT,vtx);
    //cout << " CalcIP(InBT,vtx) " <<  CalcIP(InBT,vtx)  << endl;
 
    //-----------------------------------
    // 2) Loop over (whole) eAli, check BT for Cuts
    // eAli_Sub
    // Loop structure:
    // Loop over plates [0..NPatterns-1] for pattern one
    // Loop over plates [0..NPatterns-1] for pattern two
    // Take only plate pairings for |PID diff|<=3
    // Loop over all BT of pattern one
    // Loop over all BT of pattern two
 
    // This doesnt yet distinguish the FIRSTPLATE, MIDDLEPLATE, LATPLATE, NUMBERPLATES
    // labelling, this will be built in later....
    //-----------------------------------
 
    Int_t npat=eAli_Sub->Npatterns();
    Int_t pat_one_bt_cnt_max,pat_two_bt_cnt_max=0;
    EdbPattern* pat_one=0;
    EdbPattern* pat_two=0;
 
    //-----------------------------------
    // ---   Loop over first patterns
    for (Int_t pat_one_cnt=0; pat_one_cnt<npat; ++pat_one_cnt) {
        pat_one=(EdbPattern*)eAli_Sub->GetPattern(pat_one_cnt);
        pat_one_bt_cnt_max=eAli_Sub->GetPattern(pat_one_cnt)->GetN();
 
        // Check if dist Z to vtx (BT) is ok:
        distZ=pat_one->Z()-InBT->Z();
        if (gEDBDEBUGLEVEL>3) {
            cout << "EdbShowerAlg_GS::FindPairs  Check if dist Z to vtx (BT) is ok:  distZ=" << distZ << endl;
        }
        if (distZ>eParaValue[3]) continue;
        if (distZ<-1000) continue;
 
 
        if (gEDBDEBUGLEVEL>2) cout << "EdbShowerAlg_GS::FindPairs   pat_one_cnt=" << pat_one_cnt << "  pat_one->Z() = " << pat_one->Z() << " pat_one_bt_cnt_max= "<< pat_one_bt_cnt_max <<endl;
 
        //-----------------------------------
        // ---   Loop over second patterns
        for (Int_t pat_two_cnt=0; pat_two_cnt<npat; ++pat_two_cnt) {
 
            // Now apply cut conditions: GS  GAMMA SEARCH Alg:
            if (TMath::Abs(pat_one_cnt-pat_two_cnt)>eParaValue[5]) continue;
 
            pat_two=(EdbPattern*)eAli_Sub->GetPattern(pat_two_cnt);
            pat_two_bt_cnt_max=eAli_Sub->GetPattern(pat_two_cnt)->GetN();
 
            distZ=pat_two->Z()-InBT->Z();
            if (distZ>eParaValue[3]) continue;
            if (distZ<-1000) continue;
 
            if (gEDBDEBUGLEVEL>2) cout << "EdbShowerAlg_GS::FindPairs   pat_two_cnt=" << pat_two_cnt << "  pat_two->Z() = " << pat_two->Z() << " pat_two_bt_cnt_max= "<< pat_two_bt_cnt_max <<endl;
 
            //-----------------------------------
            // ---   Loop over basetracks of first pattern
            for (Int_t pat_one_bt_cnt=0; pat_one_bt_cnt<pat_one_bt_cnt_max; ++pat_one_bt_cnt) {
                Segment =  (EdbSegP*)pat_one->GetSegment(pat_one_bt_cnt);
 
                IP_Pair_To_InBT_Seg =CalcIP(Segment,vtx);
                IP_Pair_To_InBT=IP_Pair_To_InBT_Seg;
                IP_Seg1_To_Vtx = CalcIP(Segment,vtx);
 
                // Check if IP to vtx (BT) is ok:
                if (IP_Pair_To_InBT>eParaValue[0]) continue;
 
                //-----------------------------------
                // ---   Loop over basetracks of second pattern
                for (Int_t pat_two_bt_cnt=0; pat_two_bt_cnt<pat_two_bt_cnt_max; ++pat_two_bt_cnt) {
                    Segment2 = (EdbSegP*)pat_two->GetSegment(pat_two_bt_cnt);
 
                    // For statistics: increase the number of total searched combinations:
                    ++NPairTriesTotal;
 
                    if (Segment2==Segment) continue;
                    if (Segment2->ID()==Segment->ID()&&Segment2->PID()==Segment->PID()) continue;
 
                    // Now apply cut conditions: GS  GAMMA SEARCH Alg  --------------------
                    // if InBT is flagged as MC InBT, take care that only BG or same MC basetracks are taken:
                    if (InBT->MCEvt()>0) if (Segment->MCEvt()>0&&Segment2->MCEvt()>0) if (Segment->MCEvt()!=Segment2->MCEvt()) continue;
                    if (InBT->MCEvt()>0) if (Segment->MCEvt()>0&&Segment2->MCEvt()>0) if (Segment->MCEvt()!=InBT->MCEvt()) continue;
                    if (InBT->MCEvt()>0) if (Segment->MCEvt()>0&&Segment2->MCEvt()>0) if (Segment2->MCEvt()!=InBT->MCEvt()) continue;
 
                    // In case of two MC events, check for e+ e- pairs
                    // Only if Parameter is set to choose different Flag() pairs:
                    if (InBT->MCEvt()>0 && eParaValue[6]==1)
                        if (Segment->MCEvt()>0&&Segment2->MCEvt()>0)
                            if ((Segment2->Flag()+Segment->Flag())!=0) continue;
 
                    // a) Check dR between tracks:
                    if (DeltaR_WithPropagation(Segment,Segment2)>eParaValue[2]) continue;
                    // b) Check dT between tracks:
                    if (DeltaThetaSingleAngles(Segment,Segment2)>eParaValue[4]) continue;
                    // c) Check dMinDist between tracks:
                    if (GetMinimumDist(Segment,Segment2)>eParaValue[1]) continue;
 
                    // d) Check if dist Z to vtx (BT) is ok:
                    distZ=Segment->Z()-InBT->Z();
                    if (distZ>eParaValue[3]) continue;
 
                    // At first:  Check for already duplicated pairings:
                    if (CheckPairDuplications(Segment->PID(),Segment->ID(),Segment2->PID(),Segment2->ID(), eSegmentPIDArray,eSegmentIDArray,eSegment2PIDArray,eSegment2IDArray,RecoShowerArrayN)) continue;
 
                    // Check if both basetracks have a vertex which is upstream
                    // of both tracks (only then the two BT are really pointing).
                    TObjArray *segments = new TObjArray(2);
                    segments->Add(Segment);
                    segments->Add(Segment2);
                    EdbVertex* vetex = CalcVertex(segments);
                    if (vetex ->Z()> TMath::Min(Segment->Z(),Segment2->Z()) ) continue;
 
                    // Set Sum values of the to BTs:
                    x_av=Segment2->X()+(Segment->X()-Segment2->X())/2.0;
                    y_av=Segment2->Y()+(Segment->Y()-Segment2->Y())/2.0;
                    z_av=Segment2->Z()+(Segment->Z()-Segment2->Z())/2.0;
                    tx_av=Segment2->TX()+(Segment->TX()-Segment2->TX())/2.0;
                    ty_av=Segment2->TY()+(Segment->TY()-Segment2->TY())/2.0;
                    Segment_Sum->SetX(x_av);
                    Segment_Sum->SetY(y_av);
                    Segment_Sum->SetTX(tx_av);
                    Segment_Sum->SetTY(ty_av);
                    Segment_Sum->SetZ(z_av);
 
                    // Check if IP to vtx (BT) is ok:
                    IP_Pair_To_InBT_Seg2    = CalcIP(Segment2,vtx);
                    if (IP_Pair_To_InBT_Seg2>eParaValue[0]) continue;
 
                    IP_Seg2_To_Vtx = CalcIP(Segment2,vtx);
                    IP_Pair_To_InBT_SegSum=CalcIP(Segment_Sum, vtx);
 
                    Float_t IP_Seg1ToVtxSeg1Seg2=0;
                    Float_t IP_Seg2ToVtxSeg1Seg2=0;
                    IP_Seg1ToVtxSeg1Seg2 = CalcIP(Segment ,vetex);
                    IP_Seg2ToVtxSeg1Seg2 = CalcIP(Segment2,vetex);
 
                    if (gEDBDEBUGLEVEL>2) {
                        cout << "EdbShowerAlg_GS::FindPairs  IP_Pair_To_InBT_Seg   = " << IP_Pair_To_InBT_Seg << endl;
                        cout << "EdbShowerAlg_GS::FindPairs  IP_Pair_To_InBT_Seg2  = " << IP_Pair_To_InBT_Seg2 << endl;
                        cout << "EdbShowerAlg_GS::FindPairs  IP_Pair_To_InBT_SegSum= " << IP_Pair_To_InBT_SegSum << endl;
                        cout << "EdbShowerAlg_GS::FindPairs  IP_Seg1ToVtxSeg1Seg2= " << IP_Seg1ToVtxSeg1Seg2 << endl;
                        cout << "EdbShowerAlg_GS::FindPairs  IP_Seg2ToVtxSeg1Seg2= " << IP_Seg2ToVtxSeg1Seg2 << endl;
                    }
 
                    // Save the segment which has smaller IP, this will be the first BT in the RecoShower
                    if ( IP_Pair_To_InBT_Seg>IP_Pair_To_InBT_Seg2 ) {
                        IP_Pair_To_InBT_SegSmaller=0; // take Segment first
                    }
                    else {
                        IP_Pair_To_InBT_SegSmaller=1; // take Segment2 first
                    }
 
                    // e) Check if this is not a possible fake doublet which is
                    //    sometimes caused by view overlap in the scanning:
                    //    Deprecated, because this is now incorporated in the cleaning routine.
                    //
                    // f) new: Last check if ePair is from e+e-Pair or from fake BG;
                    //    do a NN calculation estimate:
                    //cout << "DEBUGTEST: PRINTNICE SEGMENT AND SEGMENT2 AND VTX" << endl;
                    //Segment->PrintNice();
                    //Segment2->PrintNice();
                    //cout << "Vtx:  Print MC XYZ: " << vtx ->MCEvt()  << " " << vtx ->X()  << " " << vtx ->Y()  << " " << vtx ->Z() <<endl;
 
                    eValueGSNN_varInput=-1;
                    // Purity 1:    Input =  1.0;
                    // Purity 0.5:  Input =  0.5;
                    // Purity else: Input =  0.0;
                    if (Segment2->Flag()+Segment->Flag()==0&&TMath::Abs(Segment2->Flag())==11&&Segment->MCEvt()>0) {
                        eValueGSNN_varInput=1;
                    }
                    else if (Segment2->Flag()+Segment->Flag()!=0&&TMath::Abs(Segment2->Flag())==11&&Segment->MCEvt()>0) {
                        eValueGSNN_varInput=0.5;
                    }
                    else {
                        eValueGSNN_varInput=0;
                    }
 
                    eValueGSNN_var00=TMath::Min(IP_Seg1_To_Vtx,IP_Seg2_To_Vtx);
                    if (eRecoMode==1) eValueGSNN_var00=IP_InBT_To_Vtx;
                    if (eRecoMode==2) {
                        eValueGSNN_var00=TMath::Min(IP_Seg1ToVtxSeg1Seg2,IP_Seg2ToVtxSeg1Seg2);
                    }
 
                    // Set Tree Input Variable Values:
                    eValueGSNN_var01=GetMinimumDist(Segment,Segment2);
                    eValueGSNN_var02=DeltaR_WithPropagation(Segment,Segment2);
                    eValueGSNN_var03=InBT->Z()-vtx->Z();
                    if (eRecoMode==2) {
                        eValueGSNN_var03=TMath::Min(Segment->Z(),Segment2->Z())-vetex->Z();
                    }
                    if (eRecoMode==0) {
                        eValueGSNN_var03=TMath::Min(Segment->Z(),Segment2->Z())-vtx->Z();
                    }
                    eValueGSNN_var04=DeltaThetaSingleAngles(Segment,Segment2);
                    eValueGSNN_var05=TMath::Abs(pat_one_cnt-pat_two_cnt);
                    eValueGSNN_var06=Segment2->Flag()+Segment->Flag();
 
                    if (gEDBDEBUGLEVEL>2) {
                        cout <<"-------------   Print eValueGSNN   ------------- " << endl;
                        Segment->PrintNice();
                        Segment2->PrintNice();
                        cout <<"-------------   Print eValueGSNN_var00= " << eValueGSNN_var00 << "   ---- " << endl;
                        cout <<"-------------   Print eValueGSNN_var01= " << eValueGSNN_var01 << "   ---- " << endl;
                        cout <<"-------------   Print eValueGSNN_var02= " << eValueGSNN_var02 << "   ---- " << endl;
                        cout <<"-------------   Print eValueGSNN_var03= " << eValueGSNN_var03 << "   ---- " << endl;
                        cout <<"-------------   Print eValueGSNN_var04= " << eValueGSNN_var04 << "   ---- " << endl;
                        cout <<"-------------   Print eValueGSNN_var05= " << eValueGSNN_var05 << "   ---- " << endl;
                    }
 
                    // Severe Warning! If in this part of the code:
                    // eANNPairTree->Show(eANNPairTree->GetEntries()-1) is invoked,
                    // then the Variables of the Treeentry are loaded into the memory
                    // address and the current value is overwritten!!!!!!
                    // DO NOT DO THIS!
 
                    // Now we have to fill params values.
                    // Distinguish CASE A,C and CASE B, because of one less input neuron.
                    if (eRecoMode==0||eRecoMode==2) {
                        paramsAC[0]=eValueGSNN_var00; //dIP
                        paramsAC[1]=eValueGSNN_var01; //dMinDist
                        paramsAC[2]=eValueGSNN_var02; //dR
                        paramsAC[3]=eValueGSNN_var03; //deltaZ
                        paramsAC[4]=eValueGSNN_var04; //dT
                        paramsAC[5]=eValueGSNN_var05; //dNPL
                        //paramsAC[6]=eValueGSNN_var06; // cannot be used, uses MC info...
                    }
                    else {
                        paramsB[0]=eValueGSNN_var00; //dIP
                        paramsB[0]=eValueGSNN_var01; //dMinDist
                        paramsB[1]=eValueGSNN_var02; //dR
                        //paramsB[3]=eValueGSNN_var03;
                        //deltaZ // constant, so take out (otherwise ANN doesnt work!)
                        paramsB[2]=eValueGSNN_var04; //dT
                        paramsB[3]=eValueGSNN_var05; //dNPL
                    }
 
                    //---------
                    Double_t value=0;
                    if (gEDBDEBUGLEVEL>2) {
                        int nparams=6;
                        if (eRecoMode==1) {
                            nparams=4;
                            cout << "EdbShowerAlg_GS::FindPairs  ParamsB:  ";
                            for (int hj=0; hj<nparams; hj++) cout << "  " << paramsB[hj];
                            cout << "  - Input: " << eValueGSNN_varInput << endl;
                            cout << "EdbShowerAlg_GS::FindPairs  value=eANNPair->Evaluate(0,params) = ";
                            value=eANNPair->Evaluate(0,paramsB);
                            eValueGSNN_varOutput=value;
                            cout << value << endl;
                            cout << "EdbShowerAlg_GS::FindPairs  --------  " << endl;
                        }
                        else {
                            nparams=6;
                            cout << "EdbShowerAlg_GS::FindPairs  ParamsAC:  ";
                            for (int hj=0; hj<nparams; hj++) cout << "  " << paramsAC[hj];
                            cout << "  - Input: " << eValueGSNN_varInput << endl;
                            cout << "EdbShowerAlg_GS::FindPairs  value=eANNPair->Evaluate(0,params) = ";
                            value=eANNPair->Evaluate(0,paramsAC);
                            eValueGSNN_varOutput=value;
                            cout << value << endl;
                            cout << "EdbShowerAlg_GS::FindPairs  --------  " << endl;
                        }
                    }
 
 
                    // Evaluate now the NN.
                    if (eRecoMode==1) {
                        value=eANNPair->Evaluate(0,paramsB);
                    }
                    else {
                        value=eANNPair->Evaluate(0,paramsAC);
                    }
 
 
                    eValueGSNN_varOutput=value;
                    // Fill the Tree with NN variables later. Then it can be done
                    // also for the cleaned pairs.
                    // eANNPairTree->Fill();
 
                    // After Full Tree is Filled, check wheter last cut condition is satisfied.
                    if ( eCutModeFull ) if (value<eANNPairCut[eRecoMode]) continue;
                    if (value<eANNPairCut[eRecoMode]) continue;
 
 
//                     if (gEDBDEBUGLEVEL>2) {
                    if (value<eANNPairCut[eRecoMode]) {
                        cout << "EdbShowerAlg_GS::FindPairs  This Pair has passed all cuts w.r.t to InBT:" << endl;
                        cout << "EdbShowerAlg_GS::FindPairs  IP_Pair_To_InBT  = " << IP_Pair_To_InBT << endl;
                        cout << "EdbShowerAlg_GS::FindPairs  GetMinimumDist(Segment,Segment2)  = " << GetMinimumDist(Segment,Segment2) << endl;
                        cout << "EdbShowerAlg_GS::FindPairs  CalcIP(Segment_Sum,InBT)  = " << IP_Pair_To_InBT << endl;
                        cout << "EdbShowerAlg_GS::FindPairs  IP_Pair_To_InBT_Seg   = " << IP_Pair_To_InBT_Seg << endl;
                        cout << "EdbShowerAlg_GS::FindPairs  IP_Pair_To_InBT_Seg2  = " << IP_Pair_To_InBT_Seg2 << endl;
                        cout << "EdbShowerAlg_GS::FindPairs  GetSpatialDist(Segment_Sum,InBT)  = " << GetSpatialDist(Segment_Sum,InBT) << endl;
                        cout << "EdbShowerAlg_GS::FindPairs  IP_InBT_To_Vtx = " << IP_InBT_To_Vtx << endl;
                        cout << "EdbShowerAlg_GS::FindPairs  IP_Seg1_To_Vtx = " << IP_Seg1_To_Vtx << endl;
                        cout << "EdbShowerAlg_GS::FindPairs  IP_Seg2_To_Vtx = " << IP_Seg2_To_Vtx << endl;
                        cout << "InBT: Print InBT: " << InBT ->X()  << " " << InBT ->Y()  << " " << InBT ->Z() <<endl;
                        cout << "Vtx:  Print MC XYZ: " << vtx ->MCEvt()  << " " << vtx ->X()  << " " << vtx ->Y()  << " " << vtx ->Z() <<endl;
                    }
 
                    eSegmentPIDArray->AddAt(Segment->PID(), RecoShowerArrayN);
                    eSegmentIDArray->AddAt(Segment->ID(), RecoShowerArrayN);
                    eSegment2PIDArray->AddAt(Segment2->PID(), RecoShowerArrayN);
                    eSegment2IDArray->AddAt(Segment2->ID(), RecoShowerArrayN);
 
 
                    // Create new EdbTrackP Object for storage;
                    EdbTrackP* RecoShower = new EdbTrackP();
                    RecoShower -> AddSegment(Segment);
                    RecoShower -> AddSegment(Segment2);
                    // Set X and Y and Z values: /Take lower Z of both BTs)
                    RecoShower->SetZ(TMath::Min(Segment->Z(),Segment2->Z()));
                    RecoShower->SetX(Segment_Sum->X());
                    RecoShower->SetY(Segment_Sum->Y());
                    RecoShower->SetTX(Segment_Sum->TX());
                    RecoShower->SetTY(Segment_Sum->TY());
                    RecoShower->SetMC(InBT->MCEvt(),InBT->MCEvt());
                    RecoShower->SetID(eRecoShowerArrayN);
                    RecoShower->SetPID(Segment->PID());
 
                    // Add Shower to interim  Array:
                    RecoShowerArray->Add(RecoShower);
                    ++RecoShowerArrayN;
 
                    if (gEDBDEBUGLEVEL>2) {
                        cout << "Added shower at " << RecoShower << " to reco shower array. Until now: ";
                        cout << RecoShowerArrayN << " entries in it: Print Shower:" << endl;
                        RecoShower->PrintNice();
                    }
 
                    delete segments;
                }
            }
        }
    }  //for (Int_t pat_one_cnt=0; ...
 
    // ------------------------------------
    if (gEDBDEBUGLEVEL>2) {
        cout << "EdbShowerAlg_GS::FindPairs   For the InBT/Vtx at __" << InBT << "__, we have searched "  <<  NPairTriesTotal  << " pair combinations in the PVRec volume." << endl;
        cout << "EdbShowerAlg_GS::FindPairs   For the InBT/Vtx at __" << InBT << "__, we have found "  <<  RecoShowerArray->GetEntries()  << " compatible pairs in the PVRec volume." << endl;
    }
    // ------------------------------------
 
    // Delete unnecessary objects: important, else memory overflow!
    delete vtx;
    delete Segment_Sum;
 
    Log(3,"EdbShowerAlg_GS::FindPairs","Starting FindPairs(InBT,eAli_Sub) now....done.");
    return RecoShowerArray;
}

FindPairsPreselected() [1/2]

TObjArray * EdbShowerAlg_GS::FindPairsPreselected

(

EdbPVRec *

eAli_Sub

)

---

Small Trick.... The BTs are sorted w.r.t. z position!

---

{
    // FindPairsPreselected()
    // What is the goal of this function?
    // A simple example:
    // Consider 30 plates, Each plate has 20000 basetracks in it.
    // GammaPairing from volume would then have to search all these
    // 30*20000=600000 BTs for combinations.
    // That means that for _each_ of these 600k basetracks, we would have
    // to loop and compare with 4*20000 other basetracks to see if we
    // have BT pair combinations possible this _each_ basetrack.
    // (4 is the maximum difference number in plates according to this algo).
    // So 80000 possible combinations to check for one BT, makes in total
    // 80000*600000 bt pairs. That is too much for reasonable time.
    // ---------------
 
    Log(2,"EdbShowerAlg_GS::FindPairsPreselected","FindPairsPreselected()...");
 
    EdbSegP* InBT=NULL;
    EdbSegP* Segment=NULL;
    EdbSegP* Segment2=NULL;
    TObjArray* PreselectedArray= new TObjArray();
    TObjArray* NewPreselectedArray= new TObjArray();
    Int_t counter=0;
    Int_t newcounter=0;
    Int_t npat=eAli_Sub->Npatterns();
    Int_t pat_one_bt_cnt_max,pat_two_bt_cnt_max=0;
    EdbPattern* pat_one=0;
    EdbPattern* pat_two=0;
    EdbPVRec* pvr=eAli_Sub;
 
    cout <<  "EdbShowerAlg_GS::FindPairsPreselected   eAli_Sub->GetEntries()  " << eAli_Sub->NSeg() << endl;
 
    for (int l=0; l<eAli_Sub->Npatterns(); l++) {
        if (gEDBDEBUGLEVEL>2) cout << "Doing pattern()  " << l << endl;
 
        // EdbPattern with couples
        EdbPattern *pat = pvr->GetPattern(l);
 
        // Make a hash table in X,Theta space, with division of x microns and y rad.
        pat->FillCell(500, 500, 0.1, 0.1);
 
        // Loop over all segments.
        for (int j=0; j<pat->N(); j++) {
            EdbSegP *s1 = pat->GetSegment(j);
            // s1 is used as prediction. give your criteria as one sigma.
            // set value should be square.   (dx^2, dy^2, dz^2, dtx^2, dty^2)
            s1->SetErrors(3000,3000,100, 0.05,0.05);
            // dx=10microns, dy=10microns, dz=sqrt(10)microns, dtx=0.1rad, dty=0.1rad.
            // get corresponding tracks in pattern.
            // pat->FindCompliments() will search segments in hash table and put them in segments.
            // last 2 numbers are number of sigma to be taken.
            TObjArray segments; // array to be filled with segments.
            Int_t ncompl = pat->FindCompliments(*s1, segments, 1, 1);
            //cout << "j= " << j << " ncompl=  " << ncompl << "   segments.GetEntries()  " << segments.GetEntries() << endl;
 
            // ncompl==1 means only the BT itself was found:
            if (ncompl==1) continue;
 
            for (int k=0; k<segments.GetEntries(); k++) {
                EdbSegP *s2 = (EdbSegP *) segments.At(k);
                NewPreselectedArray->Add(s2);
            }
        }
        // for (int j=0; j<pat->N(); j++)
    } // for (int l=0; l<eAli_Sub->Npatterns(); l++)
    cout <<  "EdbShowerAlg_GS::FindPairsPreselected   Preselected: "<< NewPreselectedArray->GetEntries() << endl;
 
    // Here it might be that there are still duplicated basetracks (per plate) in the array,
    // these we have to take out:
    cout <<  "EdbShowerAlg_GS::FindPairsPreselected   Check for duplicated basetracks:" << endl;
    cout <<  "EdbShowerAlg_GS::FindPairsPreselected   This might take a while...:" << endl;
    cout <<  "EdbShowerAlg_GS::FindPairsPreselected   (status bar):" << endl;
    Int_t count=0;
    Int_t Zseparator[114];
    for (Int_t l=0; l< 114; l++) {
        Zseparator[l]=0;
    }
    EdbSegP *sl = (EdbSegP*) NewPreselectedArray->At(0);
    Zseparator[0]=0;
    Float_t ZStart=sl->Z();
    Int_t incrementor=1;
    Int_t NewPreselectedArrayN=NewPreselectedArray->GetEntries();
 
    for (Int_t l=0; l< NewPreselectedArrayN; l++) {
        EdbSegP *sl = (EdbSegP*) NewPreselectedArray->At(l);
        if (sl->Z()!=ZStart) {
            Zseparator[incrementor]=l;
            ZStart=sl->Z();
            incrementor++;
        }
    }
 
 
    if (gEDBDEBUGLEVEL>2) {
        for (Int_t l=0; l< 114; l++) {
            cout <<  "EdbShowerAlg_GS::FindPairsPreselected   incrementor = "  << incrementor << endl;
        }
    }
 
 
    for (Int_t o=0; o < incrementor-1; o++) {
 
        cout << "." ;
 
        if (gEDBDEBUGLEVEL>2) cout << "EdbShowerAlg_GS::FindPairsPreselected   Doing ZSeparator o= "  << o << " loop from BT #" << Zseparator[o] << " to BT # " << Zseparator[o+1]-1 << endl;
 
 
        for (Int_t l=Zseparator[o]; l< Zseparator[o+1]-1; l++) {
 
            Bool_t IsDuplicateThere=kFALSE;
            EdbSegP *sl = (EdbSegP*) NewPreselectedArray->At(l);
 
 
            for (Int_t k=l+1; k<  Zseparator[o+1]; k++) {
                EdbSegP *sk = (EdbSegP*) NewPreselectedArray->At(k);
 
                if (TMath::Abs(sl->X()-sk->X())<0.1 && TMath::Abs(sl->Y()-sk->Y())<0.1 ) {
                    if (gEDBDEBUGLEVEL>2) {
                        cout << "Duplicate found: " << sl << " and " << sk << endl;
                        sl->PrintNice();
                        sk->PrintNice();
                    }
                    ++count;
                    IsDuplicateThere=kTRUE;
                }
            } //for (Int_t k=l+1; k<  Zseparator[o+1]; k++)
            if (!IsDuplicateThere) PreselectedArray->Add(sl);
 
        } //    for (Int_t l=Zseparator[o]; l< Zseparator[o+1]-1; l++)
    }       // for (Int_t o=0; o < incrementor; o++) {
 
 
    if (gEDBDEBUGLEVEL>2) {
        cout << "EdbShowerAlg_GS::FindPairsPreselected   NDuplicates = " << count << endl;
        cout << "EdbShowerAlg_GS::FindPairsPreselected   NNewPreselectedArray = " << NewPreselectedArray->GetEntries()  << endl;
    }
    cout << "EdbShowerAlg_GS::FindPairsPreselected   NPreselectedArray = " << PreselectedArray->GetEntries() << endl;
 
    Log(2,"EdbShowerAlg_GS::FindPairsPreselected","FindPairsPreselected()...done.");
    return PreselectedArray;
}

FindPairsPreselected() [2/2]

Bool_t EdbShowerAlg_GS::FindPairsPreselected	(	EdbSegP *	InitiatorBT,
		EdbPVRec *	eAli_Sub
	)

{
    Log(3,"EdbShowerAlg_GS::FindPairsPreselected","Starting FindPairsPreselected(InBT,eAli_Sub) now....");
 
    EdbSegP* InBT=NULL;
    InBT=InitiatorBT;
    EdbSegP* Segment=NULL;
    EdbSegP* Segment2=NULL;
 
    Int_t npat=eAli_Sub->Npatterns();
    Int_t pat_one_bt_cnt_max,pat_two_bt_cnt_max=0;
    EdbPattern* pat_one=0;
    EdbPattern* pat_two=0;
 
    //-----------------------------------
    // ---   Loop over first patterns
    for (Int_t pat_one_cnt=0; pat_one_cnt<npat; ++pat_one_cnt) {
        pat_one=(EdbPattern*)eAli_Sub->GetPattern(pat_one_cnt);
        pat_one_bt_cnt_max=eAli_Sub->GetPattern(pat_one_cnt)->GetN();
 
        // Check if dist Z to vtx (BT) is ok:
        Float_t distZ=pat_one->Z()-InBT->Z();
        if (gEDBDEBUGLEVEL>4) {
            cout << "EdbShowerAlg_GS::FindPairsPreselected  Check if dist Z to vtx (BT) is ok:  distZ=" << distZ << endl;
        }
        if (distZ>eParaValue[3]) continue;
        if (distZ<-1000) continue;
 
        if (gEDBDEBUGLEVEL>2) cout << "EdbShowerAlg_GS::FindPairsPreselected   pat_one_cnt=" << pat_one_cnt << "  pat_one->Z() = " << pat_one->Z() << " pat_one_bt_cnt_max= "<< pat_one_bt_cnt_max <<endl;
 
        //-----------------------------------
        // ---   Loop over second patterns
        for (Int_t pat_two_cnt=0; pat_two_cnt<npat; ++pat_two_cnt) {
 
            // Now apply cut conditions: GS  GAMMA SEARCH Alg:
            if (TMath::Abs(pat_one_cnt-pat_two_cnt)>eParaValue[5]) continue;
 
            pat_two=(EdbPattern*)eAli_Sub->GetPattern(pat_two_cnt);
            pat_two_bt_cnt_max=eAli_Sub->GetPattern(pat_two_cnt)->GetN();
 
            distZ=pat_two->Z()-InBT->Z();
            if (distZ>eParaValue[3]) continue;
            if (distZ<-1000) continue;
 
            if (gEDBDEBUGLEVEL>2) cout << "EdbShowerAlg_GS::FindPairsPreselected    pat_two_cnt=" << pat_two_cnt << "  pat_two->Z() = " << pat_two->Z() << " pat_two_bt_cnt_max= "<< pat_two_bt_cnt_max <<endl;
 
            //-----------------------------------
            // ---   Loop over basetracks of first pattern
            for (Int_t pat_one_bt_cnt=0; pat_one_bt_cnt<pat_one_bt_cnt_max; ++pat_one_bt_cnt) {
                Segment =  (EdbSegP*)pat_one->GetSegment(pat_one_bt_cnt);
 
                //-----------------------------------
                // ---   Loop over basetracks of second pattern
                for (Int_t pat_two_bt_cnt=0; pat_two_bt_cnt<pat_two_bt_cnt_max; ++pat_two_bt_cnt) {
 
                    Segment2 = (EdbSegP*)pat_two->GetSegment(pat_two_bt_cnt);
                    if (Segment2==Segment) continue;
                    if (Segment2->ID()==Segment->ID()&&Segment2->PID()==Segment->PID()) continue;
 
                    // a) Check dR between tracks:
                    if (DeltaR_WithPropagation(Segment,Segment2)<eParaValue[2]) return kTRUE;
                    // b) Check dT between tracks:
                    if (DeltaThetaSingleAngles(Segment,Segment2)>eParaValue[4]) return kTRUE;
                    // c) Check dMinDist between tracks:
                    if (GetMinimumDist(Segment,Segment2)>eParaValue[1]) return kTRUE;
 
                }
            }
        }
    }  //for (Int_t pat_one_cnt=0; ...
 
    Log(3,"EdbShowerAlg_GS::FindPairsPreselected","Starting FindPairsPreselected(InBT,eAli_Sub) now....done.");
    return kFALSE;
}

GetInVtxArray()

TObjArray* EdbShowerAlg_GS::GetInVtxArray ( ) const

inline

                                                            {
        return eInVtxArray;
    }

GetInVtxArrayN()

Int_t EdbShowerAlg_GS::GetInVtxArrayN ( ) const

inline

                                                            {
        return eInVtxArrayN;
    }

Init()

void EdbShowerAlg_GS::Init

(

void

)

{
    //  Init with values according to GS Alg.
    //  Create essential objects.
    //  Init function is supposed to be executed exact once for
    //  each class instance.
    //
    //  Took over from "FindGamma.C" script I develoved before:
    //  and with recorded values from a GammaSearch best algo parameterset.
    //  See also table "tableshowerchosenparametersets"
    //  and figure "" in FWM thesis.
    //
 
    Log(2,"EdbShowerAlg_GS::EdbShowerAlg_GS","Init()");
 
    eDebug=kFALSE;
 
    eRecoMode=0;
    eCutModeFull=kTRUE;
 
    eAlgName="GS";
    eAlgValue=999; // see default.par_SHOWREC for labeling (labeling identical with ShowRec program)
 
    eFindPairsPreselected=kTRUE;
 
    eANNPairCut[0]=0.5;
    eANNPairCut[1]=0.48;
    eANNPairCut[2]=0.45;
 
    //
    // min IP cut; this cut is used for the -better- IP of both BTs to Vertex/BT
    eParaValue[0]=393;
    eParaString[0]="PARA_GS_CUT_dIP";
    //
    // min minDist.e between pair BTs
    eParaValue[1]=35;
    eParaString[1]="PARA_GS_CUT_dMin";
    //
    // min dR (with propagation) between pair BTs
    eParaValue[2]=85;
    eParaString[2]="PARA_GS_CUT_dR";
    //
    // max Z distance between pair BTs and Vertex/BT
    eParaValue[3]=25000;
    eParaString[3]="PARA_GS_CUT_dZ";
    //
    // max Angle between pair BTs
    eParaValue[4]=0.11;
    eParaString[4]="PARA_GS_CUT_dtheta";
    //
    // max plates difference between pair BTs
    eParaValue[5]=1;
    eParaString[5]="PARA_GS_CUT_PIDDIFF";
    //
    // in MC case: check for opposite flag sign
    eParaValue[6]=1;
    eParaString[6]="PARA_GS_CUT_OPPOSITEFLAG";
 
 
    // Variable to clean the found parings once more.
    // Default is kTRUE.
    eSetCleanPairs=kFALSE;
    eSetCleanPairs=kTRUE;
 
    // Create needed objects for storing Initiator Vertices/Basetracks.
    eRecoShowerArray= new TObjArray(999);
    eRecoShowerArrayN=0;
 
    eInVtxArray = new TObjArray();
    eInVtxArrayN=0;
 
    eInBTArray = new TObjArray();
    eInBTArrayN=0;
 
    eSegmentPIDArray = new TArrayI(9999);
    eSegmentIDArray = new TArrayI(9999);
    eSegment2PIDArray = new TArrayI(9999);
    eSegment2IDArray = new TArrayI(9999);
 
    if (gEDBDEBUGLEVEL>2) {
        cout << "EdbShowerAlg_GS eRecoShowerArray = " << eRecoShowerArray << endl;
        cout << "EdbShowerAlg_GS eInVtxArray = " << eInVtxArray << endl;
        cout << "EdbShowerAlg_GS eInBTArray = " << eInBTArray << endl;
    }
 
    Log(2,"EdbShowerAlg_GS::EdbShowerAlg_GS","Init()...done.");
    return;
}

Initialize()

void EdbShowerAlg_GS::Initialize ( )

virtual

Reimplemented from EdbShowerAlg.

{
    Log(2,"EdbShowerAlg_GS::EdbShowerAlg_GS","Initialize()");
    return;
}

IsPossibleFakeDoublet()

Bool_t EdbShowerAlg_GS::IsPossibleFakeDoublet	(	EdbSegP *	s1,
		EdbSegP *	s2
	)

{
    if (TMath::Abs(s1->X()-s2->X())<1) return kTRUE;  // minimum distance of 1micron
    if (TMath::Abs(s1->Y()-s2->Y())<1) return kTRUE;// minimum distance of 1micron
    if (TMath::Abs(s1->TX()-s2->TX())<0.005) return kTRUE;// minimum angle of 5mrad
    if (TMath::Abs(s1->TY()-s2->TY())<0.005) return kTRUE;// minimum angle of 5mrad
    return kFALSE;
}

ReloadANNs()

void EdbShowerAlg_GS::ReloadANNs

(

Int_t

RecoMode

)

{
    Log(2,"EdbShowerAlg_GS::ReloadANNs","ReloadANNs(Int_t RecoMode)");
 
    if (RecoMode==0) {
        eANNPair=eANNPairCaseA;
        Log(2,"EdbShowerAlg_GS::ReloadANNs","Set eANNPair to eANNPairCaseA");
        cout << "EdbShowerAlg_GS::ReloadANNs: eANNPairCaseA at " << eANNPairCaseA << endl;
    }
    else if (RecoMode==1) {
        eANNPair=eANNPairCaseB;
        Log(2,"EdbShowerAlg_GS::ReloadANNs","Set eANNPair to eANNPairCaseB");
        cout << "EdbShowerAlg_GS::ReloadANNs: eANNPairCaseB at " << eANNPairCaseB << endl;
    }
    else if (RecoMode==2) {
        eANNPair=eANNPairCaseC;
        Log(2,"EdbShowerAlg_GS::ReloadANNs","Set eANNPair to eANNPairCaseC");
        cout << "EdbShowerAlg_GS::ReloadANNs: eANNPairCaseC at " << eANNPairCaseC << endl;
    }
    else {
        cout << "ELSE  RETURN" << endl;
    }
    Log(2,"EdbShowerAlg_GS::ReloadANNs","ReloadANNs()...done.");
}

SelectHighestPInMCArray()

TObjArray * EdbShowerAlg_GS::SelectHighestPInMCArray

(

TObjArray *

BTArray

)

{
    // Helper function for MC purposes:
    // This  function inputs a EdbSegP (Initiator basetrack array) and
    // for all MC Events it looks for the highest P of the BT.
    // Mainly used for looking the first BT of a photon/electron shower.
    // Affects only MC Events, on data it should have no effect.
    // Only internal usage, should not needed by the everyday user.
 
    Log(2,"EdbShowerAlg_GS::SelectHighestPInMCArray","SelectHighestPInMCArray()");
 
    if (BTArray->GetEntries()==0) return NULL;
 
    TObjArray* HighestPBTArray = new TObjArray();
 
    Float_t maxP[999999];
    Int_t index[999999];
    for (int i=0; i<100000; ++i) {
        maxP[i]=0;
        index[i]=-1;
    }
 
    for (int i=0; i<BTArray->GetEntries(); ++i) {
        EdbSegP* s = (EdbSegP*)BTArray->At(i);
        if (s->MCEvt()<0) continue;
        if (s->P() > maxP[s->MCEvt()]) {
            maxP[s->MCEvt()] = s->P();
            index[s->MCEvt()]=i;
        }
    }
 
    for (int i=0; i<BTArray->GetEntries(); ++i) {
        EdbSegP* s = (EdbSegP*)BTArray->At(i);
        if (s->MCEvt()<0) HighestPBTArray->Add(s);
        if (s->MCEvt()>0) {
            if (s->P()< maxP[s->MCEvt()] ) {
                // do nothing, since it is not maximum for this event.
            }
            else {
                HighestPBTArray->Add(s);
            }
        }
    }
 
    if (gEDBDEBUGLEVEL>2) {
        cout << "EdbShowerAlg_GS::SelectHighestPInMCArray   BTArray->GetEntries() " << BTArray->GetEntries() << endl;
        cout << "EdbShowerAlg_GS::SelectHighestPInMCArray   HighestPBTArray->GetEntries() " << HighestPBTArray->GetEntries() << endl;
 
        cout << "EdbShowerAlg_GS::SelectHighestPInMCArray   Print BTArray: " << BTArray->GetEntries() << endl;
 
        for (int i=0; i<BTArray->GetEntries(); ++i) {
            EdbSegP* s = (EdbSegP*)BTArray->At(i);
            s->PrintNice();
        }
        cout << "EdbShowerAlg_GS::SelectHighestPInMCArray   Print HighestPBTArray: " << HighestPBTArray->GetEntries() << endl;
        for (int i=0; i<HighestPBTArray->GetEntries(); ++i) {
            EdbSegP* s = (EdbSegP*)HighestPBTArray->At(i);
            s->PrintNice();
        }
    }
 
    Log(2,"EdbShowerAlg_GS::SelectHighestPInMCArray","SelectHighestPInMCArray()...done.");
 
    return HighestPBTArray;
}

Set0()

void EdbShowerAlg_GS::Set0

(

)

{
    // Reset all internal variables.
 
    Log(2,"EdbShowerAlg_GS::EdbShowerAlg_GS","Set0()");
 
    eRecoMode=0;
    eCutModeFull=kTRUE;
 
    eInVtxArrayN=0;
    eInVtxArray->Clear();
 
    eInBTArrayN=0;
    eInBTArray->Clear();
 
    eRecoShowerArrayN=0;
    eRecoShowerArray->Clear();
 
    eAlgName="GS";
    eAlgValue=999;
    eANNPairCut[0]=0.5;
    eANNPairCut[1]=0.48;
    eANNPairCut[2]=0.45;
    eSetCleanPairs=kTRUE;
    eFindPairsPreselected=kTRUE;
 
    eParaValue[0]=393;
    eParaString[0]="PARA_GS_CUT_dIP";
    eParaValue[1]=35;
    eParaString[1]="PARA_GS_CUT_dMin";
    eParaValue[2]=85;
    eParaString[2]="PARA_GS_CUT_dR";
    eParaValue[3]=25000;
    eParaString[3]="PARA_GS_CUT_dZ";
    eParaValue[4]=0.11;
    eParaString[4]="PARA_GS_CUT_dtheta";
    eParaValue[5]=1;
    eParaString[5]="PARA_GS_CUT_PIDDIFF";
    eParaValue[6]=1;
    eParaString[6]="PARA_GS_CUT_OPPOSITEFLAG";
 
 
    eSegmentPIDArray->Reset();
    eSegmentIDArray->Reset();
    eSegment2PIDArray->Reset();
    eSegment2IDArray->Reset();
 
    Log(2,"EdbShowerAlg_GS::EdbShowerAlg_GS","Set0()...done.");
    return;
}

SetCleanPairs()

void EdbShowerAlg_GS::SetCleanPairs ( Bool_t  CleanPairs )

inline

                                                        {
        eSetCleanPairs = CleanPairs;
    }

SetCutModeFull()

void EdbShowerAlg_GS::SetCutModeFull ( Bool_t  CutModeFull )

inline

                                                          {
        eCutModeFull = CutModeFull;
    }

SetFindPairsPreselected()

void EdbShowerAlg_GS::SetFindPairsPreselected ( Bool_t  FindPairsPreselected )

inline

                                                                            {
        eFindPairsPreselected = FindPairsPreselected;
    }

SetInVtx()

void EdbShowerAlg_GS::SetInVtx

(

EdbVertex *

vtx

)

{
    Log(3,"EdbShowerAlg_GS::SetInVtx","SetInVtx()");
    if (eInVtxArrayN!=0) {
        Log(2,"EdbShowerAlg_GS::SetInVtx","SetInVtx()   WARNING! Array not empty. Clear/Reset it!");
        eInVtxArray -> Clear();
        eInVtxArrayN=0;
    }
    eInVtxArray->Add(vtx);
    ++eInVtxArrayN;
    eInVtxArraySet=kTRUE;
    cout << "EdbShowerAlg_GS::SetInVtx Added One Vertex. Now there are " << eInVtxArrayN << " InVtx stored." << endl;
    Log(3,"EdbShowerAlg_GS::SetInVtx","SetInVtx()...done");
    return;
}

SetInVtxArray()

void EdbShowerAlg_GS::SetInVtxArray ( TObjArray *  InVtxArray )

inline

                                                               {
        eInVtxArray = InVtxArray;
        eInVtxArrayN = eInVtxArray->GetEntries();
        cout << "SetInVtxArray:: " << eInVtxArrayN << "  entries set"<<endl;
    }

SetRecoMode()

void EdbShowerAlg_GS::SetRecoMode

(

Int_t

RecoMode

)

{
    // Manually Set Reco Mode A or B or C
    // Check if InVtxArray is Filled.
    // If not, do nothing, cause it does make no sense when
    // we wanna do Reco Mode and have no vertices filled.
 
    Log(2,"EdbShowerAlg_GS::SetRecoMode","SetRecoMode()...");
 
    if (eRecoMode>2) {
        Log(2,"EdbShowerAlg_GS::SetRecoMode","WARNING eRecoMode>2!!!");
        Log(2,"EdbShowerAlg_GS::SetRecoMode","Set automatically Reco Mode C (eRecoMode=2).");
        eRecoMode=2;
    }
 
    if (eInVtxArrayN==0) {
 
        // So, welches nehmen wir nun denn jetzt by default??
        // Die ROC Curve sagt eigentlich, dass CASE C besser sei...
 
        Log(2,"EdbShowerAlg_GS::SetRecoMode","No InVtxArray there. Will not set Reco Mode A.");
        Log(2,"EdbShowerAlg_GS::SetRecoMode","Set automatically Reco Mode B (eRecoMode=1).");
        eRecoMode=1;
 
        Log(2,"EdbShowerAlg_GS::SetRecoMode","No InVtxArray there. Will not set Reco Mode A.");
        Log(2,"EdbShowerAlg_GS::SetRecoMode","Set automatically Reco Mode C (eRecoMode=2).");
        eRecoMode=2;
 
    }
    else {
        eRecoMode=RecoMode;
        ReloadANNs(RecoMode);
    }
 
    cout << "DEBUG  eRecoMode  " << eRecoMode << endl;
 
    // Set the dZ cutvalue to just one plate after, since we start
    // in these cases from real Initiator Basetracks (and not vertices)
    // so we just look for pairs directly related to this InBT.
    if (eRecoMode==1 || eRecoMode == 2) {
        SetParameter(3,1800);
    }
    else {
        SetParameter(3,25000);
    }
 
    Log(2,"EdbShowerAlg_GS::SetRecoMode","SetRecoMode() eRecoMode = %d ",eRecoMode);
    Log(2,"EdbShowerAlg_GS::SetRecoMode","SetRecoMode()...done.");
    return;
}

Member Data Documentation

eANNPair

TMultiLayerPerceptron* EdbShowerAlg_GS::eANNPair

private

eANNPairCaseA

TMultiLayerPerceptron* EdbShowerAlg_GS::eANNPairCaseA

private

eANNPairCaseB

TMultiLayerPerceptron* EdbShowerAlg_GS::eANNPairCaseB

private

eANNPairCaseC

TMultiLayerPerceptron* EdbShowerAlg_GS::eANNPairCaseC

private

eANNPairCut

Float_t EdbShowerAlg_GS::eANNPairCut[3]

private

eANNPairTree

TTree* EdbShowerAlg_GS::eANNPairTree

private

eCutModeFull

Bool_t EdbShowerAlg_GS::eCutModeFull

private

eFindPairsPreselected

Bool_t EdbShowerAlg_GS::eFindPairsPreselected

private

eInVtxArray

TObjArray* EdbShowerAlg_GS::eInVtxArray

private

eInVtxArrayN

Int_t EdbShowerAlg_GS::eInVtxArrayN

private

eInVtxArraySet

Bool_t EdbShowerAlg_GS::eInVtxArraySet

private

eRecoMode

Int_t EdbShowerAlg_GS::eRecoMode

private

eSegment2IDArray

TArrayI* EdbShowerAlg_GS::eSegment2IDArray

private

eSegment2PIDArray

TArrayI* EdbShowerAlg_GS::eSegment2PIDArray

private

eSegmentIDArray

TArrayI* EdbShowerAlg_GS::eSegmentIDArray

private

eSegmentPIDArray

TArrayI* EdbShowerAlg_GS::eSegmentPIDArray

private

eSetCleanPairs

Bool_t EdbShowerAlg_GS::eSetCleanPairs

private

eValueGSNN_var00

Float_t EdbShowerAlg_GS::eValueGSNN_var00

private

eValueGSNN_var01

Float_t EdbShowerAlg_GS::eValueGSNN_var01

private

eValueGSNN_var02

Float_t EdbShowerAlg_GS::eValueGSNN_var02

private

eValueGSNN_var03

Float_t EdbShowerAlg_GS::eValueGSNN_var03

private

eValueGSNN_var04

Float_t EdbShowerAlg_GS::eValueGSNN_var04

private

eValueGSNN_var05

Float_t EdbShowerAlg_GS::eValueGSNN_var05

private

eValueGSNN_var06

Float_t EdbShowerAlg_GS::eValueGSNN_var06

private

eValueGSNN_varInput

Float_t EdbShowerAlg_GS::eValueGSNN_varInput

private

eValueGSNN_varOutput

Float_t EdbShowerAlg_GS::eValueGSNN_varOutput

private

---

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

• src/libShower/EdbShowerAlg.h
• src/libShower/EdbShowerAlg.cxx