EMCAL/AliEMCAL.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 /* History of cvs commits:
  18  *
  19  * $Log$
  20  * Revision 1.51  2007/02/24 20:42:35  pavlinov
  21  * fixed error of Geant3 parameters initialisation
  22  *
  23  * Revision 1.50  2007/02/05 10:43:25  hristov
  24  * Changes for correct initialization of Geant4 (Mihaela)
  25  *
  26  * Revision 1.49  2007/01/22 17:29:12  pavlinov
  27  * EMCAL geometry can be created independently form anything now
  28  *
  29  * Revision 1.48  2006/12/19 02:34:13  pavlinov
  30  * clean up the EMCAL name scheme : super module -> module -> tower (or cell)
  31  *
  32  * Revision 1.47  2006/12/05 17:12:03  gustavo
  33  * Updated AliEMCAL::Digits2Raw, reads first provisional RCU mapping files to make Raw data with new AliCaloAltroMapping and AliCaloRawStream
  34  *
  35  *
  36  */
  37 //_________________________________________________________________________
  38 // Base Class for EMCAL description:
  39 // This class contains material definitions
  40 // for the EMCAL - It does not place the detector in Alice
  41 //*-- Author: Yves Schutz (SUBATECH)
  42 //
  43 //*-- Additional Contributions: Sahal Yacoob (LBNL/UCT)
  44 //
  45 //////////////////////////////////////////////////////////////////////////////
  46
  47 // --- ROOT system ---
  48 class TFile;
  49 #include <TFolder.h>
  50 #include <TTree.h>
  51 #include <TVirtualMC.h>
  52 #include <TH1F.h>
  53 #include <TF1.h>
  54 #include <TRandom.h>
  55 #include <TGraph.h>
  56
  57 // --- Standard library ---
  58
  59 // --- AliRoot header files ---
  60 #include "AliMagF.h"
  61 #include "AliEMCAL.h"
  62 #include "AliRun.h"
  63 #include "AliEMCALLoader.h"
  64 #include "AliEMCALSDigitizer.h"
  65 #include "AliEMCALDigitizer.h"
  66 #include "AliEMCALDigit.h"
  67 //#include "AliAltroMapping.h"
  68 #include "AliCaloAltroMapping.h"
  69 #include "AliAltroBuffer.h"
  70 #include "AliRawReader.h"
  71 #include "AliCaloRawStream.h"
  72 #include "AliDAQ.h"
  73
  74 ClassImp(AliEMCAL)
  75 Double_t AliEMCAL::fgCapa        = 1.;        // 1pF
  76 Int_t    AliEMCAL::fgOrder       = 2 ;
  77 Double_t AliEMCAL::fgTimeMax     = 2.56E-5 ;  // each sample is over 100 ns fTimeMax/fTimeBins
  78 Double_t AliEMCAL::fgTimePeak    = 4.1E-6 ;   // 4 micro seconds
  79 Double_t AliEMCAL::fgTimeTrigger = 100E-9 ;      // 100ns, just for a reference
  80 // some digitization constants
  81 Int_t    AliEMCAL::fgThreshold = 1;
  82 Int_t    AliEMCAL::fgDDLPerSuperModule = 2;  // 2 ddls per SuperModule
  83
  84 //____________________________________________________________________________
  85 AliEMCAL::AliEMCAL()
  86   : AliDetector(),
  87     fBirkC0(0),
  88     fBirkC1(0.),
  89     fBirkC2(0.),
  90     fHighCharge(0.),
  91     fHighGain(0.),
  92     fHighLowGainFactor(0.),
  93     fLowGainOffset(0),
  94     fGeometry(0)
  95 {
  96   // Default ctor
  97   fName = "EMCAL" ;
  98   InitConstants();
  99   // Should call  AliEMCALGeometry::GetInstance(EMCAL->GetTitle(),"") for getting EMCAL geometry
 100 }
 101
 102 //____________________________________________________________________________
 103 AliEMCAL::AliEMCAL(const char* name, const char* title)
 104   : AliDetector(name,title),
 105     fBirkC0(0),
 106     fBirkC1(0.),
 107     fBirkC2(0.),
 108     fHighCharge(0.),
 109     fHighGain(0.),
 110     fHighLowGainFactor(0.),
 111     fLowGainOffset(0),
 112     fGeometry(0)
 113 {
 114   //   ctor : title is used to identify the layout
 115   InitConstants();
 116 }
 117
 118 //____________________________________________________________________________
 119 AliEMCAL::~AliEMCAL()
 120 {
 121   //dtor
 122 }
 123
 124 //____________________________________________________________________________
 125 void AliEMCAL::InitConstants()
 126 {
 127   //initialize EMCAL values
 128   fBirkC0 = 1;
 129   fBirkC1 = 0.013/1.032;
 130   fBirkC2 = 9.6e-6/(1.032 * 1.032);
 131
 132   fHighCharge        = 8.2 ;          // adjusted for a high gain range of 5.12 GeV (10 bits)
 133   fHighGain          = 6.64 ;
 134   fHighLowGainFactor = 16. ;          // adjusted for a low gain range of 82 GeV (10 bits)
 135   fLowGainOffset     = 1 ;            // offset added to the module id to distinguish high and low gain data
 136 }
 137
 138 //____________________________________________________________________________
 139 void AliEMCAL::DefineMediumParameters()
 140 {
 141   //
 142   // EMCAL cuts (Geant3)
 143   //
 144   Int_t * idtmed = fIdtmed->GetArray() - 1599 ;
 145 // --- Set decent energy thresholds for gamma and electron tracking
 146
 147   // Tracking threshold for photons and electrons in Lead
 148   Float_t cutgam=10.e-5; // 100 kev;
 149   Float_t cutele=10.e-5; // 100 kev;
 150   TString ntmp(GetTitle());
 151   ntmp.ToUpper();
 152   if(ntmp.Contains("10KEV")) {
 153     cutele = cutgam = 1.e-5;
 154   } else if(ntmp.Contains("50KEV")) {
 155     cutele = cutgam = 5.e-5;
 156   } else if(ntmp.Contains("100KEV")) {
 157     cutele = cutgam = 1.e-4;
 158   } else if(ntmp.Contains("200KEV")) {
 159     cutele = cutgam = 2.e-4;
 160   } else if(ntmp.Contains("500KEV")) {
 161     cutele = cutgam = 5.e-4;
 162   }
 163
 164   gMC->Gstpar(idtmed[1600],"CUTGAM", cutgam);
 165   gMC->Gstpar(idtmed[1600],"CUTELE", cutele); // 1MEV -> 0.1MEV; 15-aug-05
 166   gMC->Gstpar(idtmed[1600],"BCUTE",  cutgam);  // BCUTE and BCUTM start from GUTGUM
 167   gMC->Gstpar(idtmed[1600],"BCUTM",  cutgam);  // BCUTE and BCUTM start from GUTGUM
 168   // --- Generate explicitly delta rays in Lead ---
 169   gMC->Gstpar(idtmed[1600], "LOSS", 3) ;
 170   gMC->Gstpar(idtmed[1600], "DRAY", 1) ;
 171   gMC->Gstpar(idtmed[1600], "DCUTE", cutele) ;
 172   gMC->Gstpar(idtmed[1600], "DCUTM", cutele) ;
 173
 174 // --- in aluminium parts ---
 175   gMC->Gstpar(idtmed[1602],"CUTGAM", cutgam) ;
 176   gMC->Gstpar(idtmed[1602],"CUTELE", cutele) ;
 177   gMC->Gstpar(idtmed[1602],"BCUTE",  cutgam);  // BCUTE and BCUTM start from GUTGUM
 178   gMC->Gstpar(idtmed[1602],"BCUTM",  cutgam);  // BCUTE and BCUTM start from GUTGUM
 179   gMC->Gstpar(idtmed[1602], "LOSS",3.) ;
 180   gMC->Gstpar(idtmed[1602], "DRAY",1.) ;
 181   gMC->Gstpar(idtmed[1602], "DCUTE", cutele) ;
 182   gMC->Gstpar(idtmed[1602], "DCUTM", cutele) ;
 183
 184 // --- and finally thresholds for photons and electrons in the scintillator ---
 185   gMC->Gstpar(idtmed[1601],"CUTGAM", cutgam) ;
 186   gMC->Gstpar(idtmed[1601],"CUTELE", cutele) ;// 1MEV -> 0.1MEV; 15-aug-05
 187   gMC->Gstpar(idtmed[1601],"BCUTE",  cutgam);  // BCUTE and BCUTM start from GUTGUM
 188   gMC->Gstpar(idtmed[1601],"BCUTM",  cutgam);  // BCUTE and BCUTM start from GUTGUM
 189   gMC->Gstpar(idtmed[1601], "LOSS",3) ; // generate delta rays
 190   gMC->Gstpar(idtmed[1601], "DRAY",1) ;
 191   gMC->Gstpar(idtmed[1601], "DCUTE", cutele) ;
 192   gMC->Gstpar(idtmed[1601], "DCUTM", cutele) ;
 193
 194   // S steel -
 195   gMC->Gstpar(idtmed[1603],"CUTGAM", cutgam);
 196   gMC->Gstpar(idtmed[1603],"CUTELE", cutele);
 197   gMC->Gstpar(idtmed[1603],"BCUTE",  cutgam);  // BCUTE and BCUTM start from GUTGUM
 198   gMC->Gstpar(idtmed[1603],"BCUTM",  cutgam);  // BCUTE and BCUTM start from GUTGUM
 199   // --- Generate explicitly delta rays
 200   gMC->Gstpar(idtmed[1603], "LOSS",3);
 201   gMC->Gstpar(idtmed[1603], "DRAY",1);
 202   gMC->Gstpar(idtmed[1603], "DCUTE", cutele) ;
 203   gMC->Gstpar(idtmed[1603], "DCUTM", cutele) ;
 204
 205   AliEMCALGeometry* geom = GetGeometry();
 206   if(geom->GetILOSS()>=0) {
 207     for(int i=1600; i<=1603; i++) gMC->Gstpar(idtmed[i], "LOSS", geom->GetILOSS()) ;
 208   }
 209   if(geom->GetIHADR()>=0) {
 210     for(int i=1600; i<=1603; i++) gMC->Gstpar(idtmed[i], "HADR", geom->GetIHADR()) ;
 211   }
 212 }
 213
 214 //____________________________________________________________________________
 215 AliDigitizer* AliEMCAL::CreateDigitizer(AliRunDigitizer* manager) const
 216 {
 217   //create and return the digitizer
 218   return new AliEMCALDigitizer(manager);
 219 }
 220
 221 //____________________________________________________________________________
 222 void AliEMCAL::CreateMaterials()
 223 {
 224   // Definitions of materials to build EMCAL and associated tracking media.
 225   // media number in idtmed are 1599 to 1698.
 226   // --- Air ---
 227   Float_t aAir[4]={12.0107,14.0067,15.9994,39.948};
 228   Float_t zAir[4]={6.,7.,8.,18.};
 229   Float_t wAir[4]={0.000124,0.755267,0.231781,0.012827};
 230   Float_t dAir = 1.20479E-3;
 231   AliMixture(0, "Air$", aAir, zAir, dAir, 4, wAir) ;
 232
 233   // --- Lead ---
 234   AliMaterial(1, "Pb$", 207.2, 82, 11.35, 0.56, 0., 0, 0) ;
 235
 236
 237   // --- The polysterene scintillator (CH) ---
 238   Float_t aP[2] = {12.011, 1.00794} ;
 239   Float_t zP[2] = {6.0, 1.0} ;
 240   Float_t wP[2] = {1.0, 1.0} ;
 241   Float_t dP = 1.032 ;
 242
 243   AliMixture(2, "Polystyrene$", aP, zP, dP, -2, wP) ;
 244
 245   // --- Aluminium ---
 246   AliMaterial(3, "Al$", 26.98, 13., 2.7, 8.9, 999., 0, 0) ;
 247   // ---         Absorption length is ignored ^
 248
 249   // 25-aug-04 by PAI - see  PMD/AliPMDv0.cxx for STEEL definition
 250   Float_t asteel[4] = { 55.847,51.9961,58.6934,28.0855 };
 251   Float_t zsteel[4] = { 26.,24.,28.,14. };
 252   Float_t wsteel[4] = { .715,.18,.1,.005 };
 253   AliMixture(4, "STAINLESS STEEL$", asteel, zsteel, 7.88, 4, wsteel);
 254
 255   // DEFINITION OF THE TRACKING MEDIA
 256
 257   // for EMCAL: idtmed[1599->1698] equivalent to fIdtmed[0->100]
 258   Int_t   isxfld = gAlice->Field()->Integ() ;
 259   Float_t sxmgmx = gAlice->Field()->Max() ;
 260
 261   // Air                                                                         -> idtmed[1599]
 262  AliMedium(0, "Air$", 0, 0,
 263              isxfld, sxmgmx, 10.0, 1.0, 0.1, 0.1, 10.0, 0, 0) ;
 264
 265   // The Lead                                                                      -> idtmed[1600]
 266
 267   AliMedium(1, "Lead$", 1, 0,
 268              isxfld, sxmgmx, 10.0, 0.1, 0.1, 0.1, 0.1, 0, 0) ;
 269
 270  // The scintillator of the CPV made of Polystyrene scintillator                   -> idtmed[1601]
 271   float deemax = 0.1; // maximum fractional energy loss in one step (0 < DEEMAX â?¤ 1);i
 272   AliMedium(2, "Scintillator$", 2, 1,
 273             isxfld, sxmgmx, 10.0, 0.001, deemax, 0.001, 0.001, 0, 0) ;
 274
 275   // Various Aluminium parts made of Al                                            -> idtmed[1602]
 276   AliMedium(3, "Al$", 3, 0,
 277              isxfld, sxmgmx, 10.0, 0.1, 0.1, 0.001, 0.001, 0, 0) ;
 278
 279   // 25-aug-04 by PAI : see  PMD/AliPMDv0.cxx for STEEL definition                 -> idtmed[1603]
 280   AliMedium(4, "S steel$", 4, 0,
 281              isxfld, sxmgmx, 10.0, 0.1, 0.1, 0.001, 0.001, 0, 0) ;
 282
 283
 284   //set constants for Birk's Law implentation
 285   fBirkC0 =  1;
 286   fBirkC1 =  0.013/dP;
 287   fBirkC2 =  9.6e-6/(dP * dP);
 288
 289   // Call just in case of Geant3; What to do in case of Geant4 ?
 290   if(gMC->InheritsFrom("TGeant3")) DefineMediumParameters(); // Feb 20, 2007
 291 }
 292  //____________________________________________________________________________
 293 void AliEMCAL::Digits2Raw()
 294 {
 295   // convert digits of the current event to raw data
 296
 297   AliEMCALLoader * loader = dynamic_cast<AliEMCALLoader*>(fLoader) ;
 298
 299   // get the digits
 300   loader->LoadDigits("EMCAL");
 301   loader->GetEvent();
 302   TClonesArray* digits = loader->Digits() ;
 303
 304   if (!digits) {
 305     Error("Digits2Raw", "no digits found !");
 306     return;
 307   }
 308
 309   // get the digitizer
 310   loader->LoadDigitizer();
 311   AliEMCALDigitizer * digitizer = dynamic_cast<AliEMCALDigitizer *>(loader->Digitizer())  ;
 312
 313   // get the geometry
 314   AliEMCALGeometry* geom = GetGeometry();
 315   if (!geom) {
 316     AliError(Form("No geometry found !"));
 317     return;
 318   }
 319
 320   AliAltroBuffer* buffer = NULL;
 321   Int_t prevDDL = -1;
 322   Int_t adcValuesLow[fgkTimeBins];
 323   Int_t adcValuesHigh[fgkTimeBins];
 324
 325   //Load Mapping RCU files once
 326   TString path = gSystem->Getenv("ALICE_ROOT");
 327   path += "/EMCAL/mapping/RCU";
 328   TString path0 = path+"0.data";//This file will change in future
 329   TString path1 = path+"1.data";//This file will change in future
 330   AliAltroMapping * mapping[2] ; // For the moment only 2
 331   mapping[0] = new AliCaloAltroMapping(path0.Data());
 332   mapping[1] = new AliCaloAltroMapping(path1.Data());
 333
 334   // loop over digits (assume ordered digits)
 335   for (Int_t iDigit = 0; iDigit < digits->GetEntries(); iDigit++) {
 336     AliEMCALDigit* digit = dynamic_cast<AliEMCALDigit *>(digits->At(iDigit)) ;
 337     if (digit->GetAmp() < fgThreshold)
 338       continue;
 339
 340     //get cell indeces
 341     Int_t nSM = 0;
 342     Int_t nIphi = 0;
 343     Int_t nIeta = 0;
 344     Int_t iphi = 0;
 345     Int_t ieta = 0;
 346     Int_t nModule = 0;
 347     geom->GetCellIndex(digit->GetId(), nSM, nModule, nIphi, nIeta);
 348     geom->GetCellPhiEtaIndexInSModule(nSM, nModule, nIphi, nIeta,iphi, ieta) ;
 349
 350     //Check which is the RCU of the cell.
 351     Int_t iRCU = -111;
 352     //RCU0
 353     if (0<=iphi&&iphi<8) iRCU=0; // first cable row
 354     else if (8<=iphi&&iphi<16 && 0<=ieta&&ieta<24) iRCU=0; // first half;
 355     //second cable row
 356     //RCU1
 357     else if(8<=iphi&&iphi<16 && 24<=ieta&&ieta<48) iRCU=1; // second half;
 358     //second cable row
 359     else if(16<=iphi&&iphi<24) iRCU=1; // third cable row
 360
 361     //Which DDL?
 362     Int_t iDDL = fgDDLPerSuperModule* nSM + iRCU;
 363
 364     // new DDL
 365     if (iDDL != prevDDL) {
 366       // write real header and close previous file
 367
 368       if (buffer) {
 369         buffer->Flush();
 370         buffer->WriteDataHeader(kFALSE, kFALSE);
 371         delete buffer;
 372       }
 373
 374       // open new file and write dummy header
 375       TString fileName = AliDAQ::DdlFileName("EMCAL",iDDL);
 376       buffer = new AliAltroBuffer(fileName.Data(),mapping[iRCU]);
 377       buffer->WriteDataHeader(kTRUE, kFALSE);  //Dummy;
 378       prevDDL = iDDL;
 379     }
 380
 381     // out of time range signal (?)
 382     if (digit->GetTimeR() > GetRawFormatTimeMax() ) {
 383       AliInfo("Signal is out of time range.\n");
 384       buffer->FillBuffer((Int_t)digit->GetAmp());
 385       buffer->FillBuffer(GetRawFormatTimeBins() );  // time bin
 386       buffer->FillBuffer(3);          // bunch length
 387       buffer->WriteTrailer(3, ieta, iphi, nSM);  // trailer
 388       // calculate the time response function
 389     } else {
 390
 391       Double_t energy = digit->GetAmp() * digitizer->GetECAchannel() + digitizer->GetECApedestal() ;
 392
 393       Bool_t lowgain = RawSampledResponse(digit->GetTimeR(), energy, adcValuesHigh, adcValuesLow) ;
 394
 395       if (lowgain)
 396         buffer->WriteChannel(ieta, iphi, 0, GetRawFormatTimeBins(), adcValuesLow, fgThreshold);
 397       else
 398         buffer->WriteChannel(ieta,iphi, 1, GetRawFormatTimeBins(), adcValuesHigh, fgThreshold);
 399
 400     }
 401   }
 402
 403   // write real header and close last file
 404   if (buffer) {
 405     buffer->Flush();
 406     buffer->WriteDataHeader(kFALSE, kFALSE);
 407     delete buffer;
 408   }
 409   mapping[0]->Delete();
 410   mapping[1]->Delete();
 411   loader->UnloadDigits();
 412 }
 413
 414 //____________________________________________________________________________
 415 void AliEMCAL::Raw2Digits(AliRawReader* reader)
 416 {
 417   // convert raw data of the current event to digits
 418   AliEMCALGeometry * geom = GetGeometry();
 419   AliEMCALLoader * loader = dynamic_cast<AliEMCALLoader*>(fLoader) ;
 420
 421   // get the digits
 422   loader->CleanDigits(); // start from scratch
 423   loader->LoadDigits("EMCAL");
 424   TClonesArray* digits = loader->Digits() ;
 425   digits->Clear(); // yes, this is perhaps somewhat paranoid.. [clearing an extra time]
 426
 427   if (!digits) {
 428     Error("Raw2Digits", "no digits found !");
 429     return;
 430   }
 431   if (!reader) {
 432     Error("Raw2Digits", "no raw reader found !");
 433     return;
 434   }
 435
 436   // and get the digitizer too
 437   loader->LoadDigitizer();
 438   AliEMCALDigitizer * digitizer = dynamic_cast<AliEMCALDigitizer *>(loader->Digitizer())  ;
 439
 440   // Use AliAltroRawStream to read the ALTRO format.  No need to
 441   // reinvent the wheel :-)
 442   AliCaloRawStream in(reader,"EMCAL");
 443   // Select EMCAL DDL's;
 444   reader->Select("EMCAL");
 445
 446   // reading is from previously existing AliEMCALGetter.cxx
 447   // ReadRaw method
 448   Bool_t first = kTRUE ;
 449
 450   TF1 * signalF = new TF1("signal", RawResponseFunction, 0, GetRawFormatTimeMax(), 4);
 451   signalF->SetParNames("Charge", "Gain", "Amplitude", "TimeZero");
 452
 453   Bool_t lowGainFlag = kFALSE ;
 454   Int_t id =  -1;
 455   Int_t idigit = 0 ;
 456   Int_t amp = 0 ;
 457   Double_t time = 0. ;
 458   Double_t energy = 0. ;
 459
 460   TGraph * gLowGain = new TGraph(GetRawFormatTimeBins()) ;
 461   TGraph * gHighGain= new TGraph(GetRawFormatTimeBins()) ;
 462
 463   while ( in.Next() ) { // EMCAL entries loop
 464     if ( in.IsNewRow() ) {//phi
 465       if ( in.IsNewColumn() ) {//eta
 466         id =  geom->GetAbsCellIdFromCellIndexes(in.GetModule(), in.GetRow(), in.GetColumn()) ;
 467         if (!first) {
 468           FitRaw(lowGainFlag, gLowGain, gHighGain, signalF, energy, time) ;
 469
 470           if (time == 0. && energy == 0.) {
 471             amp = 0 ;
 472           }
 473           else {
 474             amp = static_cast<Int_t>( (energy - digitizer->GetECApedestal()) / digitizer->GetECAchannel() + 0.5 ) ;
 475           }
 476
 477           if (amp > 0) {
 478             new((*digits)[idigit]) AliEMCALDigit( -1, -1, id, amp, time) ;
 479             idigit++ ;
 480           }
 481           Int_t index ;
 482           for (index = 0; index < GetRawFormatTimeBins(); index++) {
 483             gLowGain->SetPoint(index, index * GetRawFormatTimeMax() / GetRawFormatTimeBins(), 0) ;
 484             gHighGain->SetPoint(index, index * GetRawFormatTimeMax() / GetRawFormatTimeBins(), 0) ;
 485           }
 486         } // not first
 487         first = kFALSE ;
 488         id =  geom->GetAbsCellIdFromCellIndexes(in.GetModule(), in.GetRow(), in.GetColumn()) ;
 489         if (in.GetModule() == GetRawFormatLowGainOffset() ) {
 490           lowGainFlag = kTRUE ;
 491         }
 492         else {
 493           lowGainFlag = kFALSE ;
 494         }
 495       } //new column-eta
 496     }// new row-phi
 497     if (lowGainFlag) {
 498       gLowGain->SetPoint(in.GetTime(),
 499                          in.GetTime()* GetRawFormatTimeMax() / GetRawFormatTimeBins(),
 500                          in.GetSignal()) ;
 501     }
 502     else {
 503       gHighGain->SetPoint(in.GetTime(),
 504                           in.GetTime() * GetRawFormatTimeMax() / GetRawFormatTimeBins(),
 505                           in.GetSignal() ) ;
 506     }
 507   } // EMCAL entries loop
 508   digits->Sort() ;
 509
 510   delete signalF ;
 511   delete gLowGain;
 512   delete gHighGain ;
 513
 514   return ;
 515 }
 516
 517 //____________________________________________________________________________
 518 void AliEMCAL::FitRaw(Bool_t lowGainFlag, TGraph * gLowGain, TGraph * gHighGain, TF1* signalF, Double_t & energy, Double_t & time)
 519 {
 520   // Fits the raw signal time distribution; from AliEMCALGetter
 521
 522   const Int_t kNoiseThreshold = 0 ;
 523   Double_t timezero1 = 0., timezero2 = 0., timemax = 0. ;
 524   Double_t signal = 0., signalmax = 0. ;
 525   energy = time = 0. ;
 526
 527   if (lowGainFlag) {
 528     timezero1 = timezero2 = signalmax = timemax = 0. ;
 529     signalF->FixParameter(0, GetRawFormatLowCharge()) ;
 530     signalF->FixParameter(1, GetRawFormatLowGain()) ;
 531     Int_t index ;
 532     for (index = 0; index < GetRawFormatTimeBins(); index++) {
 533       gLowGain->GetPoint(index, time, signal) ;
 534       if (signal > kNoiseThreshold && timezero1 == 0.)
 535         timezero1 = time ;
 536       if (signal <= kNoiseThreshold && timezero1 > 0. && timezero2 == 0.)
 537         timezero2 = time ;
 538       if (signal > signalmax) {
 539         signalmax = signal ;
 540         timemax   = time ;
 541       }
 542     }
 543     signalmax /= RawResponseFunctionMax(GetRawFormatLowCharge(),
 544                                                 GetRawFormatLowGain()) ;
 545     if ( timezero1 + GetRawFormatTimePeak() < GetRawFormatTimeMax() * 0.4 ) { // else its noise
 546       signalF->SetParameter(2, signalmax) ;
 547       signalF->SetParameter(3, timezero1) ;
 548       gLowGain->Fit(signalF, "QRON", "", 0., timezero2); //, "QRON") ;
 549       energy = signalF->GetParameter(2) ;
 550       time   = signalF->GetMaximumX() - GetRawFormatTimePeak() - GetRawFormatTimeTrigger() ;
 551     }
 552   } else {
 553     timezero1 = timezero2 = signalmax = timemax = 0. ;
 554     signalF->FixParameter(0, GetRawFormatHighCharge()) ;
 555     signalF->FixParameter(1, GetRawFormatHighGain()) ;
 556     Int_t index ;
 557     for (index = 0; index < GetRawFormatTimeBins(); index++) {
 558       gHighGain->GetPoint(index, time, signal) ;
 559       if (signal > kNoiseThreshold && timezero1 == 0.)
 560         timezero1 = time ;
 561       if (signal <= kNoiseThreshold && timezero1 > 0. && timezero2 == 0.)
 562         timezero2 = time ;
 563       if (signal > signalmax) {
 564         signalmax = signal ;
 565         timemax   = time ;
 566       }
 567     }
 568     signalmax /= RawResponseFunctionMax(GetRawFormatHighCharge(),
 569                                                 GetRawFormatHighGain()) ;;
 570     if ( timezero1 + GetRawFormatTimePeak() < GetRawFormatTimeMax() * 0.4 ) { // else its noise
 571       signalF->SetParameter(2, signalmax) ;
 572       signalF->SetParameter(3, timezero1) ;
 573       gHighGain->Fit(signalF, "QRON", "", 0., timezero2) ;
 574       energy = signalF->GetParameter(2) ;
 575       time   = signalF->GetMaximumX() - GetRawFormatTimePeak() - GetRawFormatTimeTrigger() ;
 576     }
 577   }
 578
 579   return;
 580 }
 581
 582 //____________________________________________________________________________
 583 void AliEMCAL::Hits2SDigits()
 584 {
 585 // create summable digits
 586
 587   GetGeometry();
 588   AliEMCALSDigitizer emcalDigitizer(fLoader->GetRunLoader()->GetFileName().Data()) ;
 589   emcalDigitizer.SetEventRange(0, -1) ; // do all the events
 590   emcalDigitizer.ExecuteTask() ;
 591 }
 592
 593 //____________________________________________________________________________
 594
 595 AliLoader* AliEMCAL::MakeLoader(const char* topfoldername)
 596 {
 597 //different behaviour than standard (singleton getter)
 598 // --> to be discussed and made eventually coherent
 599  fLoader = new AliEMCALLoader(GetName(),topfoldername);
 600  return fLoader;
 601 }
 602
 603 //__________________________________________________________________
 604 Double_t AliEMCAL::RawResponseFunction(Double_t *x, Double_t *par)
 605 {
 606   // Shape of the electronics raw reponse:
 607   // It is a semi-gaussian, 2nd order Gamma function of the general form
 608   // v(t) = n**n * Q * A**n / C *(t/tp)**n * exp(-n * t/tp) with
 609   // tp : peaking time par[0]
 610   // n  : order of the function
 611   // C  : integrating capacitor in the preamplifier
 612   // A  : open loop gain of the preamplifier
 613   // Q  : the total APD charge to be measured Q = C * energy
 614
 615   Double_t signal ;
 616   Double_t xx = x[0] - ( fgTimeTrigger + par[3] ) ;
 617
 618   if (xx < 0 || xx > fgTimeMax)
 619     signal = 0. ;
 620   else {
 621     Double_t fac = par[0] * TMath::Power(fgOrder, fgOrder) * TMath::Power(par[1], fgOrder) / fgCapa ;
 622     signal = fac * par[2] * TMath::Power(xx / fgTimePeak, fgOrder) * TMath::Exp(-fgOrder * (xx / fgTimePeak)) ;
 623   }
 624   return signal ;
 625 }
 626
 627 //__________________________________________________________________
 628 Double_t AliEMCAL::RawResponseFunctionMax(Double_t charge, Double_t gain)
 629 {
 630   //compute the maximum of the raw response function and return
 631   return ( charge * TMath::Power(fgOrder, fgOrder) * TMath::Power(gain, fgOrder)
 632      / ( fgCapa * TMath::Exp(fgOrder) ) );
 633
 634 }
 635 //__________________________________________________________________
 636 Bool_t AliEMCAL::RawSampledResponse(
 637 const Double_t dtime, const Double_t damp, Int_t * adcH, Int_t * adcL) const
 638 {
 639   // for a start time dtime and an amplitude damp given by digit,
 640   // calculates the raw sampled response AliEMCAL::RawResponseFunction
 641
 642   const Int_t kRawSignalOverflow = 0x3FF ;
 643   Bool_t lowGain = kFALSE ;
 644
 645   TF1 signalF("signal", RawResponseFunction, 0, GetRawFormatTimeMax(), 4);
 646
 647   for (Int_t iTime = 0; iTime < GetRawFormatTimeBins(); iTime++) {
 648     signalF.SetParameter(0, GetRawFormatHighCharge() ) ;
 649     signalF.SetParameter(1, GetRawFormatHighGain() ) ;
 650     signalF.SetParameter(2, damp) ;
 651     signalF.SetParameter(3, dtime) ;
 652     Double_t time = iTime * GetRawFormatTimeMax() / GetRawFormatTimeBins() ;
 653     Double_t signal = signalF.Eval(time) ;
 654     if ( static_cast<Int_t>(signal+0.5) > kRawSignalOverflow ){  // larger than 10 bits
 655       signal = kRawSignalOverflow ;
 656       lowGain = kTRUE ;
 657     }
 658     adcH[iTime] =  static_cast<Int_t>(signal + 0.5) ;
 659
 660     signalF.SetParameter(0, GetRawFormatLowCharge() ) ;
 661     signalF.SetParameter(1, GetRawFormatLowGain() ) ;
 662     signal = signalF.Eval(time) ;
 663     if ( static_cast<Int_t>(signal+0.5) > kRawSignalOverflow)  // larger than 10 bits
 664       signal = kRawSignalOverflow ;
 665     adcL[iTime] = static_cast<Int_t>(0.5 + signal ) ;
 666
 667   }
 668   return lowGain ;
 669 }