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