From: coppedis Date: Mon, 4 Feb 2002 09:21:58 +0000 (+0000) Subject: Merging and reconstruction code review X-Git-Url: http://git.uio.no/git/?p=u%2Fmrichter%2FAliRoot.git;a=commitdiff_plain;h=1450a7cddaa6ddb84b3fd1f13fc85727d1da8724 Merging and reconstruction code review --- diff --git a/ZDC/AliZDC.cxx b/ZDC/AliZDC.cxx index 094103e95d1..aad66cdb8ff 100644 --- a/ZDC/AliZDC.cxx +++ b/ZDC/AliZDC.cxx @@ -15,6 +15,12 @@ /* $Log$ +Revision 1.26.2.2 2001/11/12 18:41:44 hristov +All the changes from the head are merged to the release + +Revision 1.27 2001/10/21 18:27:45 hristov +Several pointers were set to zero in the default constructors to avoid memory management problems + Revision 1.26 2001/10/04 14:30:28 coppedis Event merging for ZDC @@ -94,13 +100,16 @@ Introduction of the Copyright and cvs Log #include #include #include +#include +#include // --- AliRoot header files #include "AliZDC.h" #include "AliZDCHit.h" #include "AliZDCMergedHit.h" -#include "AliZDCDigit.h" #include "AliZDCMerger.h" +#include "AliZDCDigit.h" +#include "AliZDCReco.h" #include "AliDetector.h" #include "AliCallf77.h" #include "AliConst.h" @@ -118,22 +127,22 @@ AliZDC::AliZDC() // Default constructor for the Zero Degree Calorimeter base class // - fIshunt = 1; - fNoShower = 0; - fMerger = 0; + fIshunt = 1; + fNoShower = 0; + fMerger = 0; - fHits = 0; - fNhits = 0; + fHits = 0; + fNhits = 0; - fDigits = 0; - fNdigits = 0; + fDigits = 0; + fNdigits = 0; fMergedHits = 0; - fTreeSD = 0; - fTreeMD = 0; + fTreeSD = 0; + fTreeMD = 0; fNRecPoints = 0; - fRecPoints = 0; + fRecPoints = 0; } @@ -162,7 +171,7 @@ AliZDC::AliZDC(const char *name, const char *title) fTreeMD = 0; fNRecPoints = 0; - fRecPoints = 0; + fRecPoints = 0; } //____________________________________________________________________________ @@ -250,20 +259,6 @@ void AliZDC::AddDigit(Int_t *sect, Int_t digit) // AliZDCDigit *newdigit; newdigit = new AliZDCDigit(sect, digit); - -// AliZDCDigit *curdigit; -// TClonesArray &ldigits = *fDigits; -// -// Int_t j; -// for(j=0; j sector[0] = %d, sector[1] = %d, digit = %d", // sect[0], sect[1], digit); @@ -385,15 +380,26 @@ Float_t AliZDC::ZMax(void) const MakeBranchInTree(treeD, branchname, &fDigits, kBufferSize, file) ; printf("* AliZDC::MakeBranch * Making Branch %s for Digits\n\n",branchname); +} +//_____________________________________________________________________________ + void AliZDC::MakeBranchInTreeR(TTree *treeR, const char *file) +{ + // MakeBranchInTree + const Int_t kBufferSize = 4000; + char branchname[20]; + sprintf(branchname,"%s",GetName()); + MakeBranchInTree(treeR, branchname, &fRecPoints, kBufferSize, file) ; + printf("* AliZDC::MakeBranch * Making Branch %s for RecPoints\n\n",branchname); + } //_____________________________________________________________________________ void AliZDC::Hits2SDigits() { -// printf("\n Entering AliZDC::SDigits2Digits()\n"); + printf("\n Entering AliZDC::SDigits2Digits() "); //---------------------------------------------------------------- if(!fMerger){ -// printf("\n ZDC digitization (without merging)\n"); + printf(" ZDC digitization (without merging)\n"); AliZDCMergedHit *MHit; Int_t j, sector[2]; @@ -434,7 +440,7 @@ void AliZDC::Hits2SDigits() } //---------------------------------------------------------------- else if(fMerger){ -// printf("\n ZDC merging and digitization\n"); + printf(" ZDC merging and digitization\n"); // ### Initialise merging fMerger -> InitMerging(); @@ -449,15 +455,13 @@ void AliZDC::Hits2SDigits() printf("\n ERROR -> Can't find TreeS%d in background file\n",fNEvBgr); } // Branch address - TBranch *branchSD; char branchSDname[20]; sprintf(branchSDname,"%s",GetName()); if(fTreeSD && fMergedHits){ -// printf("\n fTreeSD!=0 && fMergedHits!=0\n"); - branchSD = fTreeSD->GetBranch(branchSDname); + TBranch *branchSD = fTreeSD->GetBranch(branchSDname); if(branchSD) branchSD->SetAddress(&fMergedHits); + else if(!branchSD) MakeBranchInTreeSD(fTreeSD); } - if(!branchSD) MakeBranchInTreeSD(fTreeSD); // ### Get TCA of MergedHits from AliZDCMerger fMergedHits = fMerger->MergedHits(); @@ -483,9 +487,9 @@ void AliZDC::Hits2SDigits() //_____________________________________________________________________________ void AliZDC::SDigits2Digits() { -// printf("\n Entering AliZDC::SDigits2Digits()\n"); + //printf("\n Entering AliZDC::SDigits2Digits() "); if(!fMerger){ // Only digitization -// printf("\n ZDC digitization (without merging) \n"); + printf(" ZDC digitization (without merging) \n"); fMerger = new AliZDCMerger(); fMerger->Digitize(fNMergedhits, fMergedHits); @@ -496,13 +500,14 @@ void AliZDC::SDigits2Digits() gAlice->TreeD()->Reset(); } else if(fMerger){ // Merging and digitization -// printf("\n ZDC merging and digitization\n"); + printf(" ZDC merging and digitization\n"); fMerger->Digitize(fNMergedhits, fMergedHits); TFile *bgrFile = fMerger->BgrFile(); bgrFile->cd(); // Digits tree Int_t fNEvBgr = fMerger->EvNum(); + //printf(" fNEvBgr = %d\n",fNEvBgr); char treeDBgrName[20]; sprintf(treeDBgrName,"TreeD%d",fNEvBgr); fTreeMD = (TTree*)gDirectory->Get(treeDBgrName); // TreeH @@ -510,17 +515,15 @@ void AliZDC::SDigits2Digits() printf("\n ERROR -> Can't find TreeD%d in background file\n",fNEvBgr); } // Branch address - TBranch *branchD; char branchDname[20]; sprintf(branchDname,"%s",GetName()); if(fTreeMD && fDigits){ // printf("\n fTreeMD!=0 && fDigits!=0\n"); - branchD = fTreeMD->GetBranch(branchDname); + TBranch *branchD = fTreeMD->GetBranch(branchDname); if(branchD) branchD->SetAddress(&fDigits); + else if(!branchD) MakeBranchInTreeD(fTreeMD); } - if(!branchD) MakeBranchInTreeD(fTreeMD); -// printf("\n ### Filling Digits tree\n"); fTreeMD->Fill(); fTreeMD->Write(0,TObject::kOverwrite); } @@ -537,7 +540,194 @@ void AliZDC::Hits2Digits() //_____________________________________________________________________________ void AliZDC::Digits2Reco() { + //printf("\n Entering AliZDC::Digits2Reco() "); + Int_t fNEvBgr = fMerger->EvNum(); + //printf(" fNEvBgr = %d\n",fNEvBgr); + gAlice->GetEvent(fNEvBgr); + + AliDetector *ZDC = gAlice->GetDetector("ZDC"); + TClonesArray *ZDCdigits = ZDC->Digits(); + + char tdname[20]; + sprintf(tdname,"TreeD%d",fNEvBgr); + TTree *TD = (TTree*)gDirectory->Get(tdname); + //TTree *TD = gAlice->TreeD(); + if(TD){ + //printf(" TreeD found in gAlice object\n"); + char brname[20]; + sprintf(brname,"%s",ZDC->GetName()); + TBranch *br = TD->GetBranch(brname); + if(br) br->SetAddress(&ZDCdigits); + } + else if(!TD) printf(" ERROR -> TreeD NOT found in gAlice object\n"); + + Int_t nt = (Int_t) (TD->GetEntries()); + //printf("\n #entries in TreeD = %d\n",nt); + gAlice->ResetDigits(); + + AliZDCDigit *dig; + Int_t j, idig, ndigits, ZNraw=0, ZPraw=0, ZEMraw=0; + // --- Summing raw ADCs for each detector to obtain total light + for(j=0; jGetEvent(j); + ndigits = ZDCdigits->GetEntries(); + //printf("\n Entry #%d, ndigits = %d",j,ndigits); + ZNraw=0; + ZPraw=0; + ZEMraw=0; + // --- Loop over event digits + for(idig=0; idigUncheckedAt(idig); + if(dig->GetSector(0) == 1) ZNraw += dig->GetADCValue(); + else if(dig->GetSector(0) == 2) ZPraw += dig->GetADCValue(); + else if(dig->GetSector(0) == 3) ZEMraw += dig->GetADCValue(); + } // Digits loop + } // TreeD entries loop + printf("\n --- ZNraw = %d, ZPraw = %d, ZEMraw = %d\n",ZNraw, ZPraw, ZEMraw); + + // --- Pedestal subtraction + Int_t ZNcorr, ZPcorr, ZEMcorr, MeanPed=50; + ZNcorr = ZNraw - 5*MeanPed; + ZPcorr = ZPraw - 5*MeanPed; + ZEMcorr = ZEMraw - 2*MeanPed; + if(ZNcorr<0) ZNcorr=0; + if(ZPcorr<0) ZPcorr=0; + if(ZEMcorr<0) ZEMcorr=0; + printf("\n ZNcorr = %d, ZPcorr = %d, ZEMcorr = %d\n",ZNcorr,ZPcorr,ZEMcorr); + + // --- ADCchannel -> photoelectrons + // NB-> PM gain = 10^(5), ADC resolution = 6.4*10^(-7) + Float_t ZNphe, ZPphe, ZEMphe, ConvFactor = 0.064; + ZNphe = ZNcorr/ConvFactor; + ZPphe = ZPcorr/ConvFactor; + ZEMphe = ZEMcorr/ConvFactor; + printf("\n ZNphe = %f, ZPphe = %f, ZEMphe = %f\n",ZNphe, ZPphe, ZEMphe); + + // --- Energy calibration + // Conversion factors for hadronic ZDCs goes from phe yield to TRUE incident + // energy (conversion from GeV to TeV is included); while for EM calos + // conversion is from light yield to detected energy calculated by GEANT + // NB -> ZN and ZP conversion factors are constant since incident spectators + // have all the same energy, ZEM energy is obtained through a fit over the whole + // range of incident particle energies (obtained with full HIJING simulations) + Float_t ZNenergy, ZPenergy, ZEMenergy, ZDCenergy; + Float_t ZNphexTeV=329., ZPphexTeV=369.; + ZNenergy = ZNphe/ZNphexTeV; + ZPenergy = ZPphe/ZPphexTeV; + ZDCenergy = ZNenergy+ZPenergy; + ZEMenergy = -4.81+0.3238*ZEMphe; + if(ZEMenergy<0) ZEMenergy=0; + printf(" ZNenergy = %f TeV, ZPenergy = %f TeV, ZDCenergy = %f GeV, " + "\n ZEMenergy = %f TeV\n", ZNenergy, ZPenergy, + ZDCenergy, ZEMenergy); + + if(ZDCenergy==0) + printf("\n\n ### ATTENZIONE!!! -> ev# %d: ZNenergy = %f TeV, ZPenergy = %f TeV, ZDCenergy = %f GeV, " + " ZEMenergy = %f TeV\n\n", fNEvBgr, ZNenergy, ZPenergy, ZDCenergy, ZEMenergy); + + // --- Number of incident spectator nucleons + Int_t NDetSpecN, NDetSpecP; + NDetSpecN = (Int_t) (ZNenergy/2.760); + NDetSpecP = (Int_t) (ZPenergy/2.760); + printf("\n NDetSpecN = %d, NDetSpecP = %d\n",NDetSpecN, NDetSpecP); + + // --- Number of generated spectator nucleons and impact parameter + // Fit results for neutrons (Nspectator n true vs. EZN) + TF1 *fZNCen = new TF1("fZNCen", + "(-2.116909+sqrt(2.116909*2.116909-4*(-0.00651)*(14.556798-x)))/(2*(-0.00651))",0.,158.5); + TF1 *fZNPer = new TF1("fZNPer", + "(-34.695134-sqrt(34.695134*34.695134-4*(-0.174780)*(-1562.283443-x)))/(2*(-0.174780))",0.,158.5); + // Fit results for protons (Nspectator p true vs. EZP) + TF1 *fZPCen = new TF1("fZPCen", + "(-1.3217+sqrt(1.3217*1.3217-4*(-0.007934)*(4.742873-x)))/(2*(-0.007934))",0.,58.91); + TF1 *fZPPer = new TF1("fZPPer", + "(-15.788267-sqrt(15.788267*15.788267-4*(-0.133359)*(-383.800673-x)))/(2*(-0.133359))",0.,58.91); + // Fit results for total number of spectators (Nspectators true vs. EZDC) + TF1 *fZDCCen = new TF1("fZDCCen", + "(-1.867335+sqrt(1.867335*1.867335-4*(-0.004119)*(19.100289-x)))/(2*(-0.004119))",0.,220.4); + TF1 *fZDCPer = new TF1("fZDCPer", + "(-22.429097-sqrt(22.429097*22.429097-4*(-0.072435)*(-1482.034526-x)))/(2*(-0.072435))",0.,220.4); + // Fit results for b (b vs. EZDC) + //TF1 *fbCen = new TF1("fbCen","0.611543+0.052231*x-0.000112*x*x+0.000000374*x*x*x",0.,222.); + //TF1 *fbPer = new TF1("fbPer","16.552010-0.023866*x-0.00001*x*x",0.,222.); + TF1 *fbCen = new TF1("fbCen","0.612769+0.051929*x-0.0001074*x*x+0.0000003724*x*x*x",0.,225.); + TF1 *fbPer = new TF1("fbPer","16.6131016-0.026053*x+0.000006893*x*x",0.,225.); + // Evaluating Nspectators and b from ZEM energy + TF1 *fZEMn = new TF1("fZEMn","124.2-0.0566*x+0.000006014*x*x",0.,3500.); + TF1 *fZEMp = new TF1("fZEMp","81.3-0.03834*x+0.000004359*x*x",0.,3500.); + TF1 *fZEMsp = new TF1("fZEMsp","205.6-0.09567*x+0.00001056*x*x",0.,3500.); + TF1 *fZEMb = new TF1("fZEMb","15.8-0.02084*x+2.802e-5*x*x-2.007e-8*x*x*x+6.586e-12*x*x*x*x-8.042e-16*x*x*x*x*x",0.,3500.); + + Int_t NGenSpecN=0, NGenSpecP=0, NGenSpec=0; + Double_t ImpPar=0; + Float_t EZEMCut = 360.; // Cut value for Ezem (GeV) + if(ZEMenergy >= EZEMCut){ + NGenSpecN = (Int_t) (fZNCen->Eval(ZNenergy)); + NGenSpecP = (Int_t) (fZPCen->Eval(ZPenergy)); + NGenSpec = (Int_t) (fZDCCen->Eval(ZDCenergy)); + ImpPar = fbCen->Eval(ZDCenergy); + //printf(" fZNCen = %f, fZPCen = %f, fZDCCen = %f\n",fZNCen->Eval(ZNenergy), + // fZPCen->Eval(ZPenergy),fZDCCen->Eval(ZDCenergy)); + } + else if(ZEMenergy < EZEMCut){ + NGenSpecN = (Int_t) (fZNPer->Eval(ZNenergy)); + NGenSpecP = (Int_t) (fZPPer->Eval(ZPenergy)); + NGenSpec = (Int_t) (fZDCPer->Eval(ZDCenergy)); + ImpPar = fbPer->Eval(ZDCenergy); + //printf(" fZNPer = %f, fZPPer = %f, fZDCPer = %f\n",fZNPer->Eval(ZNenergy), + // fZPPer->Eval(ZPenergy),fZDCPer->Eval(ZDCenergy)); + } + if(ZNenergy>158.5) NGenSpecN = (Int_t) (fZEMn->Eval(ZEMenergy)); + if(ZPenergy>58.91) NGenSpecP = (Int_t) (fZEMp->Eval(ZEMenergy)); + if(ZDCenergy>220.4) NGenSpec = (Int_t)(fZEMsp->Eval(ZEMenergy)); + if(ZDCenergy>225.) ImpPar = fZEMb->Eval(ZEMenergy); + /*if(ZNenergy>158.5) NGenSpecN = -999; + if(ZPenergy>58.91) NGenSpecP = -999; + if(ZDCenergy>220.4) NGenSpec = -999; + if(ZDCenergy>225.) ImpPar = -999;*/ + + if(NGenSpecN>125) NGenSpecN=125; + else if(NGenSpecN<0) NGenSpecN=0; + if(NGenSpecP>82) NGenSpecP=82; + else if(NGenSpecP<0) NGenSpecP=0; + if(NGenSpec>207) NGenSpec=207; + else if(NGenSpec<0) NGenSpec=0; + printf(" NRecSpecN = %d, NRecSpecP = %d, NRecSpec = %d\n",NGenSpecN,NGenSpecP,NGenSpec); + + // --- Number of participants + Int_t NPart, NPartTot; + NPart = 208-NGenSpecN-NGenSpecP; + NPartTot = 208-NGenSpec; + printf(" ### NPart(ZP+ZN) = %d, NPart(ZDC) = %d, b = %f fm\n",NPart,NPartTot,ImpPar); + + // --- Writing RecPoints TCA + // Allocate the RecPoints TCA + fRecPoints = new TClonesArray("AliZDCReco",1000); + AliZDCReco *reco = new AliZDCReco(ZNenergy,ZPenergy,ZDCenergy,ZEMenergy, + NDetSpecN,NDetSpecP,NGenSpecN,NGenSpecP,NGenSpec,NPartTot,ImpPar); + new((*fRecPoints)[fNRecPoints]) AliZDCReco(*reco); + //printf(" fNRecPoints = %d \n",fNRecPoints ); + //fNRecPoints++; + //fRecPoints->Dump(); + delete reco; + + // TreeR + TTree *treeR = gAlice->TreeR(); + char tname[20]; + sprintf(tname,"TreeR%d",fNEvBgr); + if(!treeR) printf("\n ERROR -> Can't find TreeR%d in background file\n",fNEvBgr); + // Branch address + char branchRname[20]; + sprintf(branchRname,"%s",GetName()); + if(fRecPoints){ + TBranch *branchR = treeR->GetBranch(branchRname); + if(branchR) branchR->SetAddress(&fRecPoints); + else if(!branchR) MakeBranchInTreeR(treeR); + } + treeR->Fill(); + treeR->Write(tname,TObject::kOverwrite); + treeR->Reset(); } diff --git a/ZDC/AliZDC.h b/ZDC/AliZDC.h index 12e5aaf6859..982989f5ecd 100644 --- a/ZDC/AliZDC.h +++ b/ZDC/AliZDC.h @@ -31,10 +31,12 @@ public: virtual void MakeBranch(Option_t* opt, const char *file=0); virtual void MakeBranchInTreeSD(TTree *treeSD, const char *file=0); virtual void MakeBranchInTreeD(TTree *treeD, const char *file=0); + virtual void MakeBranchInTreeR(TTree *treeD, const char *file=0); virtual void Hits2SDigits(); virtual void SDigits2Digits(); virtual void Hits2Digits(); virtual void Digits2Reco(); + TClonesArray *Reconstructed() const {return fRecPoints;} virtual void SetMerger(AliZDCMerger* merger); virtual AliZDCMerger* Merger(); virtual void StepManager() {} diff --git a/ZDC/AliZDCDigit.h b/ZDC/AliZDCDigit.h index 8836730e321..ee662b44ee2 100644 --- a/ZDC/AliZDCDigit.h +++ b/ZDC/AliZDCDigit.h @@ -22,8 +22,8 @@ class AliZDCDigit : public AliDigitNew { virtual ~AliZDCDigit() {} // Getters - virtual Float_t GetSector(Int_t i) {return fSector[i];} - virtual Float_t GetADCValue() {return fADCValue;} + virtual Int_t GetSector(Int_t i) {return fSector[i];} + virtual Int_t GetADCValue() {return fADCValue;} // Operators Int_t operator == (AliZDCDigit &digit) { @@ -42,12 +42,12 @@ class AliZDCDigit : public AliDigitNew { protected: //Data members - Int_t fSector[2]; // Detecor and tower in which light is produced - Float_t fADCValue; // ADC channel value + Int_t fSector[2]; // Detector and tower in which light is produced + Int_t fADCValue; // ADC channel value // Print method virtual void Print(Option_t *) { - printf(" -> DIGIT: Detector = %d Quadrant = %d ADCCh = %f\n ", + printf(" -> DIGIT: Detector = %d Quadrant = %d ADCCh = %d\n ", fSector[0], fSector[1], fADCValue); } diff --git a/ZDC/AliZDCMerger.cxx b/ZDC/AliZDCMerger.cxx index e23e3600d0e..a1a6a1454c6 100644 --- a/ZDC/AliZDCMerger.cxx +++ b/ZDC/AliZDCMerger.cxx @@ -43,12 +43,15 @@ ClassImp(AliZDCMerger) +//int comp(const void *i,const void *j) {return *(int *)i - *(int *)j;} + //____________________________________________________________________________ AliZDCMerger::AliZDCMerger() { // Default constructor - fMerge = kDigitize; -// fMerge = kMerge; + //fMerge = kDigitize -> Only digitization + //fMerge = kMerge -> Digitization + Merging + fMerge = kMerge; fFnBgr = 0; fBgrFile = 0; fNEvBgr = 0; @@ -66,7 +69,7 @@ AliZDCMerger::AliZDCMerger() fFnSpecp = 0; fSpecpFile = 0; fNMhits = 0; - fMHits = 0; + } //____________________________________________________________________________ @@ -87,11 +90,10 @@ void AliZDCMerger::InitMerging() // Hits tree, impact parameter, num. of spectators n & p // in background (full Hijing) event Float_t b; - Int_t nspecn, nspecp; - Background(b, nspecn, nspecp); + Background(b, fSpecn, fSpecp); // Production of nuclear fragments -> num. of FREE spectators n & p - Fragmentation(b, nspecn, nspecp, fFreeSpn, fFreeSpp); + Fragmentation(b, fSpecn, fSpecp, fFreeSpn, fFreeSpp); // Extract from spectators distribution the signal events: // NFreeSpectatorN spectator n & NFreeSpectatorP spectator p @@ -119,13 +121,18 @@ void AliZDCMerger::Background(Float_t &fImpPar, Int_t &fSpecn, Int_t &fSpecp) // gAlice = new AliRun("gAlice","Alice test program"); // } // } - + + gAlice->GetEvent(fNEvBgr); AliHeader *header = gAlice->GetHeader(); AliGenEventHeader* mcHeader = header->GenEventHeader(); fImpPar = ((AliGenHijingEventHeader*) mcHeader)->ImpactParameter(); - fSpecn = ((AliGenHijingEventHeader*) mcHeader)->Spectatorsn(); - fSpecp = ((AliGenHijingEventHeader*) mcHeader)->Spectatorsp(); - printf("\n HIJING simulation - b = %f fm, Nspecn = %d, Nspecp = %d\n",fImpPar,fSpecn,fSpecp); + Int_t dSpecn = ((AliGenHijingEventHeader*) mcHeader)->Spectatorsn(); + Int_t dSpecp = ((AliGenHijingEventHeader*) mcHeader)->Spectatorsp(); + // Until there is only 1 ZDC set the # of spectators must be divided by 2!!! + fSpecn = dSpecn/2; + fSpecp = dSpecp/2; + printf("\n HIJING ev. #%d - b = %f fm, Nspecn = %d, Nspecp = %d\n", + fNEvBgr,fImpPar,fSpecn,fSpecp); } //____________________________________________________________________________ @@ -143,6 +150,7 @@ TFile* AliZDCMerger::OpenBgrFile() void AliZDCMerger::Fragmentation(Float_t fImpPar, Int_t fSpecn, Int_t fSpecp, Int_t &fFreeSpn, Int_t &fFreeSpp) { + //printf("\n Fragmentation -> fSpecn = %d, fSpecp = %d\n",fSpecn,fSpecp); Int_t j, zz[100], nn[100], nAlpha, Ztot, Ntot; AliZDCFragment *frag = new AliZDCFragment(fImpPar); for(j=0; j<=99; j++){ @@ -157,8 +165,11 @@ void AliZDCMerger::Fragmentation(Float_t fImpPar, Int_t fSpecn, Int_t fSpecp, Ztot=0; Ntot=0; frag->AttachNeutrons(zz, nn, Ztot, Ntot); - fFreeSpn = fSpecn-Ztot-2*nAlpha; - fFreeSpp = fSpecp-Ntot-2*nAlpha; + fFreeSpn = fSpecn-Ntot-2*nAlpha; + fFreeSpp = fSpecp-Ztot-2*nAlpha; + if(fFreeSpn<0) fFreeSpn=0; + if(fFreeSpp<0) fFreeSpp=0; + //printf("\n 2*nAlpha = %d, Ztot = %d, Ntot = %d\n",2*nAlpha, Ztot, Ntot); printf("\n Fragmentation -> FreeSpn = %d, FreeSpp = %d\n",fFreeSpn,fFreeSpp); } @@ -166,7 +177,7 @@ void AliZDCMerger::Fragmentation(Float_t fImpPar, Int_t fSpecn, Int_t fSpecp, void AliZDCMerger::Mixing() { -// printf("\n AliZDCMerger->Mixing\n"); + //printf("\n AliZDCMerger->Mixing\n"); // ### Background event Hits ########################################### fBgrFile->cd(); @@ -183,7 +194,6 @@ void AliZDCMerger::Mixing() fTrSDBgr = 0; // fNEvBgr = 0; // Let's suppose to have 1 full Hijing event per file - // Hits tree char treeBgrName[20]; sprintf(treeBgrName,"TreeH%d",fNEvBgr); @@ -204,7 +214,7 @@ void AliZDCMerger::Mixing() } Int_t ntracks = (Int_t) fTrHBgr->GetEntries(); -// printf("\n --- ntracks = %d\n\n", ntracks); + //printf("\n --- ntracks = %d\n\n", ntracks); Int_t itrack, nhits, ihit, j, sector[2]; AliZDCHit* zdcHit; @@ -215,16 +225,9 @@ void AliZDCMerger::Mixing() // --- Tracks loop for(itrack=0; itrackResetHits(); fTrHBgr->GetEvent(itrack); -// for(AliZDCHit* zdcHit=(AliZDCHit*)ZDC->FirstHit(-1); -// zdcHit; -// zdcHit = (AliZDCHit*)ZDC->NextHit()){ - nhits = fHitsBgr->GetEntries(); -// nhits = ZDCdp->GetNhits(); // printf(" nhits = %d \n", nhits); for(ihit=0; ihitUncheckedAt(ihit); @@ -234,7 +237,7 @@ void AliZDCMerger::Mixing() MHits[1] = zdcHit->GetXImpact(); MHits[2] = zdcHit->GetYImpact(); MHits[3] = zdcHit->GetSFlag(); - MHits[4] = zdcHit->GetLightPMQ(); + MHits[4] = zdcHit->GetLightPMQ(); MHits[5] = zdcHit->GetLightPMC(); MHits[6] = zdcHit->GetEnergy(); MHit = new AliZDCMergedHit(sector, MHits); @@ -246,11 +249,7 @@ void AliZDCMerger::Mixing() }//Hits loop } // Tracks loop -// printf("\n fNMhits (after bckg) = %d, \n",fNMhits); -// fMHits->Dump(); -// AliZDCMergedHit *pippo7 = (AliZDCMergedHit*) fMHits->At(7); -// pippo7->Dump(); - + //printf(" fNMhits (after bckg) = %d, \n",fNMhits); // ### Signal event Hits ########################################### // --- Neutrons @@ -258,7 +257,7 @@ void AliZDCMerger::Mixing() // --- Protons ExtractSignal(2); -// printf("\n fNMhits (after signal) = %d \n",fNMhits); + //printf(" fNMhits (after signal) = %d \n",fNMhits); } @@ -273,53 +272,66 @@ void AliZDCMerger::ExtractSignal(Int_t SpecType) fFnSpecn = gSystem->ExpandPathName("$ALICE/$ALICE_LEVEL/ZDC/ZNsignalntu.root"); fSpecnFile = TFile::Open(fFnSpecn,"R"); fSpecnFile->cd(); -// printf("\n --- ExtractSignal x n: file %s opened\n", fFnSpecn); + printf("\n --- ExtractSignal x n: file %s opened\n", fFnSpecn); NumEvents = fFreeSpn; } else if(SpecType == 2){ // --- Signal for spectator protons fFnSpecp = gSystem->ExpandPathName("$ALICE/$ALICE_LEVEL/ZDC/ZPsignalntu.root"); fSpecpFile = TFile::Open(fFnSpecp,"R"); fSpecpFile->cd(); -// printf("\n --- ExtractSignal x p: file %s opened\n", fFnSpecp); + printf("\n --- ExtractSignal x p: file %s opened\n", fFnSpecp); NumEvents = fFreeSpp; } -// printf("\n # of free spectator = %d\n", NumEvents); -// printf("\n fNMhits (before adding signal) = %d\n",fNMhits); + //printf("\n # of free spectators = %d\n", NumEvents); + //printf("\n fNMhits (before adding signal) = %d\n",fNMhits); TNtuple *ZDCSignal = (TNtuple*) gDirectory->Get("ZDCSignal"); Int_t nentries = (Int_t) ZDCSignal->GetEntries(); -// printf("\n # entries = %d\n", nentries); - Int_t i, j, k; + //printf("\n # entries = %d\n", nentries); AliZDCMergedHit *MHit; Float_t *entry, HitsSpec[7]; - Int_t iev, rnd, Volume[2]; - for(iev=0; ievRndm()); -// printf("\n rnd = %d\n", rnd); + Int_t pl, i, j, k, iev=0, rnd[125], Volume[2]; + for(pl=0;pl<125;pl++){ + rnd[pl] = 0; + } + for(pl=0;plRndm()); + if(rnd[pl] >= 9998) rnd[pl] = 9997; + //printf(" rnd[%d] = %d\n",pl,rnd[pl]); + } + // Sorting vector in ascending order with C function QSORT + qsort((void*)rnd,NumEvents,sizeof(Int_t),comp); + //for(pl=0;plGetEvent(i); entry = ZDCSignal->GetArgs(); - if(entry[0] == rnd){ -// printf("\n entry[0] = %f\n", entry[0]); + if(entry[0] == rnd[iev]){ for(k=0; k<2; k++) Volume[k] = (Int_t) entry[k+1]; for(j=0; j<7; j++){ HitsSpec[j] = entry[j+3]; } + //printf("\n i = %d, iev = %d, entry[0] = %f, rnd[%d] = %d ",i,iev,entry[0],iev,rnd[iev]); MHit = new AliZDCMergedHit(Volume, HitsSpec); new((*fMHits)[fNMhits++]) AliZDCMergedHit(*MHit); delete MHit; } - else if(entry[0] > rnd) break; + else if(entry[0] > rnd[iev]){ + iev++; + continue; + } } - } + }while(ievClose(); } else if(SpecType == 2){ -// printf("\n fNMhits (after p signal) = %d\n",fNMhits); + //printf("\n fNMhits (after p signal) = %d\n",fNMhits); fSpecpFile->Close(); } @@ -329,7 +341,7 @@ void AliZDCMerger::ExtractSignal(Int_t SpecType) void AliZDCMerger::Digitize(Int_t fNMhits, TClonesArray *fMHits) { -// printf("\n AliZDCMerger->Digitize()"); + printf("\n AliZDCMerger->Digitize()"); AliZDC *ZDC = (AliZDC *)gAlice->GetModule("ZDC"); // if(ZDC) printf("\n Ho trovato lo ZDC!\n"); @@ -344,21 +356,22 @@ void AliZDCMerger::Digitize(Int_t fNMhits, TClonesArray *fMHits) AliZDCMergedHit *MHit; Int_t imhit; -// printf("\n fNMHits = %d\n", fNMhits); + printf(" fNMHits = %d\n", fNMhits); // Loop over SDigits for(imhit=0; imhitUncheckedAt(imhit); sector[0] = MHit->GetSector(0); sector[1] = MHit->GetSector(1); + // tmp -> c'erano i quadranti cannati! if((sector[1]!=1) && (sector[1]!=2) && (sector[1]!=3) && (sector[1]!=4)){ -// printf("\n *** ERROR!!! sector[0] = %d, sector[1] = %d\n", -// sector[0], sector[1]); + printf("\n *** ERROR!!! sector[0] = %d, sector[1] = %d\n", + sector[0], sector[1]); sector[1] = 0; } lightQ = Int_t(MHit->GetLightPMQ()); lightC = Int_t(MHit->GetLightPMC()); -// printf("\n imhit = %d -> DET. = %d, quad = %d,PMQ = %d, PMC = %d", +// printf(" imhit = %d -> DET. = %d, quad = %d,PMQ = %d, PMC = %d\n", // imhit,sector[0], sector[1],lightQ, lightC); if(sector[0] == 1){ //ZN @@ -381,7 +394,7 @@ void AliZDCMerger::Digitize(Int_t fNMhits, TClonesArray *fMHits) sector[0] = 1; // Detector = ZN sector[1] = 0; // Common PM ADC digit = Phe2ADCch(1, 0, PMCZN); -// printf("\n\n ZN ### PMCZN = %d ADCZN = %d",PMCZN, digit); + printf("\n\n ZN ### PMCZN = %d ADCZN = %d",PMCZN, digit); PedValue = AddPedestal(); digit += PedValue; // printf(" PedValue = %d",PedValue); @@ -390,46 +403,46 @@ void AliZDCMerger::Digitize(Int_t fNMhits, TClonesArray *fMHits) for(j=0; j<4; j++){ sector[1] = j+1; // Towers PM ADCs digit = Phe2ADCch(1, j+1, PMQZN[j]); -// printf("\n PMQZN[%d] = %d phe ADCZN[%d] = %d ADCch",j,PMQZN[j],j,digit); + printf("\n PMQZN[%d] = %d phe ADCZN[%d] = %d ADCch",j,PMQZN[j],j,digit); PedValue = AddPedestal(); digit += PedValue; // printf(" PedValue = %d",PedValue); ZDC->AddDigit(sector, digit); } -// printf("\n"); + printf("\n"); // Create digits for ZP sector[0] = 2; // Detector = ZP sector[1] = 0; // Common PM ADC digit = Phe2ADCch(2, 0, PMCZP); -// printf("\n ZP --- PMCZP = %d phe ADCZP = %d ADCch",PMCZP,digit); + printf("\n ZP --- PMCZP = %d phe ADCZP = %d ADCch",PMCZP,digit); PedValue = AddPedestal(); digit += PedValue; - printf(" PedValue = %d",PedValue); +// printf(" PedValue = %d",PedValue); ZDC->AddDigit(sector, digit); for(j=0; j<4; j++){ sector[1] = j+1; // Towers PM ADCs digit = Phe2ADCch(2, j+1, PMQZP[j]); -// printf("\n PMQZP[%d] = %d phe ADCZP[%d] = %d ADCch",j,PMQZP[j],j,digit); + printf("\n PMQZP[%d] = %d phe ADCZP[%d] = %d ADCch",j,PMQZP[j],j,digit); PedValue = AddPedestal(); digit += PedValue; // printf(" PedValue = %d",PedValue); ZDC->AddDigit(sector, digit); } -// printf("\n"); + printf("\n"); // Create digits for ZEM sector[0] = 3; sector[1] = 1; // Detector = ZEM1 digit = Phe2ADCch(3, 1, PMZEM1); -// printf("\n ZEM *** PMZEM1 = %d phe ADCZEM1 = %d ADCch",PMZEM1,digit); + printf("\n ZEM *** PMZEM1 = %d phe ADCZEM1 = %d ADCch",PMZEM1,digit); PedValue = AddPedestal(); digit += PedValue; // printf(" PedValue = %d\n",PedValue); ZDC->AddDigit(sector, digit); sector[1] = 2; // Detector = ZEM2 digit = Phe2ADCch(3, 2, PMZEM2); -// printf("\n ZEM *** PMZEM2 = %d phe ADCZEM2 = %d ADCch",PMZEM2,digit); + printf("\n ZEM *** PMZEM2 = %d phe ADCZEM2 = %d ADCch\n",PMZEM2,digit); PedValue = AddPedestal(); digit += PedValue; // printf(" PedValue = %d\n",PedValue); @@ -446,9 +459,6 @@ Int_t AliZDCMerger::Phe2ADCch(Int_t Det, Int_t Quad, Int_t Light) } Int_t ADCch = 0; - // Parameters for conversion of light yield in ADC channels - Float_t fPMGain[3][5]; // PM gain - Float_t fADCRes; // ADC conversion factor Int_t j,i; for(i=0; i<3; i++){ @@ -471,7 +481,7 @@ Int_t AliZDCMerger::AddPedestal() Int_t PedValue; Float_t PedMean = 50.; - Float_t PedWidth = 10.; + Float_t PedWidth = 5.; PedValue = (Int_t) gRandom->Gaus(PedMean,PedWidth); diff --git a/ZDC/AliZDCMerger.h b/ZDC/AliZDCMerger.h index 934b953e289..b4474aa9c97 100644 --- a/ZDC/AliZDCMerger.h +++ b/ZDC/AliZDCMerger.h @@ -45,7 +45,6 @@ public: void SetBackgroundFileName(char* file) {fFnBgr = file;} void SetBackgroundEventNum(Int_t nev) {fNEvBgr = nev;} -private: //Open the background file TFile *OpenBgrFile(); @@ -75,24 +74,25 @@ protected: Int_t fNMhits; // Number of Merged hits for background TClonesArray *fMHits; // TCA for "merged" hits - -// // *** Digits -// // --- Digitization parameters setters and getters -// // PM gain -// void SetPMGain(Int_t Det, Int_t PMDet, Int_t PMGain) -// {fPMGain[Det][PMDet] = PMGain;} -// Float_t GetPMGain(Int_t Det, Int_t PMDet) -// {return fPMGain[Det][PMDet];} -// // Conversion factor from charge to ADC channels -// // F = 1.6E-19 / Resolution [Coulomb/ch] -// void SetADCRes(Int_t ADCRes) {fADCRes = ADCRes;} -// Float_t GetADCRes() {return fADCRes;} -// -// // --- Parameters for conversion of light yield in ADC channels -// Float_t fPMGain[3][5]; // PM gain -// Float_t fADCRes; // ADC conversion factor + +public: + // *** Digits + // --- Parameters for conversion of light yield in ADC channels + Float_t fPMGain[3][5]; // PM gain + Float_t fADCRes; // ADC conversion factor + // --- Digitization parameters setters and getters + // PM gain + void SetPMGain(Int_t Det, Int_t PMDet, Int_t PMGain) + {fPMGain[Det][PMDet] = PMGain;} + Float_t GetPMGain(Int_t Det, Int_t PMDet) + {return fPMGain[Det][PMDet];} + // Conversion factor from charge to ADC channels + // F = 1.6E-19 / Resolution [Coulomb/ch] + void SetADCRes(Int_t ADCRes) {fADCRes = ADCRes;} + Float_t GetADCRes() {return fADCRes;} + - ClassDef(AliZDCMerger,0) + ClassDef(AliZDCMerger,1) }; #endif diff --git a/ZDC/AliZDCReco.cxx b/ZDC/AliZDCReco.cxx new file mode 100644 index 00000000000..2296ef40af8 --- /dev/null +++ b/ZDC/AliZDCReco.cxx @@ -0,0 +1,45 @@ +/************************************************************************** + * 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. * + **************************************************************************/ + +/* +$Log +*/ + +//////////////////////////////////////////////// +// RecPoints classes for set ZDC // +//////////////////////////////////////////////// + + +#include "AliZDCReco.h" + +ClassImp(AliZDCReco) + +//_____________________________________________________________________________ +AliZDCReco::AliZDCReco(Float_t ezn, Float_t ezp, Float_t ezdc, Float_t ezem, + Int_t detspn, Int_t detspp, Int_t trspn, Int_t trspp, Int_t trsp, Int_t part, Float_t b) +{ + fZNenergy = ezn; + fZPenergy = ezp; + fZDCenergy = ezdc; + fZEMenergy = ezem; + fNDetSpecN = detspn; + fNDetSpecP = detspp; + fNTrueSpecN = trspn; + fNTrueSpecP = trspp; + fNTrueSpec = trsp; + fNPart = part; + fImpPar = b; + +} diff --git a/ZDC/AliZDCReco.h b/ZDC/AliZDCReco.h new file mode 100644 index 00000000000..d6657fd1f4d --- /dev/null +++ b/ZDC/AliZDCReco.h @@ -0,0 +1,59 @@ +#ifndef ALIZDCRECO_H +#define ALIZDCRECO_H +/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * See cxx source for full Copyright notice */ + +//////////////////////////////////////////////// +// Classe for ZDC RecPoints // +//////////////////////////////////////////////// + +#include "TObject.h" + +class AliZDCReco : public TObject { + +public: + AliZDCReco() {} + AliZDCReco(Float_t ezn, Float_t ezp, Float_t ezdc, Float_t ezem, Int_t detspn, + Int_t detspp, Int_t trspn, Int_t trspp, Int_t trsp, Int_t part, Float_t b); + AliZDCReco(AliZDCReco* oldreco) {*this=*oldreco;} + virtual ~AliZDCReco() {} + + // Getters + virtual Float_t GetZNenergy() {return fZNenergy;} + virtual Float_t GetZPenergy() {return fZPenergy;} + virtual Float_t GetZDCenergy() {return fZDCenergy;} + virtual Float_t GetZEMenergy() {return fZEMenergy;} + virtual Int_t GetNDetSpecN() {return fNDetSpecN;} + virtual Int_t GetNDetSpecP() {return fNDetSpecP;} + virtual Int_t GetNTrueSpecN() {return fNTrueSpecN;} + virtual Int_t GetNTrueSpecP() {return fNTrueSpecP;} + virtual Int_t GetNTrueSpec() {return fNTrueSpec;} + virtual Int_t GetNPart() {return fNPart;} + virtual Float_t GetImpPar() {return fImpPar;} + + // Data members + Float_t fZNenergy; // Energy detected in neutron ZDC + Float_t fZPenergy; // Energy detected in proton ZDC + Float_t fZDCenergy; // Total hadronic energy detcted in ZDCs + Float_t fZEMenergy; // Energy detected in EM ZDC + Int_t fNDetSpecN; // Number of spectator neutrons detected + Int_t fNDetSpecP; // Number of spectator protons detected + Int_t fNTrueSpecN; // Estimate of the number of spectator neutrons generated + Int_t fNTrueSpecP; // Estimate of the number of spectator protons generated + Int_t fNTrueSpec ; // Estimate of the total number of spectators + Int_t fNPart; // Estimate of the number of participants for 1 nucleus + Float_t fImpPar; // Estimate of the impact parameter + + // Print method + virtual void Print(Option_t *) const { + printf(" --- Reconstruction -> EZN = %f TeV, EZP = %f TeV, EZDC = %f TeV," + " EZEM = %f GeV \n NDetSpecN = %d, NDetSpecP = %d, Nspecn = %d," + " Nspecp = %d, Npart = %d, b = %f fm.\n ", + fZNenergy,fZPenergy,fZDCenergy,fZEMenergy,fNDetSpecN,fNDetSpecP, + fNTrueSpecN,fNTrueSpecP,fNPart,fImpPar); + } + + ClassDef(AliZDCReco,1) // RecPoints for the Zero Degree Calorimeters +}; + +#endif