+ clustersTree->Branch("ZDC", "AliZDCReco", &reco, kBufferSize);
+ // write the output tree
+ clustersTree->Fill();
+}
+
+//_____________________________________________________________________________
+void AliZDCReconstructor::ReconstructEventPbPb(TTree *clustersTree,
+ Float_t* corrADCZN1, Float_t* corrADCZP1, Float_t* corrADCZN2, Float_t* corrADCZP2,
+ Float_t* corrADCZEM1, Float_t* corrADCZEM2, Float_t* sPMRef1, Float_t* sPMRef2,
+ Bool_t isScalerOn, UInt_t* scaler,
+ Int_t* evQualityBlock, Int_t* triggerBlock, Int_t* chBlock, UInt_t puBits) const
+{
+ // ****************** Reconstruct one event ******************
+ // ---------------------- Setting reco flags for ESD
+ UInt_t rFlags[32];
+ for(Int_t ifl=0; ifl<32; ifl++) rFlags[ifl]=0;
+
+ if(evQualityBlock[0] == 1) rFlags[31] = 0x0;
+ else rFlags[31] = 0x1;
+ //
+ if(evQualityBlock[1] == 1) rFlags[30] = 0x1;
+ if(evQualityBlock[2] == 1) rFlags[29] = 0x1;
+ if(evQualityBlock[3] == 1) rFlags[28] = 0x1;
+
+ if(triggerBlock[0] == 1) rFlags[27] = 0x1;
+ if(triggerBlock[1] == 1) rFlags[26] = 0x1;
+ if(triggerBlock[2] == 1) rFlags[25] = 0x1;
+ if(triggerBlock[3] == 1) rFlags[24] = 0x1;
+
+ if(chBlock[0] == 1) rFlags[18] = 0x1;
+ if(chBlock[1] == 1) rFlags[17] = 0x1;
+ if(chBlock[2] == 1) rFlags[16] = 0x1;
+
+ rFlags[13] = puBits & 0x00000020;
+ rFlags[12] = puBits & 0x00000010;
+ rFlags[11] = puBits & 0x00000080;
+ rFlags[10] = puBits & 0x00000040;
+ rFlags[9] = puBits & 0x00000020;
+ rFlags[8] = puBits & 0x00000010;
+
+ if(corrADCZP1[0]>fSignalThreshold) rFlags[5] = 0x1;
+ if(corrADCZN1[0]>fSignalThreshold) rFlags[4] = 0x1;
+ if(corrADCZEM2[0]>fSignalThreshold) rFlags[3] = 0x1;
+ if(corrADCZEM1[0]>fSignalThreshold) rFlags[2] = 0x1;
+ if(corrADCZP2[0]>fSignalThreshold) rFlags[1] = 0x1;
+ if(corrADCZN2[0]>fSignalThreshold) rFlags[0] = 0x1;
+
+ UInt_t recoFlag = rFlags[31] << 31 | rFlags[30] << 30 | rFlags[29] << 29 | rFlags[28] << 28 |
+ rFlags[27] << 27 | rFlags[26] << 26 | rFlags[25] << 25 | rFlags[24] << 24 |
+ 0x0 << 23 | 0x0 << 22 | 0x0 << 21 | 0x0 << 20 |
+ 0x0 << 19 | rFlags[18] << 18 | rFlags[17] << 17 | rFlags[16] << 16 |
+ 0x0 << 15 | 0x0 << 14 | rFlags[13] << 13 | rFlags[12] << 12 |
+ rFlags[11] << 11 |rFlags[10] << 10 | rFlags[9] << 9 | rFlags[8] << 8 |
+ 0x0 << 7 | 0x0 << 6 | rFlags[5] << 5 | rFlags[4] << 4 |
+ rFlags[3] << 3 | rFlags[2] << 2 | rFlags[1] << 1 | rFlags[0];
+ // --------------------------------------------------
+
+
+ // ****** Retrieving calibration data
+ // --- Equalization coefficients ---------------------------------------------
+ Float_t equalCoeffZN1[5], equalCoeffZP1[5], equalCoeffZN2[5], equalCoeffZP2[5];
+ for(Int_t ji=0; ji<5; ji++){
+ equalCoeffZN1[ji] = fTowCalibData->GetZN1EqualCoeff(ji);
+ equalCoeffZP1[ji] = fTowCalibData->GetZP1EqualCoeff(ji);
+ equalCoeffZN2[ji] = fTowCalibData->GetZN2EqualCoeff(ji);
+ equalCoeffZP2[ji] = fTowCalibData->GetZP2EqualCoeff(ji);
+ }
+ // --- Energy calibration factors ------------------------------------
+ Float_t valFromOCDB[6], calibEne[6];
+ for(Int_t ij=0; ij<6; ij++){
+ valFromOCDB[ij] = fEnCalibData->GetEnCalib(ij);
+ if(ij<4){
+ if(valFromOCDB[ij]!=0) calibEne[ij] = fBeamEnergy/valFromOCDB[ij];
+ else AliWarning(" Value from OCDB for E calibration = 0 !!!\n");
+ }
+ else calibEne[ij] = valFromOCDB[ij];
+ }
+
+ // ****** Equalization of detector responses
+ Float_t equalTowZN1[10], equalTowZN2[10], equalTowZP1[10], equalTowZP2[10];
+ for(Int_t gi=0; gi<10; gi++){
+ if(gi<5){
+ equalTowZN1[gi] = corrADCZN1[gi]*equalCoeffZN1[gi];
+ equalTowZP1[gi] = corrADCZP1[gi]*equalCoeffZP1[gi];
+ equalTowZN2[gi] = corrADCZN2[gi]*equalCoeffZN2[gi];
+ equalTowZP2[gi] = corrADCZP2[gi]*equalCoeffZP2[gi];
+ }
+ else{
+ equalTowZN1[gi] = corrADCZN1[gi]*equalCoeffZN1[gi-5];
+ equalTowZP1[gi] = corrADCZP1[gi]*equalCoeffZP1[gi-5];
+ equalTowZN2[gi] = corrADCZN2[gi]*equalCoeffZN2[gi-5];
+ equalTowZP2[gi] = corrADCZP2[gi]*equalCoeffZP2[gi-5];
+ }
+ }
+
+ // ****** Summed response for hadronic calorimeter (SUMMED and then CALIBRATED!)
+ Float_t calibSumZN1[]={0,0}, calibSumZN2[]={0,0}, calibSumZP1[]={0,0}, calibSumZP2[]={0,0};
+ for(Int_t gi=0; gi<5; gi++){
+ calibSumZN1[0] += equalTowZN1[gi];
+ calibSumZP1[0] += equalTowZP1[gi];
+ calibSumZN2[0] += equalTowZN2[gi];
+ calibSumZP2[0] += equalTowZP2[gi];
+ //
+ calibSumZN1[1] += equalTowZN1[gi+5];
+ calibSumZP1[1] += equalTowZP1[gi+5];
+ calibSumZN2[1] += equalTowZN2[gi+5];
+ calibSumZP2[1] += equalTowZP2[gi+5];
+ }
+ // High gain chain
+ calibSumZN1[0] = calibSumZN1[0]*calibEne[0]*8.;
+ calibSumZP1[0] = calibSumZP1[0]*calibEne[1]*8.;
+ calibSumZN2[0] = calibSumZN2[0]*calibEne[2]*8.;
+ calibSumZP2[0] = calibSumZP2[0]*calibEne[3]*8.;
+ // Low gain chain
+ calibSumZN1[1] = calibSumZN1[1]*calibEne[0];
+ calibSumZP1[1] = calibSumZP1[1]*calibEne[1];
+ calibSumZN2[1] = calibSumZN2[1]*calibEne[2];
+ calibSumZP2[1] = calibSumZP2[1]*calibEne[3];
+ //
+ Float_t sumZEM[]={0,0}, calibZEM1[]={0,0}, calibZEM2[]={0,0};
+ calibZEM1[0] = corrADCZEM1[0]*calibEne[5]*8.;
+ calibZEM1[1] = corrADCZEM1[1]*calibEne[5];
+ calibZEM2[0] = corrADCZEM2[0]*calibEne[5]*8.;
+ calibZEM2[1] = corrADCZEM2[1]*calibEne[5];
+ for(Int_t k=0; k<2; k++) sumZEM[k] = calibZEM1[k] + calibZEM2[k];
+
+ // ****** Energy calibration of detector responses
+ Float_t calibTowZN1[10], calibTowZN2[10], calibTowZP1[10], calibTowZP2[10];
+ for(Int_t gi=0; gi<5; gi++){
+ // High gain chain
+ calibTowZN1[gi] = equalTowZN1[gi]*calibEne[0]*8.;
+ calibTowZP1[gi] = equalTowZP1[gi]*calibEne[1]*8.;
+ calibTowZN2[gi] = equalTowZN2[gi]*calibEne[2]*8.;
+ calibTowZP2[gi] = equalTowZP2[gi]*calibEne[3]*8.;
+ // Low gain chain
+ calibTowZN1[gi+5] = equalTowZN1[gi+5]*calibEne[0];
+ calibTowZP1[gi+5] = equalTowZP1[gi+5]*calibEne[1];
+ calibTowZN2[gi+5] = equalTowZN2[gi+5]*calibEne[2];
+ calibTowZP2[gi+5] = equalTowZP2[gi+5]*calibEne[3];
+ }
+
+ // ****** Number of detected spectator nucleons
+ Int_t nDetSpecNLeft=0, nDetSpecPLeft=0, nDetSpecNRight=0, nDetSpecPRight=0;
+ if(fBeamEnergy!=0){
+ nDetSpecNLeft = (Int_t) (calibSumZN1[0]/fBeamEnergy);
+ nDetSpecPLeft = (Int_t) (calibSumZP1[0]/fBeamEnergy);
+ nDetSpecNRight = (Int_t) (calibSumZN2[0]/fBeamEnergy);
+ nDetSpecPRight = (Int_t) (calibSumZP2[0]/fBeamEnergy);
+ }
+ else AliWarning(" ATTENTION!!! fBeamEnergy=0 -> N_spec will be ZERO!!! \n");
+ /*printf("\n\t AliZDCReconstructor -> nDetSpecNLeft %d, nDetSpecPLeft %d,"
+ " nDetSpecNRight %d, nDetSpecPRight %d\n",nDetSpecNLeft, nDetSpecPLeft,
+ nDetSpecNRight, nDetSpecPRight);*/
+
+ Int_t nGenSpec=0, nGenSpecA=0, nGenSpecC=0;
+ Int_t nPart=0, nPartA=0, nPartC=0;
+ Double_t b=0., bA=0., bC=0.;
+
+ if(fIsCalibrationMB == kFALSE){
+ // ****** Reconstruction parameters ------------------
+ // Ch. debug
+ //fRecoParam->Print("");
+ //
+
+ if (!fRecoParam) fRecoParam = const_cast<AliZDCRecoParam*>(GetRecoParam());
+
+ TH2F *hZDCvsZEM = fRecoParam->GethZDCvsZEM();
+ TH2F *hZDCCvsZEM = fRecoParam->GethZDCCvsZEM();
+ TH2F *hZDCAvsZEM = fRecoParam->GethZDCAvsZEM();
+ TH1D *hNpartDist = fRecoParam->GethNpartDist();
+ TH1D *hbDist = fRecoParam->GethbDist();
+ Float_t ClkCenter = fRecoParam->GetClkCenter();
+ //
+ Double_t xHighEdge = hZDCvsZEM->GetXaxis()->GetXmax();
+ Double_t origin = xHighEdge*ClkCenter;
+ // Ch. debug
+ //printf("\n\n xHighEdge %1.2f, origin %1.4f \n", xHighEdge, origin);
+ //
+ // ====> Summed ZDC info (sideA+side C)
+ TF1 *line = new TF1("line","[0]*x+[1]",0.,xHighEdge);
+ Float_t y = (calibSumZN1[0]+calibSumZP1[0]+calibSumZN2[0]+calibSumZP2[0])/1000.;
+ Float_t x = (calibZEM1[0]+calibZEM2[0])/1000.;
+ line->SetParameter(0, y/(x-origin));
+ line->SetParameter(1, -origin*y/(x-origin));
+ // Ch. debug
+ //printf(" ***************** Summed ZDC info (sideA+side C) \n");
+ //printf(" E_{ZEM} %1.4f, E_{ZDC} %1.2f, TF1: %1.2f*x + %1.2f ", x, y,y/(x-origin),-origin*y/(x-origin));
+ //
+ Double_t countPerc=0;
+ Double_t xBinCenter=0, yBinCenter=0;
+ for(Int_t nbinx=1; nbinx<=hZDCvsZEM->GetNbinsX(); nbinx++){
+ for(Int_t nbiny=1; nbiny<=hZDCvsZEM->GetNbinsY(); nbiny++){
+ xBinCenter = hZDCvsZEM->GetXaxis()->GetBinCenter(nbinx);
+ yBinCenter = hZDCvsZEM->GetYaxis()->GetBinCenter(nbiny);
+ //
+ if(line->GetParameter(0)>0){
+ if(yBinCenter < (line->GetParameter(0)*xBinCenter + line->GetParameter(1))){
+ countPerc += hZDCvsZEM->GetBinContent(nbinx,nbiny);
+ // Ch. debug
+ /*printf(" xBinCenter %1.3f, yBinCenter %1.0f, countPerc %1.0f\n",
+ xBinCenter, yBinCenter, countPerc);*/
+ }
+ }
+ else{
+ if(yBinCenter > (line->GetParameter(0)*xBinCenter + line->GetParameter(1))){
+ countPerc += hZDCvsZEM->GetBinContent(nbinx,nbiny);
+ // Ch. debug
+ /*printf(" xBinCenter %1.3f, yBinCenter %1.0f, countPerc %1.0f\n",
+ xBinCenter, yBinCenter, countPerc);*/
+ }
+ }
+ }
+ }
+ //
+ Double_t xSecPerc = 0.;
+ if(hZDCvsZEM->GetEntries()!=0){
+ xSecPerc = countPerc/hZDCvsZEM->GetEntries();
+ }
+ else{
+ AliWarning(" Histogram hZDCvsZEM from OCDB has no entries!!!");
+ }
+ // Ch. debug
+ //printf(" xSecPerc %1.4f \n", xSecPerc);