TRD/Macros/AliTRDtest.C

   1
   2 Int_t AliTRDtest()
   3 {
   4   //
   5   // Test macro for the TRD code
   6   //
   7
   8   Int_t rc = 0;
   9
  10   AliSimulation sim;
  11   sim.SetConfigFile("$(ALICE_ROOT)/TRD/Macros/AliTRDconfig.C");
  12   sim.SetLoadAlignFromCDB(0);
  13   sim.Run();
  14
  15   // Analyze the TRD hits
  16   if (rc = AliTRDanalyzeHits()) return rc;
  17
  18   // Analyze the digits
  19   //if (rc = AliTRDanalyzeDigits()) return rc;
  20
  21   return rc;
  22
  23 }
  24
  25 //_____________________________________________________________________________
  26 Int_t AliTRDanalyzeHits()
  27 {
  28   //
  29   // Analyzes the hits and fills QA-histograms
  30   //
  31
  32   Int_t rc = 0;
  33
  34   //AliRunLoader *rl = gAlice->GetRunLoader();
  35   AliRunLoader *rl = AliRunLoader::Open("TRD_test.root"
  36                                        ,AliConfig::GetDefaultEventFolderName());
  37   if (!rl) {
  38     cout << "<AliTRDanalyzeHits> No RunLoader found" << endl;
  39     rc = 1;
  40     return rc;
  41   }
  42
  43   AliLoader* loader = rl->GetLoader("TRDLoader");
  44   if (!loader) {
  45     cout << "<AliTRDanalyzeHits> No TRDLoader found" << endl;
  46     rc = 2;
  47     return rc;
  48   }
  49
  50   rl->LoadgAlice();
  51   rl->LoadHeader();
  52   rl->LoadKinematics();
  53   rl->LoadHits();
  54
  55   // Get the pointer to the TRD detector
  56   gAlice = rl->GetAliRun();
  57   AliTRD *trd = (AliTRD *) gAlice->GetDetector("TRD");
  58   if (!trd) {
  59     cout << "<AliTRDanalyzeHits> No TRD detector found" << endl;
  60     rc = 2;
  61     return rc;
  62   }
  63
  64   // Get the pointer to the geometry object
  65   AliTRDgeometry *geo;
  66   if (trd) {
  67     geo = (AliTRDgeometry *) trd->GetGeometry();
  68   }
  69   else {
  70     cout << "<AliTRDanalyzeHits> No TRD geometry found" << endl;
  71     rc = 3;
  72     return rc;
  73   }
  74
  75   AliTRDCommonParam *par = AliTRDCommonParam::Instance();
  76
  77   // Define the histograms
  78   TH1F *hQdedx  = new TH1F("hQdedx","Charge dedx-hits",100,0.0,1000.0);
  79   TH1F *hQtr    = new TH1F("hQtr"  ,"Charge TR-hits"  ,100,0.0,1000.0);
  80
  81   Float_t rmin   = geo->Rmin();
  82   Float_t rmax   = geo->Rmax();
  83   Float_t length = geo->GetChamberLength(0,2);
  84   Float_t width  = geo->GetChamberWidth(0);
  85   Int_t   ncol   = par->GetColMax(0);
  86   Int_t   nrow   = par->GetRowMax(0,2,13);
  87   Int_t   ntime  = ((Int_t) (rmax - rmin) / 22.0);
  88
  89   TH2F *hZY     = new TH2F("hZY"   ,"Y vs Z (chamber 0)", nrow,-length/2.,length/2.
  90                                                         ,ntime,      rmin,     rmax);
  91   TH2F *hXZ     = new TH2F("hXZ"   ,"Z vs X (plane 0)"  , ncol, -width/2., width/2.
  92                                                         , nrow,-length/2.,length/2.);
  93
  94   // Get the pointer hit tree
  95   TTree *hitTree = loader->TreeH();
  96   if (!hitTree) {
  97     cout << "<AliTRDanalyzeHits> No hit tree found" << endl;
  98     rc = 4;
  99     return rc;
 100   }
 101
 102   Int_t countHits = 0;
 103   Int_t nBytes    = 0;
 104
 105   // Get the number of entries in the hit tree
 106   // (Number of primary particles creating a hit somewhere)
 107   Int_t nTrack = (Int_t) hitTree->GetEntries();
 108   cout << "<AliTRDanalyzeHits> Found " << nTrack
 109        << " primary particles with hits" << endl;
 110
 111   // Loop through all entries in the tree
 112   for (Int_t iTrack = 0; iTrack < nTrack; iTrack++) {
 113
 114     gAlice->ResetHits();
 115     nBytes += hitTree->GetEvent(iTrack);
 116
 117     // Loop through the TRD hits
 118     Int_t iHit = 0;
 119     AliTRDhit *hit = (AliTRDhit *) trd->FirstHit(-1);
 120     while (hit) {
 121
 122       countHits++;
 123       iHit++;
 124
 125       Float_t x     = hit->X();
 126       Float_t y     = hit->Y();
 127       Float_t z     = hit->Z();
 128       Float_t q     = hit->GetCharge();
 129       Int_t   track = hit->Track();
 130       Int_t   det   = hit->GetDetector();
 131       Int_t   plane = geo->GetPlane(det);
 132
 133       if      (q > 0) {
 134         hQdedx->Fill(q);
 135         hZY->Fill(z,y);
 136         if (plane == 0) {
 137           hXZ->Fill(x,z);
 138         }
 139       }
 140       else if (q < 0) {
 141         hQtr->Fill(TMath::Abs(q));
 142       }
 143
 144       hit = (AliTRDhit *) trd->NextHit();
 145
 146     }
 147
 148   }
 149
 150   cout << "<AliTRDanalyzeHits> Found " << countHits << " hits in total" << endl;
 151
 152   TCanvas *cHits = new TCanvas("cHits","AliTRDanalyzeHits",50,50,600,600);
 153   cHits->Divide(2,2);
 154   cHits->cd(1);
 155   hXZ->Draw("COL");
 156   cHits->cd(2);
 157   hZY->Draw("COL");
 158   cHits->cd(3);
 159   gPad->SetLogy();
 160   hQdedx->Draw();
 161   cHits->cd(4);
 162   gPad->SetLogy();
 163   hQtr->Draw();
 164
 165   return rc;
 166
 167 }
 168
 169 //_____________________________________________________________________________
 170 Int_t AliTRDanalyzeDigits()
 171 {
 172   //
 173   // Analyzes the digits
 174   //
 175
 176   Int_t rc = 0;
 177
 178   const Int_t kNpla = AliTRDgeometry::Nplan();
 179
 180   if (!gAlice) {
 181     cout << "<AliTRDanalyzeDigits> No AliRun object found" << endl;
 182     rc = 1;
 183     return rc;
 184   }
 185
 186   AliRunLoader *rl = gAlice->GetRunLoader();
 187   if (!rl) {
 188     cout << "<AliTRDanalyzeHits> No RunLoader found" << endl;
 189     rc = 2;
 190     return rc;
 191   }
 192
 193   // Import the Trees for the event nEvent in the file
 194   rl->LoadDigits();
 195
 196   AliLoader* loader = rl->GetLoader("TRDLoader");
 197   if (!loader) {
 198     cout << "<AliTRDanalyzeHits> No TRDLoader found" << endl;
 199     rc = 3;
 200     return rc;
 201   }
 202
 203   // Get the pointer to the TRD detector
 204   AliTRD *trd = (AliTRD *) gAlice->GetDetector("TRD");
 205   if (!trd) {
 206     cout << "<AliTRDanalyzeDigits> No TRD detector found" << endl;
 207     rc = 4;
 208     return rc;
 209   }
 210
 211   // Get the parameter object
 212   AliTRDSimParam    *parSim = AliTRDSimParam::Instance();
 213   AliTRDCommonParam *parCom = AliTRDCommonParam::Instance();
 214   AliTRDcalibDB     *cal    = AliTRDcalibDB::Instance();
 215
 216   // Define the histograms
 217   Int_t adcRange = ((Int_t) parSim->GetADCoutRange());
 218   TH1F *hAmpAll   = new TH1F("hAmpAll"  ,"Amplitude of the digits (all)"
 219                             ,adcRange+1,-0.5,((Float_t) adcRange)+0.5);
 220   TH1F *hAmpEl    = new TH1F("hAmpEl"   ,"Amplitude of the digits (electrons)"
 221                             ,adcRange+1,-0.5,((Float_t) adcRange)+0.5);
 222   TH1F *hAmpPi    = new TH1F("hAmpPi"   ,"Amplitude of the digits (pions)"
 223                             ,adcRange+1,-0.5,((Float_t) adcRange)+0.5);
 224   TH1F *hAmpNoise = new TH1F("hAmpNoise","Amplitude of the digits (noise)"
 225                             ,5,-0.5,4.5);
 226
 227   // Get the pointer to the geometry object
 228   AliTRDgeometry *geo;
 229   if (trd) {
 230     geo = trd->GetGeometry();
 231   }
 232   else {
 233     cout << "<AliTRDanalyzeDigits> No TRD geometry found" << endl;
 234     rc = 5;
 235     return rc;
 236   }
 237
 238   // Create the digits manager
 239   AliTRDdigitsManager *digitsManager = new AliTRDdigitsManager();
 240   digitsManager->SetDebug(1);
 241
 242   // Read the digits from the file
 243   if (!(digitsManager->ReadDigits(loader->TreeD()))) {
 244     cout << "<AliTRDanalyzeDigits> Cannot read the digits" << endl;
 245     rc = 6;
 246     return rc;
 247   }
 248
 249   // Get the particle stack
 250   AliStack *kineStack = rl->Stack();
 251   if (!kineStack) {
 252     cout << "<AliTRDanalyzeDigits> Cannot find the KINE stack" << endl;
 253     rc = 7;
 254     return rc;
 255   }
 256
 257   Int_t countDigits = 0;
 258   Int_t iSec = 0;
 259   Int_t iCha = 2;
 260   Int_t timeMax     = cal->GetNumberOfTimeBins();
 261
 262   TProfile *hAmpTimeEl = new TProfile("hAmpTimeEl","Amplitude of the digits (electrons)"
 263                                       ,timeMax,-0.5,((Double_t) timeMax)-0.5);
 264   TProfile *hAmpTimePi = new TProfile("hAmpTimePi","Amplitude of the digits (pions)"
 265                                       ,timeMax,-0.5,((Double_t) timeMax)-0.5);
 266
 267   // Loop over all planes
 268   for (Int_t iPla = 0; iPla < kNpla; iPla++) {
 269
 270     Int_t iDet   = geo->GetDetector(iPla,iCha,iSec);
 271     Int_t rowMax = parCom->GetRowMax(iPla,iCha,iSec);
 272     Int_t colMax = parCom->GetColMax(iPla);
 273
 274     // Loop through the detector pixel
 275     for (Int_t time = 0; time < timeMax; time++) {
 276       for (Int_t  col = 0;  col <  colMax;  col++) {
 277         for (Int_t  row = 0;  row <  rowMax;  row++) {
 278
 279           AliTRDdigit *digit    = digitsManager->GetDigit(row,col,time,iDet);
 280           Int_t        amp      = digit->GetAmp();
 281           Int_t        track0   = digitsManager->GetTrack(0,row,col,time,iDet);
 282           Int_t        track1   = digitsManager->GetTrack(1,row,col,time,iDet);
 283           TParticle   *particle = 0;
 284           if (track0 > -1) {
 285             particle = (TParticle *) kineStack->Particle(track0);
 286           }
 287
 288           if (amp > 0) {
 289             countDigits++;
 290           }
 291
 292           // Total spectrum
 293           hAmpAll->Fill(amp);
 294
 295           // Noise spectrum
 296           if (track0 < 0) {
 297             hAmpNoise->Fill(amp);
 298           }
 299
 300           // Electron digit
 301           if ((particle) && (particle->GetPdgCode() ==   11) && (track1 < 0)) {
 302             hAmpEl->Fill(amp);
 303             hAmpTimeEl->Fill(time,amp);
 304           }
 305
 306           // Pion digit
 307           if ((particle) && (particle->GetPdgCode() == -211) && (track1 < 0)) {
 308             hAmpPi->Fill(amp);
 309             hAmpTimePi->Fill(time,amp);
 310           }
 311
 312           delete digit;
 313
 314         }
 315       }
 316     }
 317
 318   }
 319
 320   cout << "<AliTRDanalyzeDigits> Found " << countDigits << " digits in total" << endl;
 321
 322   TCanvas *cDigits = new TCanvas("cDigits","AliTRDanalyzeDigits",100,100,600,800);
 323   cDigits->Divide(2,3);
 324   cDigits->cd(1);
 325   gPad->SetLogy();
 326   hAmpAll->SetXTitle("Amplitude (ADC-channels)");
 327   hAmpAll->SetYTitle("Entries");
 328   hAmpAll->Draw();
 329   cDigits->cd(2);
 330   gPad->SetLogy();
 331   hAmpNoise->SetXTitle("Amplitude (ADC-channels)");
 332   hAmpNoise->SetYTitle("Entries");
 333   hAmpNoise->Draw();
 334   cDigits->cd(3);
 335   gPad->SetLogy();
 336   hAmpEl->SetXTitle("Amplitude (ADC-channels)");
 337   hAmpEl->SetYTitle("Entries");
 338   hAmpEl->Draw();
 339   cDigits->cd(4);
 340   gPad->SetLogy();
 341   hAmpPi->SetXTitle("Amplitude (ADC-channels)");
 342   hAmpPi->SetYTitle("Entries");
 343   hAmpPi->Draw();
 344   cDigits->cd(5);
 345   hAmpTimeEl->SetXTitle("Timebin number");
 346   hAmpTimeEl->SetYTitle("Mean amplitude");
 347   hAmpTimeEl->Draw("HIST");
 348   cDigits->cd(6);
 349   hAmpTimePi->SetXTitle("Timebin number");
 350   hAmpTimePi->SetYTitle("Mean amplitude");
 351   hAmpTimePi->Draw("HIST");
 352
 353   return rc;
 354
 355 }
 356