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