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->Run(1);
  17
  18   //if (gAlice) delete gAlice;
  19   //file   = (TFile *) gROOT->GetListOfFiles()->FindObject("TRD_test.root");
  20   //gAlice = (AliRun *) file->Get("gAlice");
  21
  22   // Analyze the TRD hits
  23   if (rc = AliTRDanalyzeHits()) return rc;
  24
  25   // Run the digitization
  26   if (rc = AliTRDcreateDigits()) return rc;
  27
  28 //    if (gAlice) delete gAlice;
  29 //    file   = (TFile *) gROOT->GetListOfFiles()->FindObject("TRD_test.root");
  30 //    gAlice = (AliRun *) file->Get("gAlice");
  31
  32   // Analyze the digits
  33   if (rc = AliTRDanalyzeDigits()) return rc;
  34
  35 //    // Create the cluster
  36 //    if (rc = AliTRDcreateCluster()) return rc;
  37
  38 //    if (gAlice) delete gAlice;
  39 //    file   = (TFile *) gROOT->GetListOfFiles()->FindObject("TRD_test.root");
  40 //    gAlice = (AliRun *) file->Get("gAlice");
  41
  42 //    // Analyze the cluster
  43 //    if (rc = AliTRDanalyzeCluster()) return rc;
  44
  45   //file   = (TFile *) gROOT->GetListOfFiles()->FindObject("TRD_test.root");
  46   //file->Close();
  47
  48   // Find the tracks
  49   //if (rc = AliTRDcreateTracks()) return rc;
  50
  51   return rc;
  52
  53 }
  54
  55 //_____________________________________________________________________________
  56 Int_t AliTRDanalyzeHits()
  57 {
  58   //
  59   // Analyzes the hits and fills QA-histograms
  60   //
  61
  62   Int_t rc = 0;
  63
  64   AliRunLoader *rl = gAlice->GetRunLoader();
  65   if (!rl) {
  66     cout << "<AliTRDanalyzeHits> No RunLoader found" << endl;
  67     rc = 1;
  68     return rc;
  69   }
  70
  71   // Import the Trees for the event nEvent in the file
  72   rl->GetEvent(0);
  73
  74   AliLoader* loader = rl->GetLoader("TRDLoader");
  75   if (!loader) {
  76     cout << "<AliTRDanalyzeHits> No TRDLoader found" << endl;
  77     rc = 2;
  78     return rc;
  79   }
  80
  81   // Get the pointer to the TRD detector
  82   AliTRD *trd = (AliTRD *) gAlice->GetDetector("TRD");
  83   if (!trd) {
  84     cout << "<AliTRDanalyzeHits> No TRD detector found" << endl;
  85     rc = 2;
  86     return rc;
  87   }
  88
  89   // Get the pointer to the geometry object
  90   AliTRDgeometryFull *geo;
  91   if (trd) {
  92     geo = (AliTRDgeometryFull *) trd->GetGeometry();
  93   }
  94   else {
  95     cout << "<AliTRDanalyzeHits> No TRD geometry found" << endl;
  96     rc = 3;
  97     return rc;
  98   }
  99
 100   AliTRDCommonParam *par = AliTRDCommonParam::Instance();
 101
 102   // Define the histograms
 103   TH1F *hQdedx  = new TH1F("hQdedx","Charge dedx-hits",100,0.0,1000.0);
 104   TH1F *hQtr    = new TH1F("hQtr"  ,"Charge TR-hits"  ,100,0.0,1000.0);
 105
 106   Float_t rmin   = geo->Rmin();
 107   Float_t rmax   = geo->Rmax();
 108   Float_t length = geo->GetChamberLength(0,2);
 109   Float_t width  = geo->GetChamberWidth(0);
 110   Int_t   ncol   = par->GetColMax(0);
 111   Int_t   nrow   = par->GetRowMax(0,2,13);
 112   Int_t   ntime  = ((Int_t) (rmax - rmin) / 22.0);
 113
 114   TH2F *hZY     = new TH2F("hZY"   ,"Y vs Z (chamber 0)", nrow,-length/2.,length/2.
 115                                                         ,ntime,      rmin,     rmax);
 116   TH2F *hXZ     = new TH2F("hXZ"   ,"Z vs X (plane 0)"  , ncol, -width/2., width/2.
 117                                                         , nrow,-length/2.,length/2.);
 118
 119   // Get the pointer hit tree
 120   TTree *hitTree = loader->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     // Loop through the TRD hits
 143     Int_t iHit = 0;
 144     AliTRDhit *hit = (AliTRDhit *) trd->FirstHit(-1);
 145     while (hit) {
 146
 147       countHits++;
 148       iHit++;
 149
 150       Float_t x     = hit->X();
 151       Float_t y     = hit->Y();
 152       Float_t z     = hit->Z();
 153       Float_t q     = hit->GetCharge();
 154       Int_t   track = hit->Track();
 155       Int_t   det   = hit->GetDetector();
 156       Int_t   plane = geo->GetPlane(det);
 157
 158       if      (q > 0) {
 159         hQdedx->Fill(q);
 160         hZY->Fill(z,y);
 161         if (plane == 0) {
 162           hXZ->Fill(x,z);
 163         }
 164       }
 165       else if (q < 0) {
 166         hQtr->Fill(TMath::Abs(q));
 167       }
 168
 169       hit = (AliTRDhit *) trd->NextHit();
 170
 171     }
 172
 173   }
 174
 175   cout << "<AliTRDanalyzeHits> Found " << countHits << " hits in total" << endl;
 176
 177   TCanvas *cHits = new TCanvas("cHits","AliTRDanalyzeHits",50,50,600,600);
 178   cHits->Divide(2,2);
 179   cHits->cd(1);
 180   hXZ->Draw("COL");
 181   cHits->cd(2);
 182   hZY->Draw("COL");
 183   cHits->cd(3);
 184   gPad->SetLogy();
 185   hQdedx->Draw();
 186   cHits->cd(4);
 187   gPad->SetLogy();
 188   hQtr->Draw();
 189
 190   return rc;
 191
 192 }
 193
 194 //_____________________________________________________________________________
 195 Int_t AliTRDcreateDigits()
 196 {
 197   //
 198   // Creates the digits from the hits of the slow simulator
 199   //
 200
 201   Int_t rc = 0;
 202
 203   if (!gAlice) {
 204     cout << "<AliTRDcreateDigits> No AliRun object found" << endl;
 205     rc = 1;
 206     return rc;
 207   }
 208   gAlice->GetEvent(0);
 209
 210   // Create the TRD digitzer
 211   AliTRDdigitizer *digitizer = new AliTRDdigitizer("TRDdigitizer"
 212                                                   ,"TRD digitizer class");
 213   digitizer->Open("TRD_test.root");
 214   digitizer->InitDetector();
 215   digitizer->InitOutput(0);
 216
 217   // Set the parameter
 218   digitizer->SetDebug(1);
 219
 220   // Define the parameter object
 221   // If no external parameter object is defined,
 222   // default parameter will be used
 223   AliTRDparameter *parameter = new AliTRDparameter("TRDparameter"
 224                                                   ,"TRD parameter class");
 225   digitizer->SetParameter(parameter);
 226
 227   // Create the digits
 228   if (!(digitizer->MakeDigits())) {
 229     rc = 2;
 230     return rc;
 231   }
 232
 233   // Write the digits into the input file
 234   //if (!(digitizer->MakeBranch())) {
 235   //  rc = 3;
 236   //  return rc;
 237   //}
 238
 239   // Write the digits into the input file
 240   if (!(digitizer->WriteDigits())) {
 241     rc = 4;
 242     return rc;
 243   }
 244
 245   // Save the parameter object in the AliROOT file
 246   if (!(parameter->Write())) {
 247     rc = 4;
 248     return rc;
 249   }
 250
 251   return rc;
 252
 253 }
 254
 255 //_____________________________________________________________________________
 256 Int_t AliTRDanalyzeDigits()
 257 {
 258   //
 259   // Analyzes the digits
 260   //
 261
 262   Int_t rc = 0;
 263
 264   const Int_t kNpla = AliTRDgeometry::Nplan();
 265
 266   if (!gAlice) {
 267     cout << "<AliTRDanalyzeDigits> No AliRun object found" << endl;
 268     rc = 1;
 269     return rc;
 270   }
 271
 272   AliRunLoader *rl = gAlice->GetRunLoader();
 273   if (!rl) {
 274     cout << "<AliTRDanalyzeHits> No RunLoader found" << endl;
 275     rc = 2;
 276     return rc;
 277   }
 278
 279   // Import the Trees for the event nEvent in the file
 280   rl->LoadKinematics();
 281   rl->GetEvent(0);
 282   rl->GetHeader();
 283
 284   AliLoader* loader = rl->GetLoader("TRDLoader");
 285   if (!loader) {
 286     cout << "<AliTRDanalyzeHits> No TRDLoader found" << endl;
 287     rc = 3;
 288     return rc;
 289   }
 290
 291   // Get the pointer to the TRD detector
 292   AliTRD *trd = (AliTRD *) gAlice->GetDetector("TRD");
 293   if (!trd) {
 294     cout << "<AliTRDanalyzeDigits> No TRD detector found" << endl;
 295     rc = 4;
 296     return rc;
 297   }
 298
 299   // Get the parameter object
 300   AliTRDSimParam    *parSim = AliTRDSimParam::Instance();
 301   AliTRDCommonParam *parCom = AliTRDCommonParam::Instance();
 302   AliTRDcalibDB     *cal    = AliTRDcalibDB::Instance();
 303
 304   // Define the histograms
 305   Int_t adcRange = ((Int_t) parSim->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     = cal->GetNumberOfTimeBins();
 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 = parCom->GetRowMax(iPla,iCha,iSec);
 362     Int_t colMax = parCom->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 }