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