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