+inline void AliDielectronVarManager::GetZDCRP(const AliVEvent* event, Double_t qvec[][2]) {
+
+ //
+ // Get the reaction plane from the ZDC detector for first harmonic
+ //
+ // Q{x,y} = SUM{ri(x,y)*Ei} / SUM{Ei}
+ //
+
+ const Int_t nZDCSides = 2;
+ const Int_t nZDCplanes = 3;
+ const Int_t Aside = 0, Cside = 1, ACside = 2;
+ const Int_t nZDCTowers = 4;// number of ZDCtowers
+ const Double_t ZDCTowerCenters[nZDCTowers][2] = { {-1.75, -1.75}, { 1.75, -1.75},
+ {-1.75, 1.75}, { 1.75, 1.75} };
+
+ Double_t *ZDCTEnergy[nZDCSides]; //reco E in 5 ZDC sectors - high gain chain
+ Double_t qvecNUM[nZDCplanes][2];
+ Double_t qvecDEN[nZDCplanes];
+ memset( qvecNUM, 0, sizeof(qvecNUM)); //format
+ memset(qvecDEN, 0, sizeof(qvecDEN)); //format
+
+ Double_t TPCRefMulti = 999, vtxX = 999, vtxY = 999;
+ Int_t multiBin = 0, vtxXBin = 0, vtxYBin = 0;
+ Double_t recentdim[3][3] = { { 50, 0, 2500}, //multiplicity nbin, min, max
+ { 20, 0.04, 0.08}, // vertex x nbin, min, max
+ { 20, 0.25, 0.29} }; // vertex y nbin, min, max
+
+ if(!event->GetZDCData()) return;
+ AliVZDC* aliZDC = event->GetZDCData();
+ ZDCTEnergy[Aside] = (Double_t *)aliZDC -> GetZNATowerEnergy();
+ ZDCTEnergy[Cside] = (Double_t *)aliZDC -> GetZNCTowerEnergy();
+
+ for(int j = 0; j < nZDCSides ; j++){
+ for(int k = 0; k < nZDCTowers ; k++){
+ qvecNUM[j][0] += ZDCTowerCenters[k][0]*ZDCTEnergy[j][k+1]; // zdcQ += xE
+ qvecNUM[j][1] += ZDCTowerCenters[k][1]*ZDCTEnergy[j][k+1]; // zdcQ += yE
+ qvecDEN[j] += ZDCTEnergy[j][k+1]; // zdcQsum += E
+
+ }
+ if(j == Aside){
+ qvecNUM[j][0] = -qvecNUM[j][0];
+ }
+
+ if(j == Cside){
+ qvecNUM[j][0] = -qvecNUM[j][0];
+ qvecNUM[j][1] = -qvecNUM[j][1];
+ }
+
+
+ qvecNUM[ACside][0] += qvecNUM[j][0];
+ qvecNUM[ACside][1] += qvecNUM[j][1];
+ qvecDEN[ACside] += qvecDEN[j];
+
+ }
+
+ for(int j = 0; j < nZDCplanes; j++){
+ if(qvecDEN[j] != 0){
+ qvec[j][0] = (qvecNUM[j][0] / qvecDEN[j]);
+ qvec[j][1] = (qvecNUM[j][1] / qvecDEN[j]);
+ }
+ else if(qvecDEN[j] == 0) {
+ qvec[j][0] = 999;
+ qvec[j][1] = 999;
+ }
+
+ }
+
+ if(fgZDCRecentering[0][0]){
+ const AliAODEvent* aodEv = static_cast<const AliAODEvent*>(event);
+ AliAODHeader *header = aodEv->GetHeader();
+ if(!header) return;
+ TPCRefMulti = header -> GetTPConlyRefMultiplicity();
+
+ const AliVVertex *primVtx = event->GetPrimaryVertex();
+ if(!primVtx) return;
+ vtxX = primVtx->GetX();
+ vtxY = primVtx->GetY();
+
+ multiBin = (Int_t)((TPCRefMulti-recentdim[0][1])*recentdim[0][0] / (recentdim[0][2] - recentdim[0][1])) + 1;
+ vtxXBin = (Int_t)((vtxX-recentdim[1][1])*recentdim[1][0] / (recentdim[1][2] - recentdim[1][1])) + 1;
+ vtxYBin = (Int_t)((vtxY-recentdim[2][1])*recentdim[2][0] / (recentdim[2][2] - recentdim[2][1])) + 1;
+
+ for(int j = 0; j < nZDCplanes; j++)
+ if(qvecDEN[j] != 0){
+ qvec[j][0] -= fgZDCRecentering[j][0] -> GetBinContent(multiBin, vtxXBin, vtxYBin);
+ qvec[j][1] -= fgZDCRecentering[j][1] -> GetBinContent(multiBin, vtxXBin, vtxYBin);
+ }
+ }
+
+}
+
+