// Author: M.A. Kagarlis
// Written: 21.01.99
// Revised: 05.12.05  by R. Holzmann
// Revised: 29.08.06  by R. Holzmann
// PChannel Class Header

#ifndef _PCHANNEL_H_
#define _PCHANNEL_H_

#include "TArrayI.h"
#include "TF1.h"
#include "TF2.h"
//#include "PSaid.h"
#include "PParticle.h"
#include "PDistribution.h"
#include "PStdData.h"
#include "PBulkInterface.h"
#include "PProjector.h"
//#include "PPlutoBulkDecay.h"

#define MAX_DISTRIBUTIONS 10
#define MAX_NUM_NOT_FINALIZED 5
#define MAX_BULKDECAY 5

#define STATUS_OK          0
#define STATUS_EMPTY_EVENT 80
#define STATUS_PSEUDO_EMPTY_EVENT 81
#define STATUS_REDO_CHAIN  82
#define STATUS_NOT_DECAYED 99



class PChannel : public TObject {
  
 public:

    PChannel(PParticle **, int nt=2, int mf=0, int af=0, int bf=0);
    // Channel constructor by particle array

    PChannel(int idx, PParticle **particles=NULL, int mf=0, int af=0, int bf=0);
    // Channel constructor by decay-mode index

    PChannel(int idx, PParticle &parent, int mf=0, int af=0, int bf=0);
    // Channel constructor by decay-mode index and parent reference

    ~PChannel() {
	// Channel destructor
	if (ipid.GetArray()) 
	    ipid.~TArrayI();
	// necessary, because of constr 2+3, for constr 1 also better
	if (ptcls && ptcls_extern>=0) {
	    for(int i=ptcls_extern; i<=n; i++)
		if (ptcls[i]) delete ptcls[i];
	    if (ptcls) delete [] ptcls;
	}
    }

    Bool_t Init();
    // sets on the first call all parent, siblings etc IDs. This
    // has to be done only once per PReaction

    Bool_t Reset();
    // cleans the channel, so an init will be done in the next decay()

    Bool_t SetDaughters();

    void SetSourceId(Int_t mysourceid) {
	sourceid = mysourceid;
    }

    int Decay();
    // N-body phase-space decay function based on the CERNLIB routine 
    // GENBOD. Resonance and dilepton masses (for Dalitz decays)
    // are sampled in PData. Scattering angles are sampled here, for
    // selected channels (see idChannel).

    int decay() {
	//kept for comp.
	return Decay();
    }; 

    PParticle **GetParticles() {
	// address of array of pointers to the channel particles
	return ptcls;
    }

    int *GetPids() {
	// address of array of channel-particle ids
	return ipid.GetArray();
    }   

    int GetNumPar() {  
	// number of decay products
	return n; 
    }
    
    //******for distribution handling:
    int SetDistribution(PDistribution *distribution);
    // return -1 if distribution does not match the channel
    int CheckSiblings(PParticle *p, PDistribution *dist, int flag); //check for siblings

    double GetBT() {
	PParticle *beam   = ptcls[0]->GetScattering(0);
	PParticle *target = ptcls[0]->GetScattering(1);
	
	if (beam) {
	    return (beam->KE() > target->KE()? beam->KE(): target->KE());
	}
	else return ptcls[0]->Rapidity();
    }

    //******for quasi-elastic handling:
    PChannel *GetQuasi(void) {
	return quasi_pchannel;
    };
    void ClearQuasi(void) {
	quasi_pchannel = NULL;
    };

    int GetParentSize() { 
	// 1 or 2 for elementary or quasi particles respectively
	return 1 + (ipid[0]>=1000); 
    }   

    static void SetGlobalWeight(double w) { 
	globalWeight = w; 
    } 
    static double GetGlobalWeight() { 
	return globalWeight; 
    }

    int GetDMIndex() { 
	//returns the decay mode index of this channel
	return DMIndex; 
    }

    void GetMessage();
    // status message

    void Print(const Option_t *delme="") const;
    void PrintReaction(Int_t check_key = 1) const;
    void PrintNew();
    
    void PrintReport() const; //Print a final report after sampling
    void CheckDecayKey() const;

    char const *GetName(void) const ;

//  double ds_dt(double);
    // ...> moved to PParticle (IF)
    // ds/dt(cos_th_cm) in the cm for N+N->N+Delta. With the mass of Delta
    // sampled independently in PData, the simulated events are distributed
    // as ds/dOmega (normalized). Ref: NPA459 (1986) 503

    void SetPrintTentative(bool t){ 
	// Print the tentative messages on the 1st call
	print_tentative = t;
    };

    void DisableHelicityAngle() {
	cout << "DisableHelicityAngle() is obsolete -> use DistributionManager" << endl;
    };

    //this is for envelope distributions:
    int GetNumNotFinalized() {
	return num_not_finalized;
    };
    PDistribution *GetDistributionNotFinalized(int i) {
	return distribution_not_finalized[i];
    };
    
    Bool_t AddBulk(PBulkInterface *mybulk);
    Bool_t AddPrologueBulk(PBulkInterface *mybulk); //Bulk IO before any decay

    Bool_t Do(const char *command) {
	return GetCurrentProjector()->AddCommand(command);
    }
    Bool_t Do(TH1F *f, const char *command, Int_t flag=1) {
	return GetCurrentProjector()->AddHistogram(f,command,flag);
    }
    Bool_t Do(TH2F *f, const char *command, Int_t flag=1) {
	return GetCurrentProjector()->AddHistogram(f,command,flag);
    }
    Bool_t Do(TH3F *f, const char *command, Int_t flag=1) {
	return GetCurrentProjector()->AddHistogram(f,command,flag);
    }

 private:

    int n, status, DMIndex;
    // number of products, error status
    // decay-mode index matching PData modes

    int thSrc, tcSrc, dlSrc, sourceid;

    Double_t weight_sum;   //Sum of all weights for debugging

    int decay_key;                       //! Key to data base entry
    TString decay_string; //! For parsing the decay string
    TString decay_string2; 

    // For generic distributions
    PDistribution *dist[MAX_DISTRIBUTIONS];
    Double_t dist_weight[MAX_DISTRIBUTIONS];
    Double_t dist_weight_sum[MAX_DISTRIBUTIONS];
    Double_t dist_counter[MAX_DISTRIBUTIONS];

    Int_t distribution_position;
    bool init_done, thermal_disable_didx;
    bool print_tentative;
    Long_t fEnablePattern;

    int bid,tid;
    PParticle *spectator, *participant, *quasi_composite;
    PChannel *quasi_pchannel;
    // general quasifree particle production on the deuteron  (V. Hejny 13/10/00)

    TArrayI ipid;
    // particle id array

    static double globalWeight;
    // global weight for the event (see PDecayManager, BUGBUG: never used)

    double *event_impact_param , *event_plane;  //from batch system
    double *weight_version;  //from batch system
    double *events;   //from batch system

    PParticle ** ptcls, *parent, *orig_parent, *grandparent, *grandgrandparent;
    // array of channel particles

    int ptcls_extern;
    // flag whether ptcls are coming from external source (do not delete those!)

    double ecm, w, emin; // invariant mass, weight of decay, energy threshold
    double e_cm;  // used to save parameters between calls of ds_dt()


    void IsReaction();
    // if a reaction channel identify beam, target, and decay mode

    void IdChannel();
    // identify reaction channel of known angular distribution
    // these are: p+p->p+p (elastic), N+N->N+Delta, pi+N->N+w, pi+ + P -> pi+ + P + w

    void ThermalSampling();
    // parent is midrapidity thermal source

    Int_t ReadFileInput();
    // parent is file input interface

    void MakeDilepton();
    // parent is a sampled dilepton


    int Genbod(int);
    // N-body phase-space decay of a particle in its rest frame.
    // This follows the CERNLIB routine GENBOD, in C++

    //this is for envelope distributions:
    int num_not_finalized;
    PDistribution *distribution_not_finalized[MAX_NUM_NOT_FINALIZED];

    int bulkdecay_pos, pro_bulkdecay_pos;
    
    PBulkInterface *bulk[MAX_BULKDECAY];
    PBulkInterface *pro_bulk[MAX_BULKDECAY];
    PProjector *current_projector;

    PProjector *GetCurrentProjector(void) {
	if (!bulkdecay_pos) {
	    current_projector = new PProjector();
	    AddBulk(current_projector);
	} else {
	    if ((strcmp("PProjector",bulk[bulkdecay_pos-1] -> GetName())==0)) {
		if ((bulk[bulkdecay_pos-1] -> GetPriority() > FILTER_PRIORITY) ||
		    ((PProjector *)bulk[bulkdecay_pos-1])->IsFileOpen()) {
		    current_projector = (PProjector *)bulk[bulkdecay_pos-1];
		} else {
		    current_projector = new PProjector();
		    AddBulk(current_projector);
		}
	    } else {		
		current_projector = new PProjector();
		AddBulk(current_projector);
	    }
	}
	return current_projector;
    }
    
    ClassDef(PChannel,0) //Pluto Channel Class

};

#endif // _PCHANNEL_H_