Update of macros
authorcblume <cblume@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 25 Mar 2002 20:03:21 +0000 (20:03 +0000)
committercblume <cblume@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 25 Mar 2002 20:03:21 +0000 (20:03 +0000)
TRD/AliTRDanaCluster.C [new file with mode: 0644]
TRD/AliTRDanaDigits.C [new file with mode: 0644]
TRD/AliTRDdigits2cluster.C [new file with mode: 0644]
TRD/AliTRDdisplayDigits3D.C
TRD/AliTRDdisplayTracks.C [new file with mode: 0644]
TRD/AliTRDhits2digits.C [new file with mode: 0644]
TRD/AliTRDhits2sdigits.C [new file with mode: 0644]
TRD/AliTRDmerge.C [new file with mode: 0644]
TRD/AliTRDrunSimple.C [new file with mode: 0644]
TRD/AliTRDsdigits2digits.C [new file with mode: 0644]
TRD/AliTRDtest.C

diff --git a/TRD/AliTRDanaCluster.C b/TRD/AliTRDanaCluster.C
new file mode 100644 (file)
index 0000000..a1b3315
--- /dev/null
@@ -0,0 +1,95 @@
+void AliTRDanaCluster() 
+{
+
+/////////////////////////////////////////////////////////////////////////
+//
+// Example macro for the analysis of the TRD cluster
+//
+/////////////////////////////////////////////////////////////////////////
+
+  // Dynamically link some shared libs
+  if (gClassTable->GetID("AliRun") < 0) {
+    gROOT->LoadMacro("loadlibs.C");
+    loadlibs();
+  }
+
+  // Input file name
+  Char_t *alifile = "galice.root"; 
+
+  // Event number
+  Int_t   nEvent  = 0;
+
+  // Define the histograms
+  TH1F *hCharge = new TH1F("hCharge","Cluster charge",100,0.0,1000.0);
+
+  // Connect the AliRoot file containing Geometry, Kine, Hits, and Digits
+  TFile *gafl = (TFile*) gROOT->GetListOfFiles()->FindObject(alifile);
+  if (!gafl) {
+    cout << "Open the ALIROOT-file " << alifile << endl;
+    gafl = new TFile(alifile);
+  }
+  else {
+    cout << alifile << " is already open" << endl;
+  }
+
+  // Get AliRun object from file or create it if not on file
+  gAlice = (AliRun*) gafl->Get("gAlice");
+  if (gAlice) {
+    cout << "AliRun object found on file" << endl;
+  }
+  else {
+    gAlice = new AliRun("gAlice","Alice test program");
+  }
+
+  // Import the Trees for the event nEvent in the file
+  Int_t nparticles = gAlice->GetEvent(nEvent);
+  if (nparticles <= 0) break;
+  
+  // Get the pointer to the hit-tree
+  Char_t treeName[14];
+  sprintf(treeName,"TreeR%d_TRD",nEvent);
+  TTree          *clusterTree  = gafl->Get(treeName);
+  clusterTree->Print();
+  // Get the pointer to the detector classes
+  AliTRDv1       *trd          = (AliTRDv1*) gAlice->GetDetector("TRD");
+  // Get the geometry
+  AliTRDgeometry *geo          = trd->GetGeometry();
+  // Get the pointer to the hit container
+  TObjArray      *clusterArray = trd->RecPoints();
+  // Set the branch address
+  clusterTree->GetBranch("TRDcluster")->SetAddress(&clusterArray);
+
+  Int_t nEntries = clusterTree->GetEntries();
+  cout << "nEntries = " << nEntries << endl;
+
+  // Loop through all entries in the tree
+  Int_t nbytes;
+  for (Int_t iEntry = 0; iEntry < nEntries; iEntry++) {
+
+    // Import the tree
+    nbytes += clusterTree->GetEvent(iEntry);
+
+    // Get the number of points in the detector 
+    Int_t ncluster = clusterArray->GetEntriesFast();
+
+    // Loop through all TRD digits
+    for (Int_t icluster = 0; icluster < ncluster; icluster++) {
+
+      // Get the information for this digit
+      AliTRDcluster *cluster = (AliTRDcluster *) clusterArray->UncheckedAt(icluster);
+      Int_t    detector = cluster->GetDetector();      
+      Int_t    sector   = geo->GetSector(detector);
+      Int_t    plane    = geo->GetPlane(detector);
+      Int_t    chamber  = geo->GetChamber(detector);
+      Float_t  charge   = cluster->GetQ();
+
+      hCharge->Fill(charge);
+
+    }
+
+  }
+
+  TCanvas *c1 = new TCanvas("c1","Cluster",50,50,600,400);
+  hCharge->Draw();
+
+}
diff --git a/TRD/AliTRDanaDigits.C b/TRD/AliTRDanaDigits.C
new file mode 100644 (file)
index 0000000..7e3e9d7
--- /dev/null
@@ -0,0 +1,115 @@
+void AliTRDanaDigits() 
+{
+
+/////////////////////////////////////////////////////////////////////////
+//
+// Example macro for the analysis of the TRD digits and the use
+// of the AliTRDmatrix class.
+//
+/////////////////////////////////////////////////////////////////////////
+
+  // Dynamically link some shared libs
+  if (gClassTable->GetID("AliRun") < 0) {
+    gROOT->LoadMacro("loadlibs.C");
+    loadlibs();
+  }
+
+  // Input file name
+  //Char_t *alifile = "galice.root"; 
+  Char_t *alifile = "galice_jiri.root"; 
+
+  // Event number
+  Int_t   nEvent  = 0;
+
+  // Define the objects
+  AliTRDv1       *trd;
+  AliTRDgeometry *geo;
+  AliTRDdigit    *digit;
+
+  Int_t           track;
+
+  // Connect the AliRoot file containing Geometry, Kine, Hits, and digits
+  TFile *gafl = (TFile*) gROOT->GetListOfFiles()->FindObject(alifile);
+  if (!gafl) {
+    cout << "Open the ALIROOT-file " << alifile << endl;
+    gafl = new TFile(alifile);
+  }
+  else {
+    cout << alifile << " is already open" << endl;
+  }
+
+  // Get AliRun object from file or create it if not on file
+  gAlice = (AliRun*) gafl->Get("gAlice");
+  if (gAlice)  
+    cout << "AliRun object found on file" << endl;
+  else
+    gAlice = new AliRun("gAlice","Alice test program");
+
+  // Import the Trees for the event nEvent in the file
+  Int_t nparticles = gAlice->GetEvent(nEvent);
+  if (nparticles <= 0) break;
+  
+  // Get the pointer to the detector object
+  trd = (AliTRDv1*) gAlice->GetDetector("TRD");
+
+  // Get the pointer to the geometry object
+  if (trd) {
+    geo = trd->GetGeometry();
+  }
+  else {
+    cout << "Cannot find the geometry" << endl;
+    break;
+  }
+
+  // Create the digits manager
+  AliTRDdigitsManager *digitsManager = new AliTRDdigitsManager();
+  digitsManager->SetDebug(1);
+  digitsManager->SetSDigits(kTRUE);
+
+  // Read the digits from the file
+  digitsManager->Open(alifile);
+  digitsManager->ReadDigits();
+
+  // Get the detector number
+  Int_t iDet = 423;
+  cout << " iDet = " << iDet << endl;
+
+  // Define the detector matrix for one chamber
+  const Int_t iSec = geo->GetSector(iDet);
+  const Int_t iCha = geo->GetChamber(iDet);
+  const Int_t iPla = geo->GetPlane(iDet);
+  Int_t  rowMax = geo->GetRowMax(iPla,iCha,iSec);
+  Int_t  colMax = geo->GetColMax(iPla);
+  Int_t timeMax = geo->GetTimeMax();
+  cout << "Geometry: rowMax = "  <<  rowMax
+                << " colMax = "  <<  colMax
+                << " timeMax = " << timeMax << endl;
+  AliTRDmatrix *matrix = new AliTRDmatrix(rowMax,colMax,timeMax,iSec,iCha,iPla);
+
+  // Loop through the detector pixel
+  for (Int_t time = 0; time < timeMax; time++) {
+    for (Int_t  col = 0;  col <  colMax;  col++) {
+      for (Int_t  row = 0;  row <  rowMax;  row++) {
+
+        digit = digitsManager->GetDigit(row,col,time,iDet);
+        track = digitsManager->GetTrack(0,row,col,time,iDet);
+        
+        printf("time=%d, col=%d, row=%d, amp=%f\n",time,col,row,digit->GetAmp());
+        matrix->SetSignal(row,col,time,digit->GetAmp());
+
+        delete digit;
+
+      }
+    }
+  }
+
+  // Display the detector matrix
+  matrix->Draw();
+  //matrix->DrawRow(18);
+  //matrix->DrawCol(58);
+  //matrix->DrawTime(20);
+  matrix->ProjRow();
+  matrix->ProjCol();
+  matrix->ProjTime();
+
+}
diff --git a/TRD/AliTRDdigits2cluster.C b/TRD/AliTRDdigits2cluster.C
new file mode 100644 (file)
index 0000000..cac782e
--- /dev/null
@@ -0,0 +1,47 @@
+void AliTRDdigits2cluster() 
+{
+
+///////////////////////////////////////////////////////////////////////// 
+//
+// Creates cluster from the digit information. 
+//
+///////////////////////////////////////////////////////////////////////// 
+
+  // Dynamically link some shared libs
+  if (gClassTable->GetID("AliRun") < 0) {
+    gROOT->LoadMacro("loadlibs.C");
+    loadlibs();
+    cout << "Loaded shared libraries" << endl;
+  }
+
+  // Input and output file names
+  Char_t *infile  = "galice.root";
+  Char_t *outfile = "TRDclusters.root";
+
+  // Create the clusterizer
+  AliTRDclusterizerV1 *clusterizer = 
+    new AliTRDclusterizerV1("clusterizer","Clusterizer class"); 
+
+  // Set the parameter
+  clusterizer->SetClusMaxThresh(0);
+  clusterizer->SetClusSigThresh(0);
+  //clusterizer->SetVerbose(1);
+  clusterizer->Dump();
+
+  // Open the AliRoot file 
+  clusterizer->Open(infile,0);
+  //clusterizer->Open(infile,outfile,0);
+
+  // Load the digits
+  clusterizer->ReadDigits();
+  // Find the cluster
+  clusterizer->MakeClusters();
+
+  // Write the cluster tree into file AliTRDclusters.root
+  clusterizer->WriteClusters(-1);
+
+  // Save the clusterizer class in the AliROOT file
+  // clusterizer->Write();
+
+}
index c60c103..2a28d3a 100644 (file)
@@ -1,5 +1,5 @@
 //_____________________________________________________________________________
-Int_t AliTRDdisplayDigits3D(Int_t event = 0, Int_t thresh = 0
+Int_t AliTRDdisplayDigits3D(Int_t event = 0, Int_t thresh = 2
                           , Bool_t sdigits = kFALSE) 
 {
   //  
diff --git a/TRD/AliTRDdisplayTracks.C b/TRD/AliTRDdisplayTracks.C
new file mode 100644 (file)
index 0000000..d38f1de
--- /dev/null
@@ -0,0 +1,143 @@
+void AliTRDdisplayTracks(Int_t track) 
+{
+
+  c1 = new TCanvas( "RecPoints", "RecPoints Display", 10, 10, 710, 740);
+  c1->SetFillColor(1);
+  TView *v=new TView(1);
+
+      v->SetRange(-350.,-350.,-400.,350.,350.,400.); // full
+  //  v->SetRange(0.,0.,0.,350.,350.,400.);  // top right
+  //  v->SetRange(-350.,0.,0.,0.,350.,400.); // top left
+  //  v->SetRange(0.,-350.,0.,350.,0.,400.); // bottom right 
+  //  v->SetRange(-350.,-350.,0.,0.,0.,400.); // bottom left
+  //  v->Side();
+  v->Top();
+  cerr<<"psi = "<<v->GetPsi()<<endl;
+
+  // Dynamically link some shared libs
+  if (gClassTable->GetID("AliRun") < 0) {
+    gROOT->LoadMacro("loadlibs.C");
+    loadlibs();
+    cout << "Loaded shared libraries" << endl;
+  }       
+
+
+// Load tracks
+
+   TFile *tf=TFile::Open("AliTRDtracks.root");
+   if (!tf->IsOpen()) {cerr<<"Can't open AliTRDtracks.root !\n"; return 3;}
+   TObjArray tarray(2000);
+   TTree *tracktree=(TTree*)tf->Get("TreeT");
+   TBranch *tbranch=tracktree->GetBranch("tracks");
+   Int_t nentr=tracktree->GetEntries();
+   cerr<<"Found "<<nentr<<" entries in the track tree"<<endl;
+   for (Int_t i=0; i<nentr; i++) {
+       AliTRDtrack *iotrack=new AliTRDtrack;
+       tbranch->SetAddress(&iotrack);
+       tracktree->GetEvent(i);
+       tarray.AddLast(iotrack);
+       cerr<<"got track "<<i<<": index ="<<iotrack->GetLabel()<<endl;
+   }
+   tf->Close();                 
+
+
+  // Load clusters
+  Char_t *alifile = "AliTRDclusters.root";
+  Int_t   nEvent  = 0;
+  TObjArray rparray(2000);
+  TObjArray *RecPointsArray = &rparray;
+  AliTRDtracker *Tracker = new AliTRDtracker("dummy","dummy");
+  Tracker->ReadClusters(RecPointsArray,alifile,-1); 
+  Int_t nRecPoints = RecPointsArray->GetEntriesFast();
+  cerr<<"Found "<<nRecPoints<<" rec. points"<<endl;
+
+
+  // Connect the AliRoot file containing Geometry, Kine, Hits, and Digits
+  alifile = "galice.root";
+  TFile *gafl = (TFile*) gROOT->GetListOfFiles()->FindObject(alifile);
+  if (!gafl) {
+    cout << "Open the ALIROOT-file " << alifile << endl;
+    gafl = new TFile(alifile);
+  }
+  else {
+    cout << alifile << " is already open" << endl;
+  }
+
+  // Get AliRun object from file or create it if not on file
+  gAlice = (AliRun*) gafl->Get("gAlice");
+  if (gAlice)
+    cout << "AliRun object found on file" << endl;
+  else
+    gAlice = new AliRun("gAlice","Alice test program");
+
+
+  AliTRDv1       *TRD           = (AliTRDv1*) gAlice->GetDetector("TRD");
+  AliTRDgeometry *fGeom   = TRD->GetGeometry();   
+
+  Int_t i,j,index,det,sector, ti[3];
+  Double_t x,y,z, cs,sn,tmp;
+  Float_t global[3], local[3];
+
+  // Display all TRD RecPoints
+  TPolyMarker3D *pm = new TPolyMarker3D(nRecPoints);
+
+  for (Int_t i = 0; i < nRecPoints; i++) {
+    printf("\r point %d out of %d",i,nRecPoints);
+    AliTRDrecPoint *rp = (AliTRDrecPoint *) RecPointsArray->UncheckedAt(i);    
+
+    local[0]=rp->GetLocalRow(); 
+    local[1]=rp->GetLocalCol(); 
+    local[2]=rp->GetLocalTime();
+    det=rp->GetDetector();
+
+    for (j = 0; j < 3; j++) { ti[j] = rp->GetTrackIndex(j); }
+
+    if((track < 0) ||
+       ((ti[0]==track)||(ti[1]==track)||(ti[2]==track))) {
+
+      if(fGeom->Local2Global(det,local,global)) {     
+       x=global[0]; y=global[1]; z=global[2];
+       pm->SetPoint(i,x,y,z);
+      }
+    }
+  }
+
+  pm->SetMarkerSize(1); pm->SetMarkerColor(10); pm->SetMarkerStyle(1);
+  pm->Draw();
+  
+
+  Int_t ntracks = tarray.GetEntriesFast();
+
+  for (i = 0; i < ntracks; i++) {
+    AliTRDtrack *pt = (AliTRDtrack *) tarray.UncheckedAt(i), &t=*pt;
+    Int_t nclusters = t.GetNclusters();
+    cerr<<"in track "<<i<<" found "<<nclusters<<" clusters"<<endl;
+   
+    TPolyMarker3D *pm = new TPolyMarker3D(nclusters);
+    for(j = 0; j < nclusters; j++) {
+      index = t.GetClusterIndex(j);
+      AliTRDrecPoint *rp = (AliTRDrecPoint *) RecPointsArray->UncheckedAt(index);    
+      local[0]=rp->GetLocalRow(); 
+      local[1]=rp->GetLocalCol(); 
+      local[2]=rp->GetLocalTime();
+      det=rp->GetDetector();
+
+      if(fGeom->Local2Global(det,local,global)) {     
+       x=global[0]; y=global[1]; z=global[2];
+       pm->SetPoint(j,x,y,z);
+      }
+    } 
+    pm->SetMarkerSize(1); pm->SetMarkerColor(i%6+3); pm->SetMarkerStyle(1);
+    pm->Draw();
+  }
+  
+  TGeometry *geom=(TGeometry*)gafl->Get("AliceGeom");
+  geom->Draw("same");
+
+  c1->Modified(); 
+  
+  c1->Update();
+  
+  gafl->Close();
+
+}
diff --git a/TRD/AliTRDhits2digits.C b/TRD/AliTRDhits2digits.C
new file mode 100644 (file)
index 0000000..28a3970
--- /dev/null
@@ -0,0 +1,47 @@
+void AliTRDhits2digits() 
+{
+
+/////////////////////////////////////////////////////////////////////////
+//
+// Creates digits from the hit information. 
+//
+/////////////////////////////////////////////////////////////////////////
+
+  // Dynamically link some shared libs
+  if (gClassTable->GetID("AliRun") < 0) {
+    gROOT->LoadMacro("loadlibs.C");
+    loadlibs();
+    cout << "Loaded shared libraries" << endl;
+  }
+
+  // Input (and output) file name
+  Char_t *alifile = "galice.root"; 
+
+  // Create the TRD digitzer 
+  AliTRDdigitizer *digitizer = new AliTRDdigitizer("TRDdigitizer"
+                                                  ,"TRD digitizer class");
+
+  // Set the parameter
+  digitizer->SetDebug(1);
+
+  // Open the AliRoot file
+  digitizer->Open(alifile);
+
+
+  // Define the parameter object
+  // If no external parameter object is defined, 
+  // default parameter will be used
+  AliTRDparameter *parameter = new AliTRDparameter("TRDparameter"
+                                                 ,"TRD parameter class");
+  digitizer->SetParameter(parameter);
+
+  // Create the digits
+  digitizer->MakeDigits();
+
+  // Write the digits into the input file
+  digitizer->WriteDigits();
+
+  // Save the parameter object in the AliROOT file
+  parameter->Write();
+
+}
diff --git a/TRD/AliTRDhits2sdigits.C b/TRD/AliTRDhits2sdigits.C
new file mode 100644 (file)
index 0000000..d1f9f58
--- /dev/null
@@ -0,0 +1,49 @@
+void AliTRDhits2sdigits() 
+{
+
+  /////////////////////////////////////////////////////////////////////////
+  //
+  // Creates summable digits from the hit information. 
+  //
+  /////////////////////////////////////////////////////////////////////////
+
+  // Dynamically link some shared libs
+  if (gClassTable->GetID("AliRun") < 0) {
+    gROOT->LoadMacro("loadlibs.C");
+    loadlibs();
+    cout << "Loaded shared libraries" << endl;
+  }
+
+  // Input (and output) file name
+  Char_t *alifile = "galice.root"; 
+
+  // Create the TRD digitzer 
+  AliTRDdigitizer *digitizer = new AliTRDdigitizer("TRDdigitizer"
+                                                  ,"TRD digitizer class");
+
+  // Set the parameter
+  digitizer->SetDebug(1);
+
+  // For the summable digits
+  digitizer->SetSDigits(kTRUE);
+
+  // Open the AliRoot file
+  digitizer->Open(alifile);
+
+  // Define the parameter object
+  // If no external parameter object is defined, 
+  // default parameter will be used
+  AliTRDparameter *parameter = new AliTRDparameter("TRDparameter"
+                                                 ,"TRD parameter class");
+  digitizer->SetParameter(parameter);
+
+  // Create the digits
+  digitizer->MakeDigits();
+
+  // Write the digits into the input file
+  digitizer->WriteDigits();
+
+  // Save the parameter object in the AliROOT file
+  parameter->Write();
+
+}
diff --git a/TRD/AliTRDmerge.C b/TRD/AliTRDmerge.C
new file mode 100644 (file)
index 0000000..2daa903
--- /dev/null
@@ -0,0 +1,27 @@
+void AliTRDmerge()
+{
+  /////////////////////////////////////////////////////////////////////////
+  //
+  // Merges different file with summable digits
+  //
+  /////////////////////////////////////////////////////////////////////////
+
+  AliRunDigitizer *manager = new AliRunDigitizer(2,1);
+  manager->SetInputStream(0,"galice_signal.root");
+  manager->SetInputStream(1,"galice_background.root");
+
+  AliTRDdigitizer *digitizer = new AliTRDdigitizer(manager
+                                                  ,"TRDdigitizer"
+                                                  ,"TRD digitizer class");
+
+  // Define the parameter object
+  // If no external parameter object is defined, 
+  // default parameter will be used
+  AliTRDparameter *parameter = new AliTRDparameter("TRDparameter"
+                                                 ,"TRD parameter class");
+  digitizer->SetParameter(parameter);
+
+  // Do the merging
+  manager->Exec("deb");
+
+}
diff --git a/TRD/AliTRDrunSimple.C b/TRD/AliTRDrunSimple.C
new file mode 100644 (file)
index 0000000..6658c03
--- /dev/null
@@ -0,0 +1,274 @@
+
+#if !defined(__CINT__) || defined(__MAKECINT__)
+
+#include "TH1.h"
+#include "TFile.h"
+#include "TCanvas.h"
+#include "TProfile.h"
+#include "TClonesArray.h"
+
+#include "AliModule.h"
+#include "AliFRAMEv1.h"
+
+#include "AliTRDv1.h"
+#include "AliTRDhit.h"
+#include "AliTRDsim.h"
+#include "AliTRDsimple.h"
+#include "AliTRDsimpleGen.h"
+#include "AliTRDsimpleMC.h"
+#include "AliTRDgeometry.h"
+#include "AliTRDdigitizer.h"
+#include "AliTRDdigitsManager.h"
+#include "AliTRDdataArrayI.h"
+
+#endif
+
+void AliTRDrunSimple()
+{
+
+  ///////////////////////////////////////////////////////////////////////////////
+  //
+  //  Macro to run the simplified simulator
+  //
+  ///////////////////////////////////////////////////////////////////////////////
+
+  //_____________________________________________________________________________ 
+  //
+  // Initialization part
+  //_____________________________________________________________________________ 
+  //
+
+  // The simple simulator
+  AliTRDsimple *simple = new AliTRDsimple();
+  simple->Init();
+
+  // Initialize a dummy frame so that the TRD works
+  new AliFRAMEv1("FRAME","Space Frame");             
+
+  // Initialize the TRD detector
+  AliTRDv1     *trd    = new AliTRDv1("TRD","TRD slow simulator");
+
+  // Needed for some material properties 
+  trd->CreateMaterials();
+
+  // Select the gas mixture (0: 97% Xe + 3% isobutane, 1: 85% Xe + 15% CO2)
+  trd->SetGasMix(1);
+  // Get the pointer to the geometry object
+  AliTRDgeometry *geometry = trd->GetGeometry();
+  // The number of timebins
+  geometry->SetNTimeBin(30);
+  // The additional timebins before and after the drift region
+  geometry->SetExpandTimeBin(5,15);
+    
+  // Switch on TR production
+  AliTRDsim *trdsim = trd->CreateTR();
+  trdsim->SetNFoils(100);
+  trdsim->SetFoilThick(0.0013);
+  trdsim->SetGapThick(0.0500);
+
+  // Initialize the TRD object
+  trd->Init();
+
+  // The digitizer
+  AliTRDdigitizer *digitizer = new AliTRDdigitizer("TRDdigitizer","Digitizer class");
+  digitizer->SetSimple();
+  digitizer->SetGasGain(3200.);      // GEANT + new Ermilova spectrum
+  digitizer->SetChipGain(6.1);
+  digitizer->SetNoise(1000.);
+  digitizer->SetADCinRange(1000.);
+  digitizer->SetADCoutRange(1023.);
+  digitizer->SetADCthreshold(0);
+  digitizer->InitDetector();
+  //digitizer->SetTimeResponse(0);
+  //digitizer->SetVerbose(1);
+                                    
+  // The event generator
+  AliTRDsimpleGen *generator = simple->GetGenerator();         
+  generator->SetMomentum(3.0,3.0);
+  generator->SetPdg(11);                             // Electrons 
+  //generator->SetPdg(211);                            // Pions 
+
+  //_____________________________________________________________________________ 
+  //
+  // Histograms
+  //_____________________________________________________________________________ 
+  //
+
+  Int_t timeMax  = geometry->GetTimeTotal();   
+  Int_t adcRange = ((Int_t) digitizer->GetADCoutRange()); 
+
+  TH1F     *hQ       = new TH1F("hQ"    ,"Charge per hit (all)" ,100,0.0,1000.0);
+  TH1F     *hQdedx   = new TH1F("hQdedx","Charge per hit (dedx)",100,0.0,1000.0);
+  TH1F     *hQtr     = new TH1F("hQtr  ","Charge per hit (tr)"  ,100,0.0,1000.0);
+
+  TProfile *hQX      = new TProfile("hQX"    ,"Charge per hit vs X (all)" ,35,0.0,3.5);
+  TProfile *hQXdedx  = new TProfile("hQXdedx","Charge per hit vs X (dedx)",35,0.0,3.5);
+  TProfile *hQXtr    = new TProfile("hQXtr  ","Charge per hit vs X (tr)"  ,35,0.0,3.5);
+
+  TH1F     *hNstep   = new TH1F("hNstep","Number of steps / cm"    , 151,-0.5, 150.5);
+  TH1F     *hNel     = new TH1F("hNel"  ,"Number of electrons / cm",1001,-0.5,1000.5);
+
+  TH1F     *hAmp     = new TH1F("hAmp","Amplitude of the digits"
+                                      ,adcRange+1,-0.5,((Float_t) adcRange)+0.5);         
+  TProfile *hAmpTime = new TProfile("hAmpTime","Amplitude vs timebin"
+                                             ,timeMax,-0.5,((Float_t) timeMax)-0.5);
+
+  //_____________________________________________________________________________ 
+  //
+  // Event loop
+  //_____________________________________________________________________________ 
+  //
+
+  // Number of events (tracks)
+  Int_t nEvent = 10000;
+
+  Float_t x0 = geometry->GetTime0(0) - AliTRDgeometry::DrThick(); 
+
+  TClonesArray *hitsArray = trd->Hits();
+
+  for (Int_t iEvent = 0; iEvent < nEvent; iEvent++) {
+
+    if (!(iEvent % 100) && (iEvent)) {
+      printf("Event no. %d\n",iEvent);
+    }
+
+    // Generate the hits for one track
+    simple->ProcessEvent(iEvent);
+
+    // Analyze the hits
+    Float_t nElectrons = 0.0;
+    for (Int_t iHit = 0; iHit < hitsArray->GetEntries(); iHit++) {
+
+      AliTRDhit *hit = (AliTRDhit *) hitsArray->UncheckedAt(iHit);
+      Int_t   charge = TMath::Abs(hit->GetCharge());
+      Float_t x      = hit->X() - x0;
+      hQ->Fill(charge);
+      hQX->Fill(x,charge);
+      if (hit->FromDrift() ||
+          hit->FromAmplification()) {
+        hQdedx->Fill(charge);
+        hQXdedx->Fill(x,charge);
+        nElectrons += charge;
+      }
+      if (hit->FromTRphoton()) {
+        hQtr->Fill(charge);
+        hQXtr->Fill(x,charge);
+      }
+
+    }
+
+    hNstep->Fill(((AliTRDsimpleMC *) gMC)->GetNStep() / 3.5);
+    hNel->Fill(nElectrons / 3.5);
+
+    // Digitize the hits
+    digitizer->MakeDigits();
+
+    // Analyze the digits
+    Int_t row =  8;
+    Int_t col = 48;
+    AliTRDdigitsManager *digitsManager = digitizer->Digits();
+    AliTRDdataArrayI    *digitsArray   = digitsManager->GetDigits(12);
+    for (Int_t time = 0; time < timeMax; time++) {  
+
+      Int_t amp = digitsArray->GetDataUnchecked(row,col,time);
+      if (amp > 0) {
+        hAmp->Fill((Double_t) amp);
+        hAmpTime->Fill((Double_t) time, (Double_t) amp);
+      }
+
+    }
+    
+  }
+
+  //_____________________________________________________________________________ 
+  //
+  // Plot the spectra
+  //_____________________________________________________________________________ 
+  //
+
+  TCanvas *cHit = new TCanvas("cHit","Hits",50,50,800,600);
+  cHit->Divide(3,2);
+
+  cHit->cd(1);
+  gPad->SetLogy();
+  hQ->SetLineColor(4);
+  hQ->SetXTitle("Q (a.u.)");
+  hQ->SetYTitle("Entries");
+  hQ->Draw();
+  hQdedx->SetLineColor(3);
+  hQdedx->Draw("SAME");
+  hQtr->SetLineColor(2);
+  hQtr->Draw("SAME"); 
+
+  cHit->cd(2); 
+  gPad->SetLogy();
+  hQX->SetLineColor(4);
+  hQX->SetXTitle("x (cm)");
+  hQX->SetYTitle("<Q> (a.u.)");
+  hQX->Draw("HIST");
+  hQXdedx->SetLineColor(3);
+  hQXdedx->Draw("SAMEHIST");
+
+  cHit->cd(3); 
+  gPad->SetLogy();
+  hQXtr->SetLineColor(2);
+  hQXtr->SetXTitle("x (cm)");
+  hQXtr->SetYTitle("<Q> (a.u.)");
+  hQXtr->Draw("HIST"); 
+
+  cHit->cd(4);
+  hNstep->SetLineColor(4);
+  hNstep->SetXTitle("N_{step}  / cm");
+  hNstep->SetYTitle("Entries");
+  hNstep->Draw();
+
+  cHit->cd(5);
+  hNel->SetLineColor(4);
+  hNel->SetXTitle("N_{electron} / cm");
+  hNel->SetYTitle("Entries");
+  hNel->Draw();
+
+  TCanvas *cDigit = new TCanvas("cDigit","Digits",50,50,600,400);
+  cDigit->Divide(2,1);
+
+  cDigit->cd(1);
+  gPad->SetLogy();
+  hAmp->SetLineColor(2);
+  hAmp->SetXTitle("ADC channels");
+  hAmp->SetYTitle("Entries");
+  hAmp->Draw();
+
+  cDigit->cd(2);
+  hAmpTime->SetLineColor(2);
+  hAmpTime->SetXTitle("Timebin");
+  hAmpTime->SetYTitle("<ADC channels>");
+  hAmpTime->Draw("HIST");
+
+  //_____________________________________________________________________________ 
+  //
+  // Save the histograms
+  //_____________________________________________________________________________ 
+  //
+
+  TFile *fileOut = new TFile("simple.root","RECREATE");
+
+  hQ->Write();
+  hQdedx->Write();
+  hQtr->Write();
+
+  hQX->Write();
+  hQXdedx->Write();
+  hQXtr->Write();
+
+  hNstep->Write();
+  hNel->Write();
+
+  hAmp->Write();
+  hAmpTime->Write();
+
+  fileOut->Close();
+
+}
diff --git a/TRD/AliTRDsdigits2digits.C b/TRD/AliTRDsdigits2digits.C
new file mode 100644 (file)
index 0000000..4f82ac9
--- /dev/null
@@ -0,0 +1,54 @@
+void AliTRDsdigits2digits()
+{
+
+  /////////////////////////////////////////////////////////////////////////
+  //
+  // Converts s-digits to normal digits
+  //
+  /////////////////////////////////////////////////////////////////////////
+
+  // Dynamically link some shared libs
+  if (gClassTable->GetID("AliRun") < 0) {
+    gROOT->LoadMacro("loadlibs.C");
+    loadlibs();
+    cout << "Loaded shared libraries" << endl;
+  }
+
+  Char_t *fileName = "galice.root";
+
+  // Create the TRD digits merger
+  AliTRDdigitizer *digitizer = new AliTRDdigitizer("TRDdigitizer"
+                                                  ,"TRD digitizer class");  
+
+  // Set the parameter
+  digitizer->SetDebug(1);
+
+  // Initialize the geometry 
+  digitizer->Open(fileName);
+
+  // Define the parameter object
+  // If no external parameter object is defined, 
+  // default parameter will be used
+  AliTRDparameter *parameter = new AliTRDparameter("TRDparameter"
+                                                 ,"TRD parameter class");
+  digitizer->SetParameter(parameter);
+
+  // Create the digits manager for the input s-digits
+  AliTRDdigitsManager *sdigitsManager = new AliTRDdigitsManager();
+  sdigitsManager->SetDebug(1);
+  sdigitsManager->SetSDigits(kTRUE);
+  sdigitsManager->ReadDigits();
+
+  // Add the s-digits to the input list 
+  digitizer->AddSDigitsManager(sdigitsManager);
+
+  // Convert the s-digits to normal digits
+  digitizer->SDigits2Digits();
+
+  // Store the digits
+  digitizer->WriteDigits();
+
+  // Save the parameter object in the AliROOT file
+  parameter->Write();
+
+}
index d7a07b2..cff2b54 100644 (file)
@@ -18,11 +18,9 @@ Int_t AliTRDtest()
   gAlice = (AliRun *) file->Get("gAlice");
 
   // Analyze the TRD hits
-  gROOT->LoadMacro("$(ALICE_ROOT)/TRD/AliTRDanalyzeHits.C");
   if (rc = AliTRDanalyzeHits()) return rc;
 
   // Run the digitization
-  gROOT->LoadMacro("$(ALICE_ROOT)/TRD/AliTRDcreateDigits.C");
   if (rc = AliTRDcreateDigits()) return rc;
 
   if (gAlice) delete gAlice;
@@ -30,11 +28,9 @@ Int_t AliTRDtest()
   gAlice = (AliRun *) file->Get("gAlice");
 
   // Analyze the digits
-  gROOT->LoadMacro("$(ALICE_ROOT)/TRD/AliTRDanalyzeDigits.C");
   if (rc = AliTRDanalyzeDigits()) return rc;
 
   // Create the cluster
-  gROOT->LoadMacro("$(ALICE_ROOT)/TRD/AliTRDcreateCluster.C");
   if (rc = AliTRDcreateCluster()) return rc;
 
   if (gAlice) delete gAlice;
@@ -42,9 +38,686 @@ Int_t AliTRDtest()
   gAlice = (AliRun *) file->Get("gAlice");
 
   // Analyze the digits
-  gROOT->LoadMacro("$(ALICE_ROOT)/TRD/AliTRDanalyzeCluster.C");
   if (rc = AliTRDanalyzeCluster()) return rc;
 
   return rc;
 
 }
+
+//_____________________________________________________________________________
+Int_t AliTRDanalyzeHits()
+{
+  //
+  // Analyzes the hits and fills QA-histograms 
+  //
+
+  Int_t rc = 0;
+
+  if (!gAlice) {
+    cout << "<AliTRDanalyzeHits> No AliRun object found" << endl;
+    rc = 1;
+    return rc;
+  }
+  gAlice->GetEvent(0);
+
+  // Get the pointer to the TRD detector 
+  AliTRD *trd = (AliTRD *) gAlice->GetDetector("TRD");
+  if (!trd) {
+    cout << "<AliTRDanalyzeHits> No TRD detector found" << endl;
+    rc = 2;
+    return rc;
+  }
+
+  // Get the pointer to the geometry object
+  AliTRDgeometryFull *geo;
+  if (trd) {
+    geo = (AliTRDgeometryFull *) trd->GetGeometry();
+  }
+  else {
+    cout << "<AliTRDanalyzeHits> No TRD geometry found" << endl;
+    rc = 3;
+    return rc;
+  }
+
+  AliTRDparameter *par = new AliTRDparameter("TRDparameter","TRD parameter class");
+
+  // Define the histograms
+  TH1F *hQdedx  = new TH1F("hQdedx","Charge dedx-hits",100,0.0,1000.0);
+  TH1F *hQtr    = new TH1F("hQtr"  ,"Charge TR-hits"  ,100,0.0,1000.0);
+
+  Float_t rmin   = geo->Rmin();
+  Float_t rmax   = geo->Rmax();
+  Float_t length = geo->GetChamberLength(0,2);
+  Float_t width  = geo->GetChamberWidth(0);
+  Int_t   ncol   = par->GetColMax(0);
+  Int_t   nrow   = par->GetRowMax(0,2,13);
+  Int_t   ntime  = ((Int_t) (rmax - rmin) / par->GetTimeBinSize());
+
+  TH2F *hZY     = new TH2F("hZY"   ,"Y vs Z (chamber 0)", nrow,-length/2.,length/2.
+                                                        ,ntime,      rmin,     rmax);
+  TH2F *hXZ     = new TH2F("hXZ"   ,"Z vs X (plane 0)"  , ncol, -width/2., width/2.
+                                                        , nrow,-length/2.,length/2.);
+
+  TH1F *hNprimP = new TH1F("hNprimP","Number of primary electrons per cm (MIP pion)"
+                                    ,50,0.0,100.0);
+  TH1F *hNprimE = new TH1F("hNprimE","Number of primary electrons per cm (3GeV electron)"
+                                    ,50,0.0,100.0);
+  TH1F *hNtotP  = new TH1F("hNtotP" ,"Number of electrons per cm (MIP pion)"
+                                    ,50,0.0,1000.0);
+  TH1F *hNtotE  = new TH1F("hNtotE" ,"Number of electrons per cm (3GeV electron)"
+                                    ,50,0.0,1000.0);
+
+  // Get the pointer hit tree
+  TTree *hitTree = gAlice->TreeH();  
+  if (!hitTree) {
+    cout << "<AliTRDanalyzeHits> No hit tree found" << endl;
+    rc = 4;
+    return rc;
+  }
+
+  Int_t countHits = 0;
+  Int_t nBytes    = 0;
+
+  // Get the number of entries in the hit tree
+  // (Number of primary particles creating a hit somewhere)
+  Int_t nTrack = (Int_t) hitTree->GetEntries();
+  cout << "<AliTRDanalyzeHits> Found " << nTrack 
+       << " primary particles with hits" << endl;
+
+  // Loop through all entries in the tree
+  for (Int_t iTrack = 0; iTrack < nTrack; iTrack++) {
+
+    gAlice->ResetHits();
+    nBytes += hitTree->GetEvent(iTrack);
+
+    Int_t nPrimE = 0;
+    Int_t nPrimP = 0;
+    Int_t nTotE  = 0;
+    Int_t nTotP  = 0;    
+
+    // Loop through the TRD hits  
+    Int_t iHit = 0;
+    AliTRDhit *hit = (AliTRDhit *) trd->FirstHit(-1);
+    while (hit) {
+
+      countHits++;
+      iHit++;
+
+      Float_t x     = hit->X();
+      Float_t y     = hit->Y();
+      Float_t z     = hit->Z();
+      Float_t q     = hit->GetCharge();
+      Int_t   track = hit->Track();
+      Int_t   det   = hit->GetDetector();
+      Int_t   plane = geo->GetPlane(det);
+
+      if      (q > 0) {
+        hQdedx->Fill(q);
+        hZY->Fill(z,y);
+        if (plane == 0) {
+          hXZ->Fill(x,z);
+       }
+      }
+      else if (q < 0) {
+        hQtr->Fill(TMath::Abs(q));
+      }
+
+      TParticle *part = gAlice->Particle(track);
+
+      if ((plane == 0) && (q > 0)) {
+
+        // 3 GeV electrons
+        if ((part->GetPdgCode() ==   11) && 
+            (part->P()          >   2.9)) {
+          nPrimE++;
+          nTotE += ((Int_t) q);
+        }
+
+        // MIP pions
+        if ((part->GetPdgCode() == -211) &&
+            (part->P()          >   0.5)) {
+          nPrimP++;
+          nTotP += ((Int_t) q);
+        }
+
+      }
+
+      hit = (AliTRDhit *) trd->NextHit();         
+
+    }
+
+    if (nPrimE > 0) hNprimE->Fill(((Double_t) nPrimE)/3.7);
+    if (nPrimP > 0) hNprimP->Fill(((Double_t) nPrimP)/3.7);
+    if (nTotE  > 0) hNtotE->Fill(((Double_t) nTotE)/3.7);
+    if (nTotP  > 0) hNtotP->Fill(((Double_t) nTotP)/3.7);
+
+  }
+
+  cout << "<AliTRDanalyzeHits> Found " << countHits << " hits in total" << endl;
+
+  TCanvas *cHits = new TCanvas("cHits","AliTRDanalyzeHits",50,50,600,600);
+  cHits->Divide(2,2);
+  cHits->cd(1);
+  hXZ->Draw("COL");
+  cHits->cd(2);
+  hZY->Draw("COL");
+  cHits->cd(3);
+  gPad->SetLogy();
+  hQdedx->Draw();
+  cHits->cd(4);
+  gPad->SetLogy();
+  hQtr->Draw();
+
+  TCanvas *cNel = new TCanvas("cNel","Ionization",50,50,600,600);
+  cNel->Divide(2,2);
+  cNel->cd(1);
+  hNprimE->SetStats();
+  hNprimE->Draw();
+  cNel->cd(2);
+  hNprimP->SetStats();
+  hNprimP->Draw();
+  cNel->cd(3);
+  hNtotE->SetStats();
+  hNtotE->Draw();
+  cNel->cd(4);
+  hNtotP->SetStats();
+  hNtotP->Draw();
+
+  TFile *fout = new TFile("TRD_ionization.root","RECREATE");
+  hNprimE->Write();
+  hNprimP->Write();
+  hNtotE->Write();
+  hNtotP->Write();
+  fout->Close(); 
+
+  return rc;
+
+}
+
+//_____________________________________________________________________________
+Int_t AliTRDcreateDigits()
+{
+  //
+  // Creates the digits from the hits of the slow simulator
+  //
+
+  Int_t rc = 0;
+
+  if (!gAlice) {
+    cout << "<AliTRDcreateDigits> No AliRun object found" << endl;
+    rc = 1;
+    return rc;
+  }
+  gAlice->GetEvent(0);
+
+  // Create the TRD digitzer 
+  AliTRDdigitizer *digitizer = new AliTRDdigitizer("TRDdigitizer"
+                                                  ,"TRD digitizer class");
+  digitizer->InitDetector();
+
+  // Set the parameter
+  digitizer->SetDebug(1);
+
+  // Define the parameter object
+  // If no external parameter object is defined, 
+  // default parameter will be used
+  AliTRDparameter *parameter = new AliTRDparameter("TRDparameter"
+                                                 ,"TRD parameter class");
+  digitizer->SetParameter(parameter);
+
+  // Create the digits
+  if (!(digitizer->MakeDigits())) {
+    rc = 2;
+    return rc;
+  }
+
+  // Write the digits into the input file
+  if (!(digitizer->MakeBranch())) {
+    rc = 3;
+    return rc;
+  }
+
+  // Write the digits into the input file
+  if (!(digitizer->WriteDigits())) {
+    rc = 4;
+    return rc;
+  }
+
+  // Save the parameter object in the AliROOT file
+  if (!(parameter->Write())) {
+    rc = 4;
+    return rc;
+  }
+
+  return rc;
+
+}
+
+//_____________________________________________________________________________
+Int_t AliTRDanalyzeDigits()
+{
+  //
+  // Analyzes the digits
+  //
+
+  Int_t rc = 0;
+
+  const Int_t kNpla = AliTRDgeometry::Nplan();
+
+  if (!gAlice) {
+    cout << "<AliTRDanalyzeDigits> No AliRun object found" << endl;
+    rc = 1;
+    return rc;
+  }
+  Int_t nPart = gAlice->GetEvent(0);
+
+  // Get the pointer to the TRD detector 
+  AliTRD *trd = (AliTRD *) gAlice->GetDetector("TRD");
+  if (!trd) {
+    cout << "<AliTRDanalyzeDigits> No TRD detector found" << endl;
+    rc = 2;
+    return rc;
+  }
+
+  // Get the digitizer object
+  TFile *file = (TFile *) gROOT->GetListOfFiles()->FindObject("TRD_test.root");
+  AliTRDparameter *parameter = (AliTRDparameter *) file->Get("TRDparameter");
+  if (!parameter) {
+    cout << "<AliTRDanalyzeDigits> No parameter object found" << endl;
+    rc = 3;
+    return rc;
+  }
+
+  // Define the histograms
+  Int_t adcRange = ((Int_t) parameter->GetADCoutRange());
+  TH1F *hAmpAll   = new TH1F("hAmpAll"  ,"Amplitude of the digits (all)"
+                            ,adcRange+1,-0.5,((Float_t) adcRange)+0.5);
+  TH1F *hAmpEl    = new TH1F("hAmpEl"   ,"Amplitude of the digits (electrons)"
+                            ,adcRange+1,-0.5,((Float_t) adcRange)+0.5);
+  TH1F *hAmpPi    = new TH1F("hAmpPi"   ,"Amplitude of the digits (pions)"
+                            ,adcRange+1,-0.5,((Float_t) adcRange)+0.5);
+  TH1F *hAmpNoise = new TH1F("hAmpNoise","Amplitude of the digits (noise)"
+                            ,5,-0.5,4.5);
+
+  // Get the pointer to the geometry object
+  AliTRDgeometry *geo;
+  if (trd) {
+    geo = trd->GetGeometry();
+  }
+  else {
+    cout << "<AliTRDanalyzeDigits> No TRD geometry found" << endl;
+    rc = 4;
+    return rc;
+  }
+
+  // Create the digits manager
+  AliTRDdigitsManager *digitsManager = new AliTRDdigitsManager();
+
+  digitsManager->Open("TRD_test.root");
+
+  // Read the digits from the file
+  if (!(digitsManager->ReadDigits())) {
+    cout << "<AliTRDanalyzeDigits> Cannot read the digits" << endl;
+    rc = 5;
+    return rc;
+  }
+
+  AliTRDmatrix *matrix;
+
+  Int_t countDigits = 0;
+  Int_t iSec        = trd->GetSensSector();
+  Int_t iCha        = trd->GetSensChamber();
+  Int_t timeMax     = geo->GetTimeTotal();
+
+  TProfile *hAmpTimeEl = new TProfile("hAmpTimeEl","Amplitude of the digits (electrons)"
+                                     ,timeMax,-0.5,((Double_t) timeMax)-0.5);
+  TProfile *hAmpTimePi = new TProfile("hAmpTimePi","Amplitude of the digits (pions)"
+                                     ,timeMax,-0.5,((Double_t) timeMax)-0.5);
+
+  // Loop over all planes
+  for (Int_t iPla = 0; iPla < kNpla; iPla++) {
+
+    Int_t iDet   = geo->GetDetector(iPla,iCha,iSec);
+    Int_t rowMax = parameter->GetRowMax(iPla,iCha,iSec);
+    Int_t colMax = parameter->GetColMax(iPla);
+  
+    if (iPla == 0) {
+      matrix = new AliTRDmatrix(rowMax,colMax,timeMax,iSec,iCha,iPla);
+    }
+
+    // Loop through the detector pixel
+    for (Int_t time = 0; time < timeMax; time++) {
+      for (Int_t  col = 0;  col <  colMax;  col++) {
+        for (Int_t  row = 0;  row <  rowMax;  row++) {
+
+          AliTRDdigit *digit    = digitsManager->GetDigit(row,col,time,iDet);
+          Int_t        amp      = digit->GetAmp();
+          Int_t        track0   = digitsManager->GetTrack(0,row,col,time,iDet);
+          Int_t        track1   = digitsManager->GetTrack(1,row,col,time,iDet);
+          TParticle   *particle = 0;
+          if (track0 > -1) {
+            particle = gAlice->Particle(track0);
+         }
+
+          if (amp > 0) {
+            countDigits++;
+            if (iPla == 0) {
+              matrix->SetSignal(row,col,time,amp);
+           }
+         }
+
+         // Total spectrum
+          hAmpAll->Fill(amp);
+
+         // Noise spectrum
+          if (track0 < 0) {
+            hAmpNoise->Fill(amp);
+         }          
+
+         // Electron digit
+          if ((particle) && (particle->GetPdgCode() ==   11) && (track1 < 0)) {
+            hAmpEl->Fill(amp);
+            hAmpTimeEl->Fill(time,amp);
+         }
+
+          // Pion digit
+          if ((particle) && (particle->GetPdgCode() == -211) && (track1 < 0)) {
+            hAmpPi->Fill(amp);
+            hAmpTimePi->Fill(time,amp);
+         }
+
+          delete digit;
+
+        }
+      }
+    }
+
+  }
+
+  cout << "<AliTRDanalyzeDigits> Found " << countDigits << " digits in total" << endl;
+
+  // Display the detector matrix
+  matrix->Draw();
+
+  TCanvas *cDigits = new TCanvas("cDigits","AliTRDanalyzeDigits",50,50,600,800);
+  cDigits->Divide(2,3);
+  cDigits->cd(1);
+  gPad->SetLogy();
+  hAmpAll->SetXTitle("Amplitude (ADC-channels)");
+  hAmpAll->SetYTitle("Entries");
+  hAmpAll->Draw();
+  cDigits->cd(2);
+  gPad->SetLogy();
+  hAmpNoise->SetXTitle("Amplitude (ADC-channels)");
+  hAmpNoise->SetYTitle("Entries");
+  hAmpNoise->Draw();
+  cDigits->cd(3);
+  gPad->SetLogy();
+  hAmpEl->SetXTitle("Amplitude (ADC-channels)");
+  hAmpEl->SetYTitle("Entries");
+  hAmpEl->Draw();
+  cDigits->cd(4);
+  gPad->SetLogy();
+  hAmpPi->SetXTitle("Amplitude (ADC-channels)");
+  hAmpPi->SetYTitle("Entries");
+  hAmpPi->Draw();
+  cDigits->cd(5);
+  hAmpTimeEl->SetXTitle("Timebin number");
+  hAmpTimeEl->SetYTitle("Mean amplitude");
+  hAmpTimeEl->Draw("HIST");
+  cDigits->cd(6);
+  hAmpTimePi->SetXTitle("Timebin number");
+  hAmpTimePi->SetYTitle("Mean amplitude");
+  hAmpTimePi->Draw("HIST");
+
+  TFile *fileOut = new TFile("digits_test.root","RECREATE");
+  hAmpAll->Write();
+  hAmpNoise->Write();
+  hAmpEl->Write();
+  hAmpPi->Write();
+  hAmpTimeEl->Write();
+  hAmpTimePi->Write();
+  fileOut->Close();
+
+  return rc;
+
+}
+
+//_____________________________________________________________________________
+Int_t AliTRDcreateCluster()
+{
+  //
+  // Creates the cluster from the digits
+  //
+
+  Int_t rc = 0;
+
+  // Create the clusterizer
+  AliTRDclusterizerV1 *clusterizer =  new AliTRDclusterizerV1("TRDclusterizer"
+                                                             ,"TRD clusterizer class");
+  clusterizer->SetVerbose(1);
+
+  // Define the parameter object
+  // If no external parameter object is defined, 
+  // default parameter will be used
+  AliTRDparameter *parameter = new AliTRDparameter("TRDparameter"
+                                                 ,"TRD parameter class");
+  parameter->SetClusMaxThresh(0);
+  parameter->SetClusSigThresh(0);
+  clusterizer->SetParameter(parameter);
+  // Open the file
+  if (!(clusterizer->Open("TRD_test.root",0))) {
+    rc = 1;
+    return rc;
+  }    
+
+  // Load the digits
+  if (!(clusterizer->ReadDigits())) {
+    rc = 2;
+    return rc;
+  }    
+
+  // Find the cluster
+  if (!(clusterizer->MakeClusters())) {
+    rc = 3;
+    return rc;
+  }
+
+  // Write the cluster tree into the file 
+  if (!(clusterizer->WriteClusters(-1))) {
+    rc = 4;
+    return rc;
+  }
+
+  // Save the clusterizer class in the file
+  if (!(clusterizer->Write())) {
+    rc = 5;
+    return rc;
+  }
+
+  return rc;
+
+}
+
+//_____________________________________________________________________________
+Int_t AliTRDanalyzeCluster()
+{
+  //
+  // Analyzes the cluster
+  //
+
+  Int_t rc = 0;
+
+  if (!gAlice) {
+    cout << "<AliTRDanalyzeCluster> No AliRun object found" << endl;
+    rc = 1;
+    return rc;
+  }
+  gAlice->GetEvent(0);
+
+  // Get the pointer to the TRD detector 
+  AliTRD *trd = (AliTRD *) gAlice->GetDetector("TRD");
+  if (!trd) {
+    cout << "<AliTRDanalyzeCluster> No TRD detector found" << endl;
+    rc = 2;
+    return rc;
+  }
+
+  // Define the histograms
+  TH1F *hClusAll   = new TH1F("hClusAll"  ,"Amplitude of the cluster (all)"     
+                                          ,501,-0.5,500.5);
+  TH1F *hClusNoise = new TH1F("hClusNoise","Amplitude of the cluster (noise)"   
+                                          ,  5,-0.5,  4.5);
+  TH1F *hClusEl    = new TH1F("hClusEl"   ,"Amplitude of the cluster (electron)"
+                                          ,501,-0.5,500.5);
+  TH1F *hClusPi    = new TH1F("hClusPi"   ,"Amplitude of the cluster (pion)"    
+                                          ,501,-0.5,500.5);
+
+  // Get the pointer to the geometry object
+  AliTRDgeometry *geo;
+  if (trd) {
+    geo = trd->GetGeometry();
+  }
+  else {
+    cout << "<AliTRDanalyzeCluster> No TRD geometry found" << endl;
+    rc = 3;
+    return rc;
+  }
+
+  // Get the pointer to the hit-tree
+  TFile *file = (TFile *) gROOT->GetListOfFiles()->FindObject("TRD_test.root");
+  TTree *clusterTree = (TTree *) file->Get("TreeR0_TRD");
+  if (!(clusterTree)) {
+    cout << "<AliTRDanalyzeCluster> No tree with clusters found" << endl;
+    rc = 4;
+    return rc;
+  }
+
+  // Get the pointer to the hit container
+  TObjArray *clusterArray = trd->RecPoints();
+  if (!(clusterArray)) {
+    cout << "<AliTRDanalyzeCluster> No clusterArray found" << endl;
+    rc = 5;
+    return rc;
+  }
+
+  // Set the branch address
+  clusterTree->GetBranch("TRDcluster")->SetAddress(&clusterArray);
+  Int_t nEntries = clusterTree->GetEntries();
+  cout << "<AliTRDanalyzeCluster> Number of entries in the cluster tree = " 
+       << nEntries 
+       << endl;
+
+  Int_t countCluster = 0;
+  Int_t countOverlap = 0;
+
+  // Loop through all entries in the tree
+  Int_t nbytes;
+  for (Int_t iEntry = 0; iEntry < nEntries; iEntry++) {
+
+    // Import the tree
+    nbytes += clusterTree->GetEvent(iEntry);
+
+    // Get the number of points in the detector 
+    Int_t nCluster = clusterArray->GetEntriesFast();
+
+    // Loop through all TRD digits
+    for (Int_t iCluster = 0; iCluster < nCluster; iCluster++) {
+
+      // Get the information for this digit
+      AliTRDcluster *cluster = (AliTRDcluster *) clusterArray->UncheckedAt(iCluster);
+      Int_t    detector = cluster->GetDetector();      
+      Int_t    sector   = geo->GetSector(detector);
+      Int_t    plane    = geo->GetPlane(detector);
+      Int_t    chamber  = geo->GetChamber(detector);
+      Float_t  energy   = cluster->GetQ();
+      Int_t    track0   = cluster->GetTrackIndex(0);
+      Int_t    track1   = cluster->GetTrackIndex(1);
+      Int_t    track2   = cluster->GetTrackIndex(2);
+      TParticle *particle = 0;
+      if (track0 > -1) {
+        particle = gAlice->Particle(track0);
+      }
+
+      countCluster++;
+      if (!cluster->Isolated()) countOverlap++;
+
+      // Total spectrum
+      hClusAll->Fill(energy);
+
+      if (cluster->Isolated()) {
+
+        // Noise spectrum
+        if (track0 < 0) {
+          hClusNoise->Fill(energy);
+        }          
+
+        // Electron cluster
+        if ((particle) && (particle->GetPdgCode() ==   11) && (track1 < 0)) {
+          hClusEl->Fill(energy);
+        }
+
+        // Pion cluster
+        if ((particle) && (particle->GetPdgCode() == -211) && (track1 < 0)) {
+          hClusPi->Fill(energy);
+        }
+
+      }
+
+    }
+
+  }
+
+  cout << "<AliTRDanalyzeCluster> Found " << countCluster << " cluster in total"    << endl;
+  cout << "<AliTRDanalyzeCluster> Found " << countOverlap << " overlapping cluster" << endl;
+  cout << endl;
+
+  TCanvas *cCluster = new TCanvas("cCluster","AliTRDanalyzeCluster",50,50,600,600);
+  cCluster->Divide(2,2);
+
+  TF1 *fun;
+  cCluster->cd(1);
+  gPad->SetLogy();
+  hClusAll->Fit("landau","0");
+  fun = (TF1 *) hClusAll->GetListOfFunctions()->First();
+  Float_t meanAll = fun->GetParameter(1);
+  hClusAll->Draw();
+  fun->SetLineColor(2);
+  fun->Draw("SAME");
+  
+  //cCluster->cd(2);
+  //gPad->SetLogy();
+  Float_t meanNoise = hClusNoise->GetMean();
+  //hClusNoise->Draw();
+
+  cCluster->cd(3);
+  gPad->SetLogy();
+  hClusEl->Fit("landau","0");
+  fun = (TF1 *) hClusEl->GetListOfFunctions()->First();
+  fun->SetLineColor(2);
+  Float_t meanEl = fun->GetParameter(1);
+  hClusEl->Draw();
+  fun->Draw("SAME");
+
+  cCluster->cd(4);
+  gPad->SetLogy();
+  hClusPi->Fit("landau","0");
+  fun = (TF1 *) hClusPi->GetListOfFunctions()->First();
+  fun->SetLineColor(2);
+  Float_t meanPi = fun->GetParameter(1);
+  hClusPi->Draw();
+  fun->Draw("SAME");
+
+  cout << endl;
+  cout << "##################################################################" << endl;
+  cout << "    Mean all       = " << meanAll   << endl;
+  cout << "    Mean noise     = " << meanNoise << endl;
+  cout << "    Mean electrons = " << meanEl    << endl;
+  cout << "    Mean pions     = " << meanPi    << endl;
+  cout << "##################################################################" << endl;
+  cout << endl;
+
+  return rc;
+
+}