[u/mrichter/AliRoot.git] / EMCAL / AliEMCAL.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$ */
//_________________________________________________________________________
// Base Class for EMCAL description:
// This class contains material definitions    
// for the EMCAL - It does not place the detector in Alice
//*-- Author: Yves Schutz (SUBATECH) 
//
//*-- Additional Contributions: Sahal Yacoob (LBNL/UCT)
//                            : Alexei Pavlinov (WSU) 
//
//////////////////////////////////////////////////////////////////////////////

// --- ROOT system ---
class TFile;
#include <TFolder.h> 
#include <TGeoGlobalMagField.h>
#include <TGraph.h> 
#include <TH1F.h> 
#include <TRandom.h> 
#include <TTree.h>
#include <TVirtualMC.h> 

// --- Standard library ---

// --- AliRoot header files ---
#include "AliMagF.h"
#include "AliLog.h"
#include "AliEMCAL.h"
#include "AliRun.h"
#include "AliRunLoader.h"
#include "AliCDBManager.h"
#include "AliEMCALLoader.h"
#include "AliEMCALSDigitizer.h"
#include "AliEMCALDigitizer.h"
#include "AliEMCALDigit.h"
#include "AliEMCALRawUtils.h"
#include "AliCDBManager.h"
#include "AliCDBEntry.h"
#include "AliEMCALRawUtils.h"
#include "AliRawReader.h"
#include "AliEMCALTriggerData.h"
#include "AliEMCALRecParam.h"
#include "AliRawEventHeaderBase.h"

ClassImp(AliEMCAL)

//for embedding
AliEMCALRawUtils*        AliEMCAL::fgRawUtils   = 0;   // EMCAL raw utilities class

//____________________________________________________________________________
AliEMCAL::AliEMCAL()
  : AliDetector(),
    fBirkC0(0),
    fBirkC1(0.),
    fBirkC2(0.),
    fGeometry(0), 
    fCheckRunNumberAndGeoVersion(kTRUE),
    fTriggerData(0x0)
{
  // Default ctor 
  fName = "EMCAL" ;
  InitConstants();
  
  // Should call  AliEMCALGeometry::GetInstance(EMCAL->GetTitle(),"") for getting EMCAL geometry
}

//____________________________________________________________________________
AliEMCAL::AliEMCAL(const char* name, const char* title)
  : AliDetector(name,title),
    fBirkC0(0),
    fBirkC1(0.),
    fBirkC2(0.),
    fGeometry(0), 
    fCheckRunNumberAndGeoVersion(kTRUE),
    fTriggerData(0x0)
{
  //   ctor : title is used to identify the layout
  InitConstants();
  
}

//____________________________________________________________________________
AliEMCAL::~AliEMCAL()
{
  //dtor
  delete fgRawUtils;
  delete fTriggerData;
}

//____________________________________________________________________________
void AliEMCAL::InitConstants()
{
  //initialize EMCAL values
  fBirkC0 = 1;
  fBirkC1 = 0.013/1.032;
  fBirkC2 = 9.6e-6/(1.032 * 1.032);
}

//Not needed, modify $ALICE_ROOT/data/galice.cuts instead.
//Load the modified one in the configuration file with SetTransPar
// //____________________________________________________________________________
// void AliEMCAL::DefineMediumParameters()
// {
//   //
//   // EMCAL cuts (Geant3) 
//   // 
//   Int_t * idtmed = fIdtmed->GetArray() - 1599 ; 
// // --- Set decent energy thresholds for gamma and electron tracking

//   // Tracking threshold for photons and electrons in Lead 
//   Float_t cutgam=10.e-5; // 100 kev;
//   Float_t cutele=10.e-5; // 100 kev;
//   TString ntmp(GetTitle()); 
//   ntmp.ToUpper();
//   if(ntmp.Contains("10KEV")) {
//     cutele = cutgam = 1.e-5;
//   } else if(ntmp.Contains("50KEV")) {
//     cutele = cutgam = 5.e-5;
//   } else if(ntmp.Contains("100KEV")) {
//     cutele = cutgam = 1.e-4;
//   } else if(ntmp.Contains("200KEV")) {
//     cutele = cutgam = 2.e-4;
//   } else if(ntmp.Contains("500KEV")) {
//     cutele = cutgam = 5.e-4;
//   }

//   gMC->Gstpar(idtmed[1600],"CUTGAM", cutgam);
//   gMC->Gstpar(idtmed[1600],"CUTELE", cutele); // 1MEV -> 0.1MEV; 15-aug-05
//   gMC->Gstpar(idtmed[1600],"BCUTE",  cutgam);  // BCUTE and BCUTM start from GUTGUM
//   gMC->Gstpar(idtmed[1600],"BCUTM",  cutgam);  // BCUTE and BCUTM start from GUTGUM
//   // --- Generate explicitly delta rays in Lead ---
//   gMC->Gstpar(idtmed[1600], "LOSS", 3) ;
//   gMC->Gstpar(idtmed[1600], "DRAY", 1) ;
//   gMC->Gstpar(idtmed[1600], "DCUTE", cutele) ;
//   gMC->Gstpar(idtmed[1600], "DCUTM", cutele) ;

// // --- in aluminium parts ---
//   gMC->Gstpar(idtmed[1602],"CUTGAM", cutgam) ;
//   gMC->Gstpar(idtmed[1602],"CUTELE", cutele) ;
//   gMC->Gstpar(idtmed[1602],"BCUTE",  cutgam);  // BCUTE and BCUTM start from GUTGUM
//   gMC->Gstpar(idtmed[1602],"BCUTM",  cutgam);  // BCUTE and BCUTM start from GUTGUM
//   gMC->Gstpar(idtmed[1602], "LOSS",3.) ;
//   gMC->Gstpar(idtmed[1602], "DRAY",1.) ;
//   gMC->Gstpar(idtmed[1602], "DCUTE", cutele) ;
//   gMC->Gstpar(idtmed[1602], "DCUTM", cutele) ;

// // --- and finally thresholds for photons and electrons in the scintillator ---
//   gMC->Gstpar(idtmed[1601],"CUTGAM", cutgam) ;
//   gMC->Gstpar(idtmed[1601],"CUTELE", cutele) ;// 1MEV -> 0.1MEV; 15-aug-05
//   gMC->Gstpar(idtmed[1601],"BCUTE",  cutgam);  // BCUTE and BCUTM start from GUTGUM
//   gMC->Gstpar(idtmed[1601],"BCUTM",  cutgam);  // BCUTE and BCUTM start from GUTGUM
//   gMC->Gstpar(idtmed[1601], "LOSS",3) ; // generate delta rays 
//   gMC->Gstpar(idtmed[1601], "DRAY",1) ;
//   gMC->Gstpar(idtmed[1601], "DCUTE", cutele) ;
//   gMC->Gstpar(idtmed[1601], "DCUTM", cutele) ;

//   // S steel - 
//   gMC->Gstpar(idtmed[1603],"CUTGAM", cutgam);
//   gMC->Gstpar(idtmed[1603],"CUTELE", cutele);
//   gMC->Gstpar(idtmed[1603],"BCUTE",  cutgam);  // BCUTE and BCUTM start from GUTGUM
//   gMC->Gstpar(idtmed[1603],"BCUTM",  cutgam);  // BCUTE and BCUTM start from GUTGUM
//   // --- Generate explicitly delta rays 
//   gMC->Gstpar(idtmed[1603], "LOSS",3);
//   gMC->Gstpar(idtmed[1603], "DRAY",1);
//   gMC->Gstpar(idtmed[1603], "DCUTE", cutele) ;
//   gMC->Gstpar(idtmed[1603], "DCUTM", cutele) ;

//   AliEMCALGeometry* geom = GetGeometry();
//   if(geom->GetILOSS()>=0) {
//     for(int i=1600; i<=1603; i++) gMC->Gstpar(idtmed[i], "LOSS", geom->GetILOSS()) ; 
//   } 
//   if(geom->GetIHADR()>=0) {
//     for(int i=1600; i<=1603; i++) gMC->Gstpar(idtmed[i], "HADR", geom->GetIHADR()) ; 
//   }
// }

//____________________________________________________________________________
AliDigitizer* AliEMCAL::CreateDigitizer(AliDigitizationInput* digInput) const
{
  //create and return the digitizer
  return new AliEMCALDigitizer(digInput);
}

//____________________________________________________________________________
void AliEMCAL::CreateMaterials()
{
  // Definitions of materials to build EMCAL and associated tracking media.
  // media number in idtmed are 1599 to 1698.
  // --- Air ---               
  Float_t aAir[4]={12.0107,14.0067,15.9994,39.948};
  Float_t zAir[4]={6.,7.,8.,18.};
  Float_t wAir[4]={0.000124,0.755267,0.231781,0.012827};
  Float_t dAir = 1.20479E-3;
  AliMixture(0, "Air$", aAir, zAir, dAir, 4, wAir) ;

  // --- Lead ---                                                                     
  AliMaterial(1, "Pb$", 207.2, 82, 11.35, 0.56, 0., 0, 0) ;


  // --- The polysterene scintillator (CH) ---
  Float_t aP[2] = {12.011, 1.00794} ;
  Float_t zP[2] = {6.0, 1.0} ;
  Float_t wP[2] = {1.0, 1.0} ;
  Float_t dP = 1.032 ;

  AliMixture(2, "Polystyrene$", aP, zP, dP, -2, wP) ;

  // --- Aluminium ---
  AliMaterial(3, "Al$", 26.98, 13., 2.7, 8.9, 999., 0, 0) ;
  // ---         Absorption length is ignored ^

  // 25-aug-04 by PAI - see  PMD/AliPMDv0.cxx for STEEL definition
  Float_t asteel[4] = { 55.847,51.9961,58.6934,28.0855 };
  Float_t zsteel[4] = { 26.,24.,28.,14. };
  Float_t wsteel[4] = { .715,.18,.1,.005 };
  AliMixture(4, "STAINLESS STEEL$", asteel, zsteel, 7.88, 4, wsteel);

  // Oct 26,2010 : Multipurpose Copy Paper UNV-21200), weiht 75 g/m**2. 
  // *Cellulose C6H10O5
  //    Component C  A=12.01   Z=6.    W=6./21.
  //    Component H  A=1.      Z=1.    W=10./21.
  //    Component O  A=16.     Z=8.    W=5./21.
  Float_t apaper[3] = { 12.01, 1.0, 16.0};
  Float_t zpaper[3] = {  6.0,  1.0,  8.0};
  Float_t wpaper[3] = {6./21., 10./21., 5./21.};
  AliMixture(5, "BondPaper$", apaper, zpaper, 0.75, 3, wpaper);
 
  // DEFINITION OF THE TRACKING MEDIA
  // Look to the $ALICE_ROOT/data/galice.cuts for particular values
  // of cuts.
  // Don't forget to add a new tracking medium with non-default cuts

  // for EMCAL: idtmed[1599->1698] equivalent to fIdtmed[0->100]
  Int_t   isxfld = ((AliMagF*)TGeoGlobalMagField::Instance()->GetField())->Integ() ;
  Float_t sxmgmx = ((AliMagF*)TGeoGlobalMagField::Instance()->GetField())->Max() ;

  // Air                                                                         -> idtmed[1599]
 AliMedium(0, "Air$", 0, 0,
	     isxfld, sxmgmx, 10.0, 1.0, 0.1, 0.1, 10.0, 0, 0) ;

  // The Lead                                                                      -> idtmed[1600]
 
  AliMedium(1, "Lead$", 1, 0,
	     isxfld, sxmgmx, 10.0, 0.1, 0.1, 0.1, 0.1, 0, 0) ;

 // The scintillator of the CPV made of Polystyrene scintillator                   -> idtmed[1601]
  float deemax = 0.1; // maximum fractional energy loss in one step (0 < DEEMAX < deemax )
  AliMedium(2, "Scintillator$", 2, 1,
            isxfld, sxmgmx, 10.0, 0.001, deemax, 0.001, 0.001, 0, 0) ;

  // Various Aluminium parts made of Al                                            -> idtmed[1602]
  AliMedium(3, "Al$", 3, 0,
             isxfld, sxmgmx, 10.0, 0.1, 0.1, 0.001, 0.001, 0, 0) ;

  // 25-aug-04 by PAI : see  PMD/AliPMDv0.cxx for STEEL definition                 -> idtmed[1603]
  AliMedium(4, "S steel$", 4, 0, 
             isxfld, sxmgmx, 10.0, 0.01, 0.1, 0.001, 0.001, 0, 0) ;

  // Oct 26,2010; Nov 24,2010                                                      -> idtmed[1604]
  deemax = 0.01;
  AliMedium(5, "Paper$", 5, 0, 
             isxfld, sxmgmx, 10.0, deemax, 0.1, 0.001, 0.001, 0, 0) ;


  //set constants for Birk's Law implentation
  fBirkC0 =  1;
  fBirkC1 =  0.013/dP;
  fBirkC2 =  9.6e-6/(dP * dP);

}

//____________________________________________________________________________
void  AliEMCAL::Init()
{ 
  // Init
  //Not needed, modify $ALICE_ROOT/data/galice.cuts instead.
  //Load the modified one in the configuration file with SetTransPar
  //DefineMediumParameters(); 
}     

//____________________________________________________________________________
void AliEMCAL::Digits2Raw() {

  static AliEMCALRawUtils rawUtils;
  rawUtils.Digits2Raw();

}
//____________________________________________________________________________
void AliEMCAL::Hits2SDigits()  
{ 
// create summable digits

  GetGeometry();
  AliEMCALSDigitizer emcalDigitizer(fLoader->GetRunLoader()->GetFileName().Data()) ;
  emcalDigitizer.SetEventRange(0, -1) ; // do all the events
  emcalDigitizer.Digitize() ;
}

//______________________________________________________________________
Bool_t AliEMCAL::Raw2SDigits(AliRawReader* rawReader){
  
  // Conversion from raw data to EMCAL sdigits. 
  // Does the same as AliEMCALReconstructor::ConvertDigits()
  // Needed to embed real data and simulation
  // Works on a single-event basis
  
  rawReader->Reset() ; 
  
  //Get/create the sdigits tree and array
	AliRunLoader *rl = AliRunLoader::Instance();
	AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL")); 
  
  if(!emcalLoader){
    AliFatal("NULL loader");
    return kFALSE;
  }
  
  emcalLoader->GetEvent();
  emcalLoader->LoadSDigits("UPDATE");

  TTree * treeS = emcalLoader->TreeS();
  if ( !treeS ) { 
    emcalLoader->MakeSDigitsContainer();
    treeS = emcalLoader->TreeS();
  }
  
  if(!emcalLoader->SDigits()) {
    AliFatal("No sdigits array available\n");
    return kFALSE;
  }
  
  TClonesArray * sdigits = emcalLoader->SDigits();
  sdigits->Clear("C");  
  
  //Trigger sdigits
  if(!fTriggerData)fTriggerData = new AliEMCALTriggerData();
  fTriggerData->SetMode(1);	
  TClonesArray *digitsTrg = new TClonesArray("AliEMCALTriggerRawDigit", 32 * 96);    
  Int_t bufsize = 32000;
  treeS->Branch("EMTRG", &digitsTrg, bufsize);
  
  
  //Only physics events
  if (rawReader->GetType()== AliRawEventHeaderBase::kPhysicsEvent) {

    if(!fgRawUtils)  fgRawUtils   = new AliEMCALRawUtils;
    //must be done here because, in constructor, option is not yet known
    fgRawUtils->SetOption(GetOption());
    
    // Set parameters from OCDB to raw utils
    AliEMCALRecParam* recpar = emcalLoader->ReconstructionParameters(0);
    // fgRawUtils->SetRawFormatHighLowGainFactor(recpar->GetHighLowGainFactor());
    // fgRawUtils->SetRawFormatOrder(recpar->GetOrderParameter());
    // fgRawUtils->SetRawFormatTau(recpar->GetTau());
    fgRawUtils->SetNoiseThreshold(recpar->GetNoiseThreshold());
    fgRawUtils->SetNPedSamples(recpar->GetNPedSamples());
    fgRawUtils->SetRemoveBadChannels(recpar->GetRemoveBadChannels());
    fgRawUtils->SetFittingAlgorithm(recpar->GetFittingAlgorithm());
    fgRawUtils->SetFALTROUsage(recpar->UseFALTRO());
    // fgRawUtils->SetTimeMin(recpar->GetTimeMin());
    // fgRawUtils->SetTimeMax(recpar->GetTimeMax());

    //Fit
    fgRawUtils->Raw2Digits(rawReader,sdigits,emcalLoader->PedestalData(),digitsTrg,fTriggerData);
    
  }//skip calibration event
  else{
    AliDebug(1," Calibration Event, skip!");
  }
  
  //Final arrangements of the array, set all sdigits as embedded
  sdigits->Sort() ;
  for (Int_t iSDigit = 0 ; iSDigit < sdigits->GetEntriesFast() ; iSDigit++) { 
    AliEMCALDigit * sdigit = dynamic_cast<AliEMCALDigit *>(sdigits->At(iSDigit)) ;
    if(sdigit){
      sdigit->SetIndexInList(iSDigit) ;
      sdigit->SetType(AliEMCALDigit::kEmbedded);
    }
    else {
      AliFatal("sdigit is NULL!");
    }
  }	
  
  AliDebug(1,Form("Embedded sdigits entries %d \n",sdigits->GetEntriesFast()));
  
  //Write array, clean arrays, unload ..
  
  Int_t bufferSize = 32000 ;    
  TBranch * sdigitsBranch = treeS->GetBranch("EMCAL");
  if (sdigitsBranch)
    sdigitsBranch->SetAddress(&sdigits);
  else
    treeS->Branch("EMCAL",&sdigits,bufferSize);
  
  treeS->Fill();
  emcalLoader->WriteSDigits("OVERWRITE");
  emcalLoader->UnloadSDigits();

  digitsTrg->Delete();
  delete digitsTrg;
    
  return kTRUE;
  
}

//____________________________________________________________________________

AliLoader* AliEMCAL::MakeLoader(const char* topfoldername)
{
//different behaviour than standard (singleton getter)
// --> to be discussed and made eventually coherent
 fLoader = new AliEMCALLoader(GetName(),topfoldername);
 return fLoader;
}

//____________________________________________________________________________

AliEMCALGeometry* AliEMCAL::GetGeometry() const
{
  //Initializes and returns geometry
  
  // Pass the transpor model name (Geant3, Geant4, Fluka) and title to the geometry
  TString mcname   = "";
  TString mctitle  = "";
  if(gMC){
    mcname  = gMC->GetName()  ;
    mctitle = gMC->GetTitle() ;
  }
  
  //Check if run number and requested geometry correspond to the same geometry as
  //in real data taking. To prevent errors in official simulation productions
  if(!(AliEMCALGeometry::GetInstance()))
  {
    // Check the transport model name and option, set sampling fraction depending on it
    
    if(!fCheckRunNumberAndGeoVersion){// Set geometry with the name used in the configuration file
      return AliEMCALGeometry::GetInstance(GetTitle(),"EMCAL",mcname,mctitle) ;
    }
    else{//Check run number and version and set the corresponding one.
      //Get run number
      //AliRunLoader *rl = AliRunLoader::Instance();
      //Int_t runNumber = rl->GetRunNumber();
      
      AliCDBManager* man = AliCDBManager::Instance();
      Int_t runNumber = man->GetRun();
      
      //Instanciate geometry depending on the run number
      TString geoName(GetTitle());
      if(runNumber > 104064 && runNumber <= 140000 ){//2009-2010 runs
        //First year geometry, 4 SM.
        
        if(!geoName.Contains("FIRSTYEARV1")){ 
          AliInfo(Form("*** ATTENTION *** \n \t Specified geometry name <<%s>> for run %d is not considered! \n \t In use <<EMCAL_FIRSTYEARV1>>, check run number and year \n ",
                       geoName.Data(),runNumber)); }
        else {
          AliDebug(1,"Initialized geometry with name <<EMCAL_FIRSTYEARV1>>");}
        
        return AliEMCALGeometry::GetInstance("EMCAL_FIRSTYEARV1","EMCAL",mcname,mctitle) ;// Set geometry with the name used in the configuration file
      }
      else{ //Default geometry
        //Complete EMCAL geometry, 10 SM.
  
        if(!geoName.Contains("COMPLETEV1")){
          AliInfo(Form("*** ATTENTION *** \n \t Specified geometry name <<%s>>  for run %d is  not considered! \n \t In use <<EMCAL_COMPLETEV1>>, check run number and year \n ",
                       geoName.Data(),runNumber));}
        else {
          AliDebug(1,"Initialized geometry with name <<EMCAL_COMPLETEV1>>");}

        return AliEMCALGeometry::GetInstance("EMCAL_COMPLETEV1","EMCAL",mcname,mctitle) ;// Set geometry with the name used in the configuration file
      }
    }
  }// Init geometry for the first time
  
  
  return AliEMCALGeometry::GetInstance();
    
}
Commit	Line	Data
	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	// Base Class for EMCAL description:
	19	// This class contains material definitions
	20	// for the EMCAL - It does not place the detector in Alice
	21	//*-- Author: Yves Schutz (SUBATECH)
	22	//
	23	//*-- Additional Contributions: Sahal Yacoob (LBNL/UCT)
	24	// : Alexei Pavlinov (WSU)
	25	//
	26	//////////////////////////////////////////////////////////////////////////////
	27
	28	// --- ROOT system ---
	29	class TFile;
	30	#include <TFolder.h>
	31	#include <TGeoGlobalMagField.h>
	32	#include <TGraph.h>
	33	#include <TH1F.h>
	34	#include <TRandom.h>
	35	#include <TTree.h>
	36	#include <TVirtualMC.h>
	37
	38	// --- Standard library ---
	39
	40	// --- AliRoot header files ---
	41	#include "AliMagF.h"
	42	#include "AliLog.h"
	43	#include "AliEMCAL.h"
	44	#include "AliRun.h"
	45	#include "AliRunLoader.h"
	46	#include "AliCDBManager.h"
	47	#include "AliEMCALLoader.h"
	48	#include "AliEMCALSDigitizer.h"
	49	#include "AliEMCALDigitizer.h"
	50	#include "AliEMCALDigit.h"
	51	#include "AliEMCALRawUtils.h"
	52	#include "AliCDBManager.h"
	53	#include "AliCDBEntry.h"
	54	#include "AliEMCALRawUtils.h"
	55	#include "AliRawReader.h"
	56	#include "AliEMCALTriggerData.h"
	57	#include "AliEMCALRecParam.h"
	58	#include "AliRawEventHeaderBase.h"
	59
	60	ClassImp(AliEMCAL)
	61
	62	//for embedding
	63	AliEMCALRawUtils* AliEMCAL::fgRawUtils = 0; // EMCAL raw utilities class
	64
	65	//____________________________________________________________________________
	66	AliEMCAL::AliEMCAL()
	67	: AliDetector(),
	68	fBirkC0(0),
	69	fBirkC1(0.),
	70	fBirkC2(0.),
	71	fGeometry(0),
	72	fCheckRunNumberAndGeoVersion(kTRUE),
	73	fTriggerData(0x0)
	74	{
	75	// Default ctor
	76	fName = "EMCAL" ;
	77	InitConstants();
	78
	79	// Should call AliEMCALGeometry::GetInstance(EMCAL->GetTitle(),"") for getting EMCAL geometry
	80	}
	81
	82	//____________________________________________________________________________
	83	AliEMCAL::AliEMCAL(const char* name, const char* title)
	84	: AliDetector(name,title),
	85	fBirkC0(0),
	86	fBirkC1(0.),
	87	fBirkC2(0.),
	88	fGeometry(0),
	89	fCheckRunNumberAndGeoVersion(kTRUE),
	90	fTriggerData(0x0)
	91	{
	92	// ctor : title is used to identify the layout
	93	InitConstants();
	94
	95	}
	96
	97	//____________________________________________________________________________
	98	AliEMCAL::~AliEMCAL()
	99	{
	100	//dtor
	101	delete fgRawUtils;
	102	delete fTriggerData;
	103	}
	104
	105	//____________________________________________________________________________
	106	void AliEMCAL::InitConstants()
	107	{
	108	//initialize EMCAL values
	109	fBirkC0 = 1;
	110	fBirkC1 = 0.013/1.032;
	111	fBirkC2 = 9.6e-6/(1.032 * 1.032);
	112	}
	113
	114	//Not needed, modify $ALICE_ROOT/data/galice.cuts instead.
	115	//Load the modified one in the configuration file with SetTransPar
	116	// //____________________________________________________________________________
	117	// void AliEMCAL::DefineMediumParameters()
	118	// {
	119	// //
	120	// // EMCAL cuts (Geant3)
	121	// //
	122	// Int_t * idtmed = fIdtmed->GetArray() - 1599 ;
	123	// // --- Set decent energy thresholds for gamma and electron tracking
	124
	125	// // Tracking threshold for photons and electrons in Lead
	126	// Float_t cutgam=10.e-5; // 100 kev;
	127	// Float_t cutele=10.e-5; // 100 kev;
	128	// TString ntmp(GetTitle());
	129	// ntmp.ToUpper();
	130	// if(ntmp.Contains("10KEV")) {
	131	// cutele = cutgam = 1.e-5;
	132	// } else if(ntmp.Contains("50KEV")) {
	133	// cutele = cutgam = 5.e-5;
	134	// } else if(ntmp.Contains("100KEV")) {
	135	// cutele = cutgam = 1.e-4;
	136	// } else if(ntmp.Contains("200KEV")) {
	137	// cutele = cutgam = 2.e-4;
	138	// } else if(ntmp.Contains("500KEV")) {
	139	// cutele = cutgam = 5.e-4;
	140	// }
	141
	142	// gMC->Gstpar(idtmed[1600],"CUTGAM", cutgam);
	143	// gMC->Gstpar(idtmed[1600],"CUTELE", cutele); // 1MEV -> 0.1MEV; 15-aug-05
	144	// gMC->Gstpar(idtmed[1600],"BCUTE", cutgam); // BCUTE and BCUTM start from GUTGUM
	145	// gMC->Gstpar(idtmed[1600],"BCUTM", cutgam); // BCUTE and BCUTM start from GUTGUM
	146	// // --- Generate explicitly delta rays in Lead ---
	147	// gMC->Gstpar(idtmed[1600], "LOSS", 3) ;
	148	// gMC->Gstpar(idtmed[1600], "DRAY", 1) ;
	149	// gMC->Gstpar(idtmed[1600], "DCUTE", cutele) ;
	150	// gMC->Gstpar(idtmed[1600], "DCUTM", cutele) ;
	151
	152	// // --- in aluminium parts ---
	153	// gMC->Gstpar(idtmed[1602],"CUTGAM", cutgam) ;
	154	// gMC->Gstpar(idtmed[1602],"CUTELE", cutele) ;
	155	// gMC->Gstpar(idtmed[1602],"BCUTE", cutgam); // BCUTE and BCUTM start from GUTGUM
	156	// gMC->Gstpar(idtmed[1602],"BCUTM", cutgam); // BCUTE and BCUTM start from GUTGUM
	157	// gMC->Gstpar(idtmed[1602], "LOSS",3.) ;
	158	// gMC->Gstpar(idtmed[1602], "DRAY",1.) ;
	159	// gMC->Gstpar(idtmed[1602], "DCUTE", cutele) ;
	160	// gMC->Gstpar(idtmed[1602], "DCUTM", cutele) ;
	161
	162	// // --- and finally thresholds for photons and electrons in the scintillator ---
	163	// gMC->Gstpar(idtmed[1601],"CUTGAM", cutgam) ;
	164	// gMC->Gstpar(idtmed[1601],"CUTELE", cutele) ;// 1MEV -> 0.1MEV; 15-aug-05
	165	// gMC->Gstpar(idtmed[1601],"BCUTE", cutgam); // BCUTE and BCUTM start from GUTGUM
	166	// gMC->Gstpar(idtmed[1601],"BCUTM", cutgam); // BCUTE and BCUTM start from GUTGUM
	167	// gMC->Gstpar(idtmed[1601], "LOSS",3) ; // generate delta rays
	168	// gMC->Gstpar(idtmed[1601], "DRAY",1) ;
	169	// gMC->Gstpar(idtmed[1601], "DCUTE", cutele) ;
	170	// gMC->Gstpar(idtmed[1601], "DCUTM", cutele) ;
	171
	172	// // S steel -
	173	// gMC->Gstpar(idtmed[1603],"CUTGAM", cutgam);
	174	// gMC->Gstpar(idtmed[1603],"CUTELE", cutele);
	175	// gMC->Gstpar(idtmed[1603],"BCUTE", cutgam); // BCUTE and BCUTM start from GUTGUM
	176	// gMC->Gstpar(idtmed[1603],"BCUTM", cutgam); // BCUTE and BCUTM start from GUTGUM
	177	// // --- Generate explicitly delta rays
	178	// gMC->Gstpar(idtmed[1603], "LOSS",3);
	179	// gMC->Gstpar(idtmed[1603], "DRAY",1);
	180	// gMC->Gstpar(idtmed[1603], "DCUTE", cutele) ;
	181	// gMC->Gstpar(idtmed[1603], "DCUTM", cutele) ;
	182
	183	// AliEMCALGeometry* geom = GetGeometry();
	184	// if(geom->GetILOSS()>=0) {
	185	// for(int i=1600; i<=1603; i++) gMC->Gstpar(idtmed[i], "LOSS", geom->GetILOSS()) ;
	186	// }
	187	// if(geom->GetIHADR()>=0) {
	188	// for(int i=1600; i<=1603; i++) gMC->Gstpar(idtmed[i], "HADR", geom->GetIHADR()) ;
	189	// }
	190	// }
	191
	192	//____________________________________________________________________________
	193	AliDigitizer* AliEMCAL::CreateDigitizer(AliDigitizationInput* digInput) const
	194	{
	195	//create and return the digitizer
	196	return new AliEMCALDigitizer(digInput);
	197	}
	198
	199	//____________________________________________________________________________
	200	void AliEMCAL::CreateMaterials()
	201	{
	202	// Definitions of materials to build EMCAL and associated tracking media.
	203	// media number in idtmed are 1599 to 1698.
	204	// --- Air ---
	205	Float_t aAir[4]={12.0107,14.0067,15.9994,39.948};
	206	Float_t zAir[4]={6.,7.,8.,18.};
	207	Float_t wAir[4]={0.000124,0.755267,0.231781,0.012827};
	208	Float_t dAir = 1.20479E-3;
	209	AliMixture(0, "Air$", aAir, zAir, dAir, 4, wAir) ;
	210
	211	// --- Lead ---
	212	AliMaterial(1, "Pb$", 207.2, 82, 11.35, 0.56, 0., 0, 0) ;
	213
	214
	215	// --- The polysterene scintillator (CH) ---
	216	Float_t aP[2] = {12.011, 1.00794} ;
	217	Float_t zP[2] = {6.0, 1.0} ;
	218	Float_t wP[2] = {1.0, 1.0} ;
	219	Float_t dP = 1.032 ;
	220
	221	AliMixture(2, "Polystyrene$", aP, zP, dP, -2, wP) ;
	222
	223	// --- Aluminium ---
	224	AliMaterial(3, "Al$", 26.98, 13., 2.7, 8.9, 999., 0, 0) ;
	225	// --- Absorption length is ignored ^
	226
	227	// 25-aug-04 by PAI - see PMD/AliPMDv0.cxx for STEEL definition
	228	Float_t asteel[4] = { 55.847,51.9961,58.6934,28.0855 };
	229	Float_t zsteel[4] = { 26.,24.,28.,14. };
	230	Float_t wsteel[4] = { .715,.18,.1,.005 };
	231	AliMixture(4, "STAINLESS STEEL$", asteel, zsteel, 7.88, 4, wsteel);
	232
	233	// Oct 26,2010 : Multipurpose Copy Paper UNV-21200), weiht 75 g/m**2.
	234	// *Cellulose C6H10O5
	235	// Component C A=12.01 Z=6. W=6./21.
	236	// Component H A=1. Z=1. W=10./21.
	237	// Component O A=16. Z=8. W=5./21.
	238	Float_t apaper[3] = { 12.01, 1.0, 16.0};
	239	Float_t zpaper[3] = { 6.0, 1.0, 8.0};
	240	Float_t wpaper[3] = {6./21., 10./21., 5./21.};
	241	AliMixture(5, "BondPaper$", apaper, zpaper, 0.75, 3, wpaper);
	242
	243	// DEFINITION OF THE TRACKING MEDIA
	244	// Look to the $ALICE_ROOT/data/galice.cuts for particular values
	245	// of cuts.
	246	// Don't forget to add a new tracking medium with non-default cuts
	247
	248	// for EMCAL: idtmed[1599->1698] equivalent to fIdtmed[0->100]
	249	Int_t isxfld = ((AliMagF*)TGeoGlobalMagField::Instance()->GetField())->Integ() ;
	250	Float_t sxmgmx = ((AliMagF*)TGeoGlobalMagField::Instance()->GetField())->Max() ;
	251
	252	// Air -> idtmed[1599]
	253	AliMedium(0, "Air$", 0, 0,
	254	isxfld, sxmgmx, 10.0, 1.0, 0.1, 0.1, 10.0, 0, 0) ;
	255
	256	// The Lead -> idtmed[1600]
	257
	258	AliMedium(1, "Lead$", 1, 0,
	259	isxfld, sxmgmx, 10.0, 0.1, 0.1, 0.1, 0.1, 0, 0) ;
	260
	261	// The scintillator of the CPV made of Polystyrene scintillator -> idtmed[1601]
	262	float deemax = 0.1; // maximum fractional energy loss in one step (0 < DEEMAX < deemax )
	263	AliMedium(2, "Scintillator$", 2, 1,
	264	isxfld, sxmgmx, 10.0, 0.001, deemax, 0.001, 0.001, 0, 0) ;
	265
	266	// Various Aluminium parts made of Al -> idtmed[1602]
	267	AliMedium(3, "Al$", 3, 0,
	268	isxfld, sxmgmx, 10.0, 0.1, 0.1, 0.001, 0.001, 0, 0) ;
	269
	270	// 25-aug-04 by PAI : see PMD/AliPMDv0.cxx for STEEL definition -> idtmed[1603]
	271	AliMedium(4, "S steel$", 4, 0,
	272	isxfld, sxmgmx, 10.0, 0.01, 0.1, 0.001, 0.001, 0, 0) ;
	273
	274	// Oct 26,2010; Nov 24,2010 -> idtmed[1604]
	275	deemax = 0.01;
	276	AliMedium(5, "Paper$", 5, 0,
	277	isxfld, sxmgmx, 10.0, deemax, 0.1, 0.001, 0.001, 0, 0) ;
	278
	279
	280	//set constants for Birk's Law implentation
	281	fBirkC0 = 1;
	282	fBirkC1 = 0.013/dP;
	283	fBirkC2 = 9.6e-6/(dP * dP);
	284
	285	}
	286
	287	//____________________________________________________________________________
	288	void AliEMCAL::Init()
	289	{
	290	// Init
	291	//Not needed, modify $ALICE_ROOT/data/galice.cuts instead.
	292	//Load the modified one in the configuration file with SetTransPar
	293	//DefineMediumParameters();
	294	}
	295
	296	//____________________________________________________________________________
	297	void AliEMCAL::Digits2Raw() {
	298
	299	static AliEMCALRawUtils rawUtils;
	300	rawUtils.Digits2Raw();
	301
	302	}
	303	//____________________________________________________________________________
	304	void AliEMCAL::Hits2SDigits()
	305	{
	306	// create summable digits
	307
	308	GetGeometry();
	309	AliEMCALSDigitizer emcalDigitizer(fLoader->GetRunLoader()->GetFileName().Data()) ;
	310	emcalDigitizer.SetEventRange(0, -1) ; // do all the events
	311	emcalDigitizer.Digitize() ;
	312	}
	313
	314	//______________________________________________________________________
	315	Bool_t AliEMCAL::Raw2SDigits(AliRawReader* rawReader){
	316
	317	// Conversion from raw data to EMCAL sdigits.
	318	// Does the same as AliEMCALReconstructor::ConvertDigits()
	319	// Needed to embed real data and simulation
	320	// Works on a single-event basis
	321
	322	rawReader->Reset() ;
	323
	324	//Get/create the sdigits tree and array
	325	AliRunLoader *rl = AliRunLoader::Instance();
	326	AliEMCALLoader emcalLoader = dynamic_cast<AliEMCALLoader>(rl->GetDetectorLoader("EMCAL"));
	327
	328	if(!emcalLoader){
	329	AliFatal("NULL loader");
	330	return kFALSE;
	331	}
	332
	333	emcalLoader->GetEvent();
	334	emcalLoader->LoadSDigits("UPDATE");
	335
	336	TTree * treeS = emcalLoader->TreeS();
	337	if ( !treeS ) {
	338	emcalLoader->MakeSDigitsContainer();
	339	treeS = emcalLoader->TreeS();
	340	}
	341
	342	if(!emcalLoader->SDigits()) {
	343	AliFatal("No sdigits array available\n");
	344	return kFALSE;
	345	}
	346
	347	TClonesArray * sdigits = emcalLoader->SDigits();
	348	sdigits->Clear("C");
	349
	350	//Trigger sdigits
	351	if(!fTriggerData)fTriggerData = new AliEMCALTriggerData();
	352	fTriggerData->SetMode(1);
	353	TClonesArray digitsTrg = new TClonesArray("AliEMCALTriggerRawDigit", 32 96);
	354	Int_t bufsize = 32000;
	355	treeS->Branch("EMTRG", &digitsTrg, bufsize);
	356
	357
	358	//Only physics events
	359	if (rawReader->GetType()== AliRawEventHeaderBase::kPhysicsEvent) {
	360
	361	if(!fgRawUtils) fgRawUtils = new AliEMCALRawUtils;
	362	//must be done here because, in constructor, option is not yet known
	363	fgRawUtils->SetOption(GetOption());
	364
	365	// Set parameters from OCDB to raw utils
	366	AliEMCALRecParam* recpar = emcalLoader->ReconstructionParameters(0);
	367	// fgRawUtils->SetRawFormatHighLowGainFactor(recpar->GetHighLowGainFactor());
	368	// fgRawUtils->SetRawFormatOrder(recpar->GetOrderParameter());
	369	// fgRawUtils->SetRawFormatTau(recpar->GetTau());
	370	fgRawUtils->SetNoiseThreshold(recpar->GetNoiseThreshold());
	371	fgRawUtils->SetNPedSamples(recpar->GetNPedSamples());
	372	fgRawUtils->SetRemoveBadChannels(recpar->GetRemoveBadChannels());
	373	fgRawUtils->SetFittingAlgorithm(recpar->GetFittingAlgorithm());
	374	fgRawUtils->SetFALTROUsage(recpar->UseFALTRO());
	375	// fgRawUtils->SetTimeMin(recpar->GetTimeMin());
	376	// fgRawUtils->SetTimeMax(recpar->GetTimeMax());
	377
	378	//Fit
	379	fgRawUtils->Raw2Digits(rawReader,sdigits,emcalLoader->PedestalData(),digitsTrg,fTriggerData);
	380
	381	}//skip calibration event
	382	else{
	383	AliDebug(1," Calibration Event, skip!");
	384	}
	385
	386	//Final arrangements of the array, set all sdigits as embedded
	387	sdigits->Sort() ;
	388	for (Int_t iSDigit = 0 ; iSDigit < sdigits->GetEntriesFast() ; iSDigit++) {
	389	AliEMCALDigit * sdigit = dynamic_cast<AliEMCALDigit *>(sdigits->At(iSDigit)) ;
	390	if(sdigit){
	391	sdigit->SetIndexInList(iSDigit) ;
	392	sdigit->SetType(AliEMCALDigit::kEmbedded);
	393	}
	394	else {
	395	AliFatal("sdigit is NULL!");
	396	}
	397	}
	398
	399	AliDebug(1,Form("Embedded sdigits entries %d \n",sdigits->GetEntriesFast()));
	400
	401	//Write array, clean arrays, unload ..
	402
	403	Int_t bufferSize = 32000 ;
	404	TBranch * sdigitsBranch = treeS->GetBranch("EMCAL");
	405	if (sdigitsBranch)
	406	sdigitsBranch->SetAddress(&sdigits);
	407	else
	408	treeS->Branch("EMCAL",&sdigits,bufferSize);
	409
	410	treeS->Fill();
	411	emcalLoader->WriteSDigits("OVERWRITE");
	412	emcalLoader->UnloadSDigits();
	413
	414	digitsTrg->Delete();
	415	delete digitsTrg;
	416
	417	return kTRUE;
	418
	419	}
	420
	421	//____________________________________________________________________________
	422
	423	AliLoader* AliEMCAL::MakeLoader(const char* topfoldername)
	424	{
	425	//different behaviour than standard (singleton getter)
	426	// --> to be discussed and made eventually coherent
	427	fLoader = new AliEMCALLoader(GetName(),topfoldername);
	428	return fLoader;
	429	}
	430
	431	//____________________________________________________________________________
	432
	433	AliEMCALGeometry* AliEMCAL::GetGeometry() const
	434	{
	435	//Initializes and returns geometry
	436
	437	// Pass the transpor model name (Geant3, Geant4, Fluka) and title to the geometry
	438	TString mcname = "";
	439	TString mctitle = "";
	440	if(gMC){
	441	mcname = gMC->GetName() ;
	442	mctitle = gMC->GetTitle() ;
	443	}
	444
	445	//Check if run number and requested geometry correspond to the same geometry as
	446	//in real data taking. To prevent errors in official simulation productions
	447	if(!(AliEMCALGeometry::GetInstance()))
	448	{
	449	// Check the transport model name and option, set sampling fraction depending on it
	450
	451	if(!fCheckRunNumberAndGeoVersion){// Set geometry with the name used in the configuration file
	452	return AliEMCALGeometry::GetInstance(GetTitle(),"EMCAL",mcname,mctitle) ;
	453	}
	454	else{//Check run number and version and set the corresponding one.
	455	//Get run number
	456	//AliRunLoader *rl = AliRunLoader::Instance();
	457	//Int_t runNumber = rl->GetRunNumber();
	458
	459	AliCDBManager* man = AliCDBManager::Instance();
	460	Int_t runNumber = man->GetRun();
	461
	462	//Instanciate geometry depending on the run number
	463	TString geoName(GetTitle());
	464	if(runNumber > 104064 && runNumber <= 140000 ){//2009-2010 runs
	465	//First year geometry, 4 SM.
	466
	467	if(!geoName.Contains("FIRSTYEARV1")){
	468	AliInfo(Form("* ATTENTION * \n \t Specified geometry name <<%s>> for run %d is not considered! \n \t In use <<EMCAL_FIRSTYEARV1>>, check run number and year \n ",
	469	geoName.Data(),runNumber)); }
	470	else {
	471	AliDebug(1,"Initialized geometry with name <<EMCAL_FIRSTYEARV1>>");}
	472
	473	return AliEMCALGeometry::GetInstance("EMCAL_FIRSTYEARV1","EMCAL",mcname,mctitle) ;// Set geometry with the name used in the configuration file
	474	}
	475	else{ //Default geometry
	476	//Complete EMCAL geometry, 10 SM.
	477
	478	if(!geoName.Contains("COMPLETEV1")){
	479	AliInfo(Form("* ATTENTION * \n \t Specified geometry name <<%s>> for run %d is not considered! \n \t In use <<EMCAL_COMPLETEV1>>, check run number and year \n ",
	480	geoName.Data(),runNumber));}
	481	else {
	482	AliDebug(1,"Initialized geometry with name <<EMCAL_COMPLETEV1>>");}
	483
	484	return AliEMCALGeometry::GetInstance("EMCAL_COMPLETEV1","EMCAL",mcname,mctitle) ;// Set geometry with the name used in the configuration file
	485	}
	486	}
	487	}// Init geometry for the first time
	488
	489
	490	return AliEMCALGeometry::GetInstance();
	491
	492	}