TRD/simple.C

   1
   2 #if !defined(__CINT__) || defined(__MAKECINT__)
   3
   4 #include "TH1.h"
   5 #include "TFile.h"
   6 #include "TCanvas.h"
   7 #include "TProfile.h"
   8 #include "TClonesArray.h"
   9
  10 #include "AliModule.h"
  11 #include "AliFRAMEv1.h"
  12
  13 #include "AliTRDv1.h"
  14 #include "AliTRDhit.h"
  15 #include "AliTRDsim.h"
  16 #include "AliTRDsimple.h"
  17 #include "AliTRDsimpleGen.h"
  18 #include "AliTRDsimpleMC.h"
  19 #include "AliTRDgeometry.h"
  20 #include "AliTRDdigitizer.h"
  21 #include "AliTRDdigitsManager.h"
  22 #include "AliTRDdataArrayI.h"
  23
  24 #endif
  25
  26 void simple()
  27 {
  28
  29   ///////////////////////////////////////////////////////////////////////////////
  30   //
  31   //  Macro to run the simplified simulator
  32   //
  33   ///////////////////////////////////////////////////////////////////////////////
  34
  35   //_____________________________________________________________________________
  36   //
  37   // Initialization part
  38   //_____________________________________________________________________________
  39   //
  40
  41   // The simple simulator
  42   AliTRDsimple *simple = new AliTRDsimple();
  43
  44   // Initialize a dummy frame so that the TRD works
  45   new AliFRAMEv1("FRAME","Space Frame");
  46
  47   // Initialize the TRD detector
  48   AliTRDv1     *trd    = new AliTRDv1("TRD","TRD slow simulator");
  49
  50   // Needed for some material properties
  51   trd->CreateMaterials();
  52
  53   // Use PAI model (most propable energy loss) for the relativistic rise
  54   //trd->SetPAI(1);
  55   // Use GEANT (mean energy loss) for the relativistic rise
  56   trd->SetPAI(0);
  57
  58   // Use the Ermilova spectrum for the delta electrons
  59   trd->SetErmilova(1);
  60   // Use the GEANT spectrum for the delta electrons
  61   //trd->SetErmilova(0);
  62
  63   // Select the gas mixture (0: 97% Xe + 3% isobutane, 1: 85% Xe + 15% CO2)
  64   trd->SetGasMix(1);
  65
  66   // Get the pointer to the geometry object
  67   AliTRDgeometry *geometry = trd->GetGeometry();
  68
  69   // The number of timebins
  70   geometry->SetNTimeBin(30);
  71
  72   // The additional timebins before and after the drift region
  73   geometry->SetExpandTimeBin(5,15);
  74
  75   // Switch on TR production
  76   AliTRDsim *trdsim = trd->CreateTR();
  77   trdsim->SetNFoils(100);
  78   trdsim->SetFoilThick(0.0013);
  79   trdsim->SetGapThick(0.0500);
  80
  81   // Initialize the TRD object
  82   trd->Init();
  83
  84   // The digitizer
  85   AliTRDdigitizer *digitizer = new AliTRDdigitizer("TRDdigitizer","Digitizer class");
  86   digitizer->SetSimple();
  87   digitizer->SetGasGain(3200.);      // GEANT + new Ermilova spectrum
  88   digitizer->SetChipGain(6.1);
  89   digitizer->SetNoise(1000.);
  90   digitizer->SetADCinRange(1000.);
  91   digitizer->SetADCoutRange(1023.);
  92   digitizer->SetADCthreshold(0);
  93   digitizer->InitDetector();
  94   //digitizer->SetTimeResponse(0);
  95   //digitizer->SetVerbose(1);
  96
  97   // The event generator
  98   AliTRDsimpleGen *generator = simple->GetGenerator();
  99   generator->SetMomentum(3.0,3.0);
 100   generator->SetPdg(11);                             // Electrons
 101   //generator->SetPdg(211);                            // Pions
 102
 103   //_____________________________________________________________________________
 104   //
 105   // Histograms
 106   //_____________________________________________________________________________
 107   //
 108
 109   Int_t timeMax  = geometry->GetTimeTotal();
 110   Int_t adcRange = ((Int_t) digitizer->GetADCoutRange());
 111
 112   TH1F     *hQ       = new TH1F("hQ"    ,"Charge per hit (all)" ,100,0.0,1000.0);
 113   TH1F     *hQdedx   = new TH1F("hQdedx","Charge per hit (dedx)",100,0.0,1000.0);
 114   TH1F     *hQtr     = new TH1F("hQtr  ","Charge per hit (tr)"  ,100,0.0,1000.0);
 115
 116   TProfile *hQX      = new TProfile("hQX"    ,"Charge per hit vs X (all)" ,35,0.0,3.5);
 117   TProfile *hQXdedx  = new TProfile("hQXdedx","Charge per hit vs X (dedx)",35,0.0,3.5);
 118   TProfile *hQXtr    = new TProfile("hQXtr  ","Charge per hit vs X (tr)"  ,35,0.0,3.5);
 119
 120   TH1F     *hNstep   = new TH1F("hNstep","Number of steps / cm"    , 151,-0.5, 150.5);
 121   TH1F     *hNel     = new TH1F("hNel"  ,"Number of electrons / cm",1001,-0.5,1000.5);
 122
 123   TH1F     *hAmp     = new TH1F("hAmp","Amplitude of the digits"
 124                                       ,adcRange+1,-0.5,((Float_t) adcRange)+0.5);
 125   TProfile *hAmpTime = new TProfile("hAmpTime","Amplitude vs timebin"
 126                                               ,timeMax,-0.5,((Float_t) timeMax)-0.5);
 127
 128   //_____________________________________________________________________________
 129   //
 130   // Event loop
 131   //_____________________________________________________________________________
 132   //
 133
 134   // Number of events (tracks)
 135   Int_t nEvent = 10000;
 136
 137   Float_t x0 = geometry->GetTime0(0) - AliTRDgeometry::DrThick();
 138
 139   TClonesArray *hitsArray = trd->Hits();
 140
 141   for (Int_t iEvent = 0; iEvent < nEvent; iEvent++) {
 142
 143     if (!(iEvent % 100) && (iEvent)) {
 144       printf("Event no. %d\n",iEvent);
 145     }
 146
 147     // Generate the hits for one track
 148     simple->ProcessEvent(iEvent);
 149
 150     // Analyze the hits
 151     Float_t nElectrons = 0.0;
 152     for (Int_t iHit = 0; iHit < hitsArray->GetEntries(); iHit++) {
 153
 154       AliTRDhit *hit = (AliTRDhit *) hitsArray->UncheckedAt(iHit);
 155       Int_t   charge = TMath::Abs(hit->GetCharge());
 156       Float_t x      = hit->X() - x0;
 157       hQ->Fill(charge);
 158       hQX->Fill(x,charge);
 159       if (hit->FromDrift() ||
 160           hit->FromAmplification()) {
 161         hQdedx->Fill(charge);
 162         hQXdedx->Fill(x,charge);
 163         nElectrons += charge;
 164       }
 165       if (hit->FromTRphoton()) {
 166         hQtr->Fill(charge);
 167         hQXtr->Fill(x,charge);
 168       }
 169
 170     }
 171
 172     hNstep->Fill(((AliTRDsimpleMC *) gMC)->GetNStep() / 3.5);
 173     hNel->Fill(nElectrons / 3.5);
 174
 175     // Digitize the hits
 176     digitizer->MakeDigits();
 177
 178     // Analyze the digits
 179     Int_t row =  8;
 180     Int_t col = 48;
 181     AliTRDdigitsManager *digitsManager = digitizer->Digits();
 182     AliTRDdataArrayI    *digitsArray   = digitsManager->GetDigits(12);
 183     for (Int_t time = 0; time < timeMax; time++) {
 184
 185       Int_t amp = digitsArray->GetDataUnchecked(row,col,time);
 186       if (amp > 0) {
 187         hAmp->Fill((Double_t) amp);
 188         hAmpTime->Fill((Double_t) time, (Double_t) amp);
 189       }
 190
 191     }
 192
 193   }
 194
 195   //_____________________________________________________________________________
 196   //
 197   // Plot the spectra
 198   //_____________________________________________________________________________
 199   //
 200
 201   TCanvas *cHit = new TCanvas("cHit","Hits",50,50,800,600);
 202   cHit->Divide(3,2);
 203
 204   cHit->cd(1);
 205   gPad->SetLogy();
 206   hQ->SetLineColor(4);
 207   hQ->SetXTitle("Q (a.u.)");
 208   hQ->SetYTitle("Entries");
 209   hQ->Draw();
 210   hQdedx->SetLineColor(3);
 211   hQdedx->Draw("SAME");
 212   hQtr->SetLineColor(2);
 213   hQtr->Draw("SAME");
 214
 215   cHit->cd(2);
 216   gPad->SetLogy();
 217   hQX->SetLineColor(4);
 218   hQX->SetXTitle("x (cm)");
 219   hQX->SetYTitle("<Q> (a.u.)");
 220   hQX->Draw("HIST");
 221   hQXdedx->SetLineColor(3);
 222   hQXdedx->Draw("SAMEHIST");
 223
 224   cHit->cd(3);
 225   gPad->SetLogy();
 226   hQXtr->SetLineColor(2);
 227   hQXtr->SetXTitle("x (cm)");
 228   hQXtr->SetYTitle("<Q> (a.u.)");
 229   hQXtr->Draw("HIST");
 230
 231   cHit->cd(4);
 232   hNstep->SetLineColor(4);
 233   hNstep->SetXTitle("N_{step}  / cm");
 234   hNstep->SetYTitle("Entries");
 235   hNstep->Draw();
 236
 237   cHit->cd(5);
 238   hNel->SetLineColor(4);
 239   hNel->SetXTitle("N_{electron} / cm");
 240   hNel->SetYTitle("Entries");
 241   hNel->Draw();
 242
 243   TCanvas *cDigit = new TCanvas("cDigit","Digits",50,50,600,400);
 244   cDigit->Divide(2,1);
 245
 246   cDigit->cd(1);
 247   gPad->SetLogy();
 248   hAmp->SetLineColor(2);
 249   hAmp->SetXTitle("ADC channels");
 250   hAmp->SetYTitle("Entries");
 251   hAmp->Draw();
 252
 253   cDigit->cd(2);
 254   hAmpTime->SetLineColor(2);
 255   hAmpTime->SetXTitle("Timebin");
 256   hAmpTime->SetYTitle("<ADC channels>");
 257   hAmpTime->Draw("HIST");
 258
 259   //_____________________________________________________________________________
 260   //
 261   // Save the histograms
 262   //_____________________________________________________________________________
 263   //
 264
 265   TFile *fileOut = new TFile("simple.root","RECREATE");
 266
 267   hQ->Write();
 268   hQdedx->Write();
 269   hQtr->Write();
 270
 271   hQX->Write();
 272   hQXdedx->Write();
 273   hQXtr->Write();
 274
 275   hNstep->Write();
 276   hNel->Write();
 277
 278   hAmp->Write();
 279   hAmpTime->Write();
 280
 281   fileOut->Close();
 282
 283 }