From 5591748ed1607a6ceb8a36febb7bed8147ca62dd Mon Sep 17 00:00:00 2001 From: dgangadh Date: Wed, 28 May 2014 17:01:17 +0200 Subject: [PATCH] code updates, Add Muon Corrections, KT4index bug fix --- PWGCF/FEMTOSCOPY/Chaoticity/AliFourPion.cxx | 372 +++++++++++--------- PWGCF/FEMTOSCOPY/Chaoticity/AliFourPion.h | 22 +- 2 files changed, 220 insertions(+), 174 deletions(-) diff --git a/PWGCF/FEMTOSCOPY/Chaoticity/AliFourPion.cxx b/PWGCF/FEMTOSCOPY/Chaoticity/AliFourPion.cxx index 4b4449880ca..2e4dc0bc4b9 100755 --- a/PWGCF/FEMTOSCOPY/Chaoticity/AliFourPion.cxx +++ b/PWGCF/FEMTOSCOPY/Chaoticity/AliFourPion.cxx @@ -142,7 +142,8 @@ AliAnalysisTaskSE(), fNormQPairSwitch_E1E2(), fNormQPairSwitch_E1E3(), fNormQPairSwitch_E2E3(), - fMomResC2(0x0) + fMomResC2(0x0), + fWeightmuonCorrection(0x0) { // Default constructor for(Int_t mb=0; mbData()); nameMuonIdeal->Append("_MuonIdeal"); - Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fMuonIdeal = new TH2D(nameMuonIdeal->Data(),"", 11,0.5,11.5, fQbinsQ3,0,fQupperBoundQ3); - if(mb==0 && edB==0) fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fMuonIdeal); + Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fMuonIdeal = new TH3D(nameMuonIdeal->Data(),"", 2,0.5,2.5, 11,0.5,11.5, fQbinsQ3,0,fQupperBoundQ3); + if(mb==0 && edB==0 && term<4) fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fMuonIdeal); TString *nameMuonSmeared=new TString(namePC3->Data()); nameMuonSmeared->Append("_MuonSmeared"); - Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fMuonSmeared = new TH2D(nameMuonSmeared->Data(),"", 11,0.5,11.5, fQbinsQ3,0,fQupperBoundQ3); - if(mb==0 && edB==0) fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fMuonSmeared); + Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fMuonSmeared = new TH3D(nameMuonSmeared->Data(),"", 2,0.5,2.5, 11,0.5,11.5, fQbinsQ3,0,fQupperBoundQ3); + if(mb==0 && edB==0 && term<4) fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fMuonSmeared); // TString *nameMuonPionK3=new TString(namePC3->Data()); nameMuonPionK3->Append("_MuonPionK3"); - Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fMuonPionK3 = new TH2D(nameMuonPionK3->Data(),"", 11,0.5,11.5, fQbinsQ3,0,fQupperBoundQ3); - if(mb==0 && edB==0) fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fMuonPionK3); + Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fMuonPionK3 = new TH3D(nameMuonPionK3->Data(),"", 2,0.5,2.5, 11,0.5,11.5, fQbinsQ3,0,fQupperBoundQ3); + if(mb==0 && edB==0 && term<4) fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fMuonPionK3); // TString *namePionPionK3=new TString(namePC3->Data()); namePionPionK3->Append("_PionPionK3"); - Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fPionPionK3 = new TH2D(namePionPionK3->Data(),"", 11,0.5,11.5, fQbinsQ3,0,fQupperBoundQ3); - if(mb==0 && edB==0) fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fPionPionK3); + Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fPionPionK3 = new TH3D(namePionPionK3->Data(),"", 2,0.5,2.5, 11,0.5,11.5, fQbinsQ3,0,fQupperBoundQ3); + if(mb==0 && edB==0 && term<4) fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fPionPionK3); }// MCcase // @@ -1243,7 +1249,7 @@ void AliFourPion::UserCreateOutputObjects() Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fTwoPartNorm = new TH2D(nameTwoPartNorm->Data(),"", kDENtypes,0.5,kDENtypes+0.5, fQbinsQ4,0,fQupperBoundQ4); fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fTwoPartNorm); } - + if(fMCcase==kTRUE){ // Momentum resolution correction histos TString *nameMomResIdeal=new TString(namePC4->Data()); @@ -1257,22 +1263,22 @@ void AliFourPion::UserCreateOutputObjects() // Muon correction histos TString *nameMuonIdeal=new TString(namePC4->Data()); nameMuonIdeal->Append("_MuonIdeal"); - Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fMuonIdeal = new TH2D(nameMuonIdeal->Data(),"", 11,0.5,11.5, fQbinsQ4,0,fQupperBoundQ4); - if(mb==0 && edB==0) fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fMuonIdeal); + Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fMuonIdeal = new TH3D(nameMuonIdeal->Data(),"", 2,0.5,2.5, 11,0.5,11.5, fQbinsQ4,0,fQupperBoundQ4); + if(mb==0 && edB==0 && term<12) fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fMuonIdeal); TString *nameMuonSmeared=new TString(namePC4->Data()); nameMuonSmeared->Append("_MuonSmeared"); - Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fMuonSmeared = new TH2D(nameMuonSmeared->Data(),"", 11,0.5,11.5, fQbinsQ4,0,fQupperBoundQ4); - if(mb==0 && edB==0) fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fMuonSmeared); + Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fMuonSmeared = new TH3D(nameMuonSmeared->Data(),"", 2,0.5,2.5, 11,0.5,11.5, fQbinsQ4,0,fQupperBoundQ4); + if(mb==0 && edB==0 && term<12) fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fMuonSmeared); // TString *nameMuonPionK4=new TString(namePC4->Data()); nameMuonPionK4->Append("_MuonPionK4"); - Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fMuonPionK4 = new TH2D(nameMuonPionK4->Data(),"", 11,0.5,11.5, fQbinsQ4,0,fQupperBoundQ4); - if(mb==0 && edB==0) fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fMuonPionK4); + Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fMuonPionK4 = new TH3D(nameMuonPionK4->Data(),"", 2,0.5,2.5, 11,0.5,11.5, fQbinsQ4,0,fQupperBoundQ4); + if(mb==0 && edB==0 && term<12) fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fMuonPionK4); // TString *namePionPionK4=new TString(namePC4->Data()); namePionPionK4->Append("_PionPionK4"); - Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fPionPionK4 = new TH2D(namePionPionK4->Data(),"", 11,0.5,11.5, fQbinsQ4,0,fQupperBoundQ4); - if(mb==0 && edB==0) fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fPionPionK4); + Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fPionPionK4 = new TH3D(namePionPionK4->Data(),"", 2,0.5,2.5, 11,0.5,11.5, fQbinsQ4,0,fQupperBoundQ4); + if(mb==0 && edB==0 && term<12) fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fPionPionK4); }// MCcase @@ -1352,6 +1358,9 @@ void AliFourPion::UserCreateOutputObjects() TH2D *fAvgQ23VersusQ3 = new TH2D("fAvgQ23VersusQ3","",10,0,0.1, 20,0,0.1); fOutputList->Add(fAvgQ23VersusQ3); + TH1D *fDistPionParents4 = new TH1D("fDistPionParents4","",4,0.5,4.5); + fOutputList->Add(fDistPionParents4); + //////////////////////////////////// /////////////////////////////////// @@ -1836,7 +1845,12 @@ void AliFourPion::UserExec(Option_t *) Float_t Pparent4[4]={0}; Float_t weight12=0, weight13=0, weight14=0, weight23=0, weight24=0, weight34=0; Float_t weight12Err=0, weight13Err=0, weight14Err=0, weight23Err=0, weight24Err=0, weight34Err=0; - Float_t weight12CC=0, weight13CC=0, weight14CC=0, weight23CC=0, weight24CC=0, weight34CC=0; + Float_t weight12CC[3]={0}; + Float_t weight13CC[3]={0}; + Float_t weight14CC[3]={0}; + Float_t weight23CC[3]={0}; + Float_t weight24CC[3]={0}; + Float_t weight34CC[3]={0}; Float_t weightTotal=0;//, weightTotalErr=0; Float_t qinv12MC=0, qinv13MC=0, qinv14MC=0, qinv23MC=0, qinv24MC=0, qinv34MC=0; Float_t parentQinv12=0, parentQinv13=0, parentQinv14=0, parentQinv23=0, parentQinv24=0, parentQinv34=0; @@ -1875,26 +1889,7 @@ void AliFourPion::UserExec(Option_t *) fNormQPairSwitch_E1E3[i]->Set(kMultLimitPbPb,fDefaultsCharSwitch); fNormQPairSwitch_E2E3[i]->Set(kMultLimitPbPb,fDefaultsCharSwitch); } - /*for(Int_t i=0; ifNtracks; i++) {// 1st particle for (Int_t j=i+1; j<(fEvt+en2)->fNtracks; j++) {// 2nd particle - + pVect1[0]=(fEvt+en1)->fTracks[i].fEaccepted; pVect2[0]=(fEvt+en2)->fTracks[j].fEaccepted; pVect1[1]=(fEvt+en1)->fTracks[i].fP[0]; pVect2[1]=(fEvt+en2)->fTracks[j].fP[0]; pVect1[2]=(fEvt+en1)->fTracks[i].fP[1]; pVect2[2]=(fEvt+en2)->fTracks[j].fP[1]; @@ -1987,13 +1982,6 @@ void AliFourPion::UserExec(Option_t *) ////////////////////////////////////////////////////////////////////////////// if(qinv12 <= fQcut) { - /*if(en1==0 && en2==0) {LowQPairSwitch_E0E0[i][j] = kTRUE; PairCount[0]++;} - if(en1==0 && en2==1) {LowQPairSwitch_E0E1[i][j] = kTRUE; PairCount[1]++;} - if(en1==0 && en2==2) {LowQPairSwitch_E0E2[i][j] = kTRUE; PairCount[2]++;} - if(en1==0 && en2==3) {LowQPairSwitch_E0E3[i][j] = kTRUE; PairCount[3]++;} - if(en1==1 && en2==2) {LowQPairSwitch_E1E2[i][j] = kTRUE; PairCount[4]++;} - if(en1==1 && en2==3) {LowQPairSwitch_E1E3[i][j] = kTRUE; PairCount[5]++;} - if(en1==2 && en2==3) {LowQPairSwitch_E2E3[i][j] = kTRUE; PairCount[6]++;}*/ if(en1==0 && en2==0) {fLowQPairSwitch_E0E0[i]->AddAt('1',j);} if(en1==0 && en2==1) {fLowQPairSwitch_E0E1[i]->AddAt('1',j);} if(en1==0 && en2==2) {fLowQPairSwitch_E0E2[i]->AddAt('1',j);} @@ -2004,13 +1992,6 @@ void AliFourPion::UserExec(Option_t *) if(en1==2 && en2==3) {fLowQPairSwitch_E2E3[i]->AddAt('1',j);} } if((qinv12 >= fNormQcutLow) && (qinv12 < fNormQcutHigh)) { - /*if(en1==0 && en2==0) {NormQPairSwitch_E0E0[i][j] = kTRUE; NormPairCount[0]++;} - if(en1==0 && en2==1) {NormQPairSwitch_E0E1[i][j] = kTRUE; NormPairCount[1]++;} - if(en1==0 && en2==2) {NormQPairSwitch_E0E2[i][j] = kTRUE; NormPairCount[2]++;} - if(en1==0 && en2==3) {NormQPairSwitch_E0E3[i][j] = kTRUE; NormPairCount[3]++;} - if(en1==1 && en2==2) {NormQPairSwitch_E1E2[i][j] = kTRUE; NormPairCount[4]++;} - if(en1==1 && en2==3) {NormQPairSwitch_E1E3[i][j] = kTRUE; NormPairCount[5]++;} - if(en1==2 && en2==3) {NormQPairSwitch_E2E3[i][j] = kTRUE; NormPairCount[6]++;}*/ if(en1==0 && en2==0) {fNormQPairSwitch_E0E0[i]->AddAt('1',j);} if(en1==0 && en2==1) {fNormQPairSwitch_E0E1[i]->AddAt('1',j);} if(en1==0 && en2==2) {fNormQPairSwitch_E0E2[i]->AddAt('1',j);} @@ -2059,8 +2040,6 @@ void AliFourPion::UserExec(Option_t *) ch1 = Int_t(((fEvt)->fTracks[i].fCharge + 1)/2.); for (Int_t j=i+1; j<(fEvt+en2)->fNtracks; j++) {// 2nd particle - //if(en2==0) {if(!NormQPairSwitch_E0E0[i][j]) continue;} - //else {if(!NormQPairSwitch_E0E1[i][j]) continue;} if(en2==0) {if(fNormQPairSwitch_E0E0[i]->At(j)=='0') continue;} else {if(fNormQPairSwitch_E0E1[i]->At(j)=='0') continue;} @@ -2070,16 +2049,6 @@ void AliFourPion::UserExec(Option_t *) ch2 = Int_t(((fEvt+en2)->fTracks[j].fCharge + 1)/2.); for (Int_t k=j+1; k<(fEvt+en3)->fNtracks; k++) {// 3rd particle - /*if(en3==0) { - if(!NormQPairSwitch_E0E0[i][k]) continue; - if(!NormQPairSwitch_E0E0[j][k]) continue; - }else if(en3==1){ - if(!NormQPairSwitch_E0E1[i][k]) continue; - if(!NormQPairSwitch_E0E1[j][k]) continue; - }else{ - if(!NormQPairSwitch_E0E2[i][k]) continue; - if(!NormQPairSwitch_E1E2[j][k]) continue; - }*/ if(en3==0) { if(fNormQPairSwitch_E0E0[i]->At(k)=='0') continue; if(fNormQPairSwitch_E0E0[j]->At(k)=='0') continue; @@ -2112,23 +2081,6 @@ void AliFourPion::UserExec(Option_t *) for (Int_t l=k+1; l<(fEvt+en4)->fNtracks; l++) {// 4th particle - /*if(en4==0){ - if(!NormQPairSwitch_E0E0[i][l]) continue; - if(!NormQPairSwitch_E0E0[j][l]) continue; - if(!NormQPairSwitch_E0E0[k][l]) continue; - }else if(en4==1){ - if(!NormQPairSwitch_E0E1[i][l]) continue; - if(!NormQPairSwitch_E0E1[j][l]) continue; - if(!NormQPairSwitch_E0E1[k][l]) continue; - }else if(en4==2){ - if(!NormQPairSwitch_E0E2[i][l]) continue; - if(!NormQPairSwitch_E0E2[j][l]) continue; - if(!NormQPairSwitch_E1E2[k][l]) continue; - }else{ - if(!NormQPairSwitch_E0E3[i][l]) continue; - if(!NormQPairSwitch_E1E3[j][l]) continue; - if(!NormQPairSwitch_E2E3[k][l]) continue; - }*/ if(en4==0){ if(fNormQPairSwitch_E0E0[i]->At(l)=='0') continue; if(fNormQPairSwitch_E0E0[j]->At(l)=='0') continue; @@ -2216,8 +2168,6 @@ void AliFourPion::UserExec(Option_t *) ///////////////////////////////////////////////////////////// for (Int_t j=i+1; j<(fEvt+en2)->fNtracks; j++) {// 2nd particle - //if(en2==0) {if(!LowQPairSwitch_E0E0[i][j]) continue;} - //else {if(!LowQPairSwitch_E0E1[i][j]) continue;} if(en2==0) {if(fLowQPairSwitch_E0E0[i]->At(j)=='0') continue;} else {if(fLowQPairSwitch_E0E1[i]->At(j)=='0') continue;} @@ -2372,8 +2322,10 @@ void AliFourPion::UserExec(Option_t *) Float_t WInput = 1.0; if(term==1) { WInput = MCWeight(chGroup2, Rvalue, 1.0, parentQinv12, 0.); + }else{ muonPionK12 = 1.0; pionPionK12=1.0; } + Charge1[bin1].Charge2[bin2].MB[0].EDB[0].TwoPT[term-1].fMuonSmeared->Fill(Rvalue, qinv12MC, WInput); Charge1[bin1].Charge2[bin2].MB[0].EDB[0].TwoPT[term-1].fMuonIdeal->Fill(Rvalue, parentQinv12, WInput); Charge1[bin1].Charge2[bin2].MB[0].EDB[0].TwoPT[term-1].fMuonPionK2->Fill(Rvalue, qinv12MC, muonPionK12); @@ -2406,16 +2358,6 @@ void AliFourPion::UserExec(Option_t *) ///////////////////////////////////////////////////////////// for (Int_t k=j+1; k<(fEvt+en3)->fNtracks; k++) {// 3rd particle - /*if(en3==0) { - if(!LowQPairSwitch_E0E0[i][k]) continue; - if(!LowQPairSwitch_E0E0[j][k]) continue; - }else if(en3==1){ - if(!LowQPairSwitch_E0E1[i][k]) continue; - if(!LowQPairSwitch_E0E1[j][k]) continue; - }else{ - if(!LowQPairSwitch_E0E2[i][k]) continue; - if(!LowQPairSwitch_E1E2[j][k]) continue; - }*/ if(en3==0) { if(fLowQPairSwitch_E0E0[i]->At(k)=='0') continue; if(fLowQPairSwitch_E0E0[j]->At(k)=='0') continue; @@ -2498,6 +2440,7 @@ void AliFourPion::UserExec(Option_t *) }else{ Float_t MomResCorr12=1.0, MomResCorr13=1.0, MomResCorr23=1.0; + Float_t MuonCorr12=1.0, MuonCorr13=1.0, MuonCorr23=1.0; if(!fGenerateSignal && !fMCcase) { Int_t momBin12 = fMomResC2->GetYaxis()->FindBin(qinv12); Int_t momBin13 = fMomResC2->GetYaxis()->FindBin(qinv13); @@ -2508,28 +2451,55 @@ void AliFourPion::UserExec(Option_t *) MomResCorr12 = fMomResC2->GetBinContent(rBinForTPNMomRes, momBin12); MomResCorr13 = fMomResC2->GetBinContent(rBinForTPNMomRes, momBin13); MomResCorr23 = fMomResC2->GetBinContent(rBinForTPNMomRes, momBin23); + MuonCorr12 = fWeightmuonCorrection->GetBinContent(rBinForTPNMomRes, momBin12); + MuonCorr13 = fWeightmuonCorrection->GetBinContent(rBinForTPNMomRes, momBin13); + MuonCorr23 = fWeightmuonCorrection->GetBinContent(rBinForTPNMomRes, momBin23); } - - weight12CC = ((weight12+1)*MomResCorr12 - ffcSq*FSICorr12 - (1-ffcSq)); - weight12CC /= FSICorr12*ffcSq; - weight13CC = ((weight13+1)*MomResCorr13 - ffcSq*FSICorr13 - (1-ffcSq)); - weight13CC /= FSICorr13*ffcSq; - weight23CC = ((weight23+1)*MomResCorr23 - ffcSq*FSICorr23 - (1-ffcSq)); - weight23CC /= FSICorr23*ffcSq; - - if(weight12CC < 0 || weight13CC < 0 || weight23CC < 0) {// C2^QS can never be less than unity + // no MRC, no Muon Correction + weight12CC[0] = ((weight12+1) - ffcSq*FSICorr12 - (1-ffcSq)); + weight12CC[0] /= FSICorr12*ffcSq; + weight13CC[0] = ((weight13+1) - ffcSq*FSICorr13 - (1-ffcSq)); + weight13CC[0] /= FSICorr13*ffcSq; + weight23CC[0] = ((weight23+1) - ffcSq*FSICorr23 - (1-ffcSq)); + weight23CC[0] /= FSICorr23*ffcSq; + if(weight12CC[0] > 0 && weight13CC[0] > 0 && weight23CC[0] > 0){ + Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[4].fTwoPartNorm->Fill(1, q3, sqrt(weight12CC[0]*weight13CC[0]*weight23CC[0])); + } + // no Muon Correction + weight12CC[1] = ((weight12+1)*MomResCorr12 - ffcSq*FSICorr12 - (1-ffcSq)); + weight12CC[1] /= FSICorr12*ffcSq; + weight13CC[1] = ((weight13+1)*MomResCorr13 - ffcSq*FSICorr13 - (1-ffcSq)); + weight13CC[1] /= FSICorr13*ffcSq; + weight23CC[1] = ((weight23+1)*MomResCorr23 - ffcSq*FSICorr23 - (1-ffcSq)); + weight23CC[1] /= FSICorr23*ffcSq; + if(weight12CC[1] > 0 && weight13CC[1] > 0 && weight23CC[1] > 0){ + Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[4].fTwoPartNorm->Fill(2, q3, sqrt(weight12CC[1]*weight13CC[1]*weight23CC[1])); + } + // both Corrections + weight12CC[2] = ((weight12+1)*MomResCorr12 - ffcSq*FSICorr12 - (1-ffcSq)); + weight12CC[2] /= FSICorr12*ffcSq; + weight12CC[2] *= MuonCorr12; + weight13CC[2] = ((weight13+1)*MomResCorr13 - ffcSq*FSICorr13 - (1-ffcSq)); + weight13CC[2] /= FSICorr13*ffcSq; + weight13CC[2] *= MuonCorr13; + weight23CC[2] = ((weight23+1)*MomResCorr23 - ffcSq*FSICorr23 - (1-ffcSq)); + weight23CC[2] /= FSICorr23*ffcSq; + weight23CC[2] *= MuonCorr23; + if(weight12CC[2] < 0 || weight13CC[2] < 0 || weight23CC[2] < 0) {// C2^QS can never be less than unity if(fMbin==0 && bin1==0) { - ((TH1D*)fOutputList->FindObject("fTPNRejects3pion2"))->Fill(q3, sqrt(fabs(weight12CC*weight13CC*weight23CC))); + ((TH1D*)fOutputList->FindObject("fTPNRejects3pion2"))->Fill(q3, sqrt(fabs(weight12CC[2]*weight13CC[2]*weight23CC[2]))); } }else{ GoodTripletWeights = kTRUE; ///////////////////////////////////////////////////// - weightTotal = sqrt(weight12CC*weight13CC*weight23CC); + weightTotal = sqrt(weight12CC[2]*weight13CC[2]*weight23CC[2]); ///////////////////////////////////////////////////// - Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[4].fTwoPartNorm->Fill(1, q3, weightTotal); - }// 2nd r3 den check + Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[4].fTwoPartNorm->Fill(3, q3, weightTotal); + }// 2nd r3 den check + + }// 1st r3 den check - + }// r3 den @@ -2599,16 +2569,25 @@ void AliFourPion::UserExec(Option_t *) ((TH2D*)fOutputList->FindObject("fAvgQ13VersusQ3"))->Fill(parentQ3, parentQinv13); ((TH2D*)fOutputList->FindObject("fAvgQ23VersusQ3"))->Fill(parentQ3, parentQinv23); - for(Int_t term=1; term<=5; term++){ + for(Int_t term=1; term<=4; term++){ + if(term==1) {} + else if(term==2) {if(!pionParent1 && !pionParent2) continue;} + else if(term==3) {if(!pionParent1 && !pionParent3) continue;} + else {if(!pionParent2 && !pionParent3) continue;} for(Int_t Riter=0; RiterFill(Rvalue, q3MC, WInput); - Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[0].EDB[0].ThreePT[term-1].fMuonIdeal->Fill(Rvalue, parentQ3, WInput); - Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[0].EDB[0].ThreePT[term-1].fMuonPionK3->Fill(Rvalue, q3MC, WInputFSI); - Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[0].EDB[0].ThreePT[term-1].fPionPionK3->Fill(Rvalue, parentQ3, WInputParentFSI); + Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[0].EDB[0].ThreePT[term-1].fMuonSmeared->Fill(1, Rvalue, q3MC, WInput); + Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[0].EDB[0].ThreePT[term-1].fMuonIdeal->Fill(1, Rvalue, parentQ3, WInput); + Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[0].EDB[0].ThreePT[term-1].fMuonPionK3->Fill(1, Rvalue, q3MC, WInputFSI); + Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[0].EDB[0].ThreePT[term-1].fPionPionK3->Fill(1, Rvalue, parentQ3, WInputParentFSI); + // + Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[0].EDB[0].ThreePT[term-1].fMuonSmeared->Fill(2, Rvalue, q3MC); + Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[0].EDB[0].ThreePT[term-1].fMuonIdeal->Fill(2, Rvalue, parentQ3); + Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[0].EDB[0].ThreePT[term-1].fMuonPionK3->Fill(2, Rvalue, q3MC); + Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[0].EDB[0].ThreePT[term-1].fPionPionK3->Fill(2, Rvalue, parentQ3); }// Riter }// term loop @@ -2628,28 +2607,11 @@ void AliFourPion::UserExec(Option_t *) }// 3rd particle label check }// MCcase and ENsum==6 - + - ///////////////////////////////////////////////////////////// + ///////////////////////////////////////////////////////////// for (Int_t l=k+1; l<(fEvt+en4)->fNtracks; l++) {// 4th particle - /*if(en4==0){ - if(!LowQPairSwitch_E0E0[i][l]) continue; - if(!LowQPairSwitch_E0E0[j][l]) continue; - if(!LowQPairSwitch_E0E0[k][l]) continue; - }else if(en4==1){ - if(!LowQPairSwitch_E0E1[i][l]) continue; - if(!LowQPairSwitch_E0E1[j][l]) continue; - if(!LowQPairSwitch_E0E1[k][l]) continue; - }else if(en4==2){ - if(!LowQPairSwitch_E0E2[i][l]) continue; - if(!LowQPairSwitch_E0E2[j][l]) continue; - if(!LowQPairSwitch_E1E2[k][l]) continue; - }else{ - if(!LowQPairSwitch_E0E3[i][l]) continue; - if(!LowQPairSwitch_E1E3[j][l]) continue; - if(!LowQPairSwitch_E2E3[k][l]) continue; - }*/ if(en4==0){ if(fLowQPairSwitch_E0E0[i]->At(l)=='0') continue; if(fLowQPairSwitch_E0E0[j]->At(l)=='0') continue; @@ -2722,6 +2684,7 @@ void AliFourPion::UserExec(Option_t *) GetWeight(pVect3, pVect4, weight34, weight34Err); Float_t MomResCorr14=1.0, MomResCorr24=1.0, MomResCorr34=1.0; + Float_t MuonCorr14=1.0, MuonCorr24=1.0, MuonCorr34=1.0; if(!fGenerateSignal && !fMCcase) { Int_t momBin14 = fMomResC2->GetYaxis()->FindBin(qinv14); Int_t momBin24 = fMomResC2->GetYaxis()->FindBin(qinv24); @@ -2729,28 +2692,65 @@ void AliFourPion::UserExec(Option_t *) if(momBin14 >= 20) momBin14 = 19; if(momBin24 >= 20) momBin24 = 19; if(momBin34 >= 20) momBin34 = 19; - MomResCorr14 = fMomResC2->GetBinContent(rBinForTPNMomRes, momBin14); - MomResCorr24 = fMomResC2->GetBinContent(rBinForTPNMomRes, momBin24); - MomResCorr34 = fMomResC2->GetBinContent(rBinForTPNMomRes, momBin34); + MomResCorr14 = fMomResC2->GetBinContent(rBinForTPNMomRes, momBin14); + MomResCorr24 = fMomResC2->GetBinContent(rBinForTPNMomRes, momBin24); + MomResCorr34 = fMomResC2->GetBinContent(rBinForTPNMomRes, momBin34); + MuonCorr14 = fWeightmuonCorrection->GetBinContent(rBinForTPNMomRes, momBin14); + MuonCorr24 = fWeightmuonCorrection->GetBinContent(rBinForTPNMomRes, momBin24); + MuonCorr34 = fWeightmuonCorrection->GetBinContent(rBinForTPNMomRes, momBin34); } - weight14CC = ((weight14+1)*MomResCorr14 - ffcSq*FSICorr14 - (1-ffcSq)); - weight14CC /= FSICorr14*ffcSq; - weight24CC = ((weight24+1)*MomResCorr24 - ffcSq*FSICorr24 - (1-ffcSq)); - weight24CC /= FSICorr24*ffcSq; - weight34CC = ((weight34+1)*MomResCorr34 - ffcSq*FSICorr34 - (1-ffcSq)); - weight34CC /= FSICorr34*ffcSq; - - if(weight14CC < 0 || weight24CC < 0 || weight34CC < 0) {// C2^QS can never be less than unity - if(fMbin==0 && bin1==0) { - ((TH1D*)fOutputList->FindObject("fTPNRejects4pion1"))->Fill(q4, sqrt(fabs(weight12CC*weight23CC*weight34CC*weight14CC))); + + // no MRC, no Muon Correction + weight14CC[0] = ((weight14+1) - ffcSq*FSICorr14 - (1-ffcSq)); + weight14CC[0] /= FSICorr14*ffcSq; + weight24CC[0] = ((weight24+1) - ffcSq*FSICorr24 - (1-ffcSq)); + weight24CC[0] /= FSICorr24*ffcSq; + weight34CC[0] = ((weight34+1) - ffcSq*FSICorr34 - (1-ffcSq)); + weight34CC[0] /= FSICorr34*ffcSq; + if(weight14CC[0] > 0 && weight24CC[0] > 0 && weight34CC[0] > 0 && weight12CC[0] > 0 && weight13CC[0] > 0 && weight23CC[0] > 0){ + weightTotal = sqrt(weight12CC[0]*weight13CC[0]*weight24CC[0]*weight34CC[0]); + weightTotal += sqrt(weight12CC[0]*weight14CC[0]*weight23CC[0]*weight34CC[0]); + weightTotal += sqrt(weight13CC[0]*weight14CC[0]*weight23CC[0]*weight24CC[0]); + weightTotal /= 3.; + Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[fMbin].EDB[KT4index].FourPT[12].fTwoPartNorm->Fill(1, q4, weightTotal); + } + // no Muon Correction + weight14CC[1] = ((weight14+1)*MomResCorr14 - ffcSq*FSICorr14 - (1-ffcSq)); + weight14CC[1] /= FSICorr14*ffcSq; + weight24CC[1] = ((weight24+1)*MomResCorr24 - ffcSq*FSICorr24 - (1-ffcSq)); + weight24CC[1] /= FSICorr24*ffcSq; + weight34CC[1] = ((weight34+1)*MomResCorr34 - ffcSq*FSICorr34 - (1-ffcSq)); + weight34CC[1] /= FSICorr34*ffcSq; + if(weight14CC[1] > 0 && weight24CC[1] > 0 && weight34CC[1] > 0 && weight12CC[1] > 0 && weight13CC[1] > 0 && weight23CC[1] > 0){ + weightTotal = sqrt(weight12CC[1]*weight13CC[1]*weight24CC[1]*weight34CC[1]); + weightTotal += sqrt(weight12CC[1]*weight14CC[1]*weight23CC[1]*weight34CC[1]); + weightTotal += sqrt(weight13CC[1]*weight14CC[1]*weight23CC[1]*weight24CC[1]); + weightTotal /= 3.; + Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[fMbin].EDB[KT4index].FourPT[12].fTwoPartNorm->Fill(2, q4, weightTotal); + } + // both corrections + weight14CC[2] = ((weight14+1)*MomResCorr14 - ffcSq*FSICorr14 - (1-ffcSq)); + weight14CC[2] /= FSICorr14*ffcSq; + weight14CC[2] *= MuonCorr14; + weight24CC[2] = ((weight24+1)*MomResCorr24 - ffcSq*FSICorr24 - (1-ffcSq)); + weight24CC[2] /= FSICorr24*ffcSq; + weight24CC[2] *= MuonCorr24; + weight34CC[2] = ((weight34+1)*MomResCorr34 - ffcSq*FSICorr34 - (1-ffcSq)); + weight34CC[2] /= FSICorr34*ffcSq; + weight34CC[2] *= MuonCorr34; + + if(weight14CC[2] < 0 || weight24CC[2] < 0 || weight34CC[2] < 0) {// C2^QS can never be less than unity + if(fMbin==0 && bin1==0 && KT4index==0) { + ((TH1D*)fOutputList->FindObject("fTPNRejects4pion1"))->Fill(q4, sqrt(fabs(weight12CC[2]*weight23CC[2]*weight34CC[2]*weight14CC[2]))); } }else{ ///////////////////////////////////////////////////// - weightTotal = sqrt(weight12CC*weight13CC*weight24CC*weight34CC) + sqrt(weight12CC*weight14CC*weight23CC*weight34CC); - weightTotal += sqrt(weight13CC*weight14CC*weight23CC*weight24CC); + weightTotal = sqrt(weight12CC[2]*weight13CC[2]*weight24CC[2]*weight34CC[2]); + weightTotal += sqrt(weight12CC[2]*weight14CC[2]*weight23CC[2]*weight34CC[2]); + weightTotal += sqrt(weight13CC[2]*weight14CC[2]*weight23CC[2]*weight24CC[2]); weightTotal /= 3.; ///////////////////////////////////////////////////// - Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[fMbin].EDB[KT3index].FourPT[12].fTwoPartNorm->Fill(1, q4, weightTotal); + Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[fMbin].EDB[KT4index].FourPT[12].fTwoPartNorm->Fill(3, q4, weightTotal); } } ///////////////////////////////////////////////////////////// @@ -2786,8 +2786,8 @@ void AliFourPion::UserExec(Option_t *) q4MC = sqrt(pow(q3MC,2) + pow(qinv14MC,2) + pow(qinv24MC,2) + pow(qinv34MC,2)); if(q4<0.1 && ch1==ch2 && ch1==ch3 && ch1==ch4) ((TH2D*)fOutputList->FindObject("fQ4Res"))->Fill(KT4, q4-q4MC); if(ch1==ch2 && ch1==ch3 && ch1==ch4) { - ((TH2D*)fOutputList->FindObject("DistQinvMC4pion"))->Fill(1, qinv12MC); ((TH2D*)fOutputList->FindObject("DistQinvMC4pion"))->Fill(2, qinv13MC); - ((TH2D*)fOutputList->FindObject("DistQinvMC4pion"))->Fill(3, qinv14MC); ((TH2D*)fOutputList->FindObject("DistQinvMC4pion"))->Fill(4, qinv23MC); + ((TH2D*)fOutputList->FindObject("DistQinvMC4pion"))->Fill(1, qinv12MC); ((TH2D*)fOutputList->FindObject("DistQinvMC4pion"))->Fill(2, qinv13MC); + ((TH2D*)fOutputList->FindObject("DistQinvMC4pion"))->Fill(3, qinv14MC); ((TH2D*)fOutputList->FindObject("DistQinvMC4pion"))->Fill(4, qinv23MC); ((TH2D*)fOutputList->FindObject("DistQinvMC4pion"))->Fill(5, qinv24MC); ((TH2D*)fOutputList->FindObject("DistQinvMC4pion"))->Fill(6, qinv34MC); } Int_t chGroup4[4]={ch1,ch2,ch3,ch4}; @@ -2813,19 +2813,40 @@ void AliFourPion::UserExec(Option_t *) if(parentQinv14 > 0.001 && parentQinv24 > 0.001 && parentQinv34 > 0.001 && parentQ4 < 0.5){ if(pionParent1 || pionParent2 || pionParent3 || pionParent4) {// want at least one pion-->muon + if(pionParent1) ((TH1D*)fOutputList->FindObject("fDistPionParents4"))->Fill(1); + if(pionParent2) ((TH1D*)fOutputList->FindObject("fDistPionParents4"))->Fill(2); + if(pionParent3) ((TH1D*)fOutputList->FindObject("fDistPionParents4"))->Fill(3); + if(pionParent4) ((TH1D*)fOutputList->FindObject("fDistPionParents4"))->Fill(4); Float_t parentQinvGroup4[6]={parentQinv12, parentQinv13, parentQinv14, parentQinv23, parentQinv24, parentQinv34}; Float_t parentkTGroup4[6]={0}; - for(Int_t term=1; term<=13; term++){ + for(Int_t term=1; term<=12; term++){ + if(term==1) {} + else if(term==2) {if(!pionParent1 && !pionParent2 && !pionParent3) continue;} + else if(term==3) {if(!pionParent1 && !pionParent2 && !pionParent4) continue;} + else if(term==4) {if(!pionParent1 && !pionParent3 && !pionParent4) continue;} + else if(term==5) {if(!pionParent2 && !pionParent3 && !pionParent4) continue;} + else if(term==6) {if(!pionParent1 && !pionParent2) continue;} + else if(term==7) {if(!pionParent1 && !pionParent3) continue;} + else if(term==8) {if(!pionParent1 && !pionParent4) continue;} + else if(term==9) {if(!pionParent2 && !pionParent3) continue;} + else if(term==10) {if(!pionParent2 && !pionParent4) continue;} + else if(term==11) {if(!pionParent3 && !pionParent4) continue;} + else {} for(Int_t Riter=0; RiterFill(Rvalue, q4MC, WInput); - Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[0].EDB[0].FourPT[term-1].fMuonIdeal->Fill(Rvalue, parentQ4, WInput); - Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[0].EDB[0].FourPT[term-1].fMuonPionK4->Fill(Rvalue, q4MC, WInputFSI); - Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[0].EDB[0].FourPT[term-1].fPionPionK4->Fill(Rvalue, parentQ4, WInputParentFSI); + Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[0].EDB[0].FourPT[term-1].fMuonSmeared->Fill(1, Rvalue, q4MC, WInput); + Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[0].EDB[0].FourPT[term-1].fMuonIdeal->Fill(1, Rvalue, parentQ4, WInput); + Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[0].EDB[0].FourPT[term-1].fMuonPionK4->Fill(1, Rvalue, q4MC, WInputFSI); + Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[0].EDB[0].FourPT[term-1].fPionPionK4->Fill(1, Rvalue, parentQ4, WInputParentFSI); + // + Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[0].EDB[0].FourPT[term-1].fMuonSmeared->Fill(2, Rvalue, q4MC); + Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[0].EDB[0].FourPT[term-1].fMuonIdeal->Fill(2, Rvalue, parentQ4); + Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[0].EDB[0].FourPT[term-1].fMuonPionK4->Fill(2, Rvalue, q4MC); + Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[0].EDB[0].FourPT[term-1].fPionPionK4->Fill(2, Rvalue, parentQ4); }// Riter }// term loop @@ -3393,8 +3414,8 @@ Float_t AliFourPion::MCWeight4(Int_t term, Float_t R, Float_t fcSq, Int_t c[4], else if(term==10) return ((1-fcSq) + fcSq*(1 + pow(EA24,2))*Kfactor24); else return ((1-fcSq) + fcSq*(1 + pow(EA34,2))*Kfactor34); }else if(term==12){ - Float_t C22 = ((1-fcSq) + fcSq*(1 + pow(EA12,2))*Kfactor12) - 1.0; - C22 *= ((1-fcSq) + fcSq*(1 + pow(EA34,2))*Kfactor34) - 1.0; + Float_t C22 = (1-fcSq) + fcSq*(1 + pow(EA12,2))*Kfactor12; + C22 *= (1-fcSq) + fcSq*(1 + pow(EA34,2))*Kfactor34; return C22; }else return 1.0; @@ -3520,8 +3541,8 @@ Float_t AliFourPion::MCWeightFSI4(Int_t term, Float_t R, Float_t fcSq, Int_t c[4 else if(term==10) return ((1-fcSq) + fcSq*Kfactor24); else return ((1-fcSq) + fcSq*Kfactor34); }else if(term==12){ - Float_t C22 = ((1-fcSq) + fcSq*Kfactor12) - 1.0; - C22 *= ((1-fcSq) + fcSq*Kfactor34) - 1.0; + Float_t C22 = (1-fcSq) + fcSq*Kfactor12; + C22 *= (1-fcSq) + fcSq*Kfactor34; return C22; }else return 1.0; @@ -3816,3 +3837,24 @@ void AliFourPion::SetFSIindex(Float_t R){ else fFSIindex = 9; } } +//________________________________________________________________________ +void AliFourPion::SetMuonCorrections(Bool_t legoCase, TH2D *tempMuon){ + if(legoCase){ + cout<<"LEGO call to SetMuonCorrections"<Clone(); + fWeightmuonCorrection->SetDirectory(0); + }else { + cout<<"non LEGO call to SetMuonCorrections"<IsOpen()) { + cout<<"No Muon file found"<Get("WeightmuonCorrection"); + fWeightmuonCorrection->SetDirectory(0); + // + MuonFile->Close(); + } + cout<<"Done reading Muon file"<