Add macro for simple simualtor

diff --git a/TRD/simple.C b/TRD/simple.C
new file mode 100644 (file)
index 0000000..638b543
--- /dev/null
+++ b/TRD/simple.C
@@ -0,0 +1,283 @@
+
+#if !defined(__CINT__) || defined(__MAKECINT__)
+
+#include "TH1.h"
+#include "TFile.h"
+#include "TCanvas.h"
+#include "TProfile.h"
+#include "TClonesArray.h"
+
+#include "AliModule.h"
+#include "AliFRAMEv1.h"
+
+#include "AliTRDv1.h"
+#include "AliTRDhit.h"
+#include "AliTRDsim.h"
+#include "AliTRDsimple.h"
+#include "AliTRDsimpleGen.h"
+#include "AliTRDsimpleMC.h"
+#include "AliTRDgeometry.h"
+#include "AliTRDdigitizer.h"
+#include "AliTRDdigitsManager.h"
+#include "AliTRDdataArrayI.h"
+
+#endif
+
+void simple()
+{
+
+  ///////////////////////////////////////////////////////////////////////////////
+  //
+  //  Macro to run the simplified simulator
+  //
+  ///////////////////////////////////////////////////////////////////////////////
+
+  //_____________________________________________________________________________ 
+  //
+  // Initialization part
+  //_____________________________________________________________________________ 
+  //
+
+  // The simple simulator
+  AliTRDsimple *simple = new AliTRDsimple();
+
+  // Initialize a dummy frame so that the TRD works
+  new AliFRAMEv1("FRAME","Space Frame");             
+
+  // Initialize the TRD detector
+  AliTRDv1     *trd    = new AliTRDv1("TRD","TRD slow simulator");
+
+  // Needed for some material properties 
+  trd->CreateMaterials();
+
+  // Use PAI model (most propable energy loss) for the relativistic rise
+  //trd->SetPAI(1); 
+  // Use GEANT (mean energy loss) for the relativistic rise
+  trd->SetPAI(0);
+ 
+  // Use the Ermilova spectrum for the delta electrons
+  trd->SetErmilova(1);
+  // Use the GEANT spectrum for the delta electrons
+  //trd->SetErmilova(0);
+ 
+  // Select the gas mixture (0: 97% Xe + 3% isobutane, 1: 85% Xe + 15% CO2)
+  trd->SetGasMix(1);
+ 
+  // Get the pointer to the geometry object
+  AliTRDgeometry *geometry = trd->GetGeometry();
+ 
+  // The number of timebins
+  geometry->SetNTimeBin(30);
+ 
+  // The additional timebins before and after the drift region
+  geometry->SetExpandTimeBin(5,15);
+    
+  // Switch on TR production
+  AliTRDsim *trdsim = trd->CreateTR();
+  trdsim->SetNFoils(100);
+  trdsim->SetFoilThick(0.0013);
+  trdsim->SetGapThick(0.0500);
+
+  // Initialize the TRD object
+  trd->Init();
+
+  // The digitizer
+  AliTRDdigitizer *digitizer = new AliTRDdigitizer("TRDdigitizer","Digitizer class");
+  digitizer->SetSimple();
+  digitizer->SetGasGain(3200.);      // GEANT + new Ermilova spectrum
+  digitizer->SetChipGain(6.1);
+  digitizer->SetNoise(1000.);
+  digitizer->SetADCinRange(1000.);
+  digitizer->SetADCoutRange(1023.);
+  digitizer->SetADCthreshold(0);
+  digitizer->InitDetector();
+  //digitizer->SetTimeResponse(0);
+  //digitizer->SetVerbose(1);
+                                    
+  // The event generator
+  AliTRDsimpleGen *generator = simple->GetGenerator();         
+  generator->SetMomentum(3.0,3.0);
+  generator->SetPdg(11);                             // Electrons 
+  //generator->SetPdg(211);                            // Pions 
+
+  //_____________________________________________________________________________ 
+  //
+  // Histograms
+  //_____________________________________________________________________________ 
+  //
+
+  Int_t timeMax  = geometry->GetTimeTotal();   
+  Int_t adcRange = ((Int_t) digitizer->GetADCoutRange()); 
+
+  TH1F     *hQ       = new TH1F("hQ"    ,"Charge per hit (all)" ,100,0.0,1000.0);
+  TH1F     *hQdedx   = new TH1F("hQdedx","Charge per hit (dedx)",100,0.0,1000.0);
+  TH1F     *hQtr     = new TH1F("hQtr  ","Charge per hit (tr)"  ,100,0.0,1000.0);
+
+  TProfile *hQX      = new TProfile("hQX"    ,"Charge per hit vs X (all)" ,35,0.0,3.5);
+  TProfile *hQXdedx  = new TProfile("hQXdedx","Charge per hit vs X (dedx)",35,0.0,3.5);
+  TProfile *hQXtr    = new TProfile("hQXtr  ","Charge per hit vs X (tr)"  ,35,0.0,3.5);
+
+  TH1F     *hNstep   = new TH1F("hNstep","Number of steps / cm"    , 151,-0.5, 150.5);
+  TH1F     *hNel     = new TH1F("hNel"  ,"Number of electrons / cm",1001,-0.5,1000.5);
+
+  TH1F     *hAmp     = new TH1F("hAmp","Amplitude of the digits"
+                                      ,adcRange+1,-0.5,((Float_t) adcRange)+0.5);         
+  TProfile *hAmpTime = new TProfile("hAmpTime","Amplitude vs timebin"
+                                             ,timeMax,-0.5,((Float_t) timeMax)-0.5);
+
+  //_____________________________________________________________________________ 
+  //
+  // Event loop
+  //_____________________________________________________________________________ 
+  //
+
+  // Number of events (tracks)
+  Int_t nEvent = 10000;
+
+  Float_t x0 = geometry->GetTime0(0) - AliTRDgeometry::DrThick(); 
+
+  TClonesArray *hitsArray = trd->Hits();
+
+  for (Int_t iEvent = 0; iEvent < nEvent; iEvent++) {
+
+    if (!(iEvent % 100) && (iEvent)) {
+      printf("Event no. %d\n",iEvent);
+    }
+
+    // Generate the hits for one track
+    simple->ProcessEvent(iEvent);
+
+    // Analyze the hits
+    Float_t nElectrons = 0.0;
+    for (Int_t iHit = 0; iHit < hitsArray->GetEntries(); iHit++) {
+
+      AliTRDhit *hit = (AliTRDhit *) hitsArray->UncheckedAt(iHit);
+      Int_t   charge = TMath::Abs(hit->GetCharge());
+      Float_t x      = hit->X() - x0;
+      hQ->Fill(charge);
+      hQX->Fill(x,charge);
+      if (hit->FromDrift() ||
+          hit->FromAmplification()) {
+        hQdedx->Fill(charge);
+        hQXdedx->Fill(x,charge);
+        nElectrons += charge;
+      }
+      if (hit->FromTRphoton()) {
+        hQtr->Fill(charge);
+        hQXtr->Fill(x,charge);
+      }
+
+    }
+
+    hNstep->Fill(((AliTRDsimpleMC *) gMC)->GetNStep() / 3.5);
+    hNel->Fill(nElectrons / 3.5);
+
+    // Digitize the hits
+    digitizer->MakeDigits();
+
+    // Analyze the digits
+    Int_t row =  8;
+    Int_t col = 48;
+    AliTRDdigitsManager *digitsManager = digitizer->Digits();
+    AliTRDdataArrayI    *digitsArray   = digitsManager->GetDigits(12);
+    for (Int_t time = 0; time < timeMax; time++) {  
+
+      Int_t amp = digitsArray->GetDataUnchecked(row,col,time);
+      if (amp > 0) {
+        hAmp->Fill((Double_t) amp);
+        hAmpTime->Fill((Double_t) time, (Double_t) amp);
+      }
+
+    }
+    
+  }
+
+  //_____________________________________________________________________________ 
+  //
+  // Plot the spectra
+  //_____________________________________________________________________________ 
+  //
+
+  TCanvas *cHit = new TCanvas("cHit","Hits",50,50,800,600);
+  cHit->Divide(3,2);
+
+  cHit->cd(1);
+  gPad->SetLogy();
+  hQ->SetLineColor(4);
+  hQ->SetXTitle("Q (a.u.)");
+  hQ->SetYTitle("Entries");
+  hQ->Draw();
+  hQdedx->SetLineColor(3);
+  hQdedx->Draw("SAME");
+  hQtr->SetLineColor(2);
+  hQtr->Draw("SAME"); 
+
+  cHit->cd(2); 
+  gPad->SetLogy();
+  hQX->SetLineColor(4);
+  hQX->SetXTitle("x (cm)");
+  hQX->SetYTitle("<Q> (a.u.)");
+  hQX->Draw("HIST");
+  hQXdedx->SetLineColor(3);
+  hQXdedx->Draw("SAMEHIST");
+
+  cHit->cd(3); 
+  gPad->SetLogy();
+  hQXtr->SetLineColor(2);
+  hQXtr->SetXTitle("x (cm)");
+  hQXtr->SetYTitle("<Q> (a.u.)");
+  hQXtr->Draw("HIST"); 
+
+  cHit->cd(4);
+  hNstep->SetLineColor(4);
+  hNstep->SetXTitle("N_{step}  / cm");
+  hNstep->SetYTitle("Entries");
+  hNstep->Draw();
+
+  cHit->cd(5);
+  hNel->SetLineColor(4);
+  hNel->SetXTitle("N_{electron} / cm");
+  hNel->SetYTitle("Entries");
+  hNel->Draw();
+
+  TCanvas *cDigit = new TCanvas("cDigit","Digits",50,50,600,400);
+  cDigit->Divide(2,1);
+
+  cDigit->cd(1);
+  gPad->SetLogy();
+  hAmp->SetLineColor(2);
+  hAmp->SetXTitle("ADC channels");
+  hAmp->SetYTitle("Entries");
+  hAmp->Draw();
+
+  cDigit->cd(2);
+  hAmpTime->SetLineColor(2);
+  hAmpTime->SetXTitle("Timebin");
+  hAmpTime->SetYTitle("<ADC channels>");
+  hAmpTime->Draw("HIST");
+
+  //_____________________________________________________________________________ 
+  //
+  // Save the histograms
+  //_____________________________________________________________________________ 
+  //
+
+  TFile *fileOut = new TFile("simple.root","RECREATE");
+
+  hQ->Write();
+  hQdedx->Write();
+  hQtr->Write();
+
+  hQX->Write();
+  hQXdedx->Write();
+  hQXtr->Write();
+
+  hNstep->Write();
+  hNel->Write();
+
+  hAmp->Write();
+  hAmpTime->Write();
+
+  fileOut->Close();
+
+}