VZERO/AliVZERODigitizer.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id$ */
  17
  18 ///_________________________________________________________________________
  19 ///
  20 /// This class constructs Digits out of Hits
  21 ///
  22 ///
  23
  24 // --- Standard library ---
  25
  26 // --- ROOT system ---
  27 #include <TMath.h>
  28 #include <TTree.h>
  29 #include <TMap.h>
  30 #include <TGeoManager.h>
  31 #include <TGeoPhysicalNode.h>
  32 #include <AliGeomManager.h>
  33 #include <TRandom.h>
  34
  35 // --- AliRoot header files ---
  36 #include "AliVZEROConst.h"
  37 #include "AliRun.h"
  38 #include "AliVZERO.h"
  39 #include "AliVZEROhit.h"
  40 #include "AliRunLoader.h"
  41 #include "AliLoader.h"
  42 #include "AliGRPObject.h"
  43 #include "AliRunDigitizer.h"
  44 #include "AliCDBManager.h"
  45 #include "AliCDBStorage.h"
  46 #include "AliCDBEntry.h"
  47 #include "AliVZEROCalibData.h"
  48
  49 #include "AliVZEROdigit.h"
  50 #include "AliVZERODigitizer.h"
  51
  52 ClassImp(AliVZERODigitizer)
  53
  54  AliVZERODigitizer::AliVZERODigitizer()
  55                    :AliDigitizer(),
  56                     fCalibData(GetCalibData()),
  57                     fPhotoCathodeEfficiency(0.18),
  58                     fPMVoltage(768.0),
  59                     fPMGain(TMath::Power((fPMVoltage / 112.5) ,7.04277)),
  60                     fNdigits(0),
  61                     fDigits(0),
  62                     fCollisionMode(0),
  63                     fBeamEnergy(0.)
  64
  65 {
  66   // default constructor
  67
  68 //    fNdigits = 0;
  69 //    fDigits  = 0;
  70 //
  71 //    fPhotoCathodeEfficiency =   0.18;
  72 //    fPMVoltage              =  768.0;
  73 //    fPMGain = TMath::Power((fPMVoltage / 112.5) ,7.04277);
  74
  75 //   fCalibData = GetCalibData();
  76 }
  77
  78 //____________________________________________________________________________
  79   AliVZERODigitizer::AliVZERODigitizer(AliRunDigitizer* manager)
  80                     :AliDigitizer(manager),
  81                      fCalibData(GetCalibData()),
  82                      fPhotoCathodeEfficiency(0.18),
  83                      fPMVoltage(768.0),
  84                      fPMGain(TMath::Power((fPMVoltage / 112.5) ,7.04277)),
  85                      fNdigits(0),
  86                      fDigits(0),
  87                      fCollisionMode(0),
  88                      fBeamEnergy(0.)
  89
  90 {
  91   // constructor
  92
  93 //   fNdigits = 0;
  94 //   fDigits  = 0;
  95 //
  96 //   fPhotoCathodeEfficiency =   0.18;
  97 //   fPMVoltage              =  768.0;
  98 //   fPMGain = TMath::Power( (fPMVoltage / 112.5) ,7.04277 );
  99
 100 //  fCalibData = GetCalibData();
 101
 102 }
 103
 104 //____________________________________________________________________________
 105   AliVZERODigitizer::~AliVZERODigitizer()
 106 {
 107   // destructor
 108
 109   if (fDigits) {
 110     fDigits->Delete();
 111     delete fDigits;
 112     fDigits=0;
 113   }
 114 }
 115
 116 //_____________________________________________________________________________
 117 Bool_t AliVZERODigitizer::Init()
 118 {
 119   // Initialises the digitizer
 120
 121   // Initialises the Digit array
 122   fDigits = new TClonesArray ("AliVZEROdigit", 1000);
 123
 124   //  TGeoHMatrix *im = AliGeomManager::GetMatrix("VZERO/V0C");
 125   //  im->Print();
 126
 127   GetCollisionMode();
 128   return kTRUE;
 129 }
 130
 131 //____________________________________________________________________________
 132 void AliVZERODigitizer::Exec(Option_t* /*option*/)
 133 {
 134   // Creates digits from hits
 135
 136   Int_t       map[80];    // 48 values on V0C + 32 on V0A
 137 //  Int_t       pmNumber[80];
 138   Int_t       adc[64];    // 32 PMs on V0C + 32 PMs on V0A
 139   Float_t     time[80], time_ref[80], time2[64];
 140   Float_t     adc_gain[80];
 141   Float_t     adc_pedestal[64],adc_sigma[64];
 142   fNdigits     =    0;
 143   Float_t     pmGain_smeared[64];
 144   Float_t     cPM[80];
 145
 146   // Smearing of the PM tubes intrinsic characteristics
 147
 148   for(Int_t ii=0; ii<64; ii++) {
 149       pmGain_smeared[ii] = gRandom->Gaus(fPMGain, fPMGain/5); }
 150
 151   // Retrieval of ADC gain values and pedestal information from local CDB
 152   // I use only the first 64th values of the calibration array in CDB
 153   // as I have no beam burst structure - odd or even beam burst number
 154
 155   // Reminder : We have 16 scintillating cells mounted on 8 PMs
 156   // on Ring 3 and Ring 4 in V0C -  added to produce  ADC outputs
 157   // on these rings...
 158
 159   for(Int_t i=0; i<16; i++) {
 160         adc_gain[i] = fCalibData->GetGain(i);
 161         cPM[i]      = fPhotoCathodeEfficiency * pmGain_smeared[i];
 162   }
 163
 164   for(Int_t j=16; j<48; j=j+2) {
 165         Int_t i=(j+17)/2;
 166         adc_gain[j]   = fCalibData->GetGain(i);
 167         adc_gain[j+1] = fCalibData->GetGain(i);
 168         cPM[j]        = fPhotoCathodeEfficiency * pmGain_smeared[i];
 169         cPM[j+1]      = fPhotoCathodeEfficiency * pmGain_smeared[i];
 170   }
 171
 172   for(Int_t i=48; i<80; i++){
 173           adc_gain[i] = fCalibData->GetGain(i-16);
 174         cPM[i]      = fPhotoCathodeEfficiency * pmGain_smeared[i-16];
 175   };
 176
 177   for(Int_t  i=0; i<64; i++){
 178           adc_pedestal[i] = fCalibData->GetPedestal(i);
 179           adc_sigma[i]    = fCalibData->GetSigma(i);
 180   };
 181
 182 //  for(Int_t i=0; i<64; i++) { printf(" i = %d pedestal = %f sigma = %f \n\n",
 183 //                                       i, adc_pedestal[i], adc_sigma[i] );}
 184
 185   AliRunLoader* outRunLoader =  AliRunLoader::GetRunLoader(fManager->GetOutputFolderName());
 186   if (!outRunLoader) {
 187     Error("Exec", "Can not get output Run Loader");
 188     return;
 189   }
 190
 191   AliLoader* outLoader = outRunLoader->GetLoader("VZEROLoader");
 192
 193   if (!outLoader) {
 194     Error("Exec", "Can not get output VZERO Loader");
 195     return;
 196   }
 197
 198   const char* mode = "update";
 199   if(outRunLoader->GetEventNumber() == 0) mode = "recreate";
 200   outLoader->LoadDigits(mode);
 201
 202   if (!outLoader->TreeD()) outLoader->MakeTree("D");
 203   outLoader->MakeDigitsContainer();
 204   TTree* treeD  = outLoader->TreeD();
 205   Int_t bufsize = 16000;
 206   treeD->Branch("VZERODigit", &fDigits, bufsize);
 207
 208   for (Int_t iInput = 0; iInput < fManager->GetNinputs(); iInput++) {
 209      AliRunLoader* runLoader = AliRunLoader::GetRunLoader(fManager->GetInputFolderName(iInput));
 210      AliLoader* loader = runLoader->GetLoader("VZEROLoader");
 211      if (!loader) {
 212        Error("Exec", "Can not get VZERO Loader for input %d", iInput);
 213        continue;
 214          }
 215
 216      if (!runLoader->GetAliRun()) runLoader->LoadgAlice();
 217
 218      AliVZERO* vzero = (AliVZERO*) runLoader->GetAliRun()->GetDetector("VZERO");
 219      if (!vzero) {
 220        Error("Exec", "No VZERO detector for input %d", iInput);
 221        continue;
 222          }
 223
 224      loader->LoadHits();
 225      TTree* treeH = loader->TreeH();
 226      if (!treeH) {
 227        Error("Exec", "Cannot get TreeH for input %d", iInput);
 228        continue;
 229          }
 230
 231      for(Int_t i=0; i<80; i++) {map[i] = 0; time[i] = 0.0;}
 232
 233      TClonesArray* hits = vzero->Hits();
 234
 235 //  Now makes Digits from hits
 236
 237      Int_t nTracks = (Int_t) treeH->GetEntries();
 238      for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
 239          for (Int_t i=0; i<80; i++) {time_ref[i] = 999999.0;}
 240          vzero->ResetHits();
 241          treeH->GetEvent(iTrack);
 242          Int_t nHits = hits->GetEntriesFast();
 243          for (Int_t iHit = 0; iHit < nHits; iHit++) {
 244                          AliVZEROhit* hit = (AliVZEROhit *)hits->UncheckedAt(iHit);
 245                          Int_t nPhot = hit->Nphot();
 246                          Int_t cell  = hit->Cell();
 247                          map[cell] += nPhot;
 248                          Float_t dt_scintillator = gRandom->Gaus(0,0.7);
 249                          Float_t t = dt_scintillator + 1e9*hit->Tof();
 250                          if (t > 0.0) {
 251                                  if(t < time_ref[cell]) time_ref[cell] = t;
 252                                  time[cell] = TMath::Min(t,time_ref[cell]);
 253                          }
 254          }           // hit   loop
 255      }             // track loop
 256
 257      loader->UnloadHits();
 258
 259   }               // input loop
 260
 261 // Now builds the scintillator cell response (80 cells i.e. 80 responses)
 262
 263    for (Int_t i=0; i<80; i++) {
 264         Float_t q1 = Float_t ( map[i] )* cPM[i] * kQe;
 265         Float_t noise = gRandom->Gaus(10.5,3.22);
 266         Float_t pmResponse  =  q1/kC*TMath::Power(ktheta/kthau,1/(1-ktheta/kthau))
 267         + noise*1e-3;
 268         if(fCollisionMode >0) adc_gain[i] = adc_gain[i]/70.0; // reduce dynamics in Ion Collision Mode
 269         map[i] = Int_t( pmResponse * adc_gain[i]);
 270         Float_t MIP = 1.0/fCalibData->GetMIPperADC(GetPMNumber(i));
 271         if(fCollisionMode >0) MIP=2.0;
 272 //      printf("cell = %d,  ADC = %d, TDC = %f \n",i,map[i], time[i]*10.0 );
 273         if(map[i] > (int(( MIP/2 ) + 0.5)) )
 274                   {map[i] = Int_t(gRandom->Gaus(map[i], (int(( MIP/6 ) + 0.5)) ));}
 275    }
 276
 277 // Now transforms 80 cell responses into 64 photomultiplier responses
 278 // Also adds the ADC pedestals taken out of the calibration data base
 279
 280    for (Int_t j=0; j<16; j++){
 281         adc[j]  = static_cast<Int_t>(map [j] + gRandom->Gaus(adc_pedestal[j], adc_sigma[j]));
 282         time2[j]= time[j];}
 283
 284    for (Int_t j=48; j<80; j++){
 285         adc[j-16]  = static_cast<Int_t>(map [j]
 286                                         + gRandom->Gaus(adc_pedestal[j-16],adc_sigma[j-16]));
 287         time2[j-16]= time[j]; }
 288
 289    for (Int_t j=0; j<16; j++){
 290         adc[16+j] = static_cast<Int_t>(map [16+2*j]+ map [16+2*j+1]
 291                                        + gRandom->Gaus(adc_pedestal[16+j], adc_sigma[16+j]));
 292         Float_t min_time = TMath::Min(time [16+2*j],time [16+2*j+1]);
 293         time2[16+j] = min_time;
 294         if(min_time==0.0){time2[16+j]=TMath::Max(time[16+2*j],time[16+2*j+1]);}
 295    }
 296
 297
 298 // Now add digits to the digit Tree
 299
 300    for (Int_t i=0; i<64; i++) {
 301       if(adc[i] > 0) {
 302 //           printf(" Event, cell, adc, tof = %d %d %d %d\n",
 303 //                    outRunLoader->GetEventNumber(),i, adc[i], Int_t((time2[i]*10.0) +0.5));
 304 //           multiply by 10 to have 100 ps per channel :
 305
 306                   AddDigit(i, adc[i], Int_t((time2[i]*10.0) +0.5)) ;
 307           }
 308    }
 309   treeD->Fill();
 310   outLoader->WriteDigits("OVERWRITE");
 311   outLoader->UnloadDigits();
 312   ResetDigit();
 313 }
 314
 315 //____________________________________________________________________________
 316 void AliVZERODigitizer::AddDigit(Int_t PMnumber, Int_t adc, Int_t time)
 317  {
 318
 319 // Adds Digit
 320
 321   TClonesArray &ldigits = *fDigits;
 322   Bool_t integrator;
 323   if (((Int_t) gRandom->Uniform(2))<1) integrator = kFALSE;
 324   else integrator = kTRUE;
 325
 326   new(ldigits[fNdigits++]) AliVZEROdigit(PMnumber,adc,time,0,kFALSE,kFALSE,integrator);
 327
 328 }
 329 //____________________________________________________________________________
 330 void AliVZERODigitizer::ResetDigit()
 331 {
 332
 333 // Clears Digits
 334
 335   fNdigits = 0;
 336   if (fDigits) fDigits->Delete();
 337 }
 338
 339 //____________________________________________________________________________
 340 void AliVZERODigitizer::GetCollisionMode()
 341 {
 342 // Retrieves the collision mode from GRP data
 343
 344 // Initialization of the GRP entry
 345
 346    Int_t run = AliCDBManager::Instance()->GetRun();
 347
 348 //   printf("\n ++++++ Run Number retrieved as %d \n",run);
 349
 350   AliCDBEntry*  entry = AliCDBManager::Instance()->Get("GRP/GRP/Data",run);
 351   AliGRPObject* grpData = 0x0;
 352
 353   if(entry){
 354     TMap* m = dynamic_cast<TMap*>(entry->GetObject());  // old GRP entry
 355     if(m){
 356        m->Print();
 357        grpData = new AliGRPObject();
 358        grpData->ReadValuesFromMap(m);
 359     }
 360     else{
 361        grpData = dynamic_cast<AliGRPObject*>(entry->GetObject());  // new GRP entry
 362        entry->SetOwner(0);
 363     }
 364     AliCDBManager::Instance()->UnloadFromCache("GRP/GRP/Data");
 365   }
 366
 367   if(!grpData) AliError("No GRP entry found in OCDB!");
 368
 369 // Retrieval of collision mode
 370
 371   TString beamType = grpData->GetBeamType();
 372   if(beamType==AliGRPObject::GetInvalidString()){
 373      AliError("GRP/GRP/Data entry:  missing value for the beam type !");
 374      AliError("\t VZERO cannot retrieve beam type\n");
 375      return;
 376   }
 377
 378    if( (beamType.CompareTo("P-P") ==0)  || (beamType.CompareTo("p-p") ==0) ){
 379        fCollisionMode=0;
 380   }
 381    else if( (beamType.CompareTo("Pb-Pb") ==0)  || (beamType.CompareTo("A-A") ==0) ){
 382        fCollisionMode=1;
 383    }
 384
 385   fBeamEnergy = grpData->GetBeamEnergy();
 386   if(fBeamEnergy==AliGRPObject::GetInvalidFloat()) {
 387      AliError("GRP/GRP/Data entry:  missing value for the beam energy ! Using 0");
 388      fBeamEnergy = 0.;
 389   }
 390
 391 //     printf("\n ++++++ Beam type and collision mode retrieved as %s %d @ %1.3f GeV ++++++\n\n",beamType.Data(), fCollisionMode, fBeamEnergy);
 392
 393 }
 394
 395 //____________________________________________________________________________
 396 AliVZEROCalibData* AliVZERODigitizer::GetCalibData() const
 397
 398 {
 399   AliCDBManager *man = AliCDBManager::Instance();
 400
 401   AliCDBEntry *entry=0;
 402
 403   entry = man->Get("VZERO/Calib/Data");
 404
 405 //   if(!entry){
 406 //     AliWarning("Load of calibration data from default storage failed!");
 407 //     AliWarning("Calibration data will be loaded from local storage ($ALICE_ROOT)");
 408 //     Int_t runNumber = man->GetRun();
 409 //     entry = man->GetStorage("local://$ALICE_ROOT/OCDB")
 410 //       ->Get("VZERO/Calib/Data",runNumber);
 411 //
 412 //   }
 413
 414   // Retrieval of data in directory VZERO/Calib/Data:
 415
 416
 417   AliVZEROCalibData *calibdata = 0;
 418
 419   if (entry) calibdata = (AliVZEROCalibData*) entry->GetObject();
 420   if (!calibdata)  AliFatal("No calibration data from calibration database !");
 421
 422   return calibdata;
 423
 424 }
 425
 426 //____________________________________________________________________________
 427 Int_t AliVZERODigitizer::GetPMNumber(Int_t cell) const
 428
 429 {
 430
 431   Int_t pmNumber[80] = { 0,  1,  2,  3,  4,  5,  6,  7,
 432                           8,  9, 10, 11, 12, 13, 14, 15,
 433                          16, 16, 17, 17, 18, 18, 19, 19, 20, 20, 21, 21, 22, 22, 23, 23,
 434                          24, 24, 25, 25, 26, 26, 27, 27, 28, 28, 29, 29, 30, 30, 31, 31,
 435                          32, 33, 34, 35, 36, 37, 38, 39,
 436                          40, 41, 42, 43, 44, 45, 46, 47,
 437                          48, 49, 50, 51, 52, 53, 54, 55,
 438                          56, 57, 58, 59, 60, 61, 62, 63 };
 439
 440   return pmNumber[cell];
 441 }
 442
 443