6658c03d591a42b71b15edeff49abb7c79c65632
[u/mrichter/AliRoot.git] / 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 }