[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;

  AliRunLoader *rl = gAlice->GetRunLoader();
  if (!rl) {
    cout << "<AliTRDanalyzeHits> No RunLoader found" << endl;
    rc = 1;
    return rc;
  }

  // Import the Trees for the event nEvent in the file
  rl->GetEvent(0);
  
  AliLoader* loader = rl->GetLoader("TRDLoader");
  if (!loader) {
    cout << "<AliTRDanalyzeHits> No TRDLoader found" << endl;
    rc = 2;
    return rc;
  }

  // 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 = loader->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));
      }

//        printf("Getting TParticle for track %d\n",track);
//        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->Open("TRD_test.root");
  digitizer->InitDetector();
  digitizer->InitOutput(0);

  // 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     = parameter->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)"   
                                          , 11,-0.5, 10.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->GetLabel(0);
      Int_t    track1   = cluster->GetLabel(1);
      Int_t    track2   = cluster->GetLabel(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
1f4ee546	19	//if (gAlice) delete gAlice;
	20	//file = (TFile *) gROOT->GetListOfFiles()->FindObject("TRD_test.root");
	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
1f4ee546	29	// if (gAlice) delete gAlice;
	30	// file = (TFile *) gROOT->GetListOfFiles()->FindObject("TRD_test.root");
	31	// gAlice = (AliRun *) file->Get("gAlice");
6244debe	32
1f4ee546	33	// // Analyze the digits
1f4ee546	34	// if (rc = AliTRDanalyzeDigits()) return rc;
793ff80c	35
1f4ee546	36	// // Create the cluster
1f4ee546	37	// if (rc = AliTRDcreateCluster()) return rc;
4a0fe73c	38
1f4ee546	39	// if (gAlice) delete gAlice;
	40	// file = (TFile *) gROOT->GetListOfFiles()->FindObject("TRD_test.root");
	41	// gAlice = (AliRun *) file->Get("gAlice");
4a0fe73c	42
1f4ee546	43	// // Analyze the cluster
1f4ee546	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
1f4ee546	65	AliRunLoader *rl = gAlice->GetRunLoader();
	66	if (!rl) {
	67	cout << "<AliTRDanalyzeHits> No RunLoader found" << endl;
fa148e6c	68	rc = 1;
	69	return rc;
	70	}
1f4ee546	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	}
fa148e6c	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
1f4ee546	130	TTree *hitTree = loader->TreeH();
fa148e6c	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
1f4ee546	184	// printf("Getting TParticle for track %d\n",track);
1f4ee546	185	// TParticle *part = gAlice->Particle(track);
fa148e6c	186
1f4ee546	187	// if ((plane == 0) && (q > 0)) {
fa148e6c	188
1f4ee546	189	// // 3 GeV electrons
	190	// if ((part->GetPdgCode() == 11) &&
	191	// (part->P() > 2.9)) {
	192	// nPrimE++;
	193	// nTotE += ((Int_t) q);
	194	// }
fa148e6c	195
1f4ee546	196	// // MIP pions
	197	// if ((part->GetPdgCode() == -211) &&
	198	// (part->P() > 0.5)) {
	199	// nPrimP++;
	200	// nTotP += ((Int_t) q);
	201	// }
fa148e6c	202
1f4ee546	203	// }
fa148e6c	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");
1f4ee546	276	digitizer->Open("TRD_test.root");
fa148e6c	277	digitizer->InitDetector();
1f4ee546	278	digitizer->InitOutput(0);
fa148e6c	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
1f4ee546	297	//if (!(digitizer->MakeBranch())) {
	298	// rc = 3;
	299	// return rc;
	300	//}
fa148e6c	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();
8e64dd77	377	digitsManager->SetDebug(1);
fa148e6c	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();
f6ff4cec	393	Int_t timeMax = parameter->GetTimeTotal();
fa148e6c	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)"
963b25ce	594	, 11,-0.5, 10.5);
fa148e6c	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();
f6ff4cec	658	Int_t track0 = cluster->GetLabel(0);
	659	Int_t track1 = cluster->GetLabel(1);
	660	Int_t track2 = cluster->GetLabel(2);
fa148e6c	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
963b25ce	712	cCluster->cd(2);
963b25ce	713	gPad->SetLogy();
fa148e6c	714	Float_t meanNoise = hClusNoise->GetMean();
963b25ce	715	hClusNoise->Draw();
fa148e6c	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	}
8e64dd77	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	}