[u/mrichter/AliRoot.git] / TFluka / mgdraw.cxx

#include <Riostream.h>
#include "TVirtualMCApplication.h"
#include "TVirtualMCStack.h"
#include "TFluka.h"
#include "TFlukaCodes.h"

// Fluka includes
#include "Fdimpar.h"  //(DIMPAR) 
#include "Fdblprc.h"  //(DBLPRC) 
#include "Ftrackr.h"  //(TRACKR) 
#include "Fopphst.h"  //(OPPHST) 
#include "Fflkstk.h"  //(FLKSTK) 
#include "Fltclcm.h"  //(LTCLCM) 
#include "Fpaprop.h"  //(PAPROP) 
#include "Falldlt.h"  //(ALLDLT) 
#include "Fdpdxcm.h"  //(DPDXCM) 
#include "Fflkmat.h"  //(FLKMAT) 

#include "TGeoManager.h"

#ifndef WIN32
# define mgdraw mgdraw_
#else
# define mgdraw MGDRAW
#endif

extern "C" {
void mgdraw(Int_t& icode, Int_t& mreg)
{
    TFluka* fluka =  (TFluka*) gMC;
    if (mreg == fluka->GetDummyRegion()) return;
//
//  Make sure that stack has currrent track Id
//
    Int_t trackId = -1;
    TVirtualMCStack* cppstack = fluka->GetStack();
    
    if (TRACKR.jtrack == -1) {
        trackId = OPPHST.louopp[OPPHST.lstopp];
        if (trackId == 0) {
            trackId = FLKSTK.ispark[FLKSTK.npflka][mkbmx2-1];
        }
    } else {
        trackId = TRACKR.ispusr[mkbmx2-1];
    }
    
    Int_t verbosityLevel = fluka->GetVerbosityLevel();
    if (TRACKR.jtrack < -6) {
       // from -7 to -12 = "heavy" fragment
       // assing parent id
       // id < -6 was skipped in stuprf =>   if (kpart < -6) return;
       if (verbosityLevel >= 3) {
          cout << "mgdraw: (heavy fragment) jtrack < -6 =" << TRACKR.jtrack
               << " assign parent pdg=" << fluka->PDGFromId(TRACKR.ispusr[mkbmx2 - 3]) << endl;
       }
//       TRACKR.jtrack = TRACKR.ispusr[mkbmx2 - 3];
    }
    
    cppstack->SetCurrentTrack(trackId);
//
//    
    Int_t mlttc = TRACKR.lt1trk; // LTCLCM.mlatm1;
    fluka->SetMreg(mreg, mlttc);
    fluka->SetIcode((FlukaProcessCode_t) icode);
    fluka->SetCaller(kMGDRAW);

    Int_t nodeId;
    Int_t volId   = fluka->CurrentVolID(nodeId);
    Int_t crtlttc = gGeoManager->GetCurrentNodeId()+1;

    // check region lattice consistency (debug Ernesto)
    // *****************************************************
    if(verbosityLevel>=3 && mreg != volId  && !gGeoManager->IsOutside() ) {
       cout << "  mgdraw:   track=" << trackId << " pdg=" << fluka->PDGFromId(TRACKR.jtrack)
            << " icode=" << icode << " gNstep=" << fluka->GetNstep() << endl
            << "               fluka   mreg=" << mreg << " mlttc=" << mlttc << endl
            << "               TGeo   volId=" << volId << " crtlttc=" << crtlttc << endl
            << "     common TRACKR   lt1trk=" << TRACKR.lt1trk << " lt2trk=" << TRACKR.lt2trk << endl
            << "     common LTCLCM   newlat=" << LTCLCM.newlat << " mlatld=" <<  LTCLCM.mlatld << endl
            << "                     mlatm1=" << LTCLCM.mlatm1 << " mltsen=" <<  LTCLCM.mltsen << endl
            << "                     mltsm1=" << LTCLCM.mltsm1 << " mlattc=" << LTCLCM.mlattc << endl;
        if( mlttc == crtlttc ) cout << "   *************************************************************" << endl;
    }
    // *****************************************************

    
    Int_t med   = FLKMAT.medium[mreg - 1];  // Medium
    Int_t msd   = DPDXCM.msdpdx[med  - 1];  // Iionisation model

    if (!TRACKR.ispusr[mkbmx2 - 2]) {
	if (verbosityLevel >= 3) {
	    cout << endl << "mgdraw: energy deposition for:" << trackId
		 << " icode=" << icode
		 << " pdg=" << fluka->PDGFromId(TRACKR.jtrack)
		 << " flukaid="<< TRACKR.jtrack
		 << " mreg=" << mreg
		 << " np  =" << ALLDLT.nalldl
		 << endl;
	         
	}

	if (msd > 0) {
//
// Primary ionsiations
	    Int_t nprim = ALLDLT.nalldl;
	    if (nprim >= mxalld) {
		nprim = mxalld;
		Warning("mgdraw", "nprim > mxalld, nprim: %6d  pdg: %6d mreg %6d p %13.3f step %13.3f\n", 
			ALLDLT.nalldl, 
			fluka->PDGFromId(TRACKR.jtrack), 
			mreg, 
			TRACKR.ptrack, 
			TRACKR.ctrack);
	    }
	    fluka->PrimaryIonisationStepping(nprim);
	} else {
	    // Single step
	    (TVirtualMCApplication::Instance())->Stepping();
	    fluka->SetTrackIsNew(kFALSE);
	}
	
    } else {
        //
        // Tracking is being resumed after secondary tracking
        //
        if (verbosityLevel >= 3) {
            cout << endl << "mgdraw: resuming Stepping(): " << trackId << endl;
        }

        fluka->SetTrackIsNew(kTRUE);
        fluka->SetCaller(kMGResumedTrack);
        (TVirtualMCApplication::Instance())->Stepping();

        // Reset flag and stored values
        TRACKR.ispusr[mkbmx2 - 2] = 0;
        for (Int_t i = 0; i < 9; i++) TRACKR.spausr[i] = -1.;


        if (verbosityLevel >= 3) {
            cout << endl << " !!! I am in mgdraw - first Stepping() after resume: " << icode << endl;
            cout << " Track= " << trackId << " region = " << mreg << endl;
        }

        fluka->SetTrackIsNew(kFALSE);
        fluka->SetCaller(kMGDRAW);
	if (msd == 0) {
	    (TVirtualMCApplication::Instance())->Stepping();
	} else {
	    Int_t nprim = ALLDLT.nalldl;
// Protection against nprim > mxalld
	    if (nprim >= mxalld) {
		nprim = mxalld;
		Warning("mgdraw", "nprim > mxalld, nprim: %6d  pdg: %6d mreg %6d p %13.3f step %13.3f\n", 
			ALLDLT.nalldl, 
			fluka->PDGFromId(TRACKR.jtrack), 
			mreg, 
			TRACKR.ptrack, 
			TRACKR.ctrack);
	    }
	    fluka->PrimaryIonisationStepping(nprim);
	} // primary ionisation switched on
    } // tracking resumed
} // end of mgdraw
} // end of extern "C"
Commit	Line	Data
fa3d1cc7	1	#include <Riostream.h>
8fd1d27e	2	#include "TVirtualMCApplication.h"
b1b2005d	3	#include "TVirtualMCStack.h"
fa3d1cc7	4	#include "TFluka.h"
d566901f	5	#include "TFlukaCodes.h"
efd8559a	6
	7	// Fluka includes
	8	#include "Fdimpar.h" //(DIMPAR)
	9	#include "Fdblprc.h" //(DBLPRC)
	10	#include "Ftrackr.h" //(TRACKR)
	11	#include "Fopphst.h" //(OPPHST)
	12	#include "Fflkstk.h" //(FLKSTK)
	13	#include "Fltclcm.h" //(LTCLCM)
	14	#include "Fpaprop.h" //(PAPROP)
	15	#include "Falldlt.h" //(ALLDLT)
	16	#include "Fdpdxcm.h" //(DPDXCM)
	17	#include "Fflkmat.h" //(FLKMAT)
	18
	19	#include "TGeoManager.h"
b1b2005d	20
fa3d1cc7	21	#ifndef WIN32
	22	# define mgdraw mgdraw_
	23	#else
	24	# define mgdraw MGDRAW
	25	#endif
	26
	27	extern "C" {
	28	void mgdraw(Int_t& icode, Int_t& mreg)
	29	{
b1b2005d	30	TFluka* fluka = (TFluka*) gMC;
d566901f	31	if (mreg == fluka->GetDummyRegion()) return;
b1b2005d	32	//
b1b2005d	33	// Make sure that stack has currrent track Id
57dc5a4a	34	//
3a625972	35	Int_t trackId = -1;
b1b2005d	36	TVirtualMCStack* cppstack = fluka->GetStack();
3a625972	37
3a625972	38	if (TRACKR.jtrack == -1) {
18e0cabb	39	trackId = OPPHST.louopp[OPPHST.lstopp];
	40	if (trackId == 0) {
	41	trackId = FLKSTK.ispark[FLKSTK.npflka][mkbmx2-1];
	42	}
3a625972	43	} else {
18e0cabb	44	trackId = TRACKR.ispusr[mkbmx2-1];
	45	}
	46
	47	Int_t verbosityLevel = fluka->GetVerbosityLevel();
18e0cabb	48	if (TRACKR.jtrack < -6) {
4aba9d66	49	// from -7 to -12 = "heavy" fragment
4aba9d66	50	// assing parent id
18e0cabb	51	// id < -6 was skipped in stuprf => if (kpart < -6) return;
18e0cabb	52	if (verbosityLevel >= 3) {
4aba9d66	53	cout << "mgdraw: (heavy fragment) jtrack < -6 =" << TRACKR.jtrack
18e0cabb	54	<< " assign parent pdg=" << fluka->PDGFromId(TRACKR.ispusr[mkbmx2 - 3]) << endl;
18e0cabb	55	}
ca01d0af	56	// TRACKR.jtrack = TRACKR.ispusr[mkbmx2 - 3];
3a625972	57	}
3a625972	58
b1b2005d	59	cppstack->SetCurrentTrack(trackId);
	60	//
	61	//
4aba9d66	62	Int_t mlttc = TRACKR.lt1trk; // LTCLCM.mlatm1;
d566901f	63	fluka->SetMreg(mreg, mlttc);
d566901f	64	fluka->SetIcode((FlukaProcessCode_t) icode);
d566901f	65	fluka->SetCaller(kMGDRAW);
57dc5a4a	66
4aba9d66	67	Int_t nodeId;
875111c2	68	Int_t volId = fluka->CurrentVolID(nodeId);
4aba9d66	69	Int_t crtlttc = gGeoManager->GetCurrentNodeId()+1;
4aba9d66	70
4aba9d66	71	// check region lattice consistency (debug Ernesto)
4aba9d66	72	// *****************************************************
4b1c3976	73	if(verbosityLevel>=3 && mreg != volId && !gGeoManager->IsOutside() ) {
4aba9d66	74	cout << " mgdraw: track=" << trackId << " pdg=" << fluka->PDGFromId(TRACKR.jtrack)
	75	<< " icode=" << icode << " gNstep=" << fluka->GetNstep() << endl
	76	<< " fluka mreg=" << mreg << " mlttc=" << mlttc << endl
	77	<< " TGeo volId=" << volId << " crtlttc=" << crtlttc << endl
	78	<< " common TRACKR lt1trk=" << TRACKR.lt1trk << " lt2trk=" << TRACKR.lt2trk << endl
	79	<< " common LTCLCM newlat=" << LTCLCM.newlat << " mlatld=" << LTCLCM.mlatld << endl
	80	<< " mlatm1=" << LTCLCM.mlatm1 << " mltsen=" << LTCLCM.mltsen << endl
	81	<< " mltsm1=" << LTCLCM.mltsm1 << " mlattc=" << LTCLCM.mlattc << endl;
	82	if( mlttc == crtlttc ) cout << " *************************************************************" << endl;
	83	}
	84	// *****************************************************
	85
efd8559a	86
	87	Int_t med = FLKMAT.medium[mreg - 1]; // Medium
	88	Int_t msd = DPDXCM.msdpdx[med - 1]; // Iionisation model
	89
57dc5a4a	90	if (!TRACKR.ispusr[mkbmx2 - 2]) {
875111c2	91	if (verbosityLevel >= 3) {
	92	cout << endl << "mgdraw: energy deposition for:" << trackId
	93	<< " icode=" << icode
	94	<< " pdg=" << fluka->PDGFromId(TRACKR.jtrack)
	95	<< " flukaid="<< TRACKR.jtrack
	96	<< " mreg=" << mreg
	97	<< " np =" << ALLDLT.nalldl
	98	<< endl;
	99
	100	}
efd8559a	101
	102	if (msd > 0) {
	103	//
	104	// Primary ionsiations
	105	Int_t nprim = ALLDLT.nalldl;
9c1667f7	106	if (nprim >= mxalld) {
	107	nprim = mxalld;
	108	Warning("mgdraw", "nprim > mxalld, nprim: %6d pdg: %6d mreg %6d p %13.3f step %13.3f\n",
ea262cc6	109	ALLDLT.nalldl,
	110	fluka->PDGFromId(TRACKR.jtrack),
	111	mreg,
	112	TRACKR.ptrack,
	113	TRACKR.ctrack);
875111c2	114	}
ea262cc6	115	fluka->PrimaryIonisationStepping(nprim);
875111c2	116	} else {
	117	// Single step
	118	(TVirtualMCApplication::Instance())->Stepping();
	119	fluka->SetTrackIsNew(kFALSE);
	120	}
	121
57dc5a4a	122	} else {
18e0cabb	123	//
	124	// Tracking is being resumed after secondary tracking
	125	//
	126	if (verbosityLevel >= 3) {
	127	cout << endl << "mgdraw: resuming Stepping(): " << trackId << endl;
	128	}
d566901f	129
18e0cabb	130	fluka->SetTrackIsNew(kTRUE);
	131	fluka->SetCaller(kMGResumedTrack);
	132	(TVirtualMCApplication::Instance())->Stepping();
5d80a015	133
18e0cabb	134	// Reset flag and stored values
	135	TRACKR.ispusr[mkbmx2 - 2] = 0;
	136	for (Int_t i = 0; i < 9; i++) TRACKR.spausr[i] = -1.;
5d80a015	137
57dc5a4a	138
18e0cabb	139	if (verbosityLevel >= 3) {
	140	cout << endl << " !!! I am in mgdraw - first Stepping() after resume: " << icode << endl;
	141	cout << " Track= " << trackId << " region = " << mreg << endl;
	142	}
5d80a015	143
18e0cabb	144	fluka->SetTrackIsNew(kFALSE);
18e0cabb	145	fluka->SetCaller(kMGDRAW);
ea262cc6	146	if (msd == 0) {
	147	(TVirtualMCApplication::Instance())->Stepping();
	148	} else {
	149	Int_t nprim = ALLDLT.nalldl;
	150	// Protection against nprim > mxalld
	151	if (nprim >= mxalld) {
	152	nprim = mxalld;
	153	Warning("mgdraw", "nprim > mxalld, nprim: %6d pdg: %6d mreg %6d p %13.3f step %13.3f\n",
	154	ALLDLT.nalldl,
	155	fluka->PDGFromId(TRACKR.jtrack),
	156	mreg,
	157	TRACKR.ptrack,
	158	TRACKR.ctrack);
	159	}
	160	fluka->PrimaryIonisationStepping(nprim);
	161	} // primary ionisation switched on
	162	} // tracking resumed
fa3d1cc7	163	} // end of mgdraw
	164	} // end of extern "C"
	165