[u/mrichter/AliRoot.git] / VZERO / AliVZERODigitizer.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/
 
/* $Id$ */

///_________________________________________________________________________
///
/// This class constructs Digits out of Hits
///
///

// --- Standard library ---

// --- ROOT system ---
#include <TMath.h>
#include <TTree.h>
#include <TMap.h>
#include <TGeoManager.h>
#include <TGeoPhysicalNode.h>
#include <AliGeomManager.h>
#include <TRandom.h>
#include <TF1.h>

// --- AliRoot header files ---
#include "AliVZEROConst.h"
#include "AliRun.h"
#include "AliVZERO.h"
#include "AliVZEROhit.h"
#include "AliRunLoader.h"
#include "AliLoader.h"
#include "AliGRPObject.h"
#include "AliRunDigitizer.h"
#include "AliCDBManager.h"
#include "AliCDBStorage.h"
#include "AliCDBEntry.h"
#include "AliVZEROCalibData.h"

#include "AliVZEROdigit.h"
#include "AliVZERODigitizer.h"

ClassImp(AliVZERODigitizer)

 AliVZERODigitizer::AliVZERODigitizer()
                   :AliDigitizer(),
                    fCalibData(GetCalibData()),
                    fPhotoCathodeEfficiency(0.18),
                    fPMVoltage(768.0),
                    fPMGain(TMath::Power((fPMVoltage / 112.5) ,7.04277)),
                    fNdigits(0),
                    fDigits(0),
                    fCollisionMode(0),
                    fBeamEnergy(0.),
                    fSignalShape(NULL)
   
{
  // default constructor

//    fNdigits = 0;
//    fDigits  = 0;
//   
//    fPhotoCathodeEfficiency =   0.18;
//    fPMVoltage              =  768.0;
//    fPMGain = TMath::Power((fPMVoltage / 112.5) ,7.04277); 
   
//   fCalibData = GetCalibData();

  fSignalShape = new TF1("VZEROSignalShape",this,&AliVZERODigitizer::SignalShape,0,200,6,"AliVZERODigitizer","SignalShape");
  fSignalShape->SetParameters(0,1.57345e1,-4.25603e-1,2.9,6.40982,3.69339e-01);
}

//____________________________________________________________________________ 
  AliVZERODigitizer::AliVZERODigitizer(AliRunDigitizer* manager)
                    :AliDigitizer(manager),
		     fCalibData(GetCalibData()),
                     fPhotoCathodeEfficiency(0.18),
                     fPMVoltage(768.0),
                     fPMGain(TMath::Power((fPMVoltage / 112.5) ,7.04277)),
		     fNdigits(0),
                     fDigits(0),
		     fCollisionMode(0),
                     fBeamEnergy(0.),
		     fSignalShape(NULL)
		                        
{
  // constructor
  
//   fNdigits = 0;
//   fDigits  = 0;
//   
//   fPhotoCathodeEfficiency =   0.18;
//   fPMVoltage              =  768.0;
//   fPMGain = TMath::Power( (fPMVoltage / 112.5) ,7.04277 );
  
//  fCalibData = GetCalibData();
 
  fSignalShape = new TF1("VZEROSignalShape",this,&AliVZERODigitizer::SignalShape,0,200,6,"AliVZERODigitizer","SignalShape");
  fSignalShape->SetParameters(0,1.57345e1,-4.25603e-1,2.9,6.40982,3.69339e-01);
  
}
           
//____________________________________________________________________________ 
  AliVZERODigitizer::~AliVZERODigitizer()
{
  // destructor
  
  if (fDigits) {
    fDigits->Delete();
    delete fDigits;
    fDigits=0; 
  }

  if (fSignalShape) {
    delete fSignalShape;
    fSignalShape = NULL;
  }
}

//_____________________________________________________________________________
Bool_t AliVZERODigitizer::Init()
{
  // Initialises the digitizer

  // Initialises the Digit array
  fDigits = new TClonesArray ("AliVZEROdigit", 1000);
  
  //  TGeoHMatrix *im = AliGeomManager::GetMatrix("VZERO/V0C");
  //  im->Print();

  GetCollisionMode();
  return kTRUE;
}

//____________________________________________________________________________
void AliVZERODigitizer::Exec(Option_t* /*option*/) 
{   
  // Creates digits from hits
     
  Float_t     map[80];    // 48 values on V0C + 32 on V0A
//  Int_t       pmNumber[80];
  Float_t     adc[64];    // 32 PMs on V0C + 32 PMs on V0A
  Float_t     time[80], time_ref[80], time2[64];
  Float_t     adc_gain[80]; 
  Float_t     adc_pedestal[64],adc_sigma[64];    
  fNdigits     =    0;  
  Float_t     pmGain_smeared[64];  
  Float_t     cPM[80];
  
  // Smearing of the PM tubes intrinsic characteristics
  
  for(Int_t ii=0; ii<64; ii++) {
      pmGain_smeared[ii] = gRandom->Gaus(fPMGain, fPMGain/5); }
              
  // Retrieval of ADC gain values and pedestal information from local CDB 
  // I use only the first 64th values of the calibration array in CDB 
  // as I have no beam burst structure - odd or even beam burst number
  
  // Reminder : We have 16 scintillating cells mounted on 8 PMs 
  // on Ring 3 and Ring 4 in V0C -  added to produce  ADC outputs 
  // on these rings... 
   
  for(Int_t i=0; i<16; i++) { 
	adc_gain[i] = fCalibData->GetGain(i); 
	cPM[i]      = fPhotoCathodeEfficiency * pmGain_smeared[i];
  }
  
  for(Int_t j=16; j<48; j=j+2) { 
	Int_t i=(j+17)/2;
	adc_gain[j]   = fCalibData->GetGain(i);	    
	adc_gain[j+1] = fCalibData->GetGain(i); 
	cPM[j]        = fPhotoCathodeEfficiency * pmGain_smeared[i];   
	cPM[j+1]      = fPhotoCathodeEfficiency * pmGain_smeared[i]; 
  }
	    
  for(Int_t i=48; i<80; i++){ 
	adc_gain[i] = fCalibData->GetGain(i-16); 
	cPM[i]      = fPhotoCathodeEfficiency * pmGain_smeared[i-16];
  };
  
  for(Int_t  i=0; i<64; i++){ 
	  adc_pedestal[i] = fCalibData->GetPedestal(i);
	  adc_sigma[i]    = fCalibData->GetSigma(i); 
  }; 
                                
//  for(Int_t i=0; i<64; i++) { printf(" i = %d pedestal = %f sigma = %f \n\n", 
//                                       i, adc_pedestal[i], adc_sigma[i] );} 
            
  AliRunLoader* outRunLoader =  AliRunLoader::GetRunLoader(fManager->GetOutputFolderName());    
  if (!outRunLoader) {
    Error("Exec", "Can not get output Run Loader");
    return;
  }
    
  AliLoader* outLoader = outRunLoader->GetLoader("VZEROLoader");

  if (!outLoader) {
    Error("Exec", "Can not get output VZERO Loader");
    return;
  }

  const char* mode = "update";
  if(outRunLoader->GetEventNumber() == 0) mode = "recreate";
  outLoader->LoadDigits(mode);

  if (!outLoader->TreeD()) outLoader->MakeTree("D");
  outLoader->MakeDigitsContainer();
  TTree* treeD  = outLoader->TreeD();
  Int_t bufsize = 16000;
  treeD->Branch("VZERODigit", &fDigits, bufsize); 
  
  for (Int_t iInput = 0; iInput < fManager->GetNinputs(); iInput++) {
     AliRunLoader* runLoader = AliRunLoader::GetRunLoader(fManager->GetInputFolderName(iInput));
     AliLoader* loader = runLoader->GetLoader("VZEROLoader");
     if (!loader) {
       Error("Exec", "Can not get VZERO Loader for input %d", iInput);
       continue;
	 }
      
     if (!runLoader->GetAliRun()) runLoader->LoadgAlice();

     AliVZERO* vzero = (AliVZERO*) runLoader->GetAliRun()->GetDetector("VZERO");
     if (!vzero) {
       Error("Exec", "No VZERO detector for input %d", iInput);
       continue;
	 }
      
     loader->LoadHits();
     TTree* treeH = loader->TreeH();
     if (!treeH) {
       Error("Exec", "Cannot get TreeH for input %d", iInput);
       continue; 
	 }
       
     for(Int_t i=0; i<80; i++) {map[i] = 0; time[i] = 0.0;}
     
     TClonesArray* hits = vzero->Hits();
             
//  Now makes Digits from hits
     
     Int_t nTracks = (Int_t) treeH->GetEntries();
     for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
         for (Int_t i=0; i<80; i++) {time_ref[i] = 999999.0;}   
         vzero->ResetHits();
         treeH->GetEvent(iTrack);
         Int_t nHits = hits->GetEntriesFast();
         for (Int_t iHit = 0; iHit < nHits; iHit++) {
			 AliVZEROhit* hit = (AliVZEROhit *)hits->UncheckedAt(iHit);
			 Int_t nPhot = hit->Nphot();
			 Int_t cell  = hit->Cell();                          
			 map[cell] += Float_t(nPhot);
			 Float_t dt_scintillator = gRandom->Gaus(0,0.7);
			 Float_t t = dt_scintillator + 1e9*hit->Tof();
			 if (t > 0.0) {
				 if(t < time_ref[cell]) time_ref[cell] = t;
				 time[cell] = TMath::Min(t,time_ref[cell]); 
			 }
         }           // hit   loop      
     }             // track loop

     loader->UnloadHits();

  }               // input loop

// Now builds the scintillator cell response (80 cells i.e. 80 responses)
         
   for (Int_t i=0; i<80; i++) {    
        Float_t q1 = Float_t ( map[i] )* cPM[i] * kQe;
        Float_t noise = gRandom->Gaus(10.5,3.22);
        Float_t pmResponse  =  q1/kC*TMath::Power(ktheta/kthau,1/(1-ktheta/kthau)) 
        + noise*1e-3; 	
	if(fCollisionMode >0) adc_gain[i] = adc_gain[i]/70.0; // reduce dynamics in Ion Collision Mode
        map[i] =  pmResponse * adc_gain[i];
        Float_t MIP = 1.0/fCalibData->GetMIPperADC(GetPMNumber(i));
	if(fCollisionMode >0) MIP=2.0;
//	printf("cell = %d,  ADC = %d, TDC = %f \n",i,map[i], time[i]*10.0 );
        if(map[i] > (MIP/2.) )
	          {map[i] = gRandom->Gaus(map[i], (MIP/6.) );}
   }
      
// Now transforms 80 cell responses into 64 photomultiplier responses
// Also adds the ADC pedestals taken out of the calibration data base
	
   for (Int_t j=0; j<16; j++){
        adc[j]  = map [j] + gRandom->Gaus(adc_pedestal[j], adc_sigma[j]);
	time2[j]= time[j];}
	
   for (Int_t j=48; j<80; j++){
        adc[j-16]  = map [j] + gRandom->Gaus(adc_pedestal[j-16],adc_sigma[j-16]);
	time2[j-16]= time[j]; }
	
   for (Int_t j=0; j<16; j++){
        adc[16+j] = map [16+2*j]+ map [16+2*j+1] + gRandom->Gaus(adc_pedestal[16+j], adc_sigma[16+j]);
	Float_t min_time = TMath::Min(time [16+2*j],time [16+2*j+1]);
	time2[16+j] = min_time;
	if(min_time==0.0){time2[16+j]=TMath::Max(time[16+2*j],time[16+2*j+1]);}
   }
   	

// Now add digits to the digit Tree 
        
   for (Int_t i=0; i<64; i++) {      
      if(adc[i] > 0) {
//           printf(" Event, cell, adc, tof = %d %d %d %d\n", 
//                    outRunLoader->GetEventNumber(),i, adc[i], Int_t((time2[i]*10.0) +0.5));
//           multiply by 10 to have 100 ps per channel :
		  
		  AddDigit(i, adc[i], time2[i] ) ;
	  }      
   }
  treeD->Fill();
  outLoader->WriteDigits("OVERWRITE");  
  outLoader->UnloadDigits();     
  ResetDigit();
}

//____________________________________________________________________________
void AliVZERODigitizer::AddDigit(Int_t PMnumber, Float_t adc, Float_t time) 
 { 
 
// Adds Digit 
 
  TClonesArray &ldigits = *fDigits;  
  Bool_t integrator;
  if (((Int_t) gRandom->Uniform(2))<1) integrator = kFALSE;
  else integrator = kTRUE;
	 
  new(ldigits[fNdigits++]) AliVZEROdigit(PMnumber,adc,time,0,kFALSE,kFALSE,integrator);
	 
}
//____________________________________________________________________________
void AliVZERODigitizer::ResetDigit()
{

// Clears Digits

  fNdigits = 0;
  if (fDigits) fDigits->Delete();
}

//____________________________________________________________________________
void AliVZERODigitizer::GetCollisionMode()
{
// Retrieves the collision mode from GRP data

// Initialization of the GRP entry 

   Int_t run = AliCDBManager::Instance()->GetRun();
  
//   printf("\n ++++++ Run Number retrieved as %d \n",run); 
 
  AliCDBEntry*  entry = AliCDBManager::Instance()->Get("GRP/GRP/Data",run);
  AliGRPObject* grpData = 0x0;
   
  if(entry){
    TMap* m = dynamic_cast<TMap*>(entry->GetObject());  // old GRP entry
    if(m){
       m->Print();
       grpData = new AliGRPObject();
       grpData->ReadValuesFromMap(m);
    }
    else{
       grpData = dynamic_cast<AliGRPObject*>(entry->GetObject());  // new GRP entry
       entry->SetOwner(0);
    }
    AliCDBManager::Instance()->UnloadFromCache("GRP/GRP/Data");
  }

  if(!grpData) AliError("No GRP entry found in OCDB!");

// Retrieval of collision mode

  TString beamType = grpData->GetBeamType();
  if(beamType==AliGRPObject::GetInvalidString()){
     AliError("GRP/GRP/Data entry:  missing value for the beam type !");
     AliError("\t VZERO cannot retrieve beam type\n");
     return;
  }

   if( (beamType.CompareTo("P-P") ==0)  || (beamType.CompareTo("p-p") ==0) ){
       fCollisionMode=0;
  }
   else if( (beamType.CompareTo("Pb-Pb") ==0)  || (beamType.CompareTo("A-A") ==0) ){
       fCollisionMode=1;
   }
    
  fBeamEnergy = grpData->GetBeamEnergy();
  if(fBeamEnergy==AliGRPObject::GetInvalidFloat()) {
     AliError("GRP/GRP/Data entry:  missing value for the beam energy ! Using 0");
     fBeamEnergy = 0.;
  }
  
//     printf("\n ++++++ Beam type and collision mode retrieved as %s %d @ %1.3f GeV ++++++\n\n",beamType.Data(), fCollisionMode, fBeamEnergy);

}

//____________________________________________________________________________
AliVZEROCalibData* AliVZERODigitizer::GetCalibData() const

{
  AliCDBManager *man = AliCDBManager::Instance();

  AliCDBEntry *entry=0;

  entry = man->Get("VZERO/Calib/Data");

//   if(!entry){
//     AliWarning("Load of calibration data from default storage failed!");
//     AliWarning("Calibration data will be loaded from local storage ($ALICE_ROOT)");
//     Int_t runNumber = man->GetRun();
//     entry = man->GetStorage("local://$ALICE_ROOT/OCDB")
//       ->Get("VZERO/Calib/Data",runNumber);
// 	
//   }

  // Retrieval of data in directory VZERO/Calib/Data:


  AliVZEROCalibData *calibdata = 0;

  if (entry) calibdata = (AliVZEROCalibData*) entry->GetObject();
  if (!calibdata)  AliFatal("No calibration data from calibration database !");

  return calibdata;

}

//____________________________________________________________________________
Int_t AliVZERODigitizer::GetPMNumber(Int_t cell) const

{
   
  Int_t pmNumber[80] = { 0,  1,  2,  3,  4,  5,  6,  7,
                          8,  9, 10, 11, 12, 13, 14, 15, 
			 16, 16, 17, 17, 18, 18, 19, 19, 20, 20, 21, 21, 22, 22, 23, 23, 
			 24, 24, 25, 25, 26, 26, 27, 27, 28, 28, 29, 29, 30, 30, 31, 31,
			 32, 33, 34, 35, 36, 37, 38, 39,
		         40, 41, 42, 43, 44, 45, 46, 47, 
			 48, 49, 50, 51, 52, 53, 54, 55,
			 56, 57, 58, 59, 60, 61, 62, 63 };
			      
  return pmNumber[cell];	
}

double AliVZERODigitizer::SignalShape(double *x, double *par)
{
  Double_t xx = x[0];
  if (xx <= par[0]) return 0;
  Double_t a = 1./TMath::Power((xx-par[0])/par[1],1./par[2]);
  if (xx <= par[3]) return a;
  Double_t b = 1./TMath::Power((xx-par[3])/par[4],1./par[5]);
  //  Double_t f = TMath::Min(a,b);
  Double_t f = a*b/(a+b);
  AliDebug(100,Form("x=%f func=%f",xx,f));
  return f;
}
Commit	Line	Data
9e04c3b6	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
090026bf	16	/* $Id$ */
090026bf	17
b0d2c2d3	18	///_________________________________________________________________________
	19	///
	20	/// This class constructs Digits out of Hits
	21	///
	22	///
9e04c3b6	23
9e04c3b6	24	// --- Standard library ---
9e04c3b6	25
9e04c3b6	26	// --- ROOT system ---
090026bf	27	#include <TMath.h>
9e04c3b6	28	#include <TTree.h>
fe0adf2a	29	#include <TMap.h>
8bd3e27a	30	#include <TGeoManager.h>
	31	#include <TGeoPhysicalNode.h>
	32	#include <AliGeomManager.h>
e939a978	33	#include <TRandom.h>
9672d66e	34	#include <TF1.h>
9e04c3b6	35
	36	// --- AliRoot header files ---
	37	#include "AliVZEROConst.h"
	38	#include "AliRun.h"
	39	#include "AliVZERO.h"
	40	#include "AliVZEROhit.h"
9e04c3b6	41	#include "AliRunLoader.h"
9e04c3b6	42	#include "AliLoader.h"
fe0adf2a	43	#include "AliGRPObject.h"
9e04c3b6	44	#include "AliRunDigitizer.h"
ce7090f5	45	#include "AliCDBManager.h"
	46	#include "AliCDBStorage.h"
	47	#include "AliCDBEntry.h"
	48	#include "AliVZEROCalibData.h"
	49
9e04c3b6	50	#include "AliVZEROdigit.h"
b0d2c2d3	51	#include "AliVZERODigitizer.h"
9e04c3b6	52
	53	ClassImp(AliVZERODigitizer)
	54
	55	AliVZERODigitizer::AliVZERODigitizer()
fe0adf2a	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),
9672d66e	64	fBeamEnergy(0.),
9672d66e	65	fSignalShape(NULL)
0b2bea8b	66
9e04c3b6	67	{
b0d2c2d3	68	// default constructor
b0d2c2d3	69
0b2bea8b	70	// fNdigits = 0;
	71	// fDigits = 0;
	72	//
	73	// fPhotoCathodeEfficiency = 0.18;
	74	// fPMVoltage = 768.0;
	75	// fPMGain = TMath::Power((fPMVoltage / 112.5) ,7.04277);
ce7090f5	76
0b2bea8b	77	// fCalibData = GetCalibData();
9672d66e	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);
9e04c3b6	81	}
	82
	83	//____________________________________________________________________________
	84	AliVZERODigitizer::AliVZERODigitizer(AliRunDigitizer* manager)
0b2bea8b	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),
fe0adf2a	91	fDigits(0),
fe0adf2a	92	fCollisionMode(0),
9672d66e	93	fBeamEnergy(0.),
9672d66e	94	fSignalShape(NULL)
0b2bea8b	95
9e04c3b6	96	{
	97	// constructor
	98
0b2bea8b	99	// fNdigits = 0;
	100	// fDigits = 0;
	101	//
	102	// fPhotoCathodeEfficiency = 0.18;
	103	// fPMVoltage = 768.0;
	104	// fPMGain = TMath::Power( (fPMVoltage / 112.5) ,7.04277 );
ce7090f5	105
0b2bea8b	106	// fCalibData = GetCalibData();
9672d66e	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);
ce7090f5	110
9e04c3b6	111	}
	112
	113	//____________________________________________________________________________
	114	AliVZERODigitizer::~AliVZERODigitizer()
	115	{
	116	// destructor
	117
	118	if (fDigits) {
	119	fDigits->Delete();
	120	delete fDigits;
b0d2c2d3	121	fDigits=0;
b0d2c2d3	122	}
9672d66e	123
	124	if (fSignalShape) {
	125	delete fSignalShape;
	126	fSignalShape = NULL;
	127	}
9e04c3b6	128	}
9e04c3b6	129
b0d2c2d3	130	//_____________________________________________________________________________
b0d2c2d3	131	Bool_t AliVZERODigitizer::Init()
9e04c3b6	132	{
b0d2c2d3	133	// Initialises the digitizer
9e04c3b6	134
b0d2c2d3	135	// Initialises the Digit array
9e04c3b6	136	fDigits = new TClonesArray ("AliVZEROdigit", 1000);
8bd3e27a	137
c61a7285	138	// TGeoHMatrix *im = AliGeomManager::GetMatrix("VZERO/V0C");
c61a7285	139	// im->Print();
fe0adf2a	140
fe0adf2a	141	GetCollisionMode();
b0d2c2d3	142	return kTRUE;
9e04c3b6	143	}
	144
	145	//____________________________________________________________________________
b0d2c2d3	146	void AliVZERODigitizer::Exec(Option_t* /option/)
ce7090f5	147	{
b0d2c2d3	148	// Creates digits from hits
ce7090f5	149
db0db003	150	Float_t map[80]; // 48 values on V0C + 32 on V0A
f25f4990	151	// Int_t pmNumber[80];
db0db003	152	Float_t adc[64]; // 32 PMs on V0C + 32 PMs on V0A
8bd3e27a	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];
ce7090f5	156	fNdigits = 0;
8adc9b44	157	Float_t pmGain_smeared[64];
8bd3e27a	158	Float_t cPM[80];
	159
	160	// Smearing of the PM tubes intrinsic characteristics
	161
	162	for(Int_t ii=0; ii<64; ii++) {
8adc9b44	163	pmGain_smeared[ii] = gRandom->Gaus(fPMGain, fPMGain/5); }
8bd3e27a	164
8bd3e27a	165	// Retrieval of ADC gain values and pedestal information from local CDB
ce7090f5	166	// I use only the first 64th values of the calibration array in CDB
ce7090f5	167	// as I have no beam burst structure - odd or even beam burst number
b0d2c2d3	168
ce7090f5	169	// Reminder : We have 16 scintillating cells mounted on 8 PMs
20277079	170	// on Ring 3 and Ring 4 in V0C - added to produce ADC outputs
20277079	171	// on these rings...
ce7090f5	172
8bd3e27a	173	for(Int_t i=0; i<16; i++) {
fad64858	174	adc_gain[i] = fCalibData->GetGain(i);
	175	cPM[i] = fPhotoCathodeEfficiency * pmGain_smeared[i];
	176	}
ce7090f5	177
ce7090f5	178	for(Int_t j=16; j<48; j=j+2) {
fad64858	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	}
8bd3e27a	185
8bd3e27a	186	for(Int_t i=48; i<80; i++){
db0db003	187	adc_gain[i] = fCalibData->GetGain(i-16);
fad64858	188	cPM[i] = fPhotoCathodeEfficiency * pmGain_smeared[i-16];
fad64858	189	};
ce7090f5	190
fad64858	191	for(Int_t i=0; i<64; i++){
	192	adc_pedestal[i] = fCalibData->GetPedestal(i);
	193	adc_sigma[i] = fCalibData->GetSigma(i);
	194	};
8bd3e27a	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] );}
ce7090f5	198
fad64858	199	AliRunLoader* outRunLoader = AliRunLoader::GetRunLoader(fManager->GetOutputFolderName());
b0d2c2d3	200	if (!outRunLoader) {
b0d2c2d3	201	Error("Exec", "Can not get output Run Loader");
fad64858	202	return;
fad64858	203	}
ce7090f5	204
b0d2c2d3	205	AliLoader* outLoader = outRunLoader->GetLoader("VZEROLoader");
fad64858	206
b0d2c2d3	207	if (!outLoader) {
b0d2c2d3	208	Error("Exec", "Can not get output VZERO Loader");
fad64858	209	return;
fad64858	210	}
b0d2c2d3	211
e539620f	212	const char* mode = "update";
	213	if(outRunLoader->GetEventNumber() == 0) mode = "recreate";
	214	outLoader->LoadDigits(mode);
fad64858	215
b0d2c2d3	216	if (!outLoader->TreeD()) outLoader->MakeTree("D");
b0d2c2d3	217	outLoader->MakeDigitsContainer();
ce7090f5	218	TTree* treeD = outLoader->TreeD();
b0d2c2d3	219	Int_t bufsize = 16000;
b0d2c2d3	220	treeD->Branch("VZERODigit", &fDigits, bufsize);
c299dbe4	221
b0d2c2d3	222	for (Int_t iInput = 0; iInput < fManager->GetNinputs(); iInput++) {
fad64858	223	AliRunLoader* runLoader = AliRunLoader::GetRunLoader(fManager->GetInputFolderName(iInput));
8bd3e27a	224	AliLoader* loader = runLoader->GetLoader("VZEROLoader");
	225	if (!loader) {
	226	Error("Exec", "Can not get VZERO Loader for input %d", iInput);
fad64858	227	continue;
fad64858	228	}
ce7090f5	229
8bd3e27a	230	if (!runLoader->GetAliRun()) runLoader->LoadgAlice();
b0d2c2d3	231
8bd3e27a	232	AliVZERO* vzero = (AliVZERO*) runLoader->GetAliRun()->GetDetector("VZERO");
	233	if (!vzero) {
	234	Error("Exec", "No VZERO detector for input %d", iInput);
fad64858	235	continue;
fad64858	236	}
9e04c3b6	237
8bd3e27a	238	loader->LoadHits();
	239	TTree* treeH = loader->TreeH();
	240	if (!treeH) {
	241	Error("Exec", "Cannot get TreeH for input %d", iInput);
fad64858	242	continue;
fad64858	243	}
c299dbe4	244
8bd3e27a	245	for(Int_t i=0; i<80; i++) {map[i] = 0; time[i] = 0.0;}
fe0adf2a	246
8bd3e27a	247	TClonesArray* hits = vzero->Hits();
9e04c3b6	248
b0d2c2d3	249	// Now makes Digits from hits
fe0adf2a	250
8bd3e27a	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++) {
fad64858	258	AliVZEROhit* hit = (AliVZEROhit *)hits->UncheckedAt(iHit);
	259	Int_t nPhot = hit->Nphot();
	260	Int_t cell = hit->Cell();
db0db003	261	map[cell] += Float_t(nPhot);
fad64858	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	}
8bd3e27a	268	} // hit loop
8bd3e27a	269	} // track loop
b0d2c2d3	270
8bd3e27a	271	loader->UnloadHits();
b0d2c2d3	272
b0d2c2d3	273	} // input loop
20277079	274
20277079	275	// Now builds the scintillator cell response (80 cells i.e. 80 responses)
841137ce	276
20277079	277	for (Int_t i=0; i<80; i++) {
8adc9b44	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))
8bd3e27a	281	+ noise*1e-3;
fe0adf2a	282	if(fCollisionMode >0) adc_gain[i] = adc_gain[i]/70.0; // reduce dynamics in Ion Collision Mode
db0db003	283	map[i] = pmResponse * adc_gain[i];
fe0adf2a	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 );
db0db003	287	if(map[i] > (MIP/2.) )
db0db003	288	{map[i] = gRandom->Gaus(map[i], (MIP/6.) );}
8adc9b44	289	}
c299dbe4	290
8bd3e27a	291	// Now transforms 80 cell responses into 64 photomultiplier responses
8bd3e27a	292	// Also adds the ADC pedestals taken out of the calibration data base
20277079	293
8bd3e27a	294	for (Int_t j=0; j<16; j++){
db0db003	295	adc[j] = map [j] + gRandom->Gaus(adc_pedestal[j], adc_sigma[j]);
20277079	296	time2[j]= time[j];}
	297
	298	for (Int_t j=48; j<80; j++){
db0db003	299	adc[j-16] = map [j] + gRandom->Gaus(adc_pedestal[j-16],adc_sigma[j-16]);
8bd3e27a	300	time2[j-16]= time[j]; }
20277079	301
20277079	302	for (Int_t j=0; j<16; j++){
db0db003	303	adc[16+j] = map [16+2j]+ map [16+2j+1] + gRandom->Gaus(adc_pedestal[16+j], adc_sigma[16+j]);
c299dbe4	304	Float_t min_time = TMath::Min(time [16+2j],time [16+2j+1]);
c299dbe4	305	time2[16+j] = min_time;
8bd3e27a	306	if(min_time==0.0){time2[16+j]=TMath::Max(time[16+2j],time[16+2j+1]);}
c299dbe4	307	}
	308
	309
20277079	310	// Now add digits to the digit Tree
	311
	312	for (Int_t i=0; i<64; i++) {
	313	if(adc[i] > 0) {
fad64858	314	// printf(" Event, cell, adc, tof = %d %d %d %d\n",
fad64858	315	// outRunLoader->GetEventNumber(),i, adc[i], Int_t((time2[i]*10.0) +0.5));
8bd3e27a	316	// multiply by 10 to have 100 ps per channel :
fad64858	317
4682ffee	318	AddDigit(i, adc[i], time2[i] ) ;
fad64858	319	}
20277079	320	}
b0d2c2d3	321	treeD->Fill();
	322	outLoader->WriteDigits("OVERWRITE");
	323	outLoader->UnloadDigits();
	324	ResetDigit();
9e04c3b6	325	}
	326
	327	//____________________________________________________________________________
db0db003	328	void AliVZERODigitizer::AddDigit(Int_t PMnumber, Float_t adc, Float_t time)
9e04c3b6	329	{
	330
	331	// Adds Digit
	332
	333	TClonesArray &ldigits = *fDigits;
fad64858	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
9e04c3b6	340	}
	341	//____________________________________________________________________________
	342	void AliVZERODigitizer::ResetDigit()
	343	{
fe0adf2a	344
9e04c3b6	345	// Clears Digits
fe0adf2a	346
9e04c3b6	347	fNdigits = 0;
b0d2c2d3	348	if (fDigits) fDigits->Delete();
9e04c3b6	349	}
ce7090f5	350
fe0adf2a	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
ce7090f5	407	//____________________________________________________________________________
	408	AliVZEROCalibData* AliVZERODigitizer::GetCalibData() const
	409
	410	{
c0b82b5a	411	AliCDBManager *man = AliCDBManager::Instance();
ce7090f5	412
c0b82b5a	413	AliCDBEntry *entry=0;
ce7090f5	414
c0b82b5a	415	entry = man->Get("VZERO/Calib/Data");
c0b82b5a	416
df791951	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();
162637e4	421	// entry = man->GetStorage("local://$ALICE_ROOT/OCDB")
df791951	422	// ->Get("VZERO/Calib/Data",runNumber);
	423	//
	424	// }
c0b82b5a	425
c0b82b5a	426	// Retrieval of data in directory VZERO/Calib/Data:
ce7090f5	427
ce7090f5	428
c0b82b5a	429	AliVZEROCalibData *calibdata = 0;
ce7090f5	430
c0b82b5a	431	if (entry) calibdata = (AliVZEROCalibData*) entry->GetObject();
df791951	432	if (!calibdata) AliFatal("No calibration data from calibration database !");
ce7090f5	433
c0b82b5a	434	return calibdata;
ce7090f5	435
	436	}
	437
8adc9b44	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
9672d66e	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	}