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