[u/mrichter/AliRoot.git] / TPC / Upgrade / AliToyMCEventGeneratorSimple.cxx

#include <iostream>

#include <TDatabasePDG.h>
#include <TRandom.h>
#include <TF1.h>
#include <TStopwatch.h>
#include <TFile.h>
#include <TTree.h>
#include <TRandom3.h>

#include <AliESDtrackCuts.h>
#include <AliESDtrack.h>
#include <AliTPCLaserTrack.h>

#include "AliToyMCEvent.h"

#include "AliToyMCEventGeneratorSimple.h"

/*


*/

ClassImp(AliToyMCEventGeneratorSimple);


AliToyMCEventGeneratorSimple::AliToyMCEventGeneratorSimple()
  :AliToyMCEventGenerator()
  ,fVertexMean(0.)
  ,fVertexSigma(7.)
  ,fNtracks(20)
  ,fInputFileNameESD("")
  ,fESDCuts(0x0)
  ,fInputIndex(0)
  ,fESDEvent(0x0)
  ,fESDTree(0x0)
  ,fInputFile(0x0)
  ,fHPt(0x0)
  ,fHEta(0x0)
  ,fHMult(0x0)
  ,fHistosSet(kFALSE)
  ,fParamFile(0x0)
{
  //
  
}
//________________________________________________________________
AliToyMCEventGeneratorSimple::AliToyMCEventGeneratorSimple(const AliToyMCEventGeneratorSimple &gen)
  :AliToyMCEventGenerator(gen)
  ,fVertexMean(gen.fVertexMean)
  ,fVertexSigma(gen.fVertexSigma)
  ,fNtracks(gen.fNtracks)
  ,fInputFileNameESD(gen.fInputFileNameESD)
  ,fESDCuts(gen.fESDCuts)
  ,fInputIndex(gen.fInputIndex)
  ,fESDEvent(gen.fESDEvent)
  ,fESDTree(gen.fESDTree)
  ,fInputFile(gen.fInputFile)
  ,fHPt(gen.fHPt)
  ,fHEta(gen.fHEta)
  ,fHMult(gen.fHMult)
  ,fHistosSet(gen.fHistosSet)
  ,fParamFile(0x0)
{
}
//________________________________________________________________
AliToyMCEventGeneratorSimple::~AliToyMCEventGeneratorSimple()
{
}
 //________________________________________________________________
AliToyMCEventGeneratorSimple& AliToyMCEventGeneratorSimple::operator = (const AliToyMCEventGeneratorSimple &gen)
{
  //assignment operator
  if (&gen == this) return *this;
  new (this) AliToyMCEventGeneratorSimple(gen);

  return *this;
}
//________________________________________________________________
void AliToyMCEventGeneratorSimple::SetParametersToyGen(const Char_t* parfilename/*="files/params.root*/, Double_t vertexMean/*=0*/, Double_t vertexSigma/*=7.*/) {
  fVertexMean = vertexMean;
  fVertexSigma = vertexSigma;
  fParamFile = new TFile(parfilename, "read");
  fHPt = (TH1F*) fParamFile->Get("hPt");
  fHEta = (TH1F*) fParamFile->Get("hEta"); 
  fHMult = (TH1I*) fParamFile->Get("hMult") ;
  fHistosSet = kTRUE;
 
  
}
//________________________________________________________________
AliToyMCEvent* AliToyMCEventGeneratorSimple::Generate(Double_t time)
{
  //
  //
  //
  
  AliToyMCEvent *retEvent = new AliToyMCEvent();
  retEvent->SetT0(time);
  retEvent->SetX(0);
  retEvent->SetX(0);
  Double_t vertexZ=0;
  //limit vertex to +- 10cm
  while ( TMath::Abs(vertexZ=gRandom->Gaus(fVertexMean,fVertexSigma))>10. );
  retEvent->SetZ(vertexZ);

  Double_t etaCuts=0.9;
  Double_t mass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
  static TF1 fpt("fpt",Form("x*(1+(sqrt(x*x+%f^2)-%f)/([0]*[1]))^(-[0])",mass,mass),0.3,100);
  if (fpt.GetParameter(0)<1){
    printf("Set Parameters\n");
    fpt.SetParameters(7.24,0.120);
    fpt.SetNpx(200);
  }
//   Int_t nTracks = 400; //TODO: draw from experim dist
//   Int_t nTracks = 20; //TODO: draw from experim dist
   Int_t nTracksLocal = fNtracks;
  if(fHistosSet) {
    nTracksLocal = fHMult->GetRandom(); 
    if(nTracksLocal > fNtracks) nTracksLocal = fNtracks;
  }
  std::cout << "Generating " << nTracksLocal << " tracks " << std::endl;

  for(Int_t iTrack=0; iTrack<nTracksLocal; iTrack++){
    Double_t phi = gRandom->Uniform(0.0, 2*TMath::Pi());
    Double_t eta = 0.;
    Double_t pt = 0.;
    
    if(fHistosSet) {//draw from histograms if set
      eta = fHEta->GetRandom();
      pt = fHPt->GetRandom();
    }
    else {
      eta = gRandom->Uniform(-etaCuts, etaCuts);
      pt = fpt.GetRandom(); // momentum for f1
    }
    Short_t sign=1;
    if(gRandom->Rndm() < 0.5){
      sign =1;
    }else{
      sign=-1;
    }
    
    Double_t theta = 2*TMath::ATan(TMath::Exp(-eta))-TMath::Pi()/2.;
    Double_t pxyz[3];
    pxyz[0]=pt*TMath::Cos(phi);
    pxyz[1]=pt*TMath::Sin(phi);
    pxyz[2]=pt*TMath::Tan(theta);
    Double_t vertex[3]={0,0,retEvent->GetZ()};
    Double_t cv[21]={0};
    AliToyMCTrack *tempTrack = retEvent->AddTrack(vertex,pxyz,cv,sign);
    // use unique ID for track number
    // this will be used in DistortTrack track to set the cluster label
    // in one simulation the track id should be unique for performance studies
    tempTrack->SetUniqueID(fCurrentTrack++);
    DistortTrack(*tempTrack, time);
  }

  return retEvent;
}

//________________________________________________________________
void AliToyMCEventGeneratorSimple::RunSimulation(const Int_t nevents/*=10*/, const Int_t ntracks/*=400*/, const Int_t rate/*=50*/)
{
  //
  // run simple simulation with equal event spacing
  // rate in kHz
  //

  if (!ConnectOutputFile()) return;
  //initialise the space charge. Should be done after the tree was set up
  InitSpaceCharge();

  // number of tracks to simulate per interaction
  fNtracks=ntracks;
  // within one simulation the track count should be unique for effeciency studies
  // don't use the track ID 0 since the label -0 is not different from 0
  fCurrentTrack=1;
  
  Double_t eventTime=0.;
  const Double_t eventSpacing=1./rate/1e3; //50kHz equally spaced
  TStopwatch s;
  for (Int_t ievent=0; ievent<nevents; ++ievent){
    printf("Generating event %3d (%.3g)\n",ievent,eventTime);
    fEvent = Generate(eventTime);
    FillTree();
    delete fEvent;
    fEvent=0x0;
    eventTime+=eventSpacing;
  }
  s.Stop();
  s.Print();
  
  CloseOutputFile();
}

//________________________________________________________________
void AliToyMCEventGeneratorSimple::RunSimulationLaser(const Int_t nevents/*=1*/)
{
  //
  // run simple simulation with equal event spacing
  //
  
  if (!ConnectOutputFile()) return;
  //initialise the space charge. Should be done after the tree was set up
  InitSpaceCharge();
  
  // within one simulation the track count should be unique for effeciency studies
  fCurrentTrack=1;
  
  Double_t eventTime=0.;
  const Double_t eventSpacing=1./10.; //laser is running at 10Hz equally spaced
  TStopwatch s;
  for (Int_t ievent=0; ievent<nevents; ++ievent){
    printf("Generating event %3d (%.3g)\n",ievent,eventTime);
    fEvent = GenerateLaser(eventTime);
    FillTree();
    delete fEvent;
    fEvent=0x0;
    eventTime+=eventSpacing;
  }
  s.Stop();
  s.Print();
  
  CloseOutputFile();
}

//________________________________________________________________
AliToyMCEvent* AliToyMCEventGeneratorSimple::GenerateESD(AliESDEvent &esdEvent, Double_t time) {

 
  AliToyMCEvent *retEvent = new AliToyMCEvent();
  retEvent->SetT0(time);
  retEvent->SetX(esdEvent.GetPrimaryVertex()->GetX());
  retEvent->SetY(esdEvent.GetPrimaryVertex()->GetY());
  retEvent->SetZ(esdEvent.GetPrimaryVertex()->GetZ());


  if(!fNtracks) fNtracks =  esdEvent.GetNumberOfTracks();
  Int_t nUsedTracks = 0;
  for(Int_t iTrack=0; iTrack<esdEvent.GetNumberOfTracks(); iTrack++){

    AliESDtrack *part = esdEvent.GetTrack(iTrack);
    if(!part) continue;
    if (!fESDCuts->AcceptTrack(/*(AliESDtrack*)*/part))continue;
    if(part->Pt() < 0.3) continue; //avoid tracks that fail to propagate through entire volume
    Double_t pxyz[3];
    pxyz[0]=part->Px();
    pxyz[1]=part->Py();
    pxyz[2]=part->Pz();
    Double_t vertex[3]={retEvent->GetX(),retEvent->GetY(),retEvent->GetZ()};
    
    Double_t cv[21]={0};
    Int_t sign = part->Charge();
    AliToyMCTrack *tempTrack = retEvent->AddTrack(vertex,pxyz,cv,sign);
    // use unique ID for track number
    // this will be used in DistortTrack track to set the cluster label
    // in one simulation the track id should be unique for performance studies
    tempTrack->SetUniqueID(fCurrentTrack++);
    DistortTrack(*tempTrack, time);
    nUsedTracks++;
    
    if(nUsedTracks >= fNtracks) break;
  }

  return retEvent;
}
//________________________________________________________________
void AliToyMCEventGeneratorSimple::RunSimulationESD(const Int_t nevents/*=10*/, const Int_t ntracks/*=400*/)
{
  //
  // run simulation using esd input with equal event spacing
  //

  if (!ConnectOutputFile()) return;
  TFile f(Form("%s%s",fInputFileNameESD.Data(),"#AliESDs.root"));
  if(!f.IsOpen()) {
    std::cout << "file "<<fInputFileNameESD.Data() <<" could not be opened" << std::endl;
    return;
  }
  TTree* esdTree = (TTree*) f.Get("esdTree");
  if(!esdTree) {
    std::cout <<"no esdTree in file " << std::endl;
    return;
  }
  InitSpaceCharge();
  AliESDEvent* esdEvent = new AliESDEvent();
 
 
  esdEvent->ReadFromTree(esdTree);

  fNtracks = ntracks;
  // within one simulation the track count should be unique for effeciency studies
  fCurrentTrack=1;
  
  Double_t eventTime=0.;
  const Double_t eventSpacing=1./50e3; //50kHz equally spaced
  TStopwatch s;
  
  fESDCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2011(kTRUE,1);
  gRandom->SetSeed();
  Int_t nUsedEvents = 0;
  for (Int_t ievent=0; ievent<esdTree->GetEntries(); ++ievent){
    printf("Generating event %3d (%.3g)\n",nUsedEvents,eventTime);
    esdTree->GetEvent(ievent);
    if(esdEvent->GetNumberOfTracks()==0) {
      std::cout << " tracks == 0" << std::endl;
      continue;
    }
   
    fEvent = GenerateESD(*esdEvent, eventTime);
    if(fEvent->GetNumberOfTracks() >=10) {
      nUsedEvents++;
      FillTree();
      eventTime+=eventSpacing;
    }
    delete fEvent;
    fEvent=0x0;
    
    if(nUsedEvents>=nevents) break;
  }
  s.Stop();
  s.Print();
  
  CloseOutputFile();
}

//________________________________________________________________
void AliToyMCEventGeneratorSimple::RunSimulationBunchTrain(const Int_t nevents/*=10*/, const Int_t ntracks/*=400*/)
{
  //
  // run simple simulation with equal event spacing
  //

  //Parameters for bunc
  const Double_t abortGap = 3e-6; //
  const Double_t collFreq = 50e3;
  const Double_t bSpacing = 50e-9; //bunch spacing
  const Int_t nTrainBunches = 48;
  const Int_t nTrains = 12;
  const Double_t revFreq = 1.11e4; //revolution frequency
  const Double_t collProb = collFreq/(nTrainBunches*nTrains*revFreq);
  const Double_t trainLength = bSpacing*(nTrainBunches-1);
  const Double_t totTrainLength = nTrains*trainLength;
  const Double_t trainSpacing = (1./revFreq - abortGap - totTrainLength)/(nTrains-1); 
  Bool_t equalSpacing = kFALSE;

  TRandom3 *rand = new TRandom3();
  //rand->SetSeed();

  if (!ConnectOutputFile()) return;
  //initialise the space charge. Should be done after the tree was set up
  InitSpaceCharge();
  
  fNtracks=ntracks;
  // within one simulation the track count should be unique for effeciency studies
  fCurrentTrack=1;
  
  Double_t eventTime=0.;
  //  const Double_t eventSpacing=1./50e3; //50kHz equally spaced
  TStopwatch s;
  Int_t nGeneratedEvents = 0;
  Int_t bunchCounter = 0;
  Int_t trainCounter = 0;

  //for (Int_t ievent=0; ievent<nevents; ++ievent){
  while (nGeneratedEvents<nevents){
    //  std::cout <<trainCounter << " " << bunchCounter << " "<< "eventTime " << eventTime << std::endl;
    
    if(equalSpacing)  {
      printf("Generating event %3d (%.3g)\n",nGeneratedEvents,eventTime);
      fEvent = Generate(eventTime);
      nGeneratedEvents++;
      FillTree();
      delete fEvent;
      fEvent=0x0;
      eventTime+=1./collFreq;
    }
    else{
      Int_t nCollsInCrossing = rand -> Poisson(collProb);
      for(Int_t iColl = 0; iColl<nCollsInCrossing; iColl++){
	printf("Generating event %3d (%.3g)\n",nGeneratedEvents,eventTime);
	fEvent = Generate(eventTime);
	nGeneratedEvents++;
	FillTree();
	delete fEvent;
	fEvent=0x0;

      }
      bunchCounter++;

      if(bunchCounter>=nTrainBunches){
	
	trainCounter++;
	if(trainCounter>=nTrains){
	  eventTime+=abortGap;
	  trainCounter=0;
	}
	else eventTime+=trainSpacing;

	bunchCounter=0;
      }
      else eventTime+= bSpacing;
      
    }


  }
  s.Stop();
  s.Print();
  delete rand;
  CloseOutputFile();
}


//________________________________________________________________
Int_t AliToyMCEventGeneratorSimple::OpenInputAndGetMaxEvents(const Int_t type, const Int_t nevents) {

  
  if(type==0) return nevents;
  
  if(type==1) {
    
    fInputFile = new TFile(Form("%s%s",fInputFileNameESD.Data(),"#AliESDs.root"),"read");
    if(!fInputFile->IsOpen()) {
      std::cout << "file "<<fInputFileNameESD.Data() <<" could not be opened" << std::endl;
      return 0;
    }
    fESDTree = (TTree*) fInputFile->Get("esdTree");
    if(!fESDTree) {
      std::cout <<"no esdTree in file " << std::endl;
      return 0;
    }
    fESDEvent = new AliESDEvent();
    fESDEvent->ReadFromTree(fESDTree);
    fESDCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2011(kTRUE,1);

    fInputIndex = 0;

    return fESDTree->GetEntries();
    gRandom->SetSeed();
   }

 
  std::cout << " no available input type (toymc, esd) specified" << std::endl;
  return 0;
 

}


//________________________________________________________________
void AliToyMCEventGeneratorSimple::RunSimulation2(const Bool_t equalspacing, const Int_t type, const Int_t nevents, const Int_t ntracks) {

  //type==0 simple toy
  //type==1 esd input

  Int_t nMaxEvents = OpenInputAndGetMaxEvents(type,nevents);
  if (!ConnectOutputFile()) return;
  
  InitSpaceCharge();


  fNtracks = ntracks;
  // within one simulation the track count should be unique for effeciency studies
  fCurrentTrack=1;
  
  Double_t eventTime=0.;
  TStopwatch s;
  // const Double_t eventSpacing=1./50e3; //50kHz equally spaced

  Int_t nGeneratedEvents = 0;
  
  for (Int_t ievent=0; ievent<nMaxEvents; ++ievent){
    printf("Generating event %3d (%.3g)\n",ievent,eventTime);
    
    Int_t nColls = 0;
    Double_t spacing = 0;
    GetNGeneratedEventsAndSpacing(equalspacing, nColls, spacing);

    for(Int_t iColl = 0; iColl<nColls; iColl++){
      if(type==0)    fEvent = Generate(eventTime);
      else if(type==1) fEvent = GenerateESD2(eventTime);
      if(fEvent) {
	FillTree();
	delete fEvent;
	fEvent=0x0;
	nGeneratedEvents++;
      }
    }
    eventTime+=spacing;
    if(nGeneratedEvents >= nevents) break;
  }
  s.Stop();
  s.Print();
  
  CloseOutputFile();
  CloseInputFile();


}
//________________________________________________________________
AliToyMCEvent* AliToyMCEventGeneratorSimple::GenerateESD2(Double_t time) {

  //test that enough tracks will pass cuts (and that there is tracks at all)
  Bool_t testEvent = kTRUE;
  while(fESDTree->GetEvent(fInputIndex) && testEvent) {  
    
    Int_t nPassedCuts = 0;
    for(Int_t iTrack = 0; iTrack<fESDEvent->GetNumberOfTracks(); iTrack++){
      if(fESDCuts->AcceptTrack(fESDEvent->GetTrack(iTrack)) && (fESDEvent->GetTrack(iTrack)->Pt() < 0.3) ) nPassedCuts++;
    }
    if(nPassedCuts<fNtracks || (fNtracks>100 && nPassedCuts <100) ) fInputIndex++; //100 is a perhaps bit arbitrary, do we want the events were only a small number of tracks are accepted?
    else testEvent = kFALSE;
  }

  if(fInputIndex>=fESDTree->GetEntries()) return 0;
  
  
  AliToyMCEvent *retEvent = new AliToyMCEvent();
  retEvent->SetT0(time);
  retEvent->SetX(fESDEvent->GetPrimaryVertex()->GetX());
  retEvent->SetY(fESDEvent->GetPrimaryVertex()->GetY());
  retEvent->SetZ(fESDEvent->GetPrimaryVertex()->GetZ());


  if(!fNtracks) fNtracks =  fESDEvent->GetNumberOfTracks();
  Int_t nUsedTracks = 0;
  for(Int_t iTrack=0; iTrack<fESDEvent->GetNumberOfTracks(); iTrack++){

    AliESDtrack *part = fESDEvent->GetTrack(iTrack);
    if(!part) continue;
    if (!fESDCuts->AcceptTrack(/*(AliESDtrack*)*/part))continue;
    if(part->Pt() < 0.3) continue; //avoid tracks that fail to propagate through entire volume

    Double_t pxyz[3];
    pxyz[0]=part->Px();
    pxyz[1]=part->Py();
    pxyz[2]=part->Pz();
    Double_t vertex[3]={retEvent->GetX(),retEvent->GetY(),retEvent->GetZ()};
    Double_t cv[21]={0};
    Int_t sign = part->Charge();

    AliToyMCTrack *tempTrack = retEvent->AddTrack(vertex,pxyz,cv,sign);
    // use unique ID for track number
    // this will be used in DistortTrack track to set the cluster label
    // in one simulation the track id should be unique for performance studies
    tempTrack->SetUniqueID(fCurrentTrack++);
    DistortTrack(*tempTrack, time);
    
    
    if(nUsedTracks >= fNtracks) break;
  }
  fInputIndex++;
 
  return retEvent;
}

//________________________________________________________________
AliToyMCEvent* AliToyMCEventGeneratorSimple::GenerateLaser(Double_t time)
{
  //
  // Generate an Event with laser tracks
  //
  
  AliToyMCEvent *retEvent = new AliToyMCEvent();
  retEvent->SetEventType(AliToyMCEvent::kLaser);
  
  retEvent->SetT0(time);
  retEvent->SetX(0);
  retEvent->SetX(0);
  retEvent->SetZ(0);
  
  AliTPCLaserTrack::LoadTracks();
  TObjArray *arr=AliTPCLaserTrack::GetTracks();

  //since we have a laser track force no material budges
  Bool_t materialBudget=GetUseMaterialBudget();
  SetUseMaterialBudget(kFALSE);

  //the laser tracks have partially large inclination angles over the pads
  // -> relax the propagation contraint and switch off error messages
  SetIsLaser(kTRUE);
  Int_t debug=AliLog::GetDebugLevel("","AliToyMCEventGeneratorSimple");
  AliLog::SetClassDebugLevel("AliToyMCEventGeneratorSimple",-3);
  
  for (Int_t iTrack=0; iTrack<arr->GetEntriesFast(); ++iTrack){
    AliExternalTrackParam *track=(AliExternalTrackParam*)arr->At(iTrack);
    AliToyMCTrack *tempTrack = retEvent->AddTrack(AliToyMCTrack(*track));
    // for laser only TPC clusters exist
    tempTrack->SetUniqueID(fCurrentTrack++);
    MakeTPCClusters(*tempTrack, time);
  }

  SetIsLaser(kFALSE);
  AliLog::SetClassDebugLevel("AliToyMCEventGeneratorSimple",debug);
  SetUseMaterialBudget(materialBudget);
  return retEvent;
}

//________________________________________________________________
void AliToyMCEventGeneratorSimple::GetNGeneratedEventsAndSpacing(const Bool_t equalSpacing, Int_t &ngen, Double_t &spacing)
{

  static Int_t bunchCounter = 0;
  static Int_t trainCounter = 0;

  if(equalSpacing) {
    ngen =1;
    spacing = 1./50e3; //50kHz equally spaced
    return;
  }
  else if(!equalSpacing) {
      //parameters for bunch train
    const Double_t abortGap = 3e-6; //
    const Double_t collFreq = 50e3;
    const Double_t bSpacing = 50e-9; //bunch spacing
    const Int_t nTrainBunches = 48;
    const Int_t nTrains = 12;
    const Double_t revFreq = 1.11e4; //revolution frequency
    const Double_t collProb = collFreq/(nTrainBunches*nTrains*revFreq);
    const Double_t trainLength = bSpacing*(nTrainBunches-1);
    const Double_t totTrainLength = nTrains*trainLength;
    const Double_t trainSpacing = (1./revFreq - abortGap - totTrainLength)/(nTrains-1); 
    Double_t time = 0;
    Int_t nCollsInCrossing = 0;
    while(nCollsInCrossing ==0){
      
      
    bunchCounter++;
    
    if(bunchCounter>=nTrainBunches){
      
      trainCounter++;
      if(trainCounter>=nTrains){
	time+=abortGap;
	trainCounter=0;
      }
      else time+=trainSpacing;

      bunchCounter=0;
    }
    else time+= bSpacing;


    nCollsInCrossing = gRandom -> Poisson(collProb);
    //std::cout << " nCollsInCrossing " << nCollsInCrossing <<  std::endl;
    }
    ngen = nCollsInCrossing;
    if(nCollsInCrossing > 1)std::cout << " nCollsInCrossing " << nCollsInCrossing <<  std::endl;
    spacing = time;
    return;

  }

}

//________________________________________________________________
Bool_t AliToyMCEventGeneratorSimple::CloseInputFile()
{
  //
  // close the input file
  //
  if (!fInputFile) return kFALSE;
  fInputFile->Close();
  delete fInputFile;
  fInputFile=0x0;

  return kTRUE;
}
Commit	Line	Data
402a4145	1	#include <iostream>
	2
	3	#include <TDatabasePDG.h>
	4	#include <TRandom.h>
	5	#include <TF1.h>
	6	#include <TStopwatch.h>
	7	#include <TFile.h>
	8	#include <TTree.h>
	9	#include <TRandom3.h>
	10
	11	#include <AliESDtrackCuts.h>
	12	#include <AliESDtrack.h>
	13	#include <AliTPCLaserTrack.h>
	14
	15	#include "AliToyMCEvent.h"
	16
	17	#include "AliToyMCEventGeneratorSimple.h"
	18
	19	/*
	20
	21
	22	*/
	23
	24	ClassImp(AliToyMCEventGeneratorSimple);
	25
	26
	27	AliToyMCEventGeneratorSimple::AliToyMCEventGeneratorSimple()
	28	:AliToyMCEventGenerator()
	29	,fVertexMean(0.)
	30	,fVertexSigma(7.)
	31	,fNtracks(20)
	32	,fInputFileNameESD("")
	33	,fESDCuts(0x0)
	34	,fInputIndex(0)
	35	,fESDEvent(0x0)
	36	,fESDTree(0x0)
	37	,fInputFile(0x0)
	38	,fHPt(0x0)
	39	,fHEta(0x0)
	40	,fHMult(0x0)
	41	,fHistosSet(kFALSE)
	42	,fParamFile(0x0)
	43	{
	44	//
	45
	46	}
	47	//________________________________________________________________
	48	AliToyMCEventGeneratorSimple::AliToyMCEventGeneratorSimple(const AliToyMCEventGeneratorSimple &gen)
	49	:AliToyMCEventGenerator(gen)
	50	,fVertexMean(gen.fVertexMean)
	51	,fVertexSigma(gen.fVertexSigma)
	52	,fNtracks(gen.fNtracks)
	53	,fInputFileNameESD(gen.fInputFileNameESD)
	54	,fESDCuts(gen.fESDCuts)
	55	,fInputIndex(gen.fInputIndex)
	56	,fESDEvent(gen.fESDEvent)
	57	,fESDTree(gen.fESDTree)
	58	,fInputFile(gen.fInputFile)
	59	,fHPt(gen.fHPt)
	60	,fHEta(gen.fHEta)
	61	,fHMult(gen.fHMult)
	62	,fHistosSet(gen.fHistosSet)
	63	,fParamFile(0x0)
	64	{
65	}
66	//________________________________________________________________
67	AliToyMCEventGeneratorSimple::~AliToyMCEventGeneratorSimple()
68	{
69	}
70	//________________________________________________________________
71	AliToyMCEventGeneratorSimple& AliToyMCEventGeneratorSimple::operator = (const AliToyMCEventGeneratorSimple &gen)
72	{
73	//assignment operator
74	if (&gen == this) return *this;
75	new (this) AliToyMCEventGeneratorSimple(gen);
76
77	return *this;
78	}
79	//________________________________________________________________
80	void AliToyMCEventGeneratorSimple::SetParametersToyGen(const Char_t* parfilename/="files/params.root/, Double_t vertexMean/=0/, Double_t vertexSigma/=7./) {
81	fVertexMean = vertexMean;
82	fVertexSigma = vertexSigma;
83	fParamFile = new TFile(parfilename, "read");
84	fHPt = (TH1F*) fParamFile->Get("hPt");
85	fHEta = (TH1F*) fParamFile->Get("hEta");
86	fHMult = (TH1I*) fParamFile->Get("hMult") ;
87	fHistosSet = kTRUE;
88
89
90	}
91	//________________________________________________________________
92	AliToyMCEvent* AliToyMCEventGeneratorSimple::Generate(Double_t time)
93	{
94	//
95	//
96	//
97
98	AliToyMCEvent *retEvent = new AliToyMCEvent();
99	retEvent->SetT0(time);
100	retEvent->SetX(0);
101	retEvent->SetX(0);
102	Double_t vertexZ=0;
103	//limit vertex to +- 10cm
104	while ( TMath::Abs(vertexZ=gRandom->Gaus(fVertexMean,fVertexSigma))>10. );
105	retEvent->SetZ(vertexZ);
106
107	Double_t etaCuts=0.9;
108	Double_t mass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
109	static TF1 fpt("fpt",Form("x(1+(sqrt(xx+%f^2)-%f)/([0]*[1]))^(-[0])",mass,mass),0.3,100);
110	if (fpt.GetParameter(0)<1){
111	printf("Set Parameters\n");
112	fpt.SetParameters(7.24,0.120);
113	fpt.SetNpx(200);
114	}
115	// Int_t nTracks = 400; //TODO: draw from experim dist
116	// Int_t nTracks = 20; //TODO: draw from experim dist
117	Int_t nTracksLocal = fNtracks;
118	if(fHistosSet) {
119	nTracksLocal = fHMult->GetRandom();
120	if(nTracksLocal > fNtracks) nTracksLocal = fNtracks;
121	}
122	std::cout << "Generating " << nTracksLocal << " tracks " << std::endl;
123
124	for(Int_t iTrack=0; iTrack<nTracksLocal; iTrack++){
125	Double_t phi = gRandom->Uniform(0.0, 2*TMath::Pi());
126	Double_t eta = 0.;
127	Double_t pt = 0.;
128
129	if(fHistosSet) {//draw from histograms if set
130	eta = fHEta->GetRandom();
131	pt = fHPt->GetRandom();
132	}
133	else {
134	eta = gRandom->Uniform(-etaCuts, etaCuts);
135	pt = fpt.GetRandom(); // momentum for f1
136	}
137	Short_t sign=1;
138	if(gRandom->Rndm() < 0.5){
139	sign =1;
140	}else{
141	sign=-1;
142	}
143
144	Double_t theta = 2*TMath::ATan(TMath::Exp(-eta))-TMath::Pi()/2.;
145	Double_t pxyz[3];
146	pxyz[0]=pt*TMath::Cos(phi);
147	pxyz[1]=pt*TMath::Sin(phi);
148	pxyz[2]=pt*TMath::Tan(theta);
149	Double_t vertex[3]={0,0,retEvent->GetZ()};
150	Double_t cv[21]={0};
151	AliToyMCTrack *tempTrack = retEvent->AddTrack(vertex,pxyz,cv,sign);
152	// use unique ID for track number
153	// this will be used in DistortTrack track to set the cluster label
154	// in one simulation the track id should be unique for performance studies
155	tempTrack->SetUniqueID(fCurrentTrack++);
156	DistortTrack(*tempTrack, time);
157	}
158
159	return retEvent;
160	}
161
162	//________________________________________________________________
163	void AliToyMCEventGeneratorSimple::RunSimulation(const Int_t nevents/=10/, const Int_t ntracks/=400/, const Int_t rate/=50/)
164	{
165	//
166	// run simple simulation with equal event spacing
167	// rate in kHz
168	//
169
170	if (!ConnectOutputFile()) return;
171	//initialise the space charge. Should be done after the tree was set up
172	InitSpaceCharge();
173
174	// number of tracks to simulate per interaction
175	fNtracks=ntracks;
176	// within one simulation the track count should be unique for effeciency studies
177	// don't use the track ID 0 since the label -0 is not different from 0
178	fCurrentTrack=1;
179
180	Double_t eventTime=0.;
181	const Double_t eventSpacing=1./rate/1e3; //50kHz equally spaced
182	TStopwatch s;
183	for (Int_t ievent=0; ievent<nevents; ++ievent){
184	printf("Generating event %3d (%.3g)\n",ievent,eventTime);
185	fEvent = Generate(eventTime);
186	FillTree();
187	delete fEvent;
188	fEvent=0x0;
189	eventTime+=eventSpacing;
190	}
191	s.Stop();
192	s.Print();
193
194	CloseOutputFile();
195	}
196
197	//________________________________________________________________
198	void AliToyMCEventGeneratorSimple::RunSimulationLaser(const Int_t nevents/=1/)
199	{
200	//
201	// run simple simulation with equal event spacing
202	//
203
204	if (!ConnectOutputFile()) return;
205	//initialise the space charge. Should be done after the tree was set up
206	InitSpaceCharge();
207
208	// within one simulation the track count should be unique for effeciency studies
209	fCurrentTrack=1;
210
211	Double_t eventTime=0.;
212	const Double_t eventSpacing=1./10.; //laser is running at 10Hz equally spaced
213	TStopwatch s;
214	for (Int_t ievent=0; ievent<nevents; ++ievent){
215	printf("Generating event %3d (%.3g)\n",ievent,eventTime);
216	fEvent = GenerateLaser(eventTime);
217	FillTree();
218	delete fEvent;
219	fEvent=0x0;
220	eventTime+=eventSpacing;
221	}
222	s.Stop();
223	s.Print();
224
225	CloseOutputFile();
226	}
227
228	//________________________________________________________________
229	AliToyMCEvent* AliToyMCEventGeneratorSimple::GenerateESD(AliESDEvent &esdEvent, Double_t time) {
230
231
232	AliToyMCEvent *retEvent = new AliToyMCEvent();
233	retEvent->SetT0(time);
234	retEvent->SetX(esdEvent.GetPrimaryVertex()->GetX());
235	retEvent->SetY(esdEvent.GetPrimaryVertex()->GetY());
236	retEvent->SetZ(esdEvent.GetPrimaryVertex()->GetZ());
237
238
239
240	if(!fNtracks) fNtracks = esdEvent.GetNumberOfTracks();
241	Int_t nUsedTracks = 0;
242	for(Int_t iTrack=0; iTrack<esdEvent.GetNumberOfTracks(); iTrack++){
243
244	AliESDtrack *part = esdEvent.GetTrack(iTrack);
245	if(!part) continue;
246	if (!fESDCuts->AcceptTrack(/(AliESDtrack)*/part))continue;
247	if(part->Pt() < 0.3) continue; //avoid tracks that fail to propagate through entire volume
248	Double_t pxyz[3];
249	pxyz[0]=part->Px();
250	pxyz[1]=part->Py();
251	pxyz[2]=part->Pz();
252	Double_t vertex[3]={retEvent->GetX(),retEvent->GetY(),retEvent->GetZ()};
253
254	Double_t cv[21]={0};
255	Int_t sign = part->Charge();
256	AliToyMCTrack *tempTrack = retEvent->AddTrack(vertex,pxyz,cv,sign);
257	// use unique ID for track number
258	// this will be used in DistortTrack track to set the cluster label
259	// in one simulation the track id should be unique for performance studies
260	tempTrack->SetUniqueID(fCurrentTrack++);
261	DistortTrack(*tempTrack, time);
262	nUsedTracks++;
263
264	if(nUsedTracks >= fNtracks) break;
265	}
266
267	return retEvent;
268	}
269	//________________________________________________________________
270	void AliToyMCEventGeneratorSimple::RunSimulationESD(const Int_t nevents/=10/, const Int_t ntracks/=400/)
271	{
272	//
273	// run simulation using esd input with equal event spacing
274	//
275
276	if (!ConnectOutputFile()) return;
277	TFile f(Form("%s%s",fInputFileNameESD.Data(),"#AliESDs.root"));
278	if(!f.IsOpen()) {
279	std::cout << "file "<<fInputFileNameESD.Data() <<" could not be opened" << std::endl;
280	return;
281	}
282	TTree* esdTree = (TTree*) f.Get("esdTree");
283	if(!esdTree) {
284	std::cout <<"no esdTree in file " << std::endl;
285	return;
286	}
287	InitSpaceCharge();
288	AliESDEvent* esdEvent = new AliESDEvent();
289
290
291	esdEvent->ReadFromTree(esdTree);
292
293	fNtracks = ntracks;
294	// within one simulation the track count should be unique for effeciency studies
295	fCurrentTrack=1;
296
297	Double_t eventTime=0.;
298	const Double_t eventSpacing=1./50e3; //50kHz equally spaced
299	TStopwatch s;
300
301	fESDCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2011(kTRUE,1);
302	gRandom->SetSeed();
303	Int_t nUsedEvents = 0;
304	for (Int_t ievent=0; ievent<esdTree->GetEntries(); ++ievent){
305	printf("Generating event %3d (%.3g)\n",nUsedEvents,eventTime);
306	esdTree->GetEvent(ievent);
307	if(esdEvent->GetNumberOfTracks()==0) {
308	std::cout << " tracks == 0" << std::endl;
309	continue;
310	}
311
312	fEvent = GenerateESD(*esdEvent, eventTime);
313	if(fEvent->GetNumberOfTracks() >=10) {
314	nUsedEvents++;
315	FillTree();
316	eventTime+=eventSpacing;
317	}
318	delete fEvent;
319	fEvent=0x0;
320
321	if(nUsedEvents>=nevents) break;
322	}
323	s.Stop();
324	s.Print();
325
326	CloseOutputFile();
327	}
328
329	//________________________________________________________________
330	void AliToyMCEventGeneratorSimple::RunSimulationBunchTrain(const Int_t nevents/=10/, const Int_t ntracks/=400/)
331	{
332	//
333	// run simple simulation with equal event spacing
334	//
335
336	//Parameters for bunc
337	const Double_t abortGap = 3e-6; //
338	const Double_t collFreq = 50e3;
339	const Double_t bSpacing = 50e-9; //bunch spacing
340	const Int_t nTrainBunches = 48;
341	const Int_t nTrains = 12;
342	const Double_t revFreq = 1.11e4; //revolution frequency
343	const Double_t collProb = collFreq/(nTrainBunchesnTrainsrevFreq);
344	const Double_t trainLength = bSpacing*(nTrainBunches-1);
345	const Double_t totTrainLength = nTrains*trainLength;
346	const Double_t trainSpacing = (1./revFreq - abortGap - totTrainLength)/(nTrains-1);
347	Bool_t equalSpacing = kFALSE;
348
349	TRandom3 *rand = new TRandom3();
350	//rand->SetSeed();
351
352	if (!ConnectOutputFile()) return;
353	//initialise the space charge. Should be done after the tree was set up
354	InitSpaceCharge();
355
356	fNtracks=ntracks;
357	// within one simulation the track count should be unique for effeciency studies
358	fCurrentTrack=1;
359
360	Double_t eventTime=0.;
361	// const Double_t eventSpacing=1./50e3; //50kHz equally spaced
362	TStopwatch s;
363	Int_t nGeneratedEvents = 0;
364	Int_t bunchCounter = 0;
365	Int_t trainCounter = 0;
366
367	//for (Int_t ievent=0; ievent<nevents; ++ievent){
368	while (nGeneratedEvents<nevents){
369	// std::cout <<trainCounter << " " << bunchCounter << " "<< "eventTime " << eventTime << std::endl;
370
371	if(equalSpacing) {
372	printf("Generating event %3d (%.3g)\n",nGeneratedEvents,eventTime);
373	fEvent = Generate(eventTime);
374	nGeneratedEvents++;
375	FillTree();
376	delete fEvent;
377	fEvent=0x0;
378	eventTime+=1./collFreq;
379	}
380	else{
381	Int_t nCollsInCrossing = rand -> Poisson(collProb);
382	for(Int_t iColl = 0; iColl<nCollsInCrossing; iColl++){
383	printf("Generating event %3d (%.3g)\n",nGeneratedEvents,eventTime);
384	fEvent = Generate(eventTime);
385	nGeneratedEvents++;
386	FillTree();
387	delete fEvent;
388	fEvent=0x0;
389
390	}
391	bunchCounter++;
392
393	if(bunchCounter>=nTrainBunches){
394
395	trainCounter++;
396	if(trainCounter>=nTrains){
397	eventTime+=abortGap;
398	trainCounter=0;
399	}
400	else eventTime+=trainSpacing;
401
402	bunchCounter=0;
403	}
404	else eventTime+= bSpacing;
405
406	}
407
408
409	}
410	s.Stop();
411	s.Print();
412	delete rand;
413	CloseOutputFile();
414	}
415
416
417
418
419
420	//________________________________________________________________
421	Int_t AliToyMCEventGeneratorSimple::OpenInputAndGetMaxEvents(const Int_t type, const Int_t nevents) {
422
423
424
425	if(type==0) return nevents;
426
427	if(type==1) {
428
429	fInputFile = new TFile(Form("%s%s",fInputFileNameESD.Data(),"#AliESDs.root"),"read");
430	if(!fInputFile->IsOpen()) {
431	std::cout << "file "<<fInputFileNameESD.Data() <<" could not be opened" << std::endl;
432	return 0;
433	}
434	fESDTree = (TTree*) fInputFile->Get("esdTree");
435	if(!fESDTree) {
436	std::cout <<"no esdTree in file " << std::endl;
437	return 0;
438	}
439	fESDEvent = new AliESDEvent();
440	fESDEvent->ReadFromTree(fESDTree);
441	fESDCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2011(kTRUE,1);
442
443	fInputIndex = 0;
444
445	return fESDTree->GetEntries();
446	gRandom->SetSeed();
447	}
448
449
450	std::cout << " no available input type (toymc, esd) specified" << std::endl;
451	return 0;
452
453
454	}
455
456
457	//________________________________________________________________
458	void AliToyMCEventGeneratorSimple::RunSimulation2(const Bool_t equalspacing, const Int_t type, const Int_t nevents, const Int_t ntracks) {
459
460	//type==0 simple toy
461	//type==1 esd input
462
463	Int_t nMaxEvents = OpenInputAndGetMaxEvents(type,nevents);
464	if (!ConnectOutputFile()) return;
465
466	InitSpaceCharge();
467
468
469	fNtracks = ntracks;
470	// within one simulation the track count should be unique for effeciency studies
471	fCurrentTrack=1;
472
473	Double_t eventTime=0.;
474	TStopwatch s;
475	// const Double_t eventSpacing=1./50e3; //50kHz equally spaced
476
477	Int_t nGeneratedEvents = 0;
478
479	for (Int_t ievent=0; ievent<nMaxEvents; ++ievent){
480	printf("Generating event %3d (%.3g)\n",ievent,eventTime);
481
482	Int_t nColls = 0;
483	Double_t spacing = 0;
484	GetNGeneratedEventsAndSpacing(equalspacing, nColls, spacing);
485
486	for(Int_t iColl = 0; iColl<nColls; iColl++){
487	if(type==0) fEvent = Generate(eventTime);
488	else if(type==1) fEvent = GenerateESD2(eventTime);
489	if(fEvent) {
490	FillTree();
491	delete fEvent;
492	fEvent=0x0;
493	nGeneratedEvents++;
494	}
495	}
496	eventTime+=spacing;
497	if(nGeneratedEvents >= nevents) break;
498	}
499	s.Stop();
500	s.Print();
501
502	CloseOutputFile();
503	CloseInputFile();
504
505
506	}
507	//________________________________________________________________
508	AliToyMCEvent* AliToyMCEventGeneratorSimple::GenerateESD2(Double_t time) {
509
510	//test that enough tracks will pass cuts (and that there is tracks at all)
511	Bool_t testEvent = kTRUE;
512	while(fESDTree->GetEvent(fInputIndex) && testEvent) {
513
514	Int_t nPassedCuts = 0;
515	for(Int_t iTrack = 0; iTrack<fESDEvent->GetNumberOfTracks(); iTrack++){
516	if(fESDCuts->AcceptTrack(fESDEvent->GetTrack(iTrack)) && (fESDEvent->GetTrack(iTrack)->Pt() < 0.3) ) nPassedCuts++;
517	}
518	if(nPassedCuts<fNtracks \|\| (fNtracks>100 && nPassedCuts <100) ) fInputIndex++; //100 is a perhaps bit arbitrary, do we want the events were only a small number of tracks are accepted?
519	else testEvent = kFALSE;
520	}
521
522	if(fInputIndex>=fESDTree->GetEntries()) return 0;
523
524
525
526	AliToyMCEvent *retEvent = new AliToyMCEvent();
527	retEvent->SetT0(time);
528	retEvent->SetX(fESDEvent->GetPrimaryVertex()->GetX());
529	retEvent->SetY(fESDEvent->GetPrimaryVertex()->GetY());
530	retEvent->SetZ(fESDEvent->GetPrimaryVertex()->GetZ());
531
532
533
534	if(!fNtracks) fNtracks = fESDEvent->GetNumberOfTracks();
535	Int_t nUsedTracks = 0;
536	for(Int_t iTrack=0; iTrack<fESDEvent->GetNumberOfTracks(); iTrack++){
537
538	AliESDtrack *part = fESDEvent->GetTrack(iTrack);
539	if(!part) continue;
540	if (!fESDCuts->AcceptTrack(/(AliESDtrack)*/part))continue;
541	if(part->Pt() < 0.3) continue; //avoid tracks that fail to propagate through entire volume
542
543	Double_t pxyz[3];
544	pxyz[0]=part->Px();
545	pxyz[1]=part->Py();
546	pxyz[2]=part->Pz();
547	Double_t vertex[3]={retEvent->GetX(),retEvent->GetY(),retEvent->GetZ()};
548	Double_t cv[21]={0};
549	Int_t sign = part->Charge();
550
551	AliToyMCTrack *tempTrack = retEvent->AddTrack(vertex,pxyz,cv,sign);
552	// use unique ID for track number
553	// this will be used in DistortTrack track to set the cluster label
554	// in one simulation the track id should be unique for performance studies
555	tempTrack->SetUniqueID(fCurrentTrack++);
556	DistortTrack(*tempTrack, time);
557
558
559	if(nUsedTracks >= fNtracks) break;
560	}
561	fInputIndex++;
562
563	return retEvent;
564	}
565
566	//________________________________________________________________
567	AliToyMCEvent* AliToyMCEventGeneratorSimple::GenerateLaser(Double_t time)
568	{
569	//
570	// Generate an Event with laser tracks
571	//
572
573	AliToyMCEvent *retEvent = new AliToyMCEvent();
574	retEvent->SetEventType(AliToyMCEvent::kLaser);
575
576	retEvent->SetT0(time);
577	retEvent->SetX(0);
578	retEvent->SetX(0);
579	retEvent->SetZ(0);
580
581	AliTPCLaserTrack::LoadTracks();
582	TObjArray *arr=AliTPCLaserTrack::GetTracks();
583
584	//since we have a laser track force no material budges
585	Bool_t materialBudget=GetUseMaterialBudget();
586	SetUseMaterialBudget(kFALSE);
587
588	//the laser tracks have partially large inclination angles over the pads
589	// -> relax the propagation contraint and switch off error messages
590	SetIsLaser(kTRUE);
591	Int_t debug=AliLog::GetDebugLevel("","AliToyMCEventGeneratorSimple");
592	AliLog::SetClassDebugLevel("AliToyMCEventGeneratorSimple",-3);
593
594	for (Int_t iTrack=0; iTrack<arr->GetEntriesFast(); ++iTrack){
595	AliExternalTrackParam track=(AliExternalTrackParam)arr->At(iTrack);
596	AliToyMCTrack tempTrack = retEvent->AddTrack(AliToyMCTrack(track));
597	// for laser only TPC clusters exist
598	tempTrack->SetUniqueID(fCurrentTrack++);
599	MakeTPCClusters(*tempTrack, time);
600	}
601
602	SetIsLaser(kFALSE);
603	AliLog::SetClassDebugLevel("AliToyMCEventGeneratorSimple",debug);
604	SetUseMaterialBudget(materialBudget);
605	return retEvent;
606	}
607
608	//________________________________________________________________
609	void AliToyMCEventGeneratorSimple::GetNGeneratedEventsAndSpacing(const Bool_t equalSpacing, Int_t &ngen, Double_t &spacing)
610	{
611
612	static Int_t bunchCounter = 0;
613	static Int_t trainCounter = 0;
614
615	if(equalSpacing) {
616	ngen =1;
617	spacing = 1./50e3; //50kHz equally spaced
618	return;
619	}
620	else if(!equalSpacing) {
621	//parameters for bunch train
622	const Double_t abortGap = 3e-6; //
623	const Double_t collFreq = 50e3;
624	const Double_t bSpacing = 50e-9; //bunch spacing
625	const Int_t nTrainBunches = 48;
626	const Int_t nTrains = 12;
627	const Double_t revFreq = 1.11e4; //revolution frequency
628	const Double_t collProb = collFreq/(nTrainBunchesnTrainsrevFreq);
629	const Double_t trainLength = bSpacing*(nTrainBunches-1);
630	const Double_t totTrainLength = nTrains*trainLength;
631	const Double_t trainSpacing = (1./revFreq - abortGap - totTrainLength)/(nTrains-1);
632	Double_t time = 0;
633	Int_t nCollsInCrossing = 0;
634	while(nCollsInCrossing ==0){
635
636
637	bunchCounter++;
638
639	if(bunchCounter>=nTrainBunches){
640
641	trainCounter++;
642	if(trainCounter>=nTrains){
643	time+=abortGap;
644	trainCounter=0;
645	}
646	else time+=trainSpacing;
647
648	bunchCounter=0;
649	}
650	else time+= bSpacing;
651
652
653	nCollsInCrossing = gRandom -> Poisson(collProb);
654	//std::cout << " nCollsInCrossing " << nCollsInCrossing << std::endl;
655	}
656	ngen = nCollsInCrossing;
657	if(nCollsInCrossing > 1)std::cout << " nCollsInCrossing " << nCollsInCrossing << std::endl;
658	spacing = time;
659	return;
660
661	}
662
663	}
664
665	//________________________________________________________________
666	Bool_t AliToyMCEventGeneratorSimple::CloseInputFile()
667	{
668	//
669	// close the input file
670	//
671	if (!fInputFile) return kFALSE;
672	fInputFile->Close();
673	delete fInputFile;
674	fInputFile=0x0;
675
676	return kTRUE;
677	}