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 }