ZDC/AliZDCMerger.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
  17 ///////////////////////////////////////////////////////////////////////////////
  18 //                                                                           //
  19 //                      ZDC event merging class                              //
  20 //                                                                           //
  21 ///////////////////////////////////////////////////////////////////////////////
  22
  23 // --- ROOT system
  24 #include <iostream.h>
  25 #include <TTree.h>
  26 #include <TFile.h>
  27 #include <TDirectory.h>
  28 #include <TNtuple.h>
  29 #include <TSystem.h>
  30
  31 // --- AliRoot header files
  32 #include "AliZDCMerger.h"
  33 #include "AliZDC.h"
  34 #include "AliZDCHit.h"
  35 #include "AliZDCMergedHit.h"
  36 #include "AliZDCDigit.h"
  37 #include "AliZDCFragment.h"
  38 #include "AliRun.h"
  39 #include "AliDetector.h"
  40 #include "AliHeader.h"
  41 #include "AliGenEventHeader.h"
  42 #include "AliGenHijingEventHeader.h"
  43
  44 ClassImp(AliZDCMerger)
  45
  46 //____________________________________________________________________________
  47 AliZDCMerger::AliZDCMerger()
  48 {
  49 // Default constructor
  50     fMerge       = kDigitize;
  51 //   fMerge       = kMerge;
  52     fFnBgr       = 0;
  53     fBgrFile     = 0;
  54     fNEvBgr      = 0;
  55     fTrHBgr      = 0;
  56     fTrSDBgr     = 0;
  57     fImpPar      = 0;
  58     fSpecn       = 0;
  59     fSpecp       = 0;
  60     fFreeSpn     = 0;
  61     fFreeSpp     = 0;
  62     fFnSpecn     = 0;
  63     fSpecnFile   = 0;
  64     fFnSpecp     = 0;
  65     fSpecpFile   = 0;
  66     fNMhits      = 0;
  67
  68 }
  69
  70 //____________________________________________________________________________
  71 AliZDCMerger::~AliZDCMerger()
  72 {
  73 // Destructor
  74     if (fBgrFile)    delete fBgrFile;
  75     if (fTrHBgr)     delete fTrHBgr;
  76     if (fTrSDBgr)    delete fTrSDBgr;
  77     if (fHitsBgr)    delete fHitsBgr;
  78     if (fSpecnFile)  delete fSpecnFile;
  79     if (fSpecpFile)  delete fSpecpFile;
  80 }
  81
  82 //____________________________________________________________________________
  83 void AliZDCMerger::InitMerging()
  84 {
  85     // Hits tree, impact parameter, num. of spectators n & p
  86     //          in background (full Hijing) event
  87     Float_t b;
  88     Int_t nspecn, nspecp;
  89     Background(b, nspecn, nspecp);
  90
  91     // Production of nuclear fragments -> num. of FREE spectators n & p
  92     Fragmentation(b, nspecn, nspecp, fFreeSpn, fFreeSpp);
  93
  94     // Extract from spectators distribution the signal events:
  95     // NFreeSpectatorN spectator n & NFreeSpectatorP spectator p
  96     Mixing();
  97 }
  98
  99 //____________________________________________________________________________
 100 void AliZDCMerger::Background(Float_t &fImpPar, Int_t &fSpecn, Int_t &fSpecp)
 101 {
 102
 103     // --- Open the background file
 104     if (fMerge && !fBgrFile) fBgrFile = OpenBgrFile();
 105
 106     // --- Read from the TreeE impact parameter (b),
 107     //     # of spectators n and p (fSpecn, fSpecp)
 108     fBgrFile->cd();
 109
 110 //    // Get AliRun object from file or create it if not on file
 111 //    gAlice = (AliRun*)fBgrFile->Get("gAlice");
 112 //    if (!gAlice) {
 113 //      gAlice = (AliRun*)fBgrFile->Get("gAlice");
 114 //      if (gAlice) printf("AliRun object found on file\n");
 115 //      if (!gAlice) {
 116 //        printf("\n create new gAlice object");
 117 //        gAlice = new AliRun("gAlice","Alice test program");
 118 //      }
 119 //    }
 120
 121     AliHeader *header = gAlice->GetHeader();
 122     AliGenEventHeader* mcHeader = header->GenEventHeader();
 123     fImpPar = ((AliGenHijingEventHeader*) mcHeader)->ImpactParameter();
 124     fSpecn  = ((AliGenHijingEventHeader*) mcHeader)->Spectatorsn();
 125     fSpecp  = ((AliGenHijingEventHeader*) mcHeader)->Spectatorsp();
 126     printf("\n  HIJING simulation - b = %f fm, Nspecn = %d, Nspecp = %d\n",fImpPar,fSpecn,fSpecp);
 127 }
 128
 129 //____________________________________________________________________________
 130 TFile* AliZDCMerger::OpenBgrFile()
 131 {
 132     // Initialise background event
 133     TFile *file = new TFile(fFnBgr,"UPDATE");
 134     printf("\n AliZDCMerger --- Background event -> %s file opened \n", fFnBgr);
 135     fHitsBgr = new TClonesArray("AliZDCHit",1000);
 136     fMHits   = new TClonesArray("AliZDCMergedHit",1000);
 137     return file;
 138 }
 139
 140 //____________________________________________________________________________
 141 void AliZDCMerger::Fragmentation(Float_t fImpPar, Int_t fSpecn, Int_t fSpecp,
 142                                  Int_t &fFreeSpn, Int_t &fFreeSpp)
 143 {
 144     Int_t j, zz[100], nn[100], nAlpha, Ztot, Ntot;
 145     AliZDCFragment *frag = new AliZDCFragment(fImpPar);
 146     for(j=0; j<=99; j++){
 147        zz[j] =0;
 148        nn[j] =0;
 149     }
 150
 151     // Fragments generation
 152     frag->GenerateIMF(zz, nAlpha);
 153
 154     // Attach neutrons
 155     Ztot=0;
 156     Ntot=0;
 157     frag->AttachNeutrons(zz, nn, Ztot, Ntot);
 158     fFreeSpn = fSpecn-Ztot-2*nAlpha;
 159     fFreeSpp = fSpecp-Ntot-2*nAlpha;
 160     printf("\n  Fragmentation -> FreeSpn = %d, FreeSpp = %d\n",fFreeSpn,fFreeSpp);
 161 }
 162
 163 //____________________________________________________________________________
 164 void AliZDCMerger::Mixing()
 165 {
 166
 167 //    printf("\n        AliZDCMerger->Mixing\n");
 168
 169     // ### Background event Hits ###########################################
 170     fBgrFile->cd();
 171 //    fBgrFile->ls();
 172
 173    AliZDC *ZDC = (AliZDC *)gAlice->GetModule("ZDC");
 174 //    if(ZDC) printf("\n        Ho trovato lo ZDC!\n");
 175
 176     // Hits tree
 177     if(fTrHBgr) delete fTrHBgr;
 178     fTrHBgr = 0;
 179     // SDigits tree
 180     if(fTrSDBgr) delete fTrSDBgr;
 181     fTrSDBgr = 0;
 182
 183 //    fNEvBgr = 0; // Let's suppose to have 1 full Hijing event per file
 184
 185     // Hits tree
 186     char treeBgrName[20];
 187     sprintf(treeBgrName,"TreeH%d",fNEvBgr);
 188     fTrHBgr = (TTree*)gDirectory->Get(treeBgrName);
 189     if(!fTrHBgr){
 190       printf("\n ERROR -> Can't find TreeH%d in background file\n",fNEvBgr);
 191     }
 192 //    fTrHBgr->Print();
 193
 194     // Branch address
 195     TBranch *branch;
 196     char branchname[20];
 197     sprintf(branchname,"%s",ZDC->GetName());
 198     if(fTrHBgr && fHitsBgr){
 199 //      printf("\n      fTrHBgr!=0 && fHitsBgr!=0\n");
 200       branch = fTrHBgr->GetBranch(branchname);
 201       if(branch) branch->SetAddress(&fHitsBgr);
 202     }
 203
 204     Int_t ntracks =  (Int_t) fTrHBgr->GetEntries();
 205 //    printf("\n        --- ntracks = %d\n\n", ntracks);
 206
 207     Int_t itrack, nhits, ihit, j, sector[2];
 208     AliZDCHit* zdcHit;
 209     AliZDCMergedHit *MHit;
 210     Float_t MHits[7];
 211     TClonesArray &sdigits = *fMHits;    // SDigits TCArray
 212     fNMhits = 0;
 213
 214     // --- Tracks loop
 215     for(itrack=0; itrack<ntracks; itrack++){
 216 //       printf("               itrack = %d", itrack);
 217 //       gAlice->ResetHits();
 218        fTrHBgr->GetEvent(itrack);
 219
 220 //       for(AliZDCHit* zdcHit=(AliZDCHit*)ZDC->FirstHit(-1);
 221 //         zdcHit;
 222 //         zdcHit = (AliZDCHit*)ZDC->NextHit()){
 223
 224        nhits = fHitsBgr->GetEntries();
 225 //       nhits = ZDCdp->GetNhits();
 226 //       printf("             nhits = %d \n", nhits);
 227        for(ihit=0; ihit<nhits; ihit++){
 228           zdcHit = (AliZDCHit*) fHitsBgr->UncheckedAt(ihit);
 229
 230             for(j=0; j<2; j++) sector[j] = zdcHit->GetVolume(j);
 231             MHits[0] = zdcHit->GetPrimKinEn();
 232             MHits[1] = zdcHit->GetXImpact();
 233             MHits[2] = zdcHit->GetYImpact();
 234             MHits[3] = zdcHit->GetSFlag();
 235             MHits[4] = zdcHit->GetLightPMQ();
 236             MHits[5] = zdcHit->GetLightPMC();
 237             MHits[6] = zdcHit->GetEnergy();
 238             MHit = new AliZDCMergedHit(sector, MHits);
 239 //          MHit->Print("");
 240             new((*fMHits)[fNMhits]) AliZDCMergedHit(*MHit);
 241             new (sdigits[fNMhits])  AliZDCMergedHit(*MHit);
 242             delete MHit;
 243             fNMhits++;
 244           }//Hits loop
 245
 246     } // Tracks loop
 247 //    printf("\n        fNMhits (after bckg) = %d, \n",fNMhits);
 248 //    fMHits->Dump();
 249 //    AliZDCMergedHit *pippo7 = (AliZDCMergedHit*) fMHits->At(7);
 250 //    pippo7->Dump();
 251
 252
 253     // ### Signal event Hits ###########################################
 254     // --- Neutrons
 255     ExtractSignal(1);
 256
 257     // --- Protons
 258     ExtractSignal(2);
 259 //    printf("\n        fNMhits (after signal) = %d \n",fNMhits);
 260
 261 }
 262
 263 //____________________________________________________________________________
 264 void AliZDCMerger::ExtractSignal(Int_t SpecType)
 265 {
 266
 267 // printf("\n   Entering in Extract Signal\n");
 268
 269  Int_t NumEvents = 0;
 270  if(SpecType == 1){             // --- Signal for spectator neutrons
 271    fFnSpecn = gSystem->ExpandPathName("$ALICE/$ALICE_LEVEL/ZDC/ZNsignalntu.root");
 272    fSpecnFile = TFile::Open(fFnSpecn,"R");
 273    fSpecnFile->cd();
 274 //   printf("\n --- ExtractSignal x n: file %s opened\n", fFnSpecn);
 275    NumEvents = fFreeSpn;
 276  }
 277  else if(SpecType == 2){        // --- Signal for spectator protons
 278    fFnSpecp = gSystem->ExpandPathName("$ALICE/$ALICE_LEVEL/ZDC/ZPsignalntu.root");
 279    fSpecpFile = TFile::Open(fFnSpecp,"R");
 280    fSpecpFile->cd();
 281 //   printf("\n --- ExtractSignal x p: file %s opened\n", fFnSpecp);
 282    NumEvents = fFreeSpp;
 283  }
 284 // printf("\n           # of free spectator = %d\n", NumEvents);
 285 // printf("\n         fNMhits (before adding signal) = %d\n",fNMhits);
 286
 287   TNtuple *ZDCSignal = (TNtuple*) gDirectory->Get("ZDCSignal");
 288   Int_t nentries = (Int_t) ZDCSignal->GetEntries();
 289 //  printf("\n   # entries = %d\n", nentries);
 290   Int_t i, j, k;
 291
 292   AliZDCMergedHit *MHit;
 293   Float_t *entry, HitsSpec[7];
 294   Int_t iev, rnd, Volume[2];
 295   for(iev=0; iev<NumEvents; iev++){
 296      rnd = (Int_t) (1000*gRandom->Rndm());
 297 //     printf("\n      rnd = %d\n", rnd);
 298      for(i=0; i<nentries; i++){
 299         ZDCSignal->GetEvent(i);
 300         entry = ZDCSignal->GetArgs();
 301         if(entry[0] == rnd){
 302 //        printf("\n   entry[0] = %f\n", entry[0]);
 303           for(k=0; k<2; k++) Volume[k] = (Int_t) entry[k+1];
 304           for(j=0; j<7; j++){
 305              HitsSpec[j] = entry[j+3];
 306           }
 307           MHit = new AliZDCMergedHit(Volume, HitsSpec);
 308           new((*fMHits)[fNMhits++]) AliZDCMergedHit(*MHit);
 309           delete MHit;
 310         }
 311         else if(entry[0] > rnd) break;
 312      }
 313   }
 314
 315   if(SpecType ==1){
 316 //    printf("\n         fNMhits (after n signal) = %d\n",fNMhits);
 317     fSpecnFile->Close();
 318   }
 319   else if(SpecType == 2){
 320 //    printf("\n         fNMhits (after p signal) = %d\n",fNMhits);
 321     fSpecpFile->Close();
 322   }
 323
 324 }
 325
 326 //____________________________________________________________________________
 327 void AliZDCMerger::Digitize(Int_t fNMhits, TClonesArray *fMHits)
 328 {
 329
 330 //  printf("\n  AliZDCMerger->Digitize()");
 331
 332   AliZDC *ZDC = (AliZDC *)gAlice->GetModule("ZDC");
 333 //  if(ZDC) printf("\n  Ho trovato lo ZDC!\n");
 334
 335   Int_t lightQ, lightC, sector[2], digit;
 336   Int_t PMCZN = 0, PMCZP = 0, PMQZN[4], PMQZP[4], PMZEM1 = 0, PMZEM2 = 0;
 337   Int_t i;
 338   for(i=0; i<4; i++){
 339      PMQZN[i] = 0;
 340      PMQZP[i] = 0;
 341   }
 342
 343   AliZDCMergedHit *MHit;
 344   Int_t imhit;
 345 //  printf("\n     fNMHits = %d\n", fNMhits);
 346   // Loop over SDigits
 347   for(imhit=0; imhit<fNMhits; imhit++){
 348
 349      MHit = (AliZDCMergedHit*) fMHits->UncheckedAt(imhit);
 350      sector[0] = MHit->GetSector(0);
 351      sector[1] = MHit->GetSector(1);
 352   if((sector[1]!=1) && (sector[1]!=2) && (sector[1]!=3) && (sector[1]!=4)){
 353 //     printf("\n  *** ERROR!!! sector[0] = %d, sector[1] = %d\n",
 354 //            sector[0], sector[1]);
 355        sector[1] = 0;
 356    }
 357      lightQ    = Int_t(MHit->GetLightPMQ());
 358      lightC    = Int_t(MHit->GetLightPMC());
 359 //     printf("\n    imhit = %d -> DET. = %d, quad = %d,PMQ = %d, PMC = %d",
 360 //      imhit,sector[0], sector[1],lightQ, lightC);
 361
 362      if(sector[0] == 1){           //ZN
 363        PMCZN = PMCZN + lightC;
 364        PMQZN[sector[1]-1] = PMQZN[sector[1]-1] + lightQ;
 365      }
 366      else if(sector[0] == 2){      //ZP
 367        PMCZP = PMCZP + lightC;
 368        PMQZP[sector[1]-1] = PMQZP[sector[1]-1] + lightQ;
 369      }
 370      else if(sector[0] == 3){      //ZEM
 371        if(sector[1] ==1) PMZEM1 = PMZEM1 + lightC;
 372        else              PMZEM2 = PMZEM2 + lightQ;
 373      }
 374   } // SDigits loop
 375
 376   // ### Digits creation ###############################################
 377   // Create digits for ZN
 378   Int_t PedValue;
 379   sector[0] = 1; // Detector = ZN
 380   sector[1] = 0; // Common PM ADC
 381   digit = Phe2ADCch(1, 0, PMCZN);
 382 //  printf("\n\n        ZN ###  PMCZN = %d      ADCZN = %d",PMCZN, digit);
 383   PedValue = AddPedestal();
 384   digit += PedValue;
 385 //  printf("    PedValue = %d",PedValue);
 386   ZDC->AddDigit(sector, digit);
 387   Int_t j;
 388   for(j=0; j<4; j++){
 389     sector[1] = j+1; // Towers PM ADCs
 390     digit = Phe2ADCch(1, j+1, PMQZN[j]);
 391 //    printf("\n                PMQZN[%d] = %d  phe     ADCZN[%d] = %d ADCch",j,PMQZN[j],j,digit);
 392     PedValue = AddPedestal();
 393     digit += PedValue;
 394 //    printf("  PedValue = %d",PedValue);
 395     ZDC->AddDigit(sector, digit);
 396   }
 397 //  printf("\n");
 398
 399   // Create digits for ZP
 400   sector[0] = 2; // Detector = ZP
 401   sector[1] = 0; // Common PM ADC
 402   digit = Phe2ADCch(2, 0, PMCZP);
 403 //  printf("\n  ZP --- PMCZP = %d       phe     ADCZP = %d ADCch",PMCZP,digit);
 404   PedValue = AddPedestal();
 405   digit += PedValue;
 406   printf("      PedValue = %d",PedValue);
 407   ZDC->AddDigit(sector, digit);
 408   for(j=0; j<4; j++){
 409     sector[1] = j+1; // Towers PM ADCs
 410     digit = Phe2ADCch(2, j+1, PMQZP[j]);
 411 //    printf("\n               PMQZP[%d] = %d   phe     ADCZP[%d] = %d ADCch",j,PMQZP[j],j,digit);
 412     PedValue = AddPedestal();
 413     digit += PedValue;
 414 //    printf("  PedValue = %d",PedValue);
 415     ZDC->AddDigit(sector, digit);
 416   }
 417 //  printf("\n");
 418
 419   // Create digits for ZEM
 420   sector[0] = 3;
 421   sector[1] = 1; // Detector = ZEM1
 422   digit  = Phe2ADCch(3, 1, PMZEM1);
 423 //  printf("\n  ZEM *** PMZEM1 = %d      phe     ADCZEM1 = %d ADCch",PMZEM1,digit);
 424   PedValue = AddPedestal();
 425   digit += PedValue;
 426 //  printf("  PedValue = %d\n",PedValue);
 427   ZDC->AddDigit(sector, digit);
 428   sector[1] = 2; // Detector = ZEM2
 429   digit  = Phe2ADCch(3, 2, PMZEM2);
 430 //  printf("\n  ZEM *** PMZEM2 = %d      phe     ADCZEM2 = %d ADCch",PMZEM2,digit);
 431   PedValue = AddPedestal();
 432   digit += PedValue;
 433 //  printf("  PedValue = %d\n",PedValue);
 434   ZDC->AddDigit(sector, digit);
 435 }
 436
 437 //_____________________________________________________________________________
 438 Int_t AliZDCMerger::Phe2ADCch(Int_t Det, Int_t Quad, Int_t Light)
 439 {
 440   // Evaluation of the ADC channel corresponding to the light yield Light
 441
 442   if(gAlice->GetDebug() > 0){
 443 //    printf("\n  Phe2ADCch -> Detector = %d, Quadrant = %d, Light = %d\n", Det, Quad, Light);
 444   }
 445
 446   Int_t ADCch = 0;
 447   // Parameters for conversion of light yield in ADC channels
 448   Float_t fPMGain[3][5];      // PM gain
 449   Float_t fADCRes;            // ADC conversion factor
 450
 451   Int_t j,i;
 452   for(i=0; i<3; i++){
 453      for(j=0; j<5; j++){
 454         fPMGain[i][j]   = 100000.;
 455      }
 456   }
 457   fADCRes   = 0.00000064; // ADC Resolution: 250 fC/ADCch
 458
 459   ADCch = (Int_t) (Light*fPMGain[Det-1][Quad]*fADCRes);
 460
 461   return ADCch;
 462 }
 463
 464 //_____________________________________________________________________________
 465 Int_t AliZDCMerger::AddPedestal()
 466 {
 467   // --- Pedestal value -> extracted from a gaussian distribution
 468   // obtained from the beam test on the ZEM prototype (Aug. 2000)
 469
 470   Int_t PedValue;
 471   Float_t PedMean  = 50.;
 472   Float_t PedWidth = 10.;
 473
 474   PedValue    = (Int_t) gRandom->Gaus(PedMean,PedWidth);
 475
 476   return PedValue;
 477 }