TRD/AliTRDtest.C

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