]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - ZDC/AliZDCDigitizer.cxx
do init of QADatamaker at one place only
[u/mrichter/AliRoot.git] / ZDC / AliZDCDigitizer.cxx
index 9f1ec32b81f6d3167318be586940d6bd7fbfd7f6..a473206546734bb3674af2fbfc5eedcd7cdc90eb 100644 (file)
 #include "AliZDCDigitizer.h"
 
 class AliCDBStorage;
-class AliZDCCalibData;
+class AliZDCPedestals;
+class AliZDCCalib;
+class AliZDCRecParam;
 
 ClassImp(AliZDCDigitizer)
 
 
 //____________________________________________________________________________
-AliZDCDigitizer::AliZDCDigitizer()
+AliZDCDigitizer::AliZDCDigitizer() :
+  fIsCalibration(0), 
+  fPedData(0), 
+  fCalibData(0),
+  fRecParam(0)
 {
 // Default constructor    
 
@@ -58,11 +64,14 @@ AliZDCDigitizer::AliZDCDigitizer()
 
 //____________________________________________________________________________
 AliZDCDigitizer::AliZDCDigitizer(AliRunDigitizer* manager):
-  AliDigitizer(manager)
+  AliDigitizer(manager),
+  fIsCalibration(0), //By default the simulation doesn't create calib. data
+  fPedData(GetPedData()), 
+  fCalibData(GetCalibData()),
+  fRecParam(GetRecParam())
 {
-
   // Get calibration data
-  fCalibData = GetCalibData(); 
+  if(fIsCalibration!=0) printf("\n\t AliZDCDigitizer -> Creating calibration data (pedestals)\n");
 
 }
 
@@ -70,10 +79,29 @@ AliZDCDigitizer::AliZDCDigitizer(AliRunDigitizer* manager):
 AliZDCDigitizer::~AliZDCDigitizer()
 {
 // Destructor
-
+// Not implemented
 }
 
 
+//____________________________________________________________________________
+AliZDCDigitizer::AliZDCDigitizer(const AliZDCDigitizer &digitizer):
+  AliDigitizer(),
+  fIsCalibration(digitizer.fIsCalibration),
+  fPedData(digitizer.fPedData),
+  fCalibData(digitizer.fCalibData),
+  fRecParam(digitizer.fRecParam)
+{
+  // Copy constructor
+
+  for(Int_t i=0; i<6; i++){
+     for(Int_t j=0; j<5; j++){
+        fPMGain[i][j]   = digitizer.fPMGain[i][j];           
+     }
+  }
+  for(Int_t i=0; i<2; i++) fADCRes[i] = digitizer.fADCRes[i];
+
+}
+
 //____________________________________________________________________________
 Bool_t AliZDCDigitizer::Init()
 {
@@ -99,7 +127,8 @@ void AliZDCDigitizer::Exec(Option_t* /*option*/)
 {
   // Execute digitization
 
-  Float_t pm[5][5]; // !!! 2nd ZDC set added (needed for trigger purposes!)
+  // ------------------------------------------------------------
+  // !!! 2nd ZDC set added 
   // *** 1st 3 arrays are digits from REAL (simulated) hits
   // *** last 2 are copied from simulated digits
   // --- pm[0][...] = light in ZN right  [C, Q1, Q2, Q3, Q4]
@@ -107,11 +136,22 @@ void AliZDCDigitizer::Exec(Option_t* /*option*/)
   // --- pm[2][...] = light in ZEM [x, 1, 2, x, x]
   // --- pm[3][...] = light in ZN left [C, Q1, Q2, Q3, Q4] ->NEW!
   // --- pm[4][...] = light in ZP left [C, Q1, Q2, Q3, Q4] ->NEW!
-  
-  for (Int_t iSector1=0; iSector1<5; iSector1++) 
-    for (Int_t iSector2=0; iSector2<5; iSector2++) {
+  // ------------------------------------------------------------
+  Float_t pm[5][5]; 
+  for(Int_t iSector1=0; iSector1<5; iSector1++) 
+    for(Int_t iSector2=0; iSector2<5; iSector2++){
       pm[iSector1][iSector2] = 0;
     }
+    
+  // ------------------------------------------------------------
+  // ### Out of time ADC added (22 channels)
+  // --- same codification as for signal PTMs (see above)
+  // ------------------------------------------------------------
+  Float_t pmoot[5][5];
+  for(Int_t iSector1=0; iSector1<5; iSector1++) 
+    for(Int_t iSector2=0; iSector2<5; iSector2++){
+      pmoot[iSector1][iSector2] = 0;
+    }
 
   // impact parameter and number of spectators
   Float_t impPar = -1;
@@ -119,28 +159,28 @@ void AliZDCDigitizer::Exec(Option_t* /*option*/)
   Int_t specP = 0;
 
   // loop over input streams
-  for (Int_t iInput = 0; iInput<fManager->GetNinputs(); iInput++) {
+  for(Int_t iInput = 0; iInput<fManager->GetNinputs(); iInput++){
 
     // get run loader and ZDC loader
     AliRunLoader* runLoader = 
       AliRunLoader::GetRunLoader(fManager->GetInputFolderName(iInput));
     AliLoader* loader = runLoader->GetLoader("ZDCLoader");
-    if (!loader) continue;
+    if(!loader) continue;
 
     // load sdigits
     loader->LoadSDigits();
     TTree* treeS = loader->TreeS();
-    if (!treeS) continue;
+    if(!treeS) continue;
     AliZDCSDigit sdigit;
     AliZDCSDigit* psdigit = &sdigit;
     treeS->SetBranchAddress("ZDC", &psdigit);
 
     // loop over sdigits
-    for (Int_t iSDigit=0; iSDigit<treeS->GetEntries(); iSDigit++) {
+    for(Int_t iSDigit=0; iSDigit<treeS->GetEntries(); iSDigit++){
       treeS->GetEntry(iSDigit);
       //
-      if (!psdigit) continue;
-      if ((sdigit.GetSector(1) < 0) || (sdigit.GetSector(1) > 4)) {
+      if(!psdigit) continue;
+      if((sdigit.GetSector(1) < 0) || (sdigit.GetSector(1) > 4)){
        AliError(Form("\nsector[0] = %d, sector[1] = %d\n", 
                       sdigit.GetSector(0), sdigit.GetSector(1)));
        continue;
@@ -150,34 +190,35 @@ void AliZDCDigitizer::Exec(Option_t* /*option*/)
       /*printf("\n\t Detector %d, Tower %d -> pm[%d][%d] = %.0f \n",
          sdigit.GetSector(0), sdigit.GetSector(1),sdigit.GetSector(0)-1,
          sdigit.GetSector(1), pm[sdigit.GetSector(0)-1][sdigit.GetSector(1)]); // Chiara debugging!
-         */
+      */
     }
 
     loader->UnloadSDigits();
 
     // get the impact parameter and the number of spectators in case of hijing
-    if (!runLoader->GetAliRun()) runLoader->LoadgAlice();
+    if(!runLoader->GetAliRun()) runLoader->LoadgAlice();
     AliHeader* header = runLoader->GetAliRun()->GetHeader();
-    if (!header) continue;
+    if(!header) continue;
     AliGenEventHeader* genHeader = header->GenEventHeader();
-    if (!genHeader) continue;
-    if (!genHeader->InheritsFrom(AliGenHijingEventHeader::Class())) continue;
+    if(!genHeader) continue;
+    if(!genHeader->InheritsFrom(AliGenHijingEventHeader::Class())) continue;
     impPar = ((AliGenHijingEventHeader*) genHeader)->ImpactParameter();
     // 
     specN = ((AliGenHijingEventHeader*) genHeader)->ProjSpectatorsn();
     specP = ((AliGenHijingEventHeader*) genHeader)->ProjSpectatorsp();
     AliDebug(2, Form("\n AliZDCDigitizer -> b = %f fm, Nspecn = %d, Nspecp = %d\n",
                      impPar, specN, specP));
-    printf("\n\t AliZDCDigitizer -> b = %f fm, Nspecn = %d, Nspecp = %d\n", impPar, specN, specP);
+    printf("\n\t AliZDCDigitizer -> b = %f fm, # generated spectator n = %d," 
+    " # generated spectator p = %d\n", impPar, specN, specP);
   }
 
   // add spectators
-  if (impPar >= 0) {
+  if(impPar >= 0) {
     Int_t freeSpecN, freeSpecP;
     Fragmentation(impPar, specN, specP, freeSpecN, freeSpecP);
-    printf("\n\t AliZDCDigitizer ---- Adding signal for %d free spectator n\n",freeSpecN);
+    printf("\n\t AliZDCDigitizer -> Adding signal for %d free spectator n\n",freeSpecN);
     SpectatorSignal(1, freeSpecN, pm);
-    printf("\t AliZDCDigitizer ---- Adding signal for %d free spectator p\n\n",freeSpecP);
+    printf("\t AliZDCDigitizer -> Adding signal for %d free spectator p\n\n",freeSpecP);
     SpectatorSignal(2, freeSpecP, pm);
   }
 
@@ -186,14 +227,14 @@ void AliZDCDigitizer::Exec(Option_t* /*option*/)
   AliRunLoader* runLoader = 
     AliRunLoader::GetRunLoader(fManager->GetOutputFolderName());
   AliLoader* loader = runLoader->GetLoader("ZDCLoader");
-  if (!loader) {
+  if(!loader) {
     AliError("no ZDC loader found");
     return;
   }
 
   // create the output tree
   const char* mode = "update";
-  if (runLoader->GetEventNumber() == 0) mode = "recreate";
+  if(runLoader->GetEventNumber() == 0) mode = "recreate";
   loader->LoadDigits(mode);
   loader->MakeTree("D");
   TTree* treeD = loader->TreeD();
@@ -203,36 +244,73 @@ void AliZDCDigitizer::Exec(Option_t* /*option*/)
   treeD->Branch("ZDC", "AliZDCDigit", &pdigit, kBufferSize);
 
   // Create digits
-  Int_t sector[2], sectorL[2];
-  Int_t digi[2], digiL[2];
-  for (sector[0]=1; sector[0]<=3; sector[0]++)
-    for (sector[1]=0; sector[1]<5; sector[1]++)  {
-        if ((sector[0]==3) && ((sector[1]<1) || (sector[1]>2))) continue;
-        for (Int_t res=0; res<2; res++){
-          digi[res] = Phe2ADCch(sector[0], sector[1], pm[sector[0]-1][sector[1]], res) 
+  Int_t sector[2];
+  Int_t digi[2], digioot[2];
+  for(sector[0]=1; sector[0]<6; sector[0]++){
+    for(sector[1]=0; sector[1]<5; sector[1]++){
+        if((sector[0]==3) && ((sector[1]<1) || (sector[1]>2))) continue;
+        for(Int_t res=0; res<2; res++){
+           digi[res] = Phe2ADCch(sector[0], sector[1], pm[sector[0]-1][sector[1]], res) 
                    + Pedestal(sector[0], sector[1], res);
        }
-       //printf("\t DIGIT added -> det = %d, quad = %d - digi[0,1] = [%d, %d]\n",
-       //      sector[0], sector[1], digi[0], digi[1]); // Chiara debugging!
-        new(pdigit) AliZDCDigit(sector, digi);
-        treeD->Fill();
+       /*printf("\t DIGIT added -> det %d quad %d - digi[0,1] = [%d, %d]\n",
+            sector[0], sector[1], digi[0], digi[1]); // Chiara debugging!
+       */
        //
-       // --- Adding digits for 2nd ZDC set (left side w.r.t. IP) ---
-       // --- they are copied from right ZDC digits
-       if(sector[0]==1 || sector[0]==2){
-          sectorL[0] = sector[0]+3;
-          sectorL[1] = sector[1];
-           for (Int_t res=0; res<2; res++){
-             digiL[res] = Phe2ADCch(sectorL[0], sectorL[1], pm[sector[0]-1][sector[1]], res) 
-                   + Pedestal(sectorL[0], sectorL[1], res);
-          }
-          //printf("\t DIGIT added -> det = %d, quad = %d - digi[0,1] = [%d, %d]\n",
-          //      sectorL[0], sectorL[1], digiL[0], digiL[1]); // Chiara debugging!
-         new(pdigit) AliZDCDigit(sectorL, digiL);
-       }
+       new(pdigit) AliZDCDigit(sector, digi);
+        treeD->Fill();
+    }
+  } // Loop over detector
+  // Adding in-time digits for 2 reference PTM signals (after signal ch.)
+  // (for the moment the ref. signal is completely invented assuming a PMgain of 5*10^4!)
+  Int_t sectorRef[2];
+  sectorRef[1] = 5;
+  Int_t sigRef[2];
+  // Reference signal are set to 100 (high gain chain) and 800 (low gain chain)
+  if(fIsCalibration==0) {sigRef[0]=100;  sigRef[1]=800;}
+  else {sigRef[0]=0;  sigRef[1]=0;} // calibration -> simulation of pedestal values
+  //
+  for(Int_t iref=0; iref<2; iref++){
+     sectorRef[0] = 3*iref+1;
+     for(Int_t res=0; res<2; res++){
+       sigRef[res] += Pedestal(sectorRef[0], sectorRef[1], res);
+     }
+     /*printf("\t RefDigit added -> det = %d, quad = %d - digi[0,1] = [%d, %d]\n",
+         sectorRef[0], sectorRef[1], sigRef[0], sigRef[1]); // Chiara debugging!
+     */
+     new(pdigit) AliZDCDigit(sectorRef, sigRef);
+     treeD->Fill();     
+  }
+  //
+  // --- Adding digits for out-of-time channels after signal digits
+  for(sector[0]=1; sector[0]<6; sector[0]++){
+    for(sector[1]=0; sector[1]<5; sector[1]++){
+        if((sector[0]==3) && ((sector[1]<1) || (sector[1]>2))) continue;
+        for(Int_t res=0; res<2; res++){
+           digioot[res] = Pedestal(sector[0], sector[1], res); // out-of-time ADCs
+       }
+       /*printf("\t DIGIToot added -> det = %d, quad = %d - digi[0,1] = [%d, %d]\n",
+            sector[0], sector[1], digioot[0], digioot[1]); // Chiara debugging!
+       */
        //
+       new(pdigit) AliZDCDigit(sector, digioot);
         treeD->Fill();
     }
+  }
+  // Adding out-of-time digits for 2 reference PTM signals (after out-of-time ch.)
+  Int_t sigRefoot[2];
+  for(Int_t iref=0; iref<2; iref++){
+     sectorRef[0] = 3*iref+1;
+     for(Int_t res=0; res<2; res++){
+       sigRefoot[res] = Pedestal(sectorRef[0], sectorRef[1], res);
+     }
+     /*printf("\t RefDigitoot added -> det = %d, quad = %d - digi[0,1] = [%d, %d]\n",
+         sectorRef[0], sectorRef[1], sigRefoot[0], sigRefoot[1]); // Chiara debugging!
+     */
+     new(pdigit) AliZDCDigit(sectorRef, sigRefoot);
+     treeD->Fill();
+  }
+  //printf("\t AliZDCDigitizer -> TreeD has %d entries\n",(Int_t) treeD->GetEntries());
 
   // write the output tree
   loader->WriteDigits("OVERWRITE");
@@ -248,7 +326,7 @@ void AliZDCDigitizer::Fragmentation(Float_t impPar, Int_t specN, Int_t specP,
 
   Int_t zz[100], nn[100];
   AliZDCFragment frag(impPar);
-  for (Int_t j=0; j<=99; j++){
+  for(Int_t j=0; j<=99; j++){
      zz[j] =0;
      nn[j] =0;
   }
@@ -263,6 +341,10 @@ void AliZDCDigitizer::Fragmentation(Float_t impPar, Int_t specN, Int_t specP,
   frag.AttachNeutrons(zz, nn, ztot, ntot);
   freeSpecN = specN-ntot-2*nAlpha;
   freeSpecP = specP-ztot-2*nAlpha;
+  // Removing deuterons
+  Int_t ndeu = (Int_t) (freeSpecN*frag.DeuteronNumber());
+  freeSpecN -= ndeu;
+  //
   if(freeSpecN<0) freeSpecN=0;
   if(freeSpecP<0) freeSpecP=0;
   AliDebug(2, Form("FreeSpn = %d, FreeSpp = %d", freeSpecN, freeSpecP));
@@ -275,12 +357,12 @@ void AliZDCDigitizer::SpectatorSignal(Int_t SpecType, Int_t numEvents,
 // add signal of the spectators
  
   TFile* file = NULL;
-  if (SpecType == 1) {         // --- Signal for spectator neutrons
-    file = TFile::Open("$ALICE/$ALICE_LEVEL/ZDC/ZNsignalntu.root");
-  } else if (SpecType == 2) {  // --- Signal for spectator protons
-    file = TFile::Open("$ALICE/$ALICE_LEVEL/ZDC/ZPsignalntu.root");
+  if(SpecType == 1) {          // --- Signal for spectator neutrons
+    file = TFile::Open("$ALICE_ROOT/ZDC/ZNsignalntu.root");
+  } else if(SpecType == 2) {   // --- Signal for spectator protons
+    file = TFile::Open("$ALICE_ROOT/ZDC/ZPsignalntu.root");
   }
-  if (!file || !file->IsOpen()) {
+  if(!file || !file->IsOpen()) {
     AliError("Opening of file failed");
     return;
   }
@@ -291,7 +373,7 @@ void AliZDCDigitizer::SpectatorSignal(Int_t SpecType, Int_t numEvents,
   Float_t *entry, hitsSpec[7];
   Int_t pl, i, j, k, iev=0, rnd[125], volume[2];
   for(pl=0;pl<125;pl++) rnd[pl] = 0;
-  if (numEvents > 125) {
+  if(numEvents > 125) {
     AliWarning(Form("numEvents (%d) is larger than 125", numEvents));
     numEvents = 125;
   }
@@ -316,15 +398,15 @@ void AliZDCDigitizer::SpectatorSignal(Int_t SpecType, Int_t numEvents,
                            volume[0], volume[1], lightQ, lightC));
          //printf("\n   Volume = (%d, %d), lightQ = %.0f, lightC = %.0f",
           //                 volume[0], volume[1], lightQ, lightC);
-         if (volume[0] < 3) {  // ZN or ZP
+         if(volume[0] != 3) {  // ZN or ZP
             pm[volume[0]-1][0] += lightC;
             pm[volume[0]-1][volume[1]] += lightQ;
            //printf("\n   pm[%d][0] = %.0f, pm[%d][%d] = %.0f\n",(volume[0]-1),pm[volume[0]-1][0],
            //  (volume[0]-1),volume[1],pm[volume[0]-1][volume[1]]);
          } 
          else { 
-            if (volume[1] == 1) pm[2][1] += lightC; // ZEM 1
-            else                pm[2][2] += lightQ; // ZEM 2
+            if(volume[1] == 1) pm[2][1] += lightC; // ZEM 1
+            else               pm[2][2] += lightQ; // ZEM 2
            //printf("\n   pm[2][1] = %.0f, pm[2][2] = %.0f\n",pm[2][1],pm[2][2]);
          }
        }
@@ -346,7 +428,7 @@ Int_t AliZDCDigitizer::Phe2ADCch(Int_t Det, Int_t Quad, Float_t Light,
 {
   // Evaluation of the ADC channel corresponding to the light yield Light
   Int_t vADCch = (Int_t) (Light * fPMGain[Det-1][Quad] * fADCRes[Res]);
-  //printf("\t Phe2ADCch -> det %d quad %d - phe %.0f  ADC %d\n", Det,Quad,Light,ADCch);
+  //printf("\t Phe2ADCch -> det %d quad %d - phe %.0f  ADC %d\n", Det,Quad,Light,vADCch);
 
   return vADCch;
 }
@@ -356,19 +438,34 @@ Int_t AliZDCDigitizer::Pedestal(Int_t Det, Int_t Quad, Int_t Res) const
 {
   // Returns a pedestal for detector det, PM quad, channel with res.
   //
-  Float_t meanPed;
-  if(Det != 3) meanPed = fCalibData->GetMeanPed(10*(Det-1)+Quad+5*Res);
-  else         meanPed = fCalibData->GetMeanPed(10*(Det-1)+Quad+1*Res);
-  
-  //printf("\t Pedestal -> det = %d, quad = %d, res = %d - Ped[%d] = %d\n",
-  //   Det, Quad, Res,10*(Det-1)+Quad+5*Res,(Int_t) meanPed); // Chiara debugging!
-  
+  Float_t pedValue;
+  // Normal run
+  if(fIsCalibration == 0){
+    Int_t index=0, kNch=24;
+    if(Quad!=5){
+      if(Det==1)       index = Quad+kNch*Res;    // ZNC
+      else if(Det==2)  index = (Quad+5)+kNch*Res;  // ZPC
+      else if(Det==3)  index = (Quad+9)+kNch*Res;  // ZEM
+      else if(Det==4)  index = (Quad+12)+kNch*Res; // ZNA
+      else if(Det==5)  index = (Quad+17)+kNch*Res; // ZPA
+    }
+    else index = (Det-1)/3+22+kNch*Res; // Reference PMs
+    //
+    Float_t meanPed = fPedData->GetMeanPed(index);
+    Float_t pedWidth = fPedData->GetMeanPedWidth(index);
+    pedValue = gRandom->Gaus(meanPed,pedWidth);
+    //
+    /*printf("\t  AliZDCDigitizer::Pedestal -> det %d quad %d res %d - Ped[%d] = %d\n",
+       Det, Quad, Res, index,(Int_t) pedValue); // Chiara debugging!
+    */
+  }
   // To create calibration object
-  /*Float_t meanPed;
-  meanPed = gRandom->Gaus((40.+10.*gRandom->Rndm()),5.);
-  */
+  else{
+    if(Res == 0) pedValue = gRandom->Gaus((35.+10.*gRandom->Rndm()),(0.5+0.2*gRandom->Rndm())); //High gain
+    else  pedValue = gRandom->Gaus((250.+100.*gRandom->Rndm()),(3.5+2.*gRandom->Rndm())); //Low gain
+  }
 
-  return (Int_t) meanPed;
+  return (Int_t) pedValue;
 }
 
 //_____________________________________________________________________________
@@ -397,15 +494,46 @@ AliCDBStorage* AliZDCDigitizer::SetStorage(const char *uri)
 }
 
 //_____________________________________________________________________________
-AliZDCCalibData* AliZDCDigitizer::GetCalibData() const
+AliZDCPedestals* AliZDCDigitizer::GetPedData() const
+{
+
+  // Getting pedestal calibration object for ZDC set
+
+  AliCDBEntry  *entry = AliCDBManager::Instance()->Get("ZDC/Calib/Pedestals");
+  if(!entry) AliFatal("No calibration data loaded!");  
+
+  AliZDCPedestals *calibdata = dynamic_cast<AliZDCPedestals*>  (entry->GetObject());
+  if(!calibdata)  AliFatal("Wrong calibration object in calibration  file!");
+
+  return calibdata;
+}
+
+//_____________________________________________________________________________
+AliZDCCalib* AliZDCDigitizer::GetCalibData() const
+{
+
+  // Getting calibration object for ZDC set
+
+  AliCDBEntry  *entry = AliCDBManager::Instance()->Get("ZDC/Calib/Calib");
+  if(!entry) AliFatal("No calibration data loaded!");  
+
+  AliZDCCalib *calibdata = dynamic_cast<AliZDCCalib*>  (entry->GetObject());
+  if(!calibdata)  AliFatal("Wrong calibration object in calibration  file!");
+
+  return calibdata;
+}
+
+//_____________________________________________________________________________
+AliZDCRecParam* AliZDCDigitizer::GetRecParam() const
 {
 
   // Getting calibration object for ZDC set
 
-  AliCDBEntry  *entry = AliCDBManager::Instance()->Get("ZDC/Calib/Data");
-  AliZDCCalibData *calibdata = (AliZDCCalibData*) entry->GetObject();
+  AliCDBEntry  *entry = AliCDBManager::Instance()->Get("ZDC/Calib/RecParam");
+  if(!entry) AliFatal("No calibration data loaded!");  
 
-  if (!calibdata)  AliWarning("No calibration data from calibration database !");
+  AliZDCRecParam *calibdata = dynamic_cast<AliZDCRecParam*>  (entry->GetObject());
+  if(!calibdata)  AliFatal("Wrong calibration object in calibration  file!");
 
   return calibdata;
 }