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