code updates, Add Muon Corrections, KT4index bug fix

diff --git a/PWGCF/FEMTOSCOPY/Chaoticity/AliFourPion.cxx b/PWGCF/FEMTOSCOPY/Chaoticity/AliFourPion.cxx
index 4b4449880cacaf9f96c61972ed44b9d23dca40c6..2e4dc0bc4b9fc6e57db9bbd95d413f4316fc84a3 100755 (executable)
--- 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; mb<fMbins; mb++){
@@ -313,7 +314,8 @@ AliFourPion::AliFourPion(const Char_t *name)
   fNormQPairSwitch_E1E2(),
   fNormQPairSwitch_E1E3(),
   fNormQPairSwitch_E2E3(),
-  fMomResC2(0x0)
+  fMomResC2(0x0),
+  fWeightmuonCorrection(0x0)
 {
   // Main constructor
   fAODcase=kTRUE;
@@ -488,7 +490,8 @@ AliFourPion::AliFourPion(const AliFourPion &obj)
     fNormQPairSwitch_E1E2(),
     fNormQPairSwitch_E1E3(),
     fNormQPairSwitch_E2E3(),
-    fMomResC2(obj.fMomResC2)
+    fMomResC2(obj.fMomResC2),
+    fWeightmuonCorrection(obj.fWeightmuonCorrection)
 {
   // Copy Constructor
   
@@ -583,6 +586,7 @@ AliFourPion &AliFourPion::operator=(const AliFourPion &obj)
   fKT3transition = obj.fKT3transition;
   fKT4transition = obj.fKT4transition;
   fMomResC2 = obj.fMomResC2;
+  fWeightmuonCorrection = obj.fWeightmuonCorrection;
   
   for(Int_t i=0; i<12; i++){
     fFSIss[i]=obj.fFSIss[i]; 
@@ -609,7 +613,8 @@ AliFourPion::~AliFourPion()
   if(fTempStruct) delete [] fTempStruct;
   if(fRandomNumber) delete fRandomNumber;
   if(fMomResC2) delete fMomResC2;
- 
+  if(fWeightmuonCorrection) delete fWeightmuonCorrection;
+
   for(Int_t j=0; j<kMultLimitPbPb; j++){
     if(fLowQPairSwitch_E0E0[j]) delete [] fLowQPairSwitch_E0E0[j];
     if(fLowQPairSwitch_E0E1[j]) delete [] fLowQPairSwitch_E0E1[j];
@@ -769,8 +774,8 @@ void AliFourPion::ParInit()
   //
   fQupperBoundWeights = 0.2;
   fQupperBoundQ2 = 2.0;
-  fQupperBoundQ3 = 0.5;
-  fQupperBoundQ4 = 0.5;
+  fQupperBoundQ3 = 0.6;
+  fQupperBoundQ4 = 0.6;
   fQbinsQ2 = fQupperBoundQ2/0.005;
   fQbinsQ3 = fQupperBoundQ3/0.005;
   fQbinsQ4 = fQupperBoundQ4/0.005;
@@ -828,6 +833,7 @@ void AliFourPion::ParInit()
     SetFSICorrelations(fLEGO);// Read in 2-particle and 3-particle FSI correlations
     if(!fTabulatePairs) SetWeightArrays(fLEGO);// Set Weight Array
     if(!fMCcase && !fTabulatePairs) SetMomResCorrections(fLEGO);// Read Momentum resolution file
+    if(!fMCcase && !fTabulatePairs) SetMuonCorrections(fLEGO);// Read Muon corrections
   }
   
   /////////////////////////////////////////////
@@ -1169,22 +1175,22 @@ void AliFourPion::UserCreateOutputObjects()
                  // Muon correction histos
                  TString *nameMuonIdeal=new TString(namePC3->Data());
                  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; i<fMultLimit; i++) {
-      for(Int_t j=0; j<fMultLimit; j++) {
-       LowQPairSwitch_E0E0[i][j] = kFALSE;
-       LowQPairSwitch_E0E1[i][j] = kFALSE;
-       LowQPairSwitch_E0E2[i][j] = kFALSE;
-       LowQPairSwitch_E0E3[i][j] = kFALSE;
-       LowQPairSwitch_E1E2[i][j] = kFALSE;
-       LowQPairSwitch_E1E3[i][j] = kFALSE;
-       LowQPairSwitch_E2E3[i][j] = kFALSE;
-       //
-       NormQPairSwitch_E0E0[i][j] = kFALSE;
-       NormQPairSwitch_E0E1[i][j] = kFALSE;
-       NormQPairSwitch_E0E2[i][j] = kFALSE;
-       NormQPairSwitch_E0E3[i][j] = kFALSE;
-       NormQPairSwitch_E1E2[i][j] = kFALSE;
-       NormQPairSwitch_E1E3[i][j] = kFALSE;
-       NormQPairSwitch_E2E3[i][j] = kFALSE;
-      }
-      }
-    */
+ 
   
   //////////////////////////////////////////
   // make low-q pair storage and normalization-pair storage
@@ -1906,7 +1901,7 @@ void AliFourPion::UserExec(Option_t *)
       for (Int_t i=0; i<(fEvt+en1)->fNtracks; 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; Riter<fRVALUES; Riter++){
                            Float_t Rvalue = 5+Riter;
                            Float_t WInput = MCWeight3(term, Rvalue, 1.0, chGroup3, parentQinvGroup3, parentkTGroup3);
                            Float_t WInputParentFSI = MCWeightFSI3(term, Rvalue, 1.0, chGroup3, parentQinvGroup3);
                            Float_t WInputFSI = MCWeightFSI3(term, Rvalue, 1.0, chGroup3, QinvMCGroup3);
-                           Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[0].EDB[0].ThreePT[term-1].fMuonSmeared->Fill(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; Riter<fRVALUES; Riter++){
                              Float_t Rvalue = 5+Riter;
                              Float_t WInput = MCWeight4(term, Rvalue, 1.0, chGroup4, parentQinvGroup4, parentkTGroup4);
                              Float_t WInputParentFSI = MCWeightFSI4(term, Rvalue, 1.0, chGroup4, parentQinvGroup4);
                              Float_t WInputFSI = MCWeightFSI4(term, Rvalue, 1.0, chGroup4, QinvMCGroup4);
-                             Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[0].EDB[0].FourPT[term-1].fMuonSmeared->Fill(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"<<endl;
+    fWeightmuonCorrection = (TH2D*)tempMuon->Clone();
+    fWeightmuonCorrection->SetDirectory(0);
+  }else {
+    cout<<"non LEGO call to SetMuonCorrections"<<endl;
+    TFile *MuonFile=new TFile("MuonCorrection.root","READ");
+    if(!MuonFile->IsOpen()) {
+      cout<<"No Muon file found"<<endl;
+      AliFatal("No Muon file found.  Kill process.");
+    }else {cout<<"Good Muon File Found!"<<endl;}
+    
+    fWeightmuonCorrection = (TH2D*)MuonFile->Get("WeightmuonCorrection");
+    fWeightmuonCorrection->SetDirectory(0);
+    //
+    MuonFile->Close();
+  }
+  cout<<"Done reading Muon file"<<endl;
+}

diff --git a/PWGCF/FEMTOSCOPY/Chaoticity/AliFourPion.h b/PWGCF/FEMTOSCOPY/Chaoticity/AliFourPion.h
index 58acd69fffec1b0292e14eaa0d557c139bb4b11b..f914136af02300b4fcae2d69110ddec23349c3bc 100755 (executable)
--- a/PWGCF/FEMTOSCOPY/Chaoticity/AliFourPion.h
+++ b/PWGCF/FEMTOSCOPY/Chaoticity/AliFourPion.h
@@ -73,6 +73,7 @@ class AliFourPion : public AliAnalysisTaskSE {
   void SetWeightArrays(Bool_t legoCase=kTRUE, TH3F *histos[AliFourPion::fKbinsT][AliFourPion::fCentBins]=0x0);
   void SetMomResCorrections(Bool_t legoCase=kTRUE, TH2D *temp2D=0x0);
   void SetFSICorrelations(Bool_t legoCase=kTRUE, TH1D *tempss[12]=0x0, TH1D *tempos[12]=0x0);
+  void SetMuonCorrections(Bool_t legoCase=kTRUE, TH2D *tempMuon=0x0);
   //
   void SetMCdecision(Bool_t mc) {fMCcase = mc;}
   void SetTabulatePairs(Bool_t tabulate) {fTabulatePairs = tabulate;}
@@ -142,10 +143,11 @@ class AliFourPion : public AliAnalysisTaskSE {
     TH2D *fIdeal; //!
     TH2D *fSmeared; //!
     //
-    TH2D *fMuonSmeared; //!
-    TH2D *fMuonIdeal; //!
-    TH2D *fMuonPionK3; //!
-    TH2D *fPionPionK3; //!
+    TH3D *fMuonSmeared; //!
+    TH3D *fMuonIdeal; //!
+    TH3D *fMuonPionK3; //!
+    TH3D *fPionPionK3; //!
+    //
     TH2D *fTwoPartNorm; //!
   };
   struct St7 {
@@ -178,10 +180,11 @@ class AliFourPion : public AliAnalysisTaskSE {
     TH2D *fIdeal; //!
     TH2D *fSmeared; //!
     //
-    TH2D *fMuonSmeared; //!
-    TH2D *fMuonIdeal; //!
-    TH2D *fMuonPionK4; //!
-    TH2D *fPionPionK4; //!
+    TH3D *fMuonSmeared; //!
+    TH3D *fMuonIdeal; //!
+    TH3D *fMuonPionK4; //!
+    TH3D *fPionPionK4; //!
+    //
     TH2D *fTwoPartNorm; //!
   };
   struct St_EDB {
@@ -332,10 +335,11 @@ class AliFourPion : public AliAnalysisTaskSE {
 
  public:
   TH2D *fMomResC2;
+  TH2D *fWeightmuonCorrection;
   TH1D *fFSIss[12];
   TH1D *fFSIos[12];
   TH3F *fNormWeight[fKbinsT][fCentBins];
-  
+ 
 
   ClassDef(AliFourPion, 1); 
 };