3a179553604e6dd280f4d8174830d1881e0cf5e2
[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   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 }