Update of the test suite
[u/mrichter/AliRoot.git] / TRD / AliTRDtest.C
CommitLineData
8e64dd77 1
793ff80c 2Int_t AliTRDtest()
3{
4 //
5 // Test macro for the TRD code
6 //
7
8 Int_t rc = 0;
9
8e64dd77 10 TFile *file;
11
793ff80c 12 // Initialize the test setup
5931c105 13 gAlice->Init("$(ALICE_ROOT)/TRD/AliTRDconfig.C");
793ff80c 14
15 // Run one event and create the hits
6244debe 16 gAlice->SetDebug(2);
793ff80c 17 gAlice->Run(1);
18
6244debe 19 if (gAlice) delete gAlice;
8e64dd77 20 file = (TFile *) gROOT->GetListOfFiles()->FindObject("TRD_test.root");
6244debe 21 gAlice = (AliRun *) file->Get("gAlice");
22
793ff80c 23 // Analyze the TRD hits
793ff80c 24 if (rc = AliTRDanalyzeHits()) return rc;
25
26 // Run the digitization
793ff80c 27 if (rc = AliTRDcreateDigits()) return rc;
28
6244debe 29 if (gAlice) delete gAlice;
8e64dd77 30 file = (TFile *) gROOT->GetListOfFiles()->FindObject("TRD_test.root");
6244debe 31 gAlice = (AliRun *) file->Get("gAlice");
32
793ff80c 33 // Analyze the digits
793ff80c 34 if (rc = AliTRDanalyzeDigits()) return rc;
35
4a0fe73c 36 // Create the cluster
4a0fe73c 37 if (rc = AliTRDcreateCluster()) return rc;
38
39 if (gAlice) delete gAlice;
8e64dd77 40 file = (TFile *) gROOT->GetListOfFiles()->FindObject("TRD_test.root");
4a0fe73c 41 gAlice = (AliRun *) file->Get("gAlice");
42
8e64dd77 43 // Analyze the cluster
4a0fe73c 44 if (rc = AliTRDanalyzeCluster()) return rc;
45
8e64dd77 46 //file = (TFile *) gROOT->GetListOfFiles()->FindObject("TRD_test.root");
47 //file->Close();
48
49 // Find the tracks
50 //if (rc = AliTRDcreateTracks()) return rc;
51
793ff80c 52 return rc;
53
54}
fa148e6c 55
56//_____________________________________________________________________________
57Int_t AliTRDanalyzeHits()
58{
59 //
60 // Analyzes the hits and fills QA-histograms
61 //
62
63 Int_t rc = 0;
64
65 if (!gAlice) {
66 cout << "<AliTRDanalyzeHits> No AliRun object found" << endl;
67 rc = 1;
68 return rc;
69 }
70 gAlice->GetEvent(0);
71
72 // Get the pointer to the TRD detector
73 AliTRD *trd = (AliTRD *) gAlice->GetDetector("TRD");
74 if (!trd) {
75 cout << "<AliTRDanalyzeHits> No TRD detector found" << endl;
76 rc = 2;
77 return rc;
78 }
79
80 // Get the pointer to the geometry object
81 AliTRDgeometryFull *geo;
82 if (trd) {
83 geo = (AliTRDgeometryFull *) trd->GetGeometry();
84 }
85 else {
86 cout << "<AliTRDanalyzeHits> No TRD geometry found" << endl;
87 rc = 3;
88 return rc;
89 }
90
91 AliTRDparameter *par = new AliTRDparameter("TRDparameter","TRD parameter class");
92
93 // Define the histograms
94 TH1F *hQdedx = new TH1F("hQdedx","Charge dedx-hits",100,0.0,1000.0);
95 TH1F *hQtr = new TH1F("hQtr" ,"Charge TR-hits" ,100,0.0,1000.0);
96
97 Float_t rmin = geo->Rmin();
98 Float_t rmax = geo->Rmax();
99 Float_t length = geo->GetChamberLength(0,2);
100 Float_t width = geo->GetChamberWidth(0);
101 Int_t ncol = par->GetColMax(0);
102 Int_t nrow = par->GetRowMax(0,2,13);
103 Int_t ntime = ((Int_t) (rmax - rmin) / par->GetTimeBinSize());
104
105 TH2F *hZY = new TH2F("hZY" ,"Y vs Z (chamber 0)", nrow,-length/2.,length/2.
106 ,ntime, rmin, rmax);
107 TH2F *hXZ = new TH2F("hXZ" ,"Z vs X (plane 0)" , ncol, -width/2., width/2.
108 , nrow,-length/2.,length/2.);
109
110 TH1F *hNprimP = new TH1F("hNprimP","Number of primary electrons per cm (MIP pion)"
111 ,50,0.0,100.0);
112 TH1F *hNprimE = new TH1F("hNprimE","Number of primary electrons per cm (3GeV electron)"
113 ,50,0.0,100.0);
114 TH1F *hNtotP = new TH1F("hNtotP" ,"Number of electrons per cm (MIP pion)"
115 ,50,0.0,1000.0);
116 TH1F *hNtotE = new TH1F("hNtotE" ,"Number of electrons per cm (3GeV electron)"
117 ,50,0.0,1000.0);
118
119 // Get the pointer hit tree
120 TTree *hitTree = gAlice->TreeH();
121 if (!hitTree) {
122 cout << "<AliTRDanalyzeHits> No hit tree found" << endl;
123 rc = 4;
124 return rc;
125 }
126
127 Int_t countHits = 0;
128 Int_t nBytes = 0;
129
130 // Get the number of entries in the hit tree
131 // (Number of primary particles creating a hit somewhere)
132 Int_t nTrack = (Int_t) hitTree->GetEntries();
133 cout << "<AliTRDanalyzeHits> Found " << nTrack
134 << " primary particles with hits" << endl;
135
136 // Loop through all entries in the tree
137 for (Int_t iTrack = 0; iTrack < nTrack; iTrack++) {
138
139 gAlice->ResetHits();
140 nBytes += hitTree->GetEvent(iTrack);
141
142 Int_t nPrimE = 0;
143 Int_t nPrimP = 0;
144 Int_t nTotE = 0;
145 Int_t nTotP = 0;
146
147 // Loop through the TRD hits
148 Int_t iHit = 0;
149 AliTRDhit *hit = (AliTRDhit *) trd->FirstHit(-1);
150 while (hit) {
151
152 countHits++;
153 iHit++;
154
155 Float_t x = hit->X();
156 Float_t y = hit->Y();
157 Float_t z = hit->Z();
158 Float_t q = hit->GetCharge();
159 Int_t track = hit->Track();
160 Int_t det = hit->GetDetector();
161 Int_t plane = geo->GetPlane(det);
162
163 if (q > 0) {
164 hQdedx->Fill(q);
165 hZY->Fill(z,y);
166 if (plane == 0) {
167 hXZ->Fill(x,z);
168 }
169 }
170 else if (q < 0) {
171 hQtr->Fill(TMath::Abs(q));
172 }
173
174 TParticle *part = gAlice->Particle(track);
175
176 if ((plane == 0) && (q > 0)) {
177
178 // 3 GeV electrons
179 if ((part->GetPdgCode() == 11) &&
180 (part->P() > 2.9)) {
181 nPrimE++;
182 nTotE += ((Int_t) q);
183 }
184
185 // MIP pions
186 if ((part->GetPdgCode() == -211) &&
187 (part->P() > 0.5)) {
188 nPrimP++;
189 nTotP += ((Int_t) q);
190 }
191
192 }
193
194 hit = (AliTRDhit *) trd->NextHit();
195
196 }
197
198 if (nPrimE > 0) hNprimE->Fill(((Double_t) nPrimE)/3.7);
199 if (nPrimP > 0) hNprimP->Fill(((Double_t) nPrimP)/3.7);
200 if (nTotE > 0) hNtotE->Fill(((Double_t) nTotE)/3.7);
201 if (nTotP > 0) hNtotP->Fill(((Double_t) nTotP)/3.7);
202
203 }
204
205 cout << "<AliTRDanalyzeHits> Found " << countHits << " hits in total" << endl;
206
207 TCanvas *cHits = new TCanvas("cHits","AliTRDanalyzeHits",50,50,600,600);
208 cHits->Divide(2,2);
209 cHits->cd(1);
210 hXZ->Draw("COL");
211 cHits->cd(2);
212 hZY->Draw("COL");
213 cHits->cd(3);
214 gPad->SetLogy();
215 hQdedx->Draw();
216 cHits->cd(4);
217 gPad->SetLogy();
218 hQtr->Draw();
219
220 TCanvas *cNel = new TCanvas("cNel","Ionization",50,50,600,600);
221 cNel->Divide(2,2);
222 cNel->cd(1);
223 hNprimE->SetStats();
224 hNprimE->Draw();
225 cNel->cd(2);
226 hNprimP->SetStats();
227 hNprimP->Draw();
228 cNel->cd(3);
229 hNtotE->SetStats();
230 hNtotE->Draw();
231 cNel->cd(4);
232 hNtotP->SetStats();
233 hNtotP->Draw();
234
235 TFile *fout = new TFile("TRD_ionization.root","RECREATE");
236 hNprimE->Write();
237 hNprimP->Write();
238 hNtotE->Write();
239 hNtotP->Write();
240 fout->Close();
241
242 return rc;
243
244}
245
246//_____________________________________________________________________________
247Int_t AliTRDcreateDigits()
248{
249 //
250 // Creates the digits from the hits of the slow simulator
251 //
252
253 Int_t rc = 0;
254
255 if (!gAlice) {
256 cout << "<AliTRDcreateDigits> No AliRun object found" << endl;
257 rc = 1;
258 return rc;
259 }
260 gAlice->GetEvent(0);
261
262 // Create the TRD digitzer
263 AliTRDdigitizer *digitizer = new AliTRDdigitizer("TRDdigitizer"
264 ,"TRD digitizer class");
265 digitizer->InitDetector();
266
267 // Set the parameter
268 digitizer->SetDebug(1);
269
270 // Define the parameter object
271 // If no external parameter object is defined,
272 // default parameter will be used
273 AliTRDparameter *parameter = new AliTRDparameter("TRDparameter"
274 ,"TRD parameter class");
275 digitizer->SetParameter(parameter);
276
277 // Create the digits
278 if (!(digitizer->MakeDigits())) {
279 rc = 2;
280 return rc;
281 }
282
283 // Write the digits into the input file
284 if (!(digitizer->MakeBranch())) {
285 rc = 3;
286 return rc;
287 }
288
289 // Write the digits into the input file
290 if (!(digitizer->WriteDigits())) {
291 rc = 4;
292 return rc;
293 }
294
295 // Save the parameter object in the AliROOT file
296 if (!(parameter->Write())) {
297 rc = 4;
298 return rc;
299 }
300
301 return rc;
302
303}
304
305//_____________________________________________________________________________
306Int_t AliTRDanalyzeDigits()
307{
308 //
309 // Analyzes the digits
310 //
311
312 Int_t rc = 0;
313
314 const Int_t kNpla = AliTRDgeometry::Nplan();
315
316 if (!gAlice) {
317 cout << "<AliTRDanalyzeDigits> No AliRun object found" << endl;
318 rc = 1;
319 return rc;
320 }
321 Int_t nPart = gAlice->GetEvent(0);
322
323 // Get the pointer to the TRD detector
324 AliTRD *trd = (AliTRD *) gAlice->GetDetector("TRD");
325 if (!trd) {
326 cout << "<AliTRDanalyzeDigits> No TRD detector found" << endl;
327 rc = 2;
328 return rc;
329 }
330
331 // Get the digitizer object
332 TFile *file = (TFile *) gROOT->GetListOfFiles()->FindObject("TRD_test.root");
333 AliTRDparameter *parameter = (AliTRDparameter *) file->Get("TRDparameter");
334 if (!parameter) {
335 cout << "<AliTRDanalyzeDigits> No parameter object found" << endl;
336 rc = 3;
337 return rc;
338 }
339
340 // Define the histograms
341 Int_t adcRange = ((Int_t) parameter->GetADCoutRange());
342 TH1F *hAmpAll = new TH1F("hAmpAll" ,"Amplitude of the digits (all)"
343 ,adcRange+1,-0.5,((Float_t) adcRange)+0.5);
344 TH1F *hAmpEl = new TH1F("hAmpEl" ,"Amplitude of the digits (electrons)"
345 ,adcRange+1,-0.5,((Float_t) adcRange)+0.5);
346 TH1F *hAmpPi = new TH1F("hAmpPi" ,"Amplitude of the digits (pions)"
347 ,adcRange+1,-0.5,((Float_t) adcRange)+0.5);
348 TH1F *hAmpNoise = new TH1F("hAmpNoise","Amplitude of the digits (noise)"
349 ,5,-0.5,4.5);
350
351 // Get the pointer to the geometry object
352 AliTRDgeometry *geo;
353 if (trd) {
354 geo = trd->GetGeometry();
355 }
356 else {
357 cout << "<AliTRDanalyzeDigits> No TRD geometry found" << endl;
358 rc = 4;
359 return rc;
360 }
361
362 // Create the digits manager
363 AliTRDdigitsManager *digitsManager = new AliTRDdigitsManager();
8e64dd77 364 digitsManager->SetDebug(1);
fa148e6c 365
366 digitsManager->Open("TRD_test.root");
367
368 // Read the digits from the file
369 if (!(digitsManager->ReadDigits())) {
370 cout << "<AliTRDanalyzeDigits> Cannot read the digits" << endl;
371 rc = 5;
372 return rc;
373 }
374
375 AliTRDmatrix *matrix;
376
377 Int_t countDigits = 0;
378 Int_t iSec = trd->GetSensSector();
379 Int_t iCha = trd->GetSensChamber();
f6ff4cec 380 Int_t timeMax = parameter->GetTimeTotal();
fa148e6c 381
382 TProfile *hAmpTimeEl = new TProfile("hAmpTimeEl","Amplitude of the digits (electrons)"
383 ,timeMax,-0.5,((Double_t) timeMax)-0.5);
384 TProfile *hAmpTimePi = new TProfile("hAmpTimePi","Amplitude of the digits (pions)"
385 ,timeMax,-0.5,((Double_t) timeMax)-0.5);
386
387 // Loop over all planes
388 for (Int_t iPla = 0; iPla < kNpla; iPla++) {
389
390 Int_t iDet = geo->GetDetector(iPla,iCha,iSec);
391 Int_t rowMax = parameter->GetRowMax(iPla,iCha,iSec);
392 Int_t colMax = parameter->GetColMax(iPla);
393
394 if (iPla == 0) {
395 matrix = new AliTRDmatrix(rowMax,colMax,timeMax,iSec,iCha,iPla);
396 }
397
398 // Loop through the detector pixel
399 for (Int_t time = 0; time < timeMax; time++) {
400 for (Int_t col = 0; col < colMax; col++) {
401 for (Int_t row = 0; row < rowMax; row++) {
402
403 AliTRDdigit *digit = digitsManager->GetDigit(row,col,time,iDet);
404 Int_t amp = digit->GetAmp();
405 Int_t track0 = digitsManager->GetTrack(0,row,col,time,iDet);
406 Int_t track1 = digitsManager->GetTrack(1,row,col,time,iDet);
407 TParticle *particle = 0;
408 if (track0 > -1) {
409 particle = gAlice->Particle(track0);
410 }
411
412 if (amp > 0) {
413 countDigits++;
414 if (iPla == 0) {
415 matrix->SetSignal(row,col,time,amp);
416 }
417 }
418
419 // Total spectrum
420 hAmpAll->Fill(amp);
421
422 // Noise spectrum
423 if (track0 < 0) {
424 hAmpNoise->Fill(amp);
425 }
426
427 // Electron digit
428 if ((particle) && (particle->GetPdgCode() == 11) && (track1 < 0)) {
429 hAmpEl->Fill(amp);
430 hAmpTimeEl->Fill(time,amp);
431 }
432
433 // Pion digit
434 if ((particle) && (particle->GetPdgCode() == -211) && (track1 < 0)) {
435 hAmpPi->Fill(amp);
436 hAmpTimePi->Fill(time,amp);
437 }
438
439 delete digit;
440
441 }
442 }
443 }
444
445 }
446
447 cout << "<AliTRDanalyzeDigits> Found " << countDigits << " digits in total" << endl;
448
449 // Display the detector matrix
450 matrix->Draw();
451
452 TCanvas *cDigits = new TCanvas("cDigits","AliTRDanalyzeDigits",50,50,600,800);
453 cDigits->Divide(2,3);
454 cDigits->cd(1);
455 gPad->SetLogy();
456 hAmpAll->SetXTitle("Amplitude (ADC-channels)");
457 hAmpAll->SetYTitle("Entries");
458 hAmpAll->Draw();
459 cDigits->cd(2);
460 gPad->SetLogy();
461 hAmpNoise->SetXTitle("Amplitude (ADC-channels)");
462 hAmpNoise->SetYTitle("Entries");
463 hAmpNoise->Draw();
464 cDigits->cd(3);
465 gPad->SetLogy();
466 hAmpEl->SetXTitle("Amplitude (ADC-channels)");
467 hAmpEl->SetYTitle("Entries");
468 hAmpEl->Draw();
469 cDigits->cd(4);
470 gPad->SetLogy();
471 hAmpPi->SetXTitle("Amplitude (ADC-channels)");
472 hAmpPi->SetYTitle("Entries");
473 hAmpPi->Draw();
474 cDigits->cd(5);
475 hAmpTimeEl->SetXTitle("Timebin number");
476 hAmpTimeEl->SetYTitle("Mean amplitude");
477 hAmpTimeEl->Draw("HIST");
478 cDigits->cd(6);
479 hAmpTimePi->SetXTitle("Timebin number");
480 hAmpTimePi->SetYTitle("Mean amplitude");
481 hAmpTimePi->Draw("HIST");
482
483 TFile *fileOut = new TFile("digits_test.root","RECREATE");
484 hAmpAll->Write();
485 hAmpNoise->Write();
486 hAmpEl->Write();
487 hAmpPi->Write();
488 hAmpTimeEl->Write();
489 hAmpTimePi->Write();
490 fileOut->Close();
491
492 return rc;
493
494}
495
496//_____________________________________________________________________________
497Int_t AliTRDcreateCluster()
498{
499 //
500 // Creates the cluster from the digits
501 //
502
503 Int_t rc = 0;
504
505 // Create the clusterizer
506 AliTRDclusterizerV1 *clusterizer = new AliTRDclusterizerV1("TRDclusterizer"
507 ,"TRD clusterizer class");
508 clusterizer->SetVerbose(1);
509
510 // Define the parameter object
511 // If no external parameter object is defined,
512 // default parameter will be used
513 AliTRDparameter *parameter = new AliTRDparameter("TRDparameter"
514 ,"TRD parameter class");
515 parameter->SetClusMaxThresh(0);
516 parameter->SetClusSigThresh(0);
517 clusterizer->SetParameter(parameter);
518
519 // Open the file
520 if (!(clusterizer->Open("TRD_test.root",0))) {
521 rc = 1;
522 return rc;
523 }
524
525 // Load the digits
526 if (!(clusterizer->ReadDigits())) {
527 rc = 2;
528 return rc;
529 }
530
531 // Find the cluster
532 if (!(clusterizer->MakeClusters())) {
533 rc = 3;
534 return rc;
535 }
536
537 // Write the cluster tree into the file
538 if (!(clusterizer->WriteClusters(-1))) {
539 rc = 4;
540 return rc;
541 }
542
543 // Save the clusterizer class in the file
544 if (!(clusterizer->Write())) {
545 rc = 5;
546 return rc;
547 }
548
549 return rc;
550
551}
552
553//_____________________________________________________________________________
554Int_t AliTRDanalyzeCluster()
555{
556 //
557 // Analyzes the cluster
558 //
559
560 Int_t rc = 0;
561
562 if (!gAlice) {
563 cout << "<AliTRDanalyzeCluster> No AliRun object found" << endl;
564 rc = 1;
565 return rc;
566 }
567 gAlice->GetEvent(0);
568
569 // Get the pointer to the TRD detector
570 AliTRD *trd = (AliTRD *) gAlice->GetDetector("TRD");
571 if (!trd) {
572 cout << "<AliTRDanalyzeCluster> No TRD detector found" << endl;
573 rc = 2;
574 return rc;
575 }
576
577 // Define the histograms
578 TH1F *hClusAll = new TH1F("hClusAll" ,"Amplitude of the cluster (all)"
579 ,501,-0.5,500.5);
580 TH1F *hClusNoise = new TH1F("hClusNoise","Amplitude of the cluster (noise)"
581 , 5,-0.5, 4.5);
582 TH1F *hClusEl = new TH1F("hClusEl" ,"Amplitude of the cluster (electron)"
583 ,501,-0.5,500.5);
584 TH1F *hClusPi = new TH1F("hClusPi" ,"Amplitude of the cluster (pion)"
585 ,501,-0.5,500.5);
586
587 // Get the pointer to the geometry object
588 AliTRDgeometry *geo;
589 if (trd) {
590 geo = trd->GetGeometry();
591 }
592 else {
593 cout << "<AliTRDanalyzeCluster> No TRD geometry found" << endl;
594 rc = 3;
595 return rc;
596 }
597
598 // Get the pointer to the hit-tree
599 TFile *file = (TFile *) gROOT->GetListOfFiles()->FindObject("TRD_test.root");
600 TTree *clusterTree = (TTree *) file->Get("TreeR0_TRD");
601 if (!(clusterTree)) {
602 cout << "<AliTRDanalyzeCluster> No tree with clusters found" << endl;
603 rc = 4;
604 return rc;
605 }
606
607 // Get the pointer to the hit container
608 TObjArray *clusterArray = trd->RecPoints();
609 if (!(clusterArray)) {
610 cout << "<AliTRDanalyzeCluster> No clusterArray found" << endl;
611 rc = 5;
612 return rc;
613 }
614
615 // Set the branch address
616 clusterTree->GetBranch("TRDcluster")->SetAddress(&clusterArray);
617 Int_t nEntries = clusterTree->GetEntries();
618 cout << "<AliTRDanalyzeCluster> Number of entries in the cluster tree = "
619 << nEntries
620 << endl;
621
622 Int_t countCluster = 0;
623 Int_t countOverlap = 0;
624
625 // Loop through all entries in the tree
626 Int_t nbytes;
627 for (Int_t iEntry = 0; iEntry < nEntries; iEntry++) {
628
629 // Import the tree
630 nbytes += clusterTree->GetEvent(iEntry);
631
632 // Get the number of points in the detector
633 Int_t nCluster = clusterArray->GetEntriesFast();
634
635 // Loop through all TRD digits
636 for (Int_t iCluster = 0; iCluster < nCluster; iCluster++) {
637
638 // Get the information for this digit
639 AliTRDcluster *cluster = (AliTRDcluster *) clusterArray->UncheckedAt(iCluster);
640 Int_t detector = cluster->GetDetector();
641 Int_t sector = geo->GetSector(detector);
642 Int_t plane = geo->GetPlane(detector);
643 Int_t chamber = geo->GetChamber(detector);
644 Float_t energy = cluster->GetQ();
f6ff4cec 645 Int_t track0 = cluster->GetLabel(0);
646 Int_t track1 = cluster->GetLabel(1);
647 Int_t track2 = cluster->GetLabel(2);
fa148e6c 648 TParticle *particle = 0;
649 if (track0 > -1) {
650 particle = gAlice->Particle(track0);
651 }
652
653 countCluster++;
654 if (!cluster->Isolated()) countOverlap++;
655
656 // Total spectrum
657 hClusAll->Fill(energy);
658
659 if (cluster->Isolated()) {
660
661 // Noise spectrum
662 if (track0 < 0) {
663 hClusNoise->Fill(energy);
664 }
665
666 // Electron cluster
667 if ((particle) && (particle->GetPdgCode() == 11) && (track1 < 0)) {
668 hClusEl->Fill(energy);
669 }
670
671 // Pion cluster
672 if ((particle) && (particle->GetPdgCode() == -211) && (track1 < 0)) {
673 hClusPi->Fill(energy);
674 }
675
676 }
677
678 }
679
680 }
681
682 cout << "<AliTRDanalyzeCluster> Found " << countCluster << " cluster in total" << endl;
683 cout << "<AliTRDanalyzeCluster> Found " << countOverlap << " overlapping cluster" << endl;
684 cout << endl;
685
686 TCanvas *cCluster = new TCanvas("cCluster","AliTRDanalyzeCluster",50,50,600,600);
687 cCluster->Divide(2,2);
688
689 TF1 *fun;
690 cCluster->cd(1);
691 gPad->SetLogy();
692 hClusAll->Fit("landau","0");
693 fun = (TF1 *) hClusAll->GetListOfFunctions()->First();
694 Float_t meanAll = fun->GetParameter(1);
695 hClusAll->Draw();
696 fun->SetLineColor(2);
697 fun->Draw("SAME");
698
699 //cCluster->cd(2);
700 //gPad->SetLogy();
701 Float_t meanNoise = hClusNoise->GetMean();
702 //hClusNoise->Draw();
703
704 cCluster->cd(3);
705 gPad->SetLogy();
706 hClusEl->Fit("landau","0");
707 fun = (TF1 *) hClusEl->GetListOfFunctions()->First();
708 fun->SetLineColor(2);
709 Float_t meanEl = fun->GetParameter(1);
710 hClusEl->Draw();
711 fun->Draw("SAME");
712
713 cCluster->cd(4);
714 gPad->SetLogy();
715 hClusPi->Fit("landau","0");
716 fun = (TF1 *) hClusPi->GetListOfFunctions()->First();
717 fun->SetLineColor(2);
718 Float_t meanPi = fun->GetParameter(1);
719 hClusPi->Draw();
720 fun->Draw("SAME");
721
722 cout << endl;
723 cout << "##################################################################" << endl;
724 cout << " Mean all = " << meanAll << endl;
725 cout << " Mean noise = " << meanNoise << endl;
726 cout << " Mean electrons = " << meanEl << endl;
727 cout << " Mean pions = " << meanPi << endl;
728 cout << "##################################################################" << endl;
729 cout << endl;
730
731 return rc;
732
733}
8e64dd77 734
735//_____________________________________________________________________________
736Int_t AliTRDcreateTracks()
737{
738 //
739 // Creates the tracks
740 //
741
742 Int_t rc = 0;
743
744 // Create the tracker
745 AliTRDtracker *tracker = new AliTRDtracker("TRDtracker","TRD tracker");
746
747 // Read in the kine tree and the cluster
748 tracker->GetEvent("TRD_test.root","TRD_test.root");
749
750 // Find the tracks
751 TH1F *hs = new TH1F("hs","hs",100,0.0,1.0);
752 TH1F *hd = new TH1F("hd","hd",100,0.0,1.0);
753 tracker->Clusters2Tracks(hs,hd);
754
755 // Store the tracks
756 tracker->WriteTracks("TRD_test.root");
757
758 return rc;
759
760}