Update of macros
[u/mrichter/AliRoot.git] / TRD / AliTRDtest.C
1 Int_t AliTRDtest() 
2 {
3   //
4   // Test macro for the TRD code
5   //
6
7   Int_t rc = 0;
8
9   // Initialize the test setup 
10   gAlice->Init("$(ALICE_ROOT)/TRD/AliTRDconfig.C");
11
12   // Run one event and create the hits
13   gAlice->SetDebug(2);
14   gAlice->Run(1);
15
16   if (gAlice) delete gAlice;
17   TFile *file = (TFile *) gROOT->GetListOfFiles()->FindObject("TRD_test.root");
18   gAlice = (AliRun *) file->Get("gAlice");
19
20   // Analyze the TRD hits
21   if (rc = AliTRDanalyzeHits()) return rc;
22
23   // Run the digitization
24   if (rc = AliTRDcreateDigits()) return rc;
25
26   if (gAlice) delete gAlice;
27   TFile *file = (TFile *) gROOT->GetListOfFiles()->FindObject("TRD_test.root");
28   gAlice = (AliRun *) file->Get("gAlice");
29
30   // Analyze the digits
31   if (rc = AliTRDanalyzeDigits()) return rc;
32
33   // Create the cluster
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
41   if (rc = AliTRDanalyzeCluster()) return rc;
42
43   return rc;
44
45 }
46
47 //_____________________________________________________________________________
48 Int_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 //_____________________________________________________________________________
238 Int_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 //_____________________________________________________________________________
297 Int_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 //_____________________________________________________________________________
487 Int_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 //_____________________________________________________________________________
544 Int_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 }