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

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

  Int_t rc = 0;

  // 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;
  TFile *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;
  TFile *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;
  TFile *file = (TFile *) gROOT->GetListOfFiles()->FindObject("TRD_test.root");
  gAlice = (AliRun *) file->Get("gAlice");

  // Analyze the digits
  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;

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