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