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