Update of the test suite
[u/mrichter/AliRoot.git] / TRD / AliTRDtest.C
1
2 Int_t AliTRDtest() 
3 {
4   //
5   // Test macro for the TRD code
6   //
7
8   Int_t rc = 0;
9
10   TFile *file;
11
12   // Initialize the test setup 
13   gAlice->Init("$(ALICE_ROOT)/TRD/AliTRDconfig.C");
14
15   // Run one event and create the hits
16   gAlice->SetDebug(2);
17   gAlice->Run(1);
18
19   if (gAlice) delete gAlice;
20   file   = (TFile *) gROOT->GetListOfFiles()->FindObject("TRD_test.root");
21   gAlice = (AliRun *) file->Get("gAlice");
22
23   // Analyze the TRD hits
24   if (rc = AliTRDanalyzeHits()) return rc;
25
26   // Run the digitization
27   if (rc = AliTRDcreateDigits()) return rc;
28
29   if (gAlice) delete gAlice;
30   file   = (TFile *) gROOT->GetListOfFiles()->FindObject("TRD_test.root");
31   gAlice = (AliRun *) file->Get("gAlice");
32
33   // Analyze the digits
34   if (rc = AliTRDanalyzeDigits()) return rc;
35
36   // Create the cluster
37   if (rc = AliTRDcreateCluster()) return rc;
38
39   if (gAlice) delete gAlice;
40   file   = (TFile *) gROOT->GetListOfFiles()->FindObject("TRD_test.root");
41   gAlice = (AliRun *) file->Get("gAlice");
42
43   // Analyze the cluster
44   if (rc = AliTRDanalyzeCluster()) return rc;
45
46   //file   = (TFile *) gROOT->GetListOfFiles()->FindObject("TRD_test.root");
47   //file->Close();
48
49   // Find the tracks
50   //if (rc = AliTRDcreateTracks()) return rc;
51
52   return rc;
53
54 }
55
56 //_____________________________________________________________________________
57 Int_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 //_____________________________________________________________________________
247 Int_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 //_____________________________________________________________________________
306 Int_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();
364   digitsManager->SetDebug(1);
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();
380   Int_t timeMax     = parameter->GetTimeTotal();
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 //_____________________________________________________________________________
497 Int_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 //_____________________________________________________________________________
554 Int_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();
645       Int_t    track0   = cluster->GetLabel(0);
646       Int_t    track1   = cluster->GetLabel(1);
647       Int_t    track2   = cluster->GetLabel(2);
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 }
734
735 //_____________________________________________________________________________
736 Int_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 }