TpcSimIonisationChamber.cxx

Go to the documentation of this file.

#include "TMath.h"
#include "TRandom.h" 

#include "HarpoConfig.h"
#include "HarpoDebug.h"
#include "TpcSimMCEvent.h"
#include "TpcSimIonisationTrack.h"
#include "TpcSimIonisationChamber.h"


#ifdef HARPO_G4
#include "DetectorConstruction.hh"
#include "PhysicsList.hh"
#include "PrimaryGeneratorAction.hh"
#include "RunAction.hh"
#include "EventAction.hh"
#endif

#include <iostream>
#include <fstream>

//using namespace std;

ClassImp(TpcSimIonisationChamber)

TpcSimIonisationChamber::TpcSimIonisationChamber()
{

  fGeom = new TpcSimGeometry();

  fIpot = 15;
  fW = 25;
  fFano = 0;

  fMfp = 1./3.641661E+01;   // MFP [1/cm]

  Double_t ipot;
  if ( ! gHConfig->Lookup("Sim.Ipot",ipot) )
    Info("Constructor","Use default Ipot %g",fIpot);
  else
    fIpot = ipot;
 
  Double_t w;
  if ( ! gHConfig->Lookup("Sim.W",w) )
    Info("Constructor","Use default W %g",fW);
  else
    fW = w;
 
  Double_t fano;
  if ( ! gHConfig->Lookup("Sim.Fano",fano) )
    Info("Constructor","Use default Fano %g",fFano);
  else
    fFano = fano;
 
  Double_t mfp;
  if ( ! gHConfig->Lookup("Sim.Mfp",mfp) )
    Info("Constructor","Use default Mfp %g",fMfp);
  else
    fMfp = mfp;

  fModuleX = 0;
  fModuleY = 0;
  fModuleZ = 0;

  Long64_t cx;
  if ( ! gHConfig->Lookup("Sim.ModuleX",cx) )
    Info("Constructor","Use default ModuleX %d",fModuleX);
  else
    fModuleX = cx;
  Long64_t cy;

  if ( ! gHConfig->Lookup("Sim.ModuleY",cy) )
    Info("Constructor","Use default ModuleY %d",fModuleY);
  else
    fModuleY = cy;

  Long64_t cz;
  if ( ! gHConfig->Lookup("Sim.ModuleZ",cz) )
    Info("Constructor","Use default ModuleZ %d",fModuleZ);
  else
    fModuleZ = cz;
 

}

TpcSimIonisationChamber::~TpcSimIonisationChamber()
{
   if(gHarpoDebug>0) 
     Info("~TpcSimIonisationChamber","Deleting");
  fGeom->Delete();
}


TpcSimIonisationTrack* TpcSimIonisationChamber::GenerateTrack(TpcSimMCTrack* tr)
{
  if(gHarpoDebug>1) 
    Info("GenerateTrack","start");

  TpcSimIonisationTrack* trackout = new TpcSimIonisationTrack();

  Double_t x = tr->GetX0(), y = tr->GetY0(), z = tr->GetZ0(), t = tr->GetT0();
  Double_t px = tr->GetPx(), py = tr->GetPy(), pz = tr->GetPz();
  Double_t p = TMath::Sqrt(px*px+py*py+pz*pz);
  Double_t d = 0;
  Double_t dMax = 10000;
  if(gHarpoDebug>1) 
    Info("GenerateTrack","%f %f %f %f   %f %f %f",x,y,z,t,px,py,pz);
  Bool_t notyet = kTRUE;
  Int_t isin = 0;
  while((notyet || isin==1) && d<dMax){
    isin = fGeom->GetVolume(x,y,z);
    if(isin==1) notyet = kFALSE;

    Double_t l = fMfp*gRandom->Exp(1);                         // *
    Double_t dx = l*px/p;
    Double_t dy = l*py/p;
    Double_t dz = l*pz/p;

    x += dx; y += dy; z += dz; d += l;

    Double_t eloss = GetEloss();
    Int_t Ne = GetNe(eloss);

    TpcSimIonisationPoint* point = new TpcSimIonisationPoint(Ne,x,y,z,t,eloss,0);
    trackout->AddPoint(point);
  }

  if(gHarpoDebug>1) 
    Info("GenerateTrack","End track %d %d %f",notyet,isin,d);
 
  return trackout;

}

//TpcSimIonisationTrack* TpcSimIonisationChamber::GenerateTrackFromText(ifstream *str, TpcSimMCEvent* mcevent)
TpcSimIonisationTrack* TpcSimIonisationChamber::GenerateTrackFromText(std::ifstream *str)
{
  if(gHarpoDebug>1) 
    Info("GenerateTrackFromText","start");


  if(!str->good()) 
    return 0;
  Char_t line[256];
  str->getline(line,256);
  if(strstr(line,"NEW TRACK") == NULL)
    return 0;

  Float_t x0 = 0, y0 = 0, z0 = 0, t0 = 0;
  Float_t px = 0, py = 0, pz = 0, e = 0;
  Int_t pid = 0, trackId = 0;
  /*  Int_t test = */ sscanf(line,"  NEW TRACK %f %f %f %f %f %f %f %f %d %d",
         &x0, &y0, &z0, &t0, 
         &px, &py, &pz, &e, 
         &pid, &trackId);

  if(gHarpoDebug>1) Info("GenerateTrackFromText"," %s",line);
  // Info("GenerateTrackFromText","%d",test);
  // Info("GenerateTrackFromText","%g %g %g %g %g %g %g %g %d %d",x0, y0, z0, t0, px, py, pz, e,pid, trackId);

  // TpcSimMCTrack* mctrack 
  //   = new TpcSimMCTrack(Double_t(x0), Double_t(y0), Double_t(z0), Double_t(t0),
  //                            Double_t(px), Double_t(py), Double_t(pz),
  //                            pid);
  // mcevent->AddTrack(mctrack);

  TpcSimIonisationTrack* trackout 
    = new TpcSimIonisationTrack(Double_t(x0), Double_t(y0), Double_t(z0), Double_t(t0),
                                Double_t(px), Double_t(py), Double_t(pz),
                                pid,
                                0, 0, trackId);
  //  TpcSimIonisationTrack* trackout = new TpcSimIonisationTrack();

  Double_t t = 2;

  while(str->good()){
    str->getline(line,256);
    if(strstr(line,"END TRACK") != NULL)
      break;
    // if(strstr(line,"BEGINNING OF EVENT") != NULL){
    //   Float_t x0, y0, z0, t0;
    //   Float_t px, py, pz, e;
    //   Int_t pid, trackId;
    //   sscanf(line,"BEGINNING OF EVENT %f %f %f %f %f %f %f %f %d %d",
    //       &x0, &y0, &z0, &t0, 
    //       &px, &py, &pz, &e, 
    //       &pid, &trackId);
    //   continue;
    // }
    Float_t x = 10, y = 10, z = 10;
    Int_t Ne;
    //    Int_t test = sscanf(line,"%f %f %f %d",&x, &y, &z, &Ne);
    Int_t cx = -1, cy = -1, cz = -1, trk = -1;
    Float_t time;
    Int_t test = sscanf(line,"%f %f %f %d %d %d %d %d %f",&x, &y, &z, &Ne, &cx, &cy, &cz, &trk, &time);
    z += fGeom->GetGasSizeZ()*0.5;
    if(test < 4) continue;

    // std::cout << cx << " " << cy << " " << cz << std::endl;
    // std::cout << "#" << fModuleX << " " << fModuleY << " " << fModuleZ << std::endl;

    if(cx != -1 && cx != fModuleX) continue;
    if(cy != -1 && cy != fModuleY) continue;
    if(cz != -1 && cz != fModuleZ) continue;
    // std::cout << "***" << cx << " " << cy << " " << cz << std::endl;

    Double_t eloss = Ne;

    if(gHarpoDebug>3) 
      Info("GenerateTrackFromText","%f %f %f %d",x,y,z,Ne);

    TpcSimIonisationPoint* point = new TpcSimIonisationPoint(Ne,x,y,z,t,eloss,0);
    trackout->AddPoint(point);
    point->Delete();
  }

  if(gHarpoDebug>1) 
    Info("GenerateTrackFromText","End track %s",line);
 
  return trackout;

}



Int_t TpcSimIonisationChamber::GetNe(Double_t Eloss)
{

  Int_t Ne;

  // ****** Number of electrons (primary + secondary) *****
  if(fFano <= 0){
    Ne = Int_t((Eloss - fIpot)/fW) + 1;         // *
  } else{
    Double_t netmp = (Eloss - fIpot)/fW;
    while(1){
      if(netmp<=0){ Ne = 0; break;}
      Double_t elosstmp = gRandom->Gaus(netmp, TMath::Sqrt(fFano*netmp));
      if(elosstmp<0) continue;
      Ne = Int_t(elosstmp);
      break;
    }
  }
  if(Ne>1000) Ne = 1000;   // Delta electron cut        *
  // ******************************************************


  return Ne;
}



Double_t TpcSimIonisationChamber::GetEloss()
{
                                    
  Double_t tmprndm = gRandom->Rndm();
  Double_t Eloss;  

  // switch(fEdist){
  // case 0: // ALICE power law
    Eloss = fIpot*TMath::Power(tmprndm+6e-04,-5./6.);
  //   break;
  // default: // from file
  //   Eloss = fStraggleFunc->Eval(tmprndm);
  //   break;
  // }

  // if(Eloss < 0){
  //   Debug(1,Form("%f | %f",tmprndm,Eloss));
  //   Eloss = 1000000;
  // }
  
  return Eloss;

}