TRD/AliTRDanalyzeHits.C

   1 Int_t AliTRDanalyzeHits()
   2 {
   3   //
   4   // Analyzes the hits and fills QA-histograms
   5   //
   6
   7   Int_t rc = 0;
   8
   9   if (!gAlice) {
  10     cout << "<AliTRDanalyzeHits> 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 << "<AliTRDanalyzeHits> No TRD detector found" << endl;
  20     rc = 2;
  21     return rc;
  22   }
  23
  24   // Get the pointer to the geometry object
  25   AliTRDgeometryFull *geo;
  26   if (trd) {
  27     geo = (AliTRDgeometryFull *) trd->GetGeometry();
  28   }
  29   else {
  30     cout << "<AliTRDanalyzeHits> No TRD geometry found" << endl;
  31     rc = 3;
  32     return rc;
  33   }
  34
  35   // Define the histograms
  36   TH1F *hQdedx  = new TH1F("hQdedx","Charge dedx-hits",100,0.0,1000.0);
  37   TH1F *hQtr    = new TH1F("hQtr"  ,"Charge TR-hits"  ,100,0.0,1000.0);
  38
  39   Float_t rmin   = geo->Rmin();
  40   Float_t rmax   = geo->Rmax();
  41   Float_t length = geo->GetChamberLength(0,2);
  42   Float_t width  = geo->GetChamberWidth(0);
  43   Int_t   ncol   = geo->GetColMax(0);
  44   Int_t   nrow   = geo->GetRowMax(0,2,13);
  45   Int_t   ntime  = ((Int_t) (rmax - rmin) / geo->GetTimeBinSize());
  46
  47   TH2F *hZY     = new TH2F("hZY"   ,"Y vs Z (chamber 0)", nrow,-length/2.,length/2.
  48                                                         ,ntime,      rmin,     rmax);
  49   TH2F *hXZ     = new TH2F("hXZ"   ,"Z vs X (plane 0)"  , ncol, -width/2., width/2.
  50                                                         , nrow,-length/2.,length/2.);
  51
  52   TH1F *hNprimP = new TH1F("hNprimP","Number of primary electrons per cm (MIP pion)"
  53                                     ,50,0.0,100.0);
  54   TH1F *hNprimE = new TH1F("hNprimE","Number of primary electrons per cm (3GeV electron)"
  55                                     ,50,0.0,100.0);
  56   TH1F *hNtotP  = new TH1F("hNtotP" ,"Number of electrons per cm (MIP pion)"
  57                                     ,50,0.0,1000.0);
  58   TH1F *hNtotE  = new TH1F("hNtotE" ,"Number of electrons per cm (3GeV electron)"
  59                                     ,50,0.0,1000.0);
  60
  61   // Get the pointer hit tree
  62   TTree *hitTree = gAlice->TreeH();
  63   if (!hitTree) {
  64     cout << "<AliTRDanalyzeHits> No hit tree found" << endl;
  65     rc = 4;
  66     return rc;
  67   }
  68
  69   Int_t countHits = 0;
  70   Int_t nBytes    = 0;
  71
  72   // Get the number of entries in the hit tree
  73   // (Number of primary particles creating a hit somewhere)
  74   Int_t nTrack = (Int_t) hitTree->GetEntries();
  75   cout << "<AliTRDanalyzeHits> Found " << nTrack
  76        << " primary particles with hits" << endl;
  77
  78   // Loop through all entries in the tree
  79   for (Int_t iTrack = 0; iTrack < nTrack; iTrack++) {
  80
  81     gAlice->ResetHits();
  82     nBytes += hitTree->GetEvent(iTrack);
  83
  84     Int_t nPrimE = 0;
  85     Int_t nPrimP = 0;
  86     Int_t nTotE  = 0;
  87     Int_t nTotP  = 0;
  88
  89     // Loop through the TRD hits
  90     Int_t iHit = 0;
  91     AliTRDhit *hit = (AliTRDhit *) trd->FirstHit(-1);
  92     while (hit) {
  93
  94       countHits++;
  95       iHit++;
  96
  97       Float_t x     = hit->X();
  98       Float_t y     = hit->Y();
  99       Float_t z     = hit->Z();
 100       Float_t q     = hit->GetCharge();
 101       Int_t   track = hit->Track();
 102       Int_t   det   = hit->GetDetector();
 103       Int_t   plane = geo->GetPlane(det);
 104
 105       if      (q > 0) {
 106         hQdedx->Fill(q);
 107         hZY->Fill(z,y);
 108         if (plane == 0) {
 109           hXZ->Fill(x,z);
 110         }
 111       }
 112       else if (q < 0) {
 113         hQtr->Fill(TMath::Abs(q));
 114       }
 115
 116       TParticle *part = gAlice->Particle(track);
 117
 118       if ((plane == 0) && (q > 0)) {
 119
 120         // 3 GeV electrons
 121         if ((part->GetPdgCode() ==   11) &&
 122             (part->P()          >   2.9)) {
 123           nPrimE++;
 124           nTotE += ((Int_t) q);
 125         }
 126
 127         // MIP pions
 128         if ((part->GetPdgCode() == -211) &&
 129             (part->P()          >   0.5)) {
 130           nPrimP++;
 131           nTotP += ((Int_t) q);
 132         }
 133
 134       }
 135
 136       hit = (AliTRDhit *) trd->NextHit();
 137
 138     }
 139
 140     if (nPrimE > 0) hNprimE->Fill(((Double_t) nPrimE)/3.7);
 141     if (nPrimP > 0) hNprimP->Fill(((Double_t) nPrimP)/3.7);
 142     if (nTotE  > 0) hNtotE->Fill(((Double_t) nTotE)/3.7);
 143     if (nTotP  > 0) hNtotP->Fill(((Double_t) nTotP)/3.7);
 144
 145   }
 146
 147   cout << "<AliTRDanalyzeHits> Found " << countHits << " hits in total" << endl;
 148
 149   TCanvas *cHits = new TCanvas("cHits","AliTRDanalyzeHits",50,50,600,600);
 150   cHits->Divide(2,2);
 151   cHits->cd(1);
 152   hXZ->Draw("COL");
 153   cHits->cd(2);
 154   hZY->Draw("COL");
 155   cHits->cd(3);
 156   gPad->SetLogy();
 157   hQdedx->Draw();
 158   cHits->cd(4);
 159   gPad->SetLogy();
 160   hQtr->Draw();
 161
 162   TCanvas *cNel = new TCanvas("cNel","Ionization",50,50,600,600);
 163   cNel->Divide(2,2);
 164   cNel->cd(1);
 165   hNprimE->SetStats();
 166   hNprimE->Draw();
 167   cNel->cd(2);
 168   hNprimP->SetStats();
 169   hNprimP->Draw();
 170   cNel->cd(3);
 171   hNtotE->SetStats();
 172   hNtotE->Draw();
 173   cNel->cd(4);
 174   hNtotP->SetStats();
 175   hNtotP->Draw();
 176
 177   TFile *fout = new TFile("TRD_ionization.root","RECREATE");
 178   hNprimE->Write();
 179   hNprimP->Write();
 180   hNtotE->Write();
 181   hNtotP->Write();
 182   fout->Close();
 183
 184   return rc;
 185
 186 }