[u/mrichter/AliRoot.git] / STEER / AliLego.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

/* $Id$ */

//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////
//                            
//  Utility class to evaluate the material budget from
//  a given radius to the surface of an arbitrary cylinder
//  along radial directions from the centre:
// 
//   - radiation length
//   - Interaction length
//   - g/cm2
// 
//  Geantinos are shot in the bins in the fNtheta bins in theta
//  and fNphi bins in phi with specified rectangular limits.
//  The statistics are accumulated per
//    fRadMin < r < fRadMax    and  <0 < z < fZMax
//
//  To activate this option, you can do:
//      Root > gAlice.RunLego();
//  or  Root > .x menu.C  then select menu item "RunLego"
//  Note that when calling gAlice->RunLego, an optional list
//  of arguments may be specified.
//
//Begin_Html
/*
<img src="picts/alilego.gif">
*/
//End_Html
//
//////////////////////////////////////////////////////////////

#include "TClonesArray.h"
#include "TH2.h"
#include "TMath.h"
#include "TString.h"
#include "TVirtualMC.h"

#include "AliDebugVolume.h"
#include "AliLego.h"
#include "AliLegoGenerator.h"
#include "AliMC.h"
#include "AliRun.h"

ClassImp(AliLego)

//_______________________________________________________________________
AliLego::AliLego():
  fGener(0),
  fTotRadl(0),
  fTotAbso(0),
  fTotGcm2(0),
  fHistRadl(0),
  fHistAbso(0),
  fHistGcm2(0),
  fHistReta(0),
  fVolumesFwd(0),
  fVolumesBwd(0),
  fStepBack(0),
  fStepsBackward(0),
  fStepsForward(0),
  fErrorCondition(0),
  fDebug(0)
{
  //
  // Default constructor
  //
}

//_______________________________________________________________________
AliLego::AliLego(const AliLego &lego):
  TNamed(lego),
  fGener(0),
  fTotRadl(0),
  fTotAbso(0),
  fTotGcm2(0),
  fHistRadl(0),
  fHistAbso(0),
  fHistGcm2(0),
  fHistReta(0),
  fVolumesFwd(0),
  fVolumesBwd(0),
  fStepBack(0),
  fStepsBackward(0),
  fStepsForward(0),
  fErrorCondition(0),
  fDebug(0)
{
  //
  // Copy constructor
  //
  lego.Copy(*this);
}


//_______________________________________________________________________
AliLego::AliLego(const char *title, Int_t ntheta, Float_t thetamin, 
                 Float_t thetamax, Int_t nphi, Float_t phimin, Float_t phimax,
                 Float_t rmin, Float_t rmax, Float_t zmax):
  TNamed("Lego Generator",title),
  fGener(0),
  fTotRadl(0),
  fTotAbso(0),
  fTotGcm2(0),
  fHistRadl(0),
  fHistAbso(0),
  fHistGcm2(0),
  fHistReta(0),
  fVolumesFwd(0),
  fVolumesBwd(0),
  fStepBack(0),
  fStepsBackward(0),
  fStepsForward(0),
  fErrorCondition(0),
  fDebug(0)
{
  //
  // specify the angular limits and the size of the rectangular box
  //
   fGener = new AliLegoGenerator(ntheta, thetamin, thetamax,
		       nphi, phimin, phimax, rmin, rmax, zmax);
   fHistRadl = new TH2F("hradl","Radiation length map",    
			ntheta,thetamin,thetamax,nphi,phimin,phimax);
   fHistAbso = new TH2F("habso","Interaction length map",  
			ntheta,thetamin,thetamax,nphi,phimin,phimax);
   fHistGcm2 = new TH2F("hgcm2","g/cm2 length map",        
			ntheta,thetamin,thetamax,nphi,phimin,phimax);
//
   fVolumesFwd     = new TClonesArray("AliDebugVolume",1000);
   fVolumesBwd     = new TClonesArray("AliDebugVolume",1000);
   fDebug          = gAlice->GetDebug();
}

//_______________________________________________________________________
AliLego::AliLego(const char *title, AliLegoGenerator* generator):
  TNamed("Lego Generator",title),
  fGener(0),
  fTotRadl(0),
  fTotAbso(0),
  fTotGcm2(0),
  fHistRadl(0),
  fHistAbso(0),
  fHistGcm2(0),
  fHistReta(0),
  fVolumesFwd(0),
  fVolumesBwd(0),
  fStepBack(0),
  fStepsBackward(0),
  fStepsForward(0),
  fErrorCondition(0),
  fDebug(0)
{
  //
  // specify the angular limits and the size of the rectangular box
  //
  fGener = generator;
  Float_t c1min, c1max, c2min, c2max;
  Int_t n1 = fGener->NCoor1();
  Int_t n2 = fGener->NCoor2();
  
  fGener->Coor1Range(c1min, c1max);
  fGener->Coor2Range(c2min, c2max);   
  
  fHistRadl = new TH2F("hradl","Radiation length map",    
                       n2, c2min, c2max, n1, c1min, c1max);
  fHistAbso = new TH2F("habso","Interaction length map",  
                       n2, c2min, c2max, n1, c1min, c1max);
  fHistGcm2 = new TH2F("hgcm2","g/cm2 length map",        
                       n2, c2min, c2max, n1, c1min, c1max);
  //
  //

  fVolumesFwd     = new TClonesArray("AliDebugVolume",1000);
  fVolumesBwd     = new TClonesArray("AliDebugVolume",1000);
  fDebug          = gAlice->GetDebug();
}

//_______________________________________________________________________
AliLego::~AliLego()
{
  //
  // Destructor
  //
  delete fHistRadl;
  delete fHistAbso;
  delete fHistGcm2;
  delete fGener;
  delete fVolumesFwd;
  delete fVolumesBwd;
}

//_______________________________________________________________________
void AliLego::BeginEvent()
{
  //
  // --- Set to 0 radiation length, absorption length and g/cm2 ---
  //
  fTotRadl = 0;
  fTotAbso = 0;
  fTotGcm2 = 0;

  if (fDebug) {
    if (fErrorCondition) DumpVolumes();
    fVolumesFwd->Delete();
    fVolumesBwd->Delete();
    fStepsForward    = 0;
    fStepsBackward   = 0;		  
    fErrorCondition  = 0;
    if (gAlice->GetMCApp()->GetCurrentTrackNumber() == 0) fStepBack = 0;
  }
}

//_______________________________________________________________________
void AliLego::FinishEvent()
{
  //
  // Finish the event and update the histos
  //
  Double_t c1, c2;
  c1 = fGener->CurCoor1();
  c2 = fGener->CurCoor2();
  fHistRadl->Fill(c2,c1,fTotRadl);
  fHistAbso->Fill(c2,c1,fTotAbso);
  fHistGcm2->Fill(c2,c1,fTotGcm2);
}

//_______________________________________________________________________
void AliLego::FinishRun()
{
  //
  // Store histograms in current Root file
  //
  fHistRadl->Write();
  fHistAbso->Write();
  fHistGcm2->Write();

  // Delete histograms from memory
  fHistRadl->Delete(); fHistRadl=0;
  fHistAbso->Delete(); fHistAbso=0;
  fHistGcm2->Delete(); fHistGcm2=0;
  //
  if (fErrorCondition) DumpVolumes();
}

//_______________________________________________________________________
void AliLego::Copy(TObject&) const
{
  //
  // Copy *this onto lego -- not implemented
  //
  Fatal("Copy","Not implemented!\n");
}

//_______________________________________________________________________
void AliLego::StepManager() 
{
  //
  // called from AliRun::Stepmanager from gustep.
  // Accumulate the 3 parameters step by step
  //
  static Float_t t;
  Float_t a,z,dens,radl,absl;
  Int_t i, id, copy;
  const char* vol;
  static Float_t vect[3], dir[3];
  
  TString tmp1, tmp2;
  copy = 1;
  id  = gMC->CurrentVolID(copy);
  vol = gMC->VolName(id);
  Float_t step  = gMC->TrackStep();
  
  TLorentzVector pos, mom; 
  gMC->TrackPosition(pos);  
  gMC->TrackMomentum(mom);
  
  Int_t status = 0;
  if (gMC->IsTrackEntering()) status = 1;
  if (gMC->IsTrackExiting())  status = 2; 
  
  if (! fStepBack) {
    //      printf("\n volume %s %d", vol, status);      
    //
    // Normal Forward stepping
    //
    if (fDebug) {
      //	  printf("\n steps fwd %d %s %d %f", fStepsForward, vol, fErrorCondition, step);	  
      
      //
      // store volume if different from previous
      //
	  
	  TClonesArray &lvols = *fVolumesFwd;
	  if (fStepsForward > 0) {
        AliDebugVolume* tmp = dynamic_cast<AliDebugVolume*>((*fVolumesFwd)[fStepsForward-1]);
        if (tmp->IsVEqual(vol, copy) && gMC->IsTrackEntering()) {
		  fStepsForward -= 2;
		  fVolumesFwd->RemoveAt(fStepsForward);
		  fVolumesFwd->RemoveAt(fStepsForward+1);		  
        }
	  }
      
	  new(lvols[fStepsForward++]) 
        AliDebugVolume(vol,copy,step,pos[0], pos[1], pos[2], status);
      
    } // Debug
    //
    // Get current material properties
    
    gMC->CurrentMaterial(a,z,dens,radl,absl);
    
    if (z < 1) return;
    
    // --- See if we have to stop now
    if (TMath::Abs(pos[2]) > fGener->ZMax()  || 
        pos[0]*pos[0] +pos[1]*pos[1] > fGener->RadMax()*fGener->RadMax()) {
      if (!gMC->IsNewTrack()) {
        // Not the first step, add past contribution
        fTotAbso += t/absl;
        fTotRadl += t/radl;
        fTotGcm2 += t*dens;
        if (fDebug) {
          //
          //  generate "mirror" particle flying back
          //
          fStepsBackward = fStepsForward;
          
          Float_t pmom[3], orig[3];
          Float_t polar[3] = {0.,0.,0.};
          Int_t ntr;
          pmom[0] = -dir[0];
          pmom[1] = -dir[1];	   
          pmom[2] = -dir[2];
          orig[0] =  vect[0];
          orig[1] =  vect[1];	   
          orig[2] =  vect[2];
          
          gAlice->GetMCApp()->PushTrack(1, gAlice->GetMCApp()->GetCurrentTrackNumber(), 
                           0, pmom, orig, polar, 0., kPNoProcess, ntr);
        } // debug
        
      } // not a new track !
      
      if (fDebug) fStepBack = 1;
      gMC->StopTrack();
      return;
    } // outside scoring region ?
    
    // --- See how long we have to go
    for(i=0;i<3;++i) {
      vect[i]=pos[i];
      dir[i]=mom[i];
    }
    
    t  = fGener->PropagateCylinder(vect,dir,fGener->RadMax(),fGener->ZMax());
    
    if(step) {
      
      fTotAbso += step/absl;
      fTotRadl += step/radl;
      fTotGcm2 += step*dens;
    }
    
  } else {
      if (fDebug) {
        //
        // Geometry debugging
        // Fly back and compare volume sequence
        //
        TClonesArray &lvols = *fVolumesBwd;
        if (fStepsBackward < fStepsForward) {
	      AliDebugVolume* tmp = dynamic_cast<AliDebugVolume*>((*fVolumesBwd)[fStepsBackward]);
	      if (tmp->IsVEqual(vol, copy) && gMC->IsTrackEntering()) {
            fStepsBackward += 2;
            fVolumesBwd->RemoveAt(fStepsBackward-1);
            fVolumesBwd->RemoveAt(fStepsBackward-2);		  
	      }
        } 
        
        fStepsBackward--;
        //	  printf("\n steps %d %s %d", fStepsBackward, vol, fErrorCondition);	  
        if (fStepsBackward < 0) {
	      gMC->StopTrack();
	      fStepBack = 0;
	      return;
        }
        
        new(lvols[fStepsBackward]) AliDebugVolume(vol,copy,step,pos[0], pos[1], pos[2], status);
        
        AliDebugVolume* tmp = dynamic_cast<AliDebugVolume*>((*fVolumesFwd)[fStepsBackward]);
        if (! (tmp->IsVEqual(vol, copy)) && (!fErrorCondition)) 
          {
            printf("\n Problem at (x,y,z): %d %f %f %f, volumes: %s %s step: %f\n", 
                   fStepsBackward, pos[0], pos[1], pos[2], tmp->GetName(), vol, step);
            fErrorCondition = 1;
          } 
      } // Debug
  } // bwd/fwd
}

//_______________________________________________________________________
void AliLego::DumpVolumes()
{
  //
  // Dump volume sequence in case of error
  //
  printf("\n Dumping Volume Sequence:");
  printf("\n ==============================================");
  
  for (Int_t i = fStepsForward-1; i>=0; i--)
    {
      AliDebugVolume* tmp1 = dynamic_cast<AliDebugVolume*>((*fVolumesFwd)[i]);
      AliDebugVolume* tmp2 = dynamic_cast<AliDebugVolume*>((*fVolumesBwd)[i]);
      if (tmp1)
        printf("\n Volume Fwd: %3d: %5s (%3d) step: %12.5e (x,y,z) (%12.5e %12.5e %12.5e) status: %9s \n"
               , i, 
               tmp1->GetName(), tmp1->CopyNumber(), tmp1->Step(), 
               tmp1->X(), tmp1->Y(), tmp1->Z(), tmp1->Status());
      if (tmp2 && i>= fStepsBackward)
        printf("\n Volume Bwd: %3d: %5s (%3d) step: %12.5e (x,y,z) (%12.5e %12.5e %12.5e) status: %9s \n"
               , i, 
               tmp2->GetName(), tmp2->CopyNumber(), tmp2->Step(), 
               tmp2->X(), tmp2->Y(), tmp2->Z(), tmp2->Status());
      
      printf("\n ............................................................................\n");
    }
  printf("\n ==============================================\n");
}
Commit	Line	Data
4c039060	1	/**************************************************************************
	2	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercial purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
	15
acd84897	16	/* $Id$ */
4c039060	17
fe4da5cc	18	//////////////////////////////////////////////////////////////
	19	//////////////////////////////////////////////////////////////
	20	//
	21	// Utility class to evaluate the material budget from
	22	// a given radius to the surface of an arbitrary cylinder
	23	// along radial directions from the centre:
	24	//
	25	// - radiation length
	26	// - Interaction length
	27	// - g/cm2
	28	//
	29	// Geantinos are shot in the bins in the fNtheta bins in theta
	30	// and fNphi bins in phi with specified rectangular limits.
	31	// The statistics are accumulated per
	32	// fRadMin < r < fRadMax and <0 < z < fZMax
	33	//
	34	// To activate this option, you can do:
	35	// Root > gAlice.RunLego();
	36	// or Root > .x menu.C then select menu item "RunLego"
	37	// Note that when calling gAlice->RunLego, an optional list
	38	// of arguments may be specified.
	39	//
	40	//Begin_Html
	41	/*
1439f98e	42	<img src="picts/alilego.gif">
fe4da5cc	43	*/
	44	//End_Html
	45	//
	46	//////////////////////////////////////////////////////////////
	47
e2afb3b6	48	#include "TClonesArray.h"
e2afb3b6	49	#include "TH2.h"
fe4da5cc	50	#include "TMath.h"
e2afb3b6	51	#include "TString.h"
5d12ce38	52	#include "TVirtualMC.h"
94de3818	53
7eb618fa	54	#include "AliDebugVolume.h"
e2afb3b6	55	#include "AliLego.h"
8918e700	56	#include "AliLegoGenerator.h"
5d12ce38	57	#include "AliMC.h"
0742d588	58	#include "AliRun.h"
7eb618fa	59
fe4da5cc	60	ClassImp(AliLego)
fe4da5cc	61
e2afb3b6	62	//_______________________________________________________________________
	63	AliLego::AliLego():
	64	fGener(0),
	65	fTotRadl(0),
	66	fTotAbso(0),
	67	fTotGcm2(0),
	68	fHistRadl(0),
	69	fHistAbso(0),
	70	fHistGcm2(0),
	71	fHistReta(0),
	72	fVolumesFwd(0),
	73	fVolumesBwd(0),
	74	fStepBack(0),
	75	fStepsBackward(0),
	76	fStepsForward(0),
	77	fErrorCondition(0),
	78	fDebug(0)
fe4da5cc	79	{
8918e700	80	//
	81	// Default constructor
	82	//
fe4da5cc	83	}
fe4da5cc	84
e2afb3b6	85	//_______________________________________________________________________
	86	AliLego::AliLego(const AliLego &lego):
	87	TNamed(lego),
	88	fGener(0),
	89	fTotRadl(0),
	90	fTotAbso(0),
	91	fTotGcm2(0),
	92	fHistRadl(0),
	93	fHistAbso(0),
	94	fHistGcm2(0),
	95	fHistReta(0),
	96	fVolumesFwd(0),
	97	fVolumesBwd(0),
	98	fStepBack(0),
	99	fStepsBackward(0),
	100	fStepsForward(0),
	101	fErrorCondition(0),
	102	fDebug(0)
	103	{
	104	//
	105	// Copy constructor
	106	//
	107	lego.Copy(*this);
	108	}
	109
	110
	111	//_______________________________________________________________________
0b27b8c7	112	AliLego::AliLego(const char *title, Int_t ntheta, Float_t thetamin,
e2afb3b6	113	Float_t thetamax, Int_t nphi, Float_t phimin, Float_t phimax,
	114	Float_t rmin, Float_t rmax, Float_t zmax):
	115	TNamed("Lego Generator",title),
	116	fGener(0),
	117	fTotRadl(0),
	118	fTotAbso(0),
	119	fTotGcm2(0),
	120	fHistRadl(0),
	121	fHistAbso(0),
	122	fHistGcm2(0),
	123	fHistReta(0),
	124	fVolumesFwd(0),
	125	fVolumesBwd(0),
	126	fStepBack(0),
	127	fStepsBackward(0),
	128	fStepsForward(0),
	129	fErrorCondition(0),
	130	fDebug(0)
fe4da5cc	131	{
8918e700	132	//
	133	// specify the angular limits and the size of the rectangular box
	134	//
0b27b8c7	135	fGener = new AliLegoGenerator(ntheta, thetamin, thetamax,
b13db077	136	nphi, phimin, phimax, rmin, rmax, zmax);
b13db077	137	fHistRadl = new TH2F("hradl","Radiation length map",
0b27b8c7	138	ntheta,thetamin,thetamax,nphi,phimin,phimax);
b13db077	139	fHistAbso = new TH2F("habso","Interaction length map",
0b27b8c7	140	ntheta,thetamin,thetamax,nphi,phimin,phimax);
b13db077	141	fHistGcm2 = new TH2F("hgcm2","g/cm2 length map",
0b27b8c7	142	ntheta,thetamin,thetamax,nphi,phimin,phimax);
7eb618fa	143	//
7eb618fa	144	fVolumesFwd = new TClonesArray("AliDebugVolume",1000);
7eb618fa	145	fVolumesBwd = new TClonesArray("AliDebugVolume",1000);
19237f5f	146	fDebug = gAlice->GetDebug();
0b27b8c7	147	}
0b27b8c7	148
e2afb3b6	149	//_______________________________________________________________________
	150	AliLego::AliLego(const char title, AliLegoGenerator generator):
	151	TNamed("Lego Generator",title),
	152	fGener(0),
	153	fTotRadl(0),
	154	fTotAbso(0),
	155	fTotGcm2(0),
	156	fHistRadl(0),
	157	fHistAbso(0),
	158	fHistGcm2(0),
	159	fHistReta(0),
	160	fVolumesFwd(0),
	161	fVolumesBwd(0),
	162	fStepBack(0),
	163	fStepsBackward(0),
	164	fStepsForward(0),
	165	fErrorCondition(0),
	166	fDebug(0)
0b27b8c7	167	{
	168	//
	169	// specify the angular limits and the size of the rectangular box
	170	//
e2afb3b6	171	fGener = generator;
	172	Float_t c1min, c1max, c2min, c2max;
	173	Int_t n1 = fGener->NCoor1();
	174	Int_t n2 = fGener->NCoor2();
	175
	176	fGener->Coor1Range(c1min, c1max);
	177	fGener->Coor2Range(c2min, c2max);
	178
	179	fHistRadl = new TH2F("hradl","Radiation length map",
	180	n2, c2min, c2max, n1, c1min, c1max);
	181	fHistAbso = new TH2F("habso","Interaction length map",
	182	n2, c2min, c2max, n1, c1min, c1max);
	183	fHistGcm2 = new TH2F("hgcm2","g/cm2 length map",
	184	n2, c2min, c2max, n1, c1min, c1max);
	185	//
	186	//
b13db077	187
e2afb3b6	188	fVolumesFwd = new TClonesArray("AliDebugVolume",1000);
	189	fVolumesBwd = new TClonesArray("AliDebugVolume",1000);
	190	fDebug = gAlice->GetDebug();
fe4da5cc	191	}
fe4da5cc	192
e2afb3b6	193	//_______________________________________________________________________
fe4da5cc	194	AliLego::~AliLego()
fe4da5cc	195	{
8918e700	196	//
	197	// Destructor
	198	//
	199	delete fHistRadl;
	200	delete fHistAbso;
	201	delete fHistGcm2;
8918e700	202	delete fGener;
e2afb3b6	203	delete fVolumesFwd;
e2afb3b6	204	delete fVolumesBwd;
fe4da5cc	205	}
fe4da5cc	206
e2afb3b6	207	//_______________________________________________________________________
dffd31ef	208	void AliLego::BeginEvent()
b13db077	209	{
8918e700	210	//
	211	// --- Set to 0 radiation length, absorption length and g/cm2 ---
	212	//
dffd31ef	213	fTotRadl = 0;
	214	fTotAbso = 0;
	215	fTotGcm2 = 0;
19237f5f	216
a74d2699	217	if (fDebug) {
e2afb3b6	218	if (fErrorCondition) DumpVolumes();
	219	fVolumesFwd->Delete();
	220	fVolumesBwd->Delete();
	221	fStepsForward = 0;
	222	fStepsBackward = 0;
	223	fErrorCondition = 0;
5d12ce38	224	if (gAlice->GetMCApp()->GetCurrentTrackNumber() == 0) fStepBack = 0;
a74d2699	225	}
dffd31ef	226	}
dffd31ef	227
e2afb3b6	228	//_______________________________________________________________________
dffd31ef	229	void AliLego::FinishEvent()
dffd31ef	230	{
8918e700	231	//
	232	// Finish the event and update the histos
	233	//
0b27b8c7	234	Double_t c1, c2;
	235	c1 = fGener->CurCoor1();
	236	c2 = fGener->CurCoor2();
10bd7535	237	fHistRadl->Fill(c2,c1,fTotRadl);
	238	fHistAbso->Fill(c2,c1,fTotAbso);
	239	fHistGcm2->Fill(c2,c1,fTotGcm2);
b13db077	240	}
b13db077	241
e2afb3b6	242	//_______________________________________________________________________
dffd31ef	243	void AliLego::FinishRun()
dffd31ef	244	{
8918e700	245	//
	246	// Store histograms in current Root file
	247	//
dffd31ef	248	fHistRadl->Write();
	249	fHistAbso->Write();
	250	fHistGcm2->Write();
dffd31ef	251
	252	// Delete histograms from memory
	253	fHistRadl->Delete(); fHistRadl=0;
	254	fHistAbso->Delete(); fHistAbso=0;
	255	fHistGcm2->Delete(); fHistGcm2=0;
a74d2699	256	//
a74d2699	257	if (fErrorCondition) DumpVolumes();
dffd31ef	258	}
dffd31ef	259
e2afb3b6	260	//_______________________________________________________________________
6c4904c2	261	void AliLego::Copy(TObject&) const
8918e700	262	{
	263	//
	264	// Copy *this onto lego -- not implemented
	265	//
	266	Fatal("Copy","Not implemented!\n");
	267	}
dffd31ef	268
e2afb3b6	269	//_______________________________________________________________________
	270	void AliLego::StepManager()
	271	{
	272	//
	273	// called from AliRun::Stepmanager from gustep.
	274	// Accumulate the 3 parameters step by step
	275	//
	276	static Float_t t;
	277	Float_t a,z,dens,radl,absl;
	278	Int_t i, id, copy;
	279	const char* vol;
	280	static Float_t vect[3], dir[3];
	281
	282	TString tmp1, tmp2;
	283	copy = 1;
	284	id = gMC->CurrentVolID(copy);
	285	vol = gMC->VolName(id);
	286	Float_t step = gMC->TrackStep();
	287
	288	TLorentzVector pos, mom;
	289	gMC->TrackPosition(pos);
	290	gMC->TrackMomentum(mom);
	291
	292	Int_t status = 0;
	293	if (gMC->IsTrackEntering()) status = 1;
	294	if (gMC->IsTrackExiting()) status = 2;
	295
7eb618fa	296	if (! fStepBack) {
e2afb3b6	297	// printf("\n volume %s %d", vol, status);
	298	//
	299	// Normal Forward stepping
	300	//
	301	if (fDebug) {
	302	// printf("\n steps fwd %d %s %d %f", fStepsForward, vol, fErrorCondition, step);
	303
	304	//
	305	// store volume if different from previous
	306	//
7eb618fa	307
a74d2699	308	TClonesArray &lvols = *fVolumesFwd;
7eb618fa	309	if (fStepsForward > 0) {
e2afb3b6	310	AliDebugVolume* tmp = dynamic_cast<AliDebugVolume>((fVolumesFwd)[fStepsForward-1]);
e2afb3b6	311	if (tmp->IsVEqual(vol, copy) && gMC->IsTrackEntering()) {
a74d2699	312	fStepsForward -= 2;
5f3dc1c9	313	fVolumesFwd->RemoveAt(fStepsForward);
5f3dc1c9	314	fVolumesFwd->RemoveAt(fStepsForward+1);
e2afb3b6	315	}
7eb618fa	316	}
e2afb3b6	317
a74d2699	318	new(lvols[fStepsForward++])
e2afb3b6	319	AliDebugVolume(vol,copy,step,pos[0], pos[1], pos[2], status);
7eb618fa	320
e2afb3b6	321	} // Debug
	322	//
	323	// Get current material properties
	324
	325	gMC->CurrentMaterial(a,z,dens,radl,absl);
	326
	327	if (z < 1) return;
	328
	329	// --- See if we have to stop now
	330	if (TMath::Abs(pos[2]) > fGener->ZMax() \|\|
	331	pos[0]pos[0] +pos[1]pos[1] > fGener->RadMax()*fGener->RadMax()) {
	332	if (!gMC->IsNewTrack()) {
	333	// Not the first step, add past contribution
	334	fTotAbso += t/absl;
	335	fTotRadl += t/radl;
	336	fTotGcm2 += t*dens;
	337	if (fDebug) {
	338	//
	339	// generate "mirror" particle flying back
	340	//
	341	fStepsBackward = fStepsForward;
	342
	343	Float_t pmom[3], orig[3];
	344	Float_t polar[3] = {0.,0.,0.};
	345	Int_t ntr;
	346	pmom[0] = -dir[0];
	347	pmom[1] = -dir[1];
	348	pmom[2] = -dir[2];
	349	orig[0] = vect[0];
	350	orig[1] = vect[1];
	351	orig[2] = vect[2];
	352
5d12ce38	353	gAlice->GetMCApp()->PushTrack(1, gAlice->GetMCApp()->GetCurrentTrackNumber(),
e2afb3b6	354	0, pmom, orig, polar, 0., kPNoProcess, ntr);
	355	} // debug
	356
	357	} // not a new track !
7eb618fa	358
e2afb3b6	359	if (fDebug) fStepBack = 1;
	360	gMC->StopTrack();
	361	return;
	362	} // outside scoring region ?
	363
	364	// --- See how long we have to go
	365	for(i=0;i<3;++i) {
	366	vect[i]=pos[i];
	367	dir[i]=mom[i];
	368	}
	369
	370	t = fGener->PropagateCylinder(vect,dir,fGener->RadMax(),fGener->ZMax());
	371
	372	if(step) {
	373
	374	fTotAbso += step/absl;
	375	fTotRadl += step/radl;
	376	fTotGcm2 += step*dens;
	377	}
	378
7eb618fa	379	} else {
7eb618fa	380	if (fDebug) {
e2afb3b6	381	//
	382	// Geometry debugging
	383	// Fly back and compare volume sequence
	384	//
	385	TClonesArray &lvols = *fVolumesBwd;
	386	if (fStepsBackward < fStepsForward) {
	387	AliDebugVolume* tmp = dynamic_cast<AliDebugVolume>((fVolumesBwd)[fStepsBackward]);
	388	if (tmp->IsVEqual(vol, copy) && gMC->IsTrackEntering()) {
	389	fStepsBackward += 2;
	390	fVolumesBwd->RemoveAt(fStepsBackward-1);
	391	fVolumesBwd->RemoveAt(fStepsBackward-2);
7eb618fa	392	}
e2afb3b6	393	}
	394
	395	fStepsBackward--;
	396	// printf("\n steps %d %s %d", fStepsBackward, vol, fErrorCondition);
	397	if (fStepsBackward < 0) {
7eb618fa	398	gMC->StopTrack();
a74d2699	399	fStepBack = 0;
7eb618fa	400	return;
e2afb3b6	401	}
	402
	403	new(lvols[fStepsBackward]) AliDebugVolume(vol,copy,step,pos[0], pos[1], pos[2], status);
	404
	405	AliDebugVolume* tmp = dynamic_cast<AliDebugVolume>((fVolumesFwd)[fStepsBackward]);
	406	if (! (tmp->IsVEqual(vol, copy)) && (!fErrorCondition))
	407	{
	408	printf("\n Problem at (x,y,z): %d %f %f %f, volumes: %s %s step: %f\n",
	409	fStepsBackward, pos[0], pos[1], pos[2], tmp->GetName(), vol, step);
	410	fErrorCondition = 1;
	411	}
7eb618fa	412	} // Debug
7eb618fa	413	} // bwd/fwd
b13db077	414	}
b13db077	415
e2afb3b6	416	//_______________________________________________________________________
a74d2699	417	void AliLego::DumpVolumes()
a74d2699	418	{
e2afb3b6	419	//
	420	// Dump volume sequence in case of error
	421	//
	422	printf("\n Dumping Volume Sequence:");
	423	printf("\n ==============================================");
	424
	425	for (Int_t i = fStepsForward-1; i>=0; i--)
a74d2699	426	{
e2afb3b6	427	AliDebugVolume* tmp1 = dynamic_cast<AliDebugVolume>((fVolumesFwd)[i]);
	428	AliDebugVolume* tmp2 = dynamic_cast<AliDebugVolume>((fVolumesBwd)[i]);
	429	if (tmp1)
	430	printf("\n Volume Fwd: %3d: %5s (%3d) step: %12.5e (x,y,z) (%12.5e %12.5e %12.5e) status: %9s \n"
	431	, i,
	432	tmp1->GetName(), tmp1->CopyNumber(), tmp1->Step(),
	433	tmp1->X(), tmp1->Y(), tmp1->Z(), tmp1->Status());
	434	if (tmp2 && i>= fStepsBackward)
	435	printf("\n Volume Bwd: %3d: %5s (%3d) step: %12.5e (x,y,z) (%12.5e %12.5e %12.5e) status: %9s \n"
	436	, i,
	437	tmp2->GetName(), tmp2->CopyNumber(), tmp2->Step(),
	438	tmp2->X(), tmp2->Y(), tmp2->Z(), tmp2->Status());
	439
	440	printf("\n ............................................................................\n");
a74d2699	441	}
e2afb3b6	442	printf("\n ==============================================\n");
a74d2699	443	}