TRD/AliTRDanalyzeCluster.C

   1 Int_t AliTRDanalyzeCluster()
   2 {
   3   //
   4   // Analyzes the cluster
   5   //
   6
   7   Int_t rc = 0;
   8
   9   if (!gAlice) {
  10     cout << "<AliTRDanalyzeCluster> No AliRun object found" << endl;
  11     rc = 1;
  12     return rc;
  13   }
  14   gAlice->GetEvent(0);
  15
  16   // Get the pointer to the TRD detector
  17   AliTRD *TRD = (AliTRD *) gAlice->GetDetector("TRD");
  18   if (!TRD) {
  19     cout << "<AliTRDanalyzeCluster> No TRD detector found" << endl;
  20     rc = 2;
  21     return rc;
  22   }
  23
  24   // Define the histograms
  25   TH1F *hClusAll   = new TH1F("hClusAll"  ,"Amplitude of the cluster (all)"     ,501,-0.5,500.5);
  26   TH1F *hClusNoise = new TH1F("hClusNoise","Amplitude of the cluster (noise)"   ,501,-0.5,500.5);
  27   TH1F *hClusEl    = new TH1F("hClusEl"   ,"Amplitude of the cluster (electron)",501,-0.5,500.5);
  28   TH1F *hClusPi    = new TH1F("hClusPi"   ,"Amplitude of the cluster (pion)"    ,501,-0.5,500.5);
  29
  30   // Get the pointer to the geometry object
  31   AliTRDgeometry *TRDgeometry;
  32   if (TRD) {
  33     TRDgeometry = TRD->GetGeometry();
  34   }
  35   else {
  36     cout << "<AliTRDanalyzeCluster> No TRD geometry found" << endl;
  37     rc = 3;
  38     return rc;
  39   }
  40
  41   // Get the pointer to the hit-tree
  42   TFile *file = (TFile *) gROOT->GetListOfFiles()->FindObject("TRD_test.root");
  43   TTree *RecTree = file->Get("ClusterTree");
  44   if (!(RecTree)) {
  45     cout << "<AliTRDanalyzeCluster> No tree with rec. points found" << endl;
  46     rc = 4;
  47     return rc;
  48   }
  49
  50   // Get the pointer to the hit container
  51   TObjArray *RecPointArray = TRD->RecPoints();
  52   if (!(RecPointArray)) {
  53     cout << "<AliTRDanalyzeCluster> No RecPointArray found" << endl;
  54     rc = 5;
  55     return rc;
  56   }
  57
  58   // Set the branch address
  59   RecTree->GetBranch("TRDrecPoints")->SetAddress(&RecPointArray);
  60   Int_t nEntries = RecTree->GetEntries();
  61   cout << "Number of entries in the rec. point tree = " << nEntries << endl;
  62
  63   Int_t countCluster = 0;
  64
  65   // Loop through all entries in the tree
  66   Int_t nbytes;
  67   for (Int_t iEntry = 0; iEntry < nEntries; iEntry++) {
  68
  69     // Import the tree
  70     nbytes += RecTree->GetEvent(iEntry);
  71
  72     // Get the number of points in the detector
  73     Int_t nRecPoint = RecPointArray->GetEntriesFast();
  74
  75     // Loop through all TRD digits
  76     for (Int_t iRecPoint = 0; iRecPoint < nRecPoint; iRecPoint++) {
  77
  78       // Get the information for this digit
  79       AliTRDrecPoint *RecPoint = (AliTRDrecPoint *) RecPointArray->UncheckedAt(iRecPoint);
  80       Int_t    detector = RecPoint->GetDetector();
  81       Int_t    sector   = TRDgeometry->GetSector(detector);
  82       Int_t    plane    = TRDgeometry->GetPlane(detector);
  83       Int_t    chamber  = TRDgeometry->GetChamber(detector);
  84       Int_t    energy   = RecPoint->GetEnergy();
  85       TVector3 pos;
  86       RecPoint->GetLocalPosition(pos);
  87       Int_t    track0   = RecPoint->GetTrackIndex(0);
  88       Int_t    track1   = RecPoint->GetTrackIndex(1);
  89       TParticle *Part = 0;
  90       if (track0 > -1) {
  91         Part = gAlice->Particle(track0);
  92       }
  93
  94       countCluster++;
  95
  96       // Total spectrum
  97       hClusAll->Fill((Float_t) energy);
  98
  99       // Noise spectrum
 100       if (track0 < 0) {
 101         hClusNoise->Fill((Float_t) energy);
 102       }
 103
 104       // Electron cluster
 105       if ((Part) && (Part->GetPdgCode() ==   11) && (track1 < 0)) {
 106         hClusEl->Fill((Float_t) energy);
 107       }
 108
 109       // Pion cluster
 110       if ((Part) && (Part->GetPdgCode() == -211) && (track1 < 0)) {
 111         hClusPi->Fill((Float_t) energy);
 112       }
 113
 114     }
 115
 116   }
 117
 118   cout << "<AliTRDanalyzeCluster> Found " << countCluster << " cluster in total" << endl;
 119
 120   TCanvas *cCluster = new TCanvas("cCluster","AliTRDanalyzeCluster",50,50,600,600);
 121   cCluster->Divide(2,2);
 122
 123   TF1 *fun;
 124   cCluster->cd(1);
 125   gPad->SetLogy();
 126   hClusAll->Fit("landau","0");
 127   fun = (TF1 *) hClusAll->GetListOfFunctions()->First();
 128   Float_t meanAll = fun->GetParameter(1);
 129   hClusAll->Draw();
 130   fun->SetLineColor(2);
 131   fun->Draw("SAME");
 132
 133   cCluster->cd(2);
 134   gPad->SetLogy();
 135   gPad->SetLogx();
 136   Float_t meanNoise = hClusNoise->GetMean();
 137   hClusNoise->Draw();
 138
 139   cCluster->cd(3);
 140   gPad->SetLogy();
 141   hClusEl->Fit("landau","0");
 142   fun = (TF1 *) hClusEl->GetListOfFunctions()->First();
 143   fun->SetLineColor(2);
 144   Float_t meanEl = fun->GetParameter(1);
 145   hClusEl->Draw();
 146   fun->Draw("SAME");
 147
 148   cCluster->cd(4);
 149   gPad->SetLogy();
 150   hClusPi->Fit("landau","0");
 151   fun = (TF1 *) hClusPi->GetListOfFunctions()->First();
 152   fun->SetLineColor(2);
 153   Float_t meanPi = fun->GetParameter(1);
 154   hClusPi->Draw();
 155   fun->Draw("SAME");
 156
 157   cout << endl;
 158   cout << "##################################################################" << endl;
 159   cout << "    Mean all       = " << meanAll   << endl;
 160   cout << "    Mean noise     = " << meanNoise << endl;
 161   cout << "    Mean electrons = " << meanEl    << endl;
 162   cout << "    Mean pions     = " << meanPi    << endl;
 163   cout << "##################################################################" << endl;
 164   cout << endl;
 165
 166   return rc;
 167
 168 }