[u/mrichter/AliRoot.git] / TRD / AliTRDtest.C


Int_t AliTRDtest() 
{
  //
  // Test macro for the TRD code
  //

  Int_t rc = 0;

  TFile *file;

  // Initialize the test setup 
  gAlice->Init("$(ALICE_ROOT)/TRD/AliTRDconfig.C");

  // Run one event and create the hits
  gAlice->SetDebug(2);
  gAlice->Run(1);

  if (gAlice) delete gAlice;
  file   = (TFile *) gROOT->GetListOfFiles()->FindObject("TRD_test.root");
  gAlice = (AliRun *) file->Get("gAlice");

  // Analyze the TRD hits
  if (rc = AliTRDanalyzeHits()) return rc;

  // Run the digitization
  if (rc = AliTRDcreateDigits()) return rc;

  if (gAlice) delete gAlice;
  file   = (TFile *) gROOT->GetListOfFiles()->FindObject("TRD_test.root");
  gAlice = (AliRun *) file->Get("gAlice");

  // Analyze the digits
  if (rc = AliTRDanalyzeDigits()) return rc;

  // Create the cluster
  if (rc = AliTRDcreateCluster()) return rc;

  if (gAlice) delete gAlice;
  file   = (TFile *) gROOT->GetListOfFiles()->FindObject("TRD_test.root");
  gAlice = (AliRun *) file->Get("gAlice");

  // Analyze the cluster
  if (rc = AliTRDanalyzeCluster()) return rc;

  //file   = (TFile *) gROOT->GetListOfFiles()->FindObject("TRD_test.root");
  //file->Close();

  // Find the tracks
  //if (rc = AliTRDcreateTracks()) 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->SetDebug(1);

  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;

}

//_____________________________________________________________________________
Int_t AliTRDcreateTracks()
{
  //
  // Creates the tracks
  //

  Int_t rc = 0;

  // Create the tracker
  AliTRDtracker *tracker = new AliTRDtracker("TRDtracker","TRD tracker");

  // Read in the kine tree and the cluster
  tracker->GetEvent("TRD_test.root","TRD_test.root");
 
  // Find the tracks
  TH1F *hs = new TH1F("hs","hs",100,0.0,1.0);
  TH1F *hd = new TH1F("hd","hd",100,0.0,1.0);
  tracker->Clusters2Tracks(hs,hd);
 
  // Store the tracks
  tracker->WriteTracks("TRD_test.root");

  return rc;                                                                    

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