added histograms
authorsnelling <snelling@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 15 Apr 2010 14:30:28 +0000 (14:30 +0000)
committersnelling <snelling@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 15 Apr 2010 14:30:28 +0000 (14:30 +0000)
PWG2/FLOW/AliFlowCommon/AliFlowAnalysisWithScalarProduct.cxx
PWG2/FLOW/AliFlowCommon/AliFlowAnalysisWithScalarProduct.h
PWG2/FLOW/AliFlowTasks/AliAnalysisTaskScalarProduct.cxx
PWG2/FLOW/AliFlowTasks/AliAnalysisTaskScalarProduct.h
PWG2/FLOW/macros/AddTaskFlow.C

index e823adf..f76ae26 100644 (file)
@@ -21,7 +21,8 @@
 #include "TMath.h"
 #include "TProfile.h"
 #include "TVector2.h"
-#include "TH1F.h"
+#include "TH1D.h"
+#include "TH2D.h"
 
 #include "AliFlowCommonConstants.h"
 #include "AliFlowEventSimple.h"
@@ -45,6 +46,7 @@ ClassImp(AliFlowAnalysisWithScalarProduct)
   AliFlowAnalysisWithScalarProduct::AliFlowAnalysisWithScalarProduct():
    fEventNumber(0),
    fDebug(kFALSE),
+   fRelDiffMsub(1.),
    fWeightsList(NULL),
    fUsePhiWeights(kFALSE),
    fPhiWeights(NULL),
@@ -54,6 +56,7 @@ ClassImp(AliFlowAnalysisWithScalarProduct)
    fHistProUQPtRP(NULL),
    fHistProUQPtPOI(NULL),
    fHistProQaQb(NULL),
+   fHistProQaQbNorm(NULL),
    fHistSumOfLinearWeights(NULL),
    fHistSumOfQuadraticWeights(NULL),
    fHistProUQQaQbPtRP(NULL),
@@ -61,7 +64,17 @@ ClassImp(AliFlowAnalysisWithScalarProduct)
    fHistProUQQaQbPtPOI(NULL),
    fHistProUQQaQbEtaPOI(NULL),
    fCommonHists(NULL),
-   fCommonHistsRes(NULL)
+   fCommonHistsRes(NULL),
+   fHistQaQb(NULL),
+   fHistQaQbNorm(NULL),
+   fHistQaQbCos(NULL),
+   fHistResolution(NULL),
+   fHistQaNorm(NULL),
+   fHistQaNormvsMa(NULL),
+   fHistQbNorm(NULL),
+   fHistQbNormvsMb(NULL),
+   fHistMavsMb(NULL)
+
 {
   // Constructor.
   fWeightsList = new TList();
@@ -163,6 +176,10 @@ void AliFlowAnalysisWithScalarProduct::Init() {
   fHistProQaQb->SetYTitle("<QaQb>");
   fHistList->Add(fHistProQaQb);
 
+  fHistProQaQbNorm = new TProfile("FlowPro_QaQbNorm_SP","FlowPro_QaQbNorm_SP", 1, 0.5, 1.5,"s");
+  fHistProQaQbNorm->SetYTitle("<QaQb/MaMb>");
+  fHistList->Add(fHistProQaQbNorm);
+
   fHistSumOfLinearWeights = new TH1D("Flow_SumOfLinearWeights_SP","Flow_SumOfLinearWeights_SP",1,-0.5, 0.5);
   fHistSumOfLinearWeights -> SetYTitle("sum (*)");
   fHistSumOfLinearWeights -> SetXTitle("sum (Ma*Mb)");
@@ -226,6 +243,52 @@ void AliFlowAnalysisWithScalarProduct::Init() {
   fCommonHistsRes = new AliFlowCommonHistResults("AliFlowCommonHistResultsSP");
   fHistList->Add(fCommonHistsRes);  
 
+  fHistQaQb = new TH1D("Flow_QaQb_SP","Flow_QaQb_SP",200,-100.,100.);
+  fHistQaQb -> SetYTitle("dN/dQaQb");
+  fHistQaQb -> SetXTitle("QaQb");
+  fHistList->Add(fHistQaQb);
+
+  fHistQaQbNorm = new TH1D("Flow_QaQbNorm_SP","Flow_QaQbNorm_SP",44,-1.1,1.1);
+  fHistQaQbNorm -> SetYTitle("dN/d(QaQb/MaMb)");
+  fHistQaQbNorm -> SetXTitle("QaQb/MaMb");
+  fHistList->Add(fHistQaQbNorm);
+
+  fHistQaQbCos = new TH1D("Flow_QaQbCos_SP","Flow_QaQbCos_SP",63,0.,TMath::Pi());
+  fHistQaQbCos -> SetYTitle("dN/d(#phi)");
+  fHistQaQbCos -> SetXTitle("#phi");
+  fHistList->Add(fHistQaQbCos);
+
+  fHistResolution = new TH1D("Flow_resolution_SP","Flow_resolution_SP",100,-1.0,1.0);
+  fHistResolution -> SetYTitle("dN/d(cos(2(#phi_a - #phi_b))");
+  fHistResolution -> SetXTitle("cos(2*(#phi_a - #phi_b))");
+  fHistList->Add(fHistResolution);
+
+  fHistQaNorm = new TH1D("Flow_QaNorm_SP","Flow_QaNorm_SP",22,0.,1.1);
+  fHistQaNorm -> SetYTitle("dN/d(|Qa/Ma|)");
+  fHistQaNorm -> SetXTitle("|Qa/Ma|");
+  fHistList->Add(fHistQaNorm);
+
+  fHistQaNormvsMa = new TH2D("Flow_QaNormvsMa_SP","Flow_QaNormvsMa_SP",100,0.,100.,22,0.,1.1);
+  fHistQaNormvsMa -> SetYTitle("|Qa/Ma|");
+  fHistQaNormvsMa -> SetXTitle("Ma");
+  fHistList->Add(fHistQaNormvsMa);
+
+  fHistQbNorm = new TH1D("Flow_QbNorm_SP","Flow_QbNorm_SP",22,0.,1.1);
+  fHistQbNorm -> SetYTitle("dN/d(|Qb/Mb|)");
+  fHistQbNorm -> SetXTitle("|Qb/Mb|");
+  fHistList->Add(fHistQbNorm);
+
+  fHistQbNormvsMb = new TH2D("Flow_QbNormvsMb_SP","Flow_QbNormvsMb_SP",100,0.,100.,22,0.,1.1);
+  fHistQbNormvsMb -> SetYTitle("|Qb/Mb|");
+  fHistQbNormvsMb -> SetXTitle("|Mb|");
+  fHistList->Add(fHistQbNormvsMb);
+
+  fHistMavsMb = new TH2D("Flow_MavsMb_SP","Flow_MavsMb_SP",100,0.,100.,100,0.,100.);
+  fHistMavsMb -> SetYTitle("Ma");
+  fHistMavsMb -> SetXTitle("Mb");
+  fHistList->Add(fHistMavsMb);
+
+
   //weights
   if(fUsePhiWeights) {
     if(!fWeightsList) {
@@ -271,153 +334,173 @@ void AliFlowAnalysisWithScalarProduct::Make(AliFlowEventSimple* anEvent) {
     Double_t dMb = vQb.GetMult();
     if (dMa != 0 && dMb != 0) {
       
-      //fill control histograms     
-      fCommonHists->FillControlHistograms(anEvent);
-
-      //get total Q vector = the sum of subevent a and subevent b
-      AliFlowVector vQ = vQa + vQb;
-      //weight the Q vectors for the subevents by the multiplicity
-      //Note: Weight Q only in the particle loop when it is clear if it should be (m-1) or M
-      Double_t dQXa = vQa.X()/dMa; 
-      Double_t dQYa = vQa.Y()/dMa;
-      vQa.Set(dQXa,dQYa);
-      
-      Double_t dQXb = vQb.X()/dMb; 
-      Double_t dQYb = vQb.Y()/dMb;
-      vQb.Set(dQXb,dQYb);
+
+      //request that the subevent multiplicities are not too different
+      Double_t dRelDiff = TMath::Abs((dMa - dMb)/(dMa + dMb));
+      if (dRelDiff < fRelDiffMsub) {
+
+       //fill control histograms     
+       fCommonHists->FillControlHistograms(anEvent);
+
+       //fill some SP control histograms
+       fHistProQaQb -> Fill(1.,vQa*vQb,1.); //Fill with weight 1 -> Weight with MaMb????
+       fHistQaQbCos ->Fill(0.5*TMath::ACos((vQa/vQa.Mod())*(vQb/vQb.Mod())));  //Acos(Qa*Qb) = 2*phi
+       fHistResolution -> Fill(TMath::Cos( vQa.Phi()- vQb.Phi() ));  //vQa.Phi() returns 2*phi
+       fHistQaQb -> Fill(vQa*vQb);
+       fHistMavsMb -> Fill(dMb,dMa);
+
+       //get total Q vector = the sum of subevent a and subevent b
+       AliFlowVector vQ = vQa + vQb;
+       //weight the Q vectors for the subevents by the multiplicity
+       //Note: Weight Q only in the particle loop when it is clear if it should be (m-1) or M
+       Double_t dQXa = vQa.X()/dMa; 
+       Double_t dQYa = vQa.Y()/dMa;
+       vQa.Set(dQXa,dQYa);
+       
+       Double_t dQXb = vQb.X()/dMb; 
+       Double_t dQYb = vQb.Y()/dMb;
+       vQb.Set(dQXb,dQYb);
         
-      //scalar product of the two subevents
-      Double_t dQaQb = (vQa*vQb);
-      fHistProQaQb -> Fill(1.,dQaQb,dMa*dMb);  //Fill (QaQb/MaMb) with weight (MaMb). 
-      //needed for the error calculation:
-      fHistSumOfLinearWeights -> Fill(0.,dMa*dMb);
-      fHistSumOfQuadraticWeights -> Fill(0.,pow(dMa*dMb,2.));
-      
-      //loop over the tracks of the event
-      AliFlowTrackSimple*   pTrack = NULL; 
-      Int_t iNumberOfTracks = anEvent->NumberOfTracks(); 
-      Double_t dMq =  vQ.GetMult();
+       //scalar product of the two subevents
+       Double_t dQaQb = (vQa*vQb);
+       fHistProQaQbNorm -> Fill(1.,dQaQb,dMa*dMb);  //Fill (QaQb/MaMb) with weight (MaMb). 
+       //needed for the error calculation:
+       fHistSumOfLinearWeights -> Fill(0.,dMa*dMb);
+       fHistSumOfQuadraticWeights -> Fill(0.,pow(dMa*dMb,2.));
+
+       //fill some SP control histograms
+       fHistQaQbNorm ->Fill(vQa*vQb);
+       fHistQaNorm ->Fill(vQa.Mod());
+       fHistQaNormvsMa->Fill(dMa,vQa.Mod());
+       fHistQbNorm ->Fill(vQb.Mod());
+       fHistQbNormvsMb->Fill(dMb,vQb.Mod());
+       
+       //loop over the tracks of the event
+       AliFlowTrackSimple*   pTrack = NULL; 
+       Int_t iNumberOfTracks = anEvent->NumberOfTracks(); 
+       Double_t dMq =  vQ.GetMult();
       
-      for (Int_t i=0;i<iNumberOfTracks;i++) 
-       {
-         pTrack = anEvent->GetTrack(i) ; 
-         if (pTrack){
-           Double_t dPhi = pTrack->Phi();
-           //set default phi weight to 1
-           Double_t dW = 1.; 
-           //phi weight of pTrack
-           if(fUsePhiWeights && fPhiWeights) {
-             Int_t iNBinsPhi = fPhiWeights->GetNbinsX();
-             dW = fPhiWeights->GetBinContent(1+(Int_t)(TMath::Floor(dPhi*iNBinsPhi/TMath::TwoPi())));  
-             //bin = 1 + value*nbins/range
-             //TMath::Floor rounds to the lower integer
-           }     
+       for (Int_t i=0;i<iNumberOfTracks;i++) 
+         {
+           pTrack = anEvent->GetTrack(i) ; 
+           if (pTrack){
+             Double_t dPhi = pTrack->Phi();
+             //set default phi weight to 1
+             Double_t dW = 1.; 
+             //phi weight of pTrack
+             if(fUsePhiWeights && fPhiWeights) {
+               Int_t iNBinsPhi = fPhiWeights->GetNbinsX();
+               dW = fPhiWeights->GetBinContent(1+(Int_t)(TMath::Floor(dPhi*iNBinsPhi/TMath::TwoPi())));  
+               //bin = 1 + value*nbins/range
+               //TMath::Floor rounds to the lower integer
+             }     
            
-           //calculate vU
-           TVector2 vU;
-           Double_t dUX = TMath::Cos(2*dPhi);
-           Double_t dUY = TMath::Sin(2*dPhi);
-           vU.Set(dUX,dUY);
-           Double_t dModulus = vU.Mod();
-           if (dModulus!=0.) vU.Set(dUX/dModulus,dUY/dModulus);  // make length 1
-           else cerr<<"dModulus is zero!"<<endl;
+             //calculate vU
+             TVector2 vU;
+             Double_t dUX = TMath::Cos(2*dPhi);
+             Double_t dUY = TMath::Sin(2*dPhi);
+             vU.Set(dUX,dUY);
+             Double_t dModulus = vU.Mod();
+             if (dModulus!=0.) vU.Set(dUX/dModulus,dUY/dModulus);  // make length 1
+             else cerr<<"dModulus is zero!"<<endl;
            
-           //redefine the Q vector and devide by its multiplicity
-           TVector2 vQm;
-           Double_t dQmX = 0.;
-           Double_t dQmY = 0.;
-           //subtract particle from the flowvector if used to define it
-           if (pTrack->InSubevent(0) || pTrack->InSubevent(1)) { 
-             dQmX = (vQ.X() - dW*dUX)/(dMq-1);
-             dQmY = (vQ.Y() - dW*dUY)/(dMq-1);
-             
-             vQm.Set(dQmX,dQmY);
+             //redefine the Q vector and devide by its multiplicity
+             TVector2 vQm;
+             Double_t dQmX = 0.;
+             Double_t dQmY = 0.;
+             //subtract particle from the flowvector if used to define it
+             if (pTrack->InSubevent(0) || pTrack->InSubevent(1)) { 
+               dQmX = (vQ.X() - dW*dUX)/(dMq-1);
+               dQmY = (vQ.Y() - dW*dUY)/(dMq-1);
+               vQm.Set(dQmX,dQmY);
              
-             //dUQ = scalar product of vU and vQm
-             Double_t dUQ = (vU * vQm);
-             Double_t dPt = pTrack->Pt();
-             Double_t dEta = pTrack->Eta();
-             //fill the profile histograms
-             if (pTrack->InRPSelection()) {
-               fHistProUQetaRP -> Fill(dEta,dUQ,(dMq-1)); //Fill (Qu/(Mq-1)) with weight (Mq-1) 
-               //needed for the error calculation:
-               fHistProUQQaQbEtaRP -> Fill(dEta,dUQ*dQaQb,(dMq-1)*dMa*dMb); //Fill [Qu/(Mq-1)]*[QaQb/MaMb] with weight (Mq-1)MaMb          
-               fHistProUQPtRP -> Fill(dPt,dUQ,(dMq-1));                     //Fill (Qu/(Mq-1)) with weight (Mq-1)
-               fHistProUQQaQbPtRP -> Fill(dPt,dUQ*dQaQb,(dMq-1)*dMa*dMb);   //Fill [Qu/(Mq-1)]*[QaQb/MaMb] with weight (Mq-1)MaMb      
-               
-               fHistSumOfWeightsEtaRP[0]->Fill(dEta,(dMq-1));        // sum of Mq-1     
-               fHistSumOfWeightsEtaRP[1]->Fill(dEta,pow((dMq-1),2.));// sum of (Mq-1)^2     
-               fHistSumOfWeightsEtaRP[2]->Fill(dEta,(dMq-1)*dMa*dMb);// sum of (Mq-1)*MaMb     
-               fHistSumOfWeightsPtRP[0]->Fill(dPt,(dMq-1));          // sum of Mq-1     
-               fHistSumOfWeightsPtRP[1]->Fill(dPt,pow((dMq-1),2.));  // sum of (Mq-1)^2     
-               fHistSumOfWeightsPtRP[2]->Fill(dPt,(dMq-1)*dMa*dMb);  // sum of (Mq-1)*MaMb     
-             }
-             if (pTrack->InPOISelection()) {
-               fHistProUQetaPOI -> Fill(dEta,dUQ,(dMq-1));//Fill (Qu/(Mq-1)) with weight (Mq-1)
-               //needed for the error calculation:
-               fHistProUQQaQbEtaPOI -> Fill(dEta,dUQ*dQaQb,(dMq-1)*dMa*dMb); //Fill [Qu/(Mq-1)]*[QaQb/MaMb] with weight (Mq-1)MaMb         
-               fHistProUQPtPOI -> Fill(dPt,dUQ,(dMq-1));                     //Fill (Qu/(Mq-1)) with weight (Mq-1)
-               fHistProUQQaQbPtPOI -> Fill(dPt,dUQ*dQaQb,(dMq-1)*dMa*dMb);   //Fill [Qu/(Mq-1)]*[QaQb/MaMb] with weight (Mq-1)MaMb         
+               //dUQ = scalar product of vU and vQm
+               Double_t dUQ = (vU * vQm);
+               Double_t dPt = pTrack->Pt();
+               Double_t dEta = pTrack->Eta();
+               //fill the profile histograms
+               if (pTrack->InRPSelection()) {
+                 fHistProUQetaRP -> Fill(dEta,dUQ,(dMq-1)); //Fill (Qu/(Mq-1)) with weight (Mq-1) 
+                 //needed for the error calculation:
+                 fHistProUQQaQbEtaRP -> Fill(dEta,dUQ*dQaQb,(dMq-1)*dMa*dMb); //Fill [Qu/(Mq-1)]*[QaQb/MaMb] with weight (Mq-1)MaMb        
+                 fHistProUQPtRP -> Fill(dPt,dUQ,(dMq-1));                     //Fill (Qu/(Mq-1)) with weight (Mq-1)
+                 fHistProUQQaQbPtRP -> Fill(dPt,dUQ*dQaQb,(dMq-1)*dMa*dMb);   //Fill [Qu/(Mq-1)]*[QaQb/MaMb] with weight (Mq-1)MaMb    
+                 
+                 fHistSumOfWeightsEtaRP[0]->Fill(dEta,(dMq-1));        // sum of Mq-1     
+                 fHistSumOfWeightsEtaRP[1]->Fill(dEta,pow((dMq-1),2.));// sum of (Mq-1)^2     
+                 fHistSumOfWeightsEtaRP[2]->Fill(dEta,(dMq-1)*dMa*dMb);// sum of (Mq-1)*MaMb     
+                 fHistSumOfWeightsPtRP[0]->Fill(dPt,(dMq-1));          // sum of Mq-1     
+                 fHistSumOfWeightsPtRP[1]->Fill(dPt,pow((dMq-1),2.));  // sum of (Mq-1)^2     
+                 fHistSumOfWeightsPtRP[2]->Fill(dPt,(dMq-1)*dMa*dMb);  // sum of (Mq-1)*MaMb     
+               }
+               if (pTrack->InPOISelection()) {
+                 fHistProUQetaPOI -> Fill(dEta,dUQ,(dMq-1));//Fill (Qu/(Mq-1)) with weight (Mq-1)
+                 //needed for the error calculation:
+                 fHistProUQQaQbEtaPOI -> Fill(dEta,dUQ*dQaQb,(dMq-1)*dMa*dMb); //Fill [Qu/(Mq-1)]*[QaQb/MaMb] with weight (Mq-1)MaMb       
+                 fHistProUQPtPOI -> Fill(dPt,dUQ,(dMq-1));                     //Fill (Qu/(Mq-1)) with weight (Mq-1)
+                 fHistProUQQaQbPtPOI -> Fill(dPt,dUQ*dQaQb,(dMq-1)*dMa*dMb);   //Fill [Qu/(Mq-1)]*[QaQb/MaMb] with weight (Mq-1)MaMb       
+                 
+                 fHistSumOfWeightsEtaPOI[0]->Fill(dEta,(dMq-1));        // sum of Mq-1     
+                 fHistSumOfWeightsEtaPOI[1]->Fill(dEta,pow((dMq-1),2.));// sum of (Mq-1)^2     
+                 fHistSumOfWeightsEtaPOI[2]->Fill(dEta,(dMq-1)*dMa*dMb);// sum of (Mq-1)*MaMb     
+                 fHistSumOfWeightsPtPOI[0]->Fill(dPt,(dMq-1));          // sum of Mq-1     
+                 fHistSumOfWeightsPtPOI[1]->Fill(dPt,pow((dMq-1),2.)); // sum of (Mq-1)^2     
+                 fHistSumOfWeightsPtPOI[2]->Fill(dPt,(dMq-1)*dMa*dMb); // sum of (Mq-1)*MaMb     
+               }  
                
-               fHistSumOfWeightsEtaPOI[0]->Fill(dEta,(dMq-1));        // sum of Mq-1     
-               fHistSumOfWeightsEtaPOI[1]->Fill(dEta,pow((dMq-1),2.));// sum of (Mq-1)^2     
-               fHistSumOfWeightsEtaPOI[2]->Fill(dEta,(dMq-1)*dMa*dMb);// sum of (Mq-1)*MaMb     
-               fHistSumOfWeightsPtPOI[0]->Fill(dPt,(dMq-1));          // sum of Mq-1     
-               fHistSumOfWeightsPtPOI[1]->Fill(dPt,pow((dMq-1),2.)); // sum of (Mq-1)^2     
-               fHistSumOfWeightsPtPOI[2]->Fill(dPt,(dMq-1)*dMa*dMb); // sum of (Mq-1)*MaMb     
-             }  
+             } else { //do not subtract the particle from the flowvector
+               dQmX = vQ.X()/dMq;
+               dQmY = vQ.Y()/dMq;
+               vQm.Set(dQmX,dQmY);
              
-           } else { //do not subtract the particle from the flowvector
-             dQmX = vQ.X()/dMq;
-             dQmY = vQ.Y()/dMq;
+               //dUQ = scalar product of vU and vQm
+               Double_t dUQ = (vU * vQm);
+               Double_t dPt = pTrack->Pt();
+               Double_t dEta = pTrack->Eta();
+               //fill the profile histograms
+               if (pTrack->InRPSelection()) {
+                 fHistProUQetaRP -> Fill(dEta,dUQ,dMq);                   //Fill (Qu/Mq) with weight Mq 
+                 //needed for the error calculation:
+                 fHistProUQQaQbEtaRP -> Fill(dEta,dUQ*dQaQb,dMq*dMa*dMb); //Fill [Qu/Mq]*[QaQb/MaMb] with weight Mq*MaMb           
+                 fHistProUQPtRP -> Fill(dPt,dUQ,dMq);                     //Fill (Qu/Mq) with weight Mq 
+                 fHistProUQQaQbPtRP -> Fill(dPt,dUQ*dQaQb,dMq*dMa*dMb);   //Fill [Qu/Mq]*[QaQb/MaMb] with weight Mq*MaMb           
+                 
+                 fHistSumOfWeightsEtaRP[0]->Fill(dEta,dMq);        // sum of Mq     
+                 fHistSumOfWeightsEtaRP[1]->Fill(dEta,pow(dMq,2.));// sum of Mq^2     
+                 fHistSumOfWeightsEtaRP[2]->Fill(dEta,dMq*dMa*dMb);// sum of Mq*MaMb     
+                 fHistSumOfWeightsPtRP[0]->Fill(dPt,dMq);          // sum of Mq     
+                 fHistSumOfWeightsPtRP[1]->Fill(dPt,pow(dMq,2.));  // sum of Mq^2     
+                 fHistSumOfWeightsPtRP[2]->Fill(dPt,dMq*dMa*dMb);  // sum of Mq*MaMb     
+               }
+               if (pTrack->InPOISelection()) {
+                 fHistProUQetaPOI -> Fill(dEta,dUQ,dMq); //Fill (Qu/Mq) with weight Mq 
+                 //needed for the error calculation:
+                 fHistProUQQaQbEtaPOI -> Fill(dEta,dUQ*dQaQb,dMq*dMa*dMb); //Fill [Qu/Mq]*[QaQb/MaMb] with weight Mq*MaMb          
+                 fHistProUQPtPOI -> Fill(dPt,dUQ,dMq);                     //Fill (Qu/Mq) with weight Mq 
+                 fHistProUQQaQbPtPOI -> Fill(dPt,dUQ*dQaQb,dMq*dMa*dMb);   //Fill [Qu/Mq]*[QaQb/MaMb] with weight Mq*MaMb          
+                 
+                 fHistSumOfWeightsEtaPOI[0]->Fill(dEta,dMq);        // sum of Mq     
+                 fHistSumOfWeightsEtaPOI[1]->Fill(dEta,pow(dMq,2.));// sum of Mq^2     
+                 fHistSumOfWeightsEtaPOI[2]->Fill(dEta,dMq*dMa*dMb);// sum of Mq*MaMb     
+                 fHistSumOfWeightsPtPOI[0]->Fill(dPt,dMq);          // sum of Mq     
+                 fHistSumOfWeightsPtPOI[1]->Fill(dPt,pow(dMq,2.));  // sum of Mq^2     
+                 fHistSumOfWeightsPtPOI[2]->Fill(dPt,dMq*dMa*dMb);  // sum of Mq*MaMb     
+               }  
+             }//track not in subevents
              
-             vQm.Set(dQmX,dQmY);
-             
-             //dUQ = scalar product of vU and vQm
-             Double_t dUQ = (vU * vQm);
-             Double_t dPt = pTrack->Pt();
-             Double_t dEta = pTrack->Eta();
-             //fill the profile histograms
-             if (pTrack->InRPSelection()) {
-               fHistProUQetaRP -> Fill(dEta,dUQ,dMq);                   //Fill (Qu/Mq) with weight Mq 
-               //needed for the error calculation:
-               fHistProUQQaQbEtaRP -> Fill(dEta,dUQ*dQaQb,dMq*dMa*dMb); //Fill [Qu/Mq]*[QaQb/MaMb] with weight Mq*MaMb     
-               fHistProUQPtRP -> Fill(dPt,dUQ,dMq);                     //Fill (Qu/Mq) with weight Mq 
-               fHistProUQQaQbPtRP -> Fill(dPt,dUQ*dQaQb,dMq*dMa*dMb);   //Fill [Qu/Mq]*[QaQb/MaMb] with weight Mq*MaMb     
-               
-               fHistSumOfWeightsEtaRP[0]->Fill(dEta,dMq);        // sum of Mq     
-               fHistSumOfWeightsEtaRP[1]->Fill(dEta,pow(dMq,2.));// sum of Mq^2     
-               fHistSumOfWeightsEtaRP[2]->Fill(dEta,dMq*dMa*dMb);// sum of Mq*MaMb     
-               fHistSumOfWeightsPtRP[0]->Fill(dPt,dMq);          // sum of Mq     
-               fHistSumOfWeightsPtRP[1]->Fill(dPt,pow(dMq,2.));  // sum of Mq^2     
-               fHistSumOfWeightsPtRP[2]->Fill(dPt,dMq*dMa*dMb);  // sum of Mq*MaMb     
-             }
-             if (pTrack->InPOISelection()) {
-               fHistProUQetaPOI -> Fill(dEta,dUQ,dMq); //Fill (Qu/Mq) with weight Mq 
-               //needed for the error calculation:
-               fHistProUQQaQbEtaPOI -> Fill(dEta,dUQ*dQaQb,dMq*dMa*dMb); //Fill [Qu/Mq]*[QaQb/MaMb] with weight Mq*MaMb            
-               fHistProUQPtPOI -> Fill(dPt,dUQ,dMq);                     //Fill (Qu/Mq) with weight Mq 
-               fHistProUQQaQbPtPOI -> Fill(dPt,dUQ*dQaQb,dMq*dMa*dMb);   //Fill [Qu/Mq]*[QaQb/MaMb] with weight Mq*MaMb            
-               
-               fHistSumOfWeightsEtaPOI[0]->Fill(dEta,dMq);        // sum of Mq     
-               fHistSumOfWeightsEtaPOI[1]->Fill(dEta,pow(dMq,2.));// sum of Mq^2     
-               fHistSumOfWeightsEtaPOI[2]->Fill(dEta,dMq*dMa*dMb);// sum of Mq*MaMb     
-               fHistSumOfWeightsPtPOI[0]->Fill(dPt,dMq);          // sum of Mq     
-               fHistSumOfWeightsPtPOI[1]->Fill(dPt,pow(dMq,2.));  // sum of Mq^2     
-               fHistSumOfWeightsPtPOI[2]->Fill(dPt,dMq*dMa*dMb);  // sum of Mq*MaMb     
-             }  
-           }//track not in subevents
+           }//track
            
-         }//track
-         
-       }//loop over tracks
-      
-      fEventNumber++;
-      //    cout<<"@@@@@ "<<fEventNumber<<" events processed"<<endl;
+         }//loop over tracks
+       
+       fEventNumber++;
+
+      } //difference Ma and Mb
 
     }// subevents not empty 
     delete [] vQarray;
+
   } //event
+
 }//end of Make()
 
 //--------------------------------------------------------------------  
@@ -432,17 +515,18 @@ void AliFlowAnalysisWithScalarProduct::GetOutputHistograms(TList *outputListHist
     AliFlowCommonHistResults *pCommonHistResults = dynamic_cast<AliFlowCommonHistResults*> 
       (outputListHistos->FindObject("AliFlowCommonHistResultsSP"));
     TProfile* pHistProQaQb     = dynamic_cast<TProfile*>(outputListHistos->FindObject("Flow_QaQb_SP"));
+    TProfile* pHistProQaQbNorm = dynamic_cast<TProfile*>(outputListHistos->FindObject("FlowPro_QaQbNorm_SP"));
+    TH1D*     pHistSumOfLinearWeights    = dynamic_cast<TH1D*>(outputListHistos->FindObject("Flow_SumOfLinearWeights_SP"));
+    TH1D*     pHistSumOfQuadraticWeights = dynamic_cast<TH1D*>(outputListHistos->FindObject("Flow_SumOfQuadraticWeights_SP"));
+
     TProfile* pHistProUQetaRP  = dynamic_cast<TProfile*>(outputListHistos->FindObject("Flow_UQetaRP_SP"));
     TProfile* pHistProUQetaPOI = dynamic_cast<TProfile*>(outputListHistos->FindObject("Flow_UQetaPOI_SP"));
     TProfile* pHistProUQPtRP   = dynamic_cast<TProfile*>(outputListHistos->FindObject("Flow_UQPtRP_SP"));
     TProfile* pHistProUQPtPOI  = dynamic_cast<TProfile*>(outputListHistos->FindObject("Flow_UQPtPOI_SP"));
-    TH1D*     pHistSumOfLinearWeights    = dynamic_cast<TH1D*>(outputListHistos->FindObject("Flow_SumOfLinearWeights_SP"));
-    TH1D*     pHistSumOfQuadraticWeights = dynamic_cast<TH1D*>(outputListHistos->FindObject("Flow_SumOfQuadraticWeights_SP"));
     TProfile* pHistProUQQaQbPtRP    = dynamic_cast<TProfile*>(outputListHistos->FindObject("Flow_UQQaQbPtRP_SP"));
     TProfile* pHistProUQQaQbEtaRP   = dynamic_cast<TProfile*>(outputListHistos->FindObject("Flow_UQQaQbEtaRP_SP"));
     TProfile* pHistProUQQaQbPtPOI   = dynamic_cast<TProfile*>(outputListHistos->FindObject("Flow_UQQaQbPtPOI_SP"));
     TProfile* pHistProUQQaQbEtaPOI  = dynamic_cast<TProfile*>(outputListHistos->FindObject("Flow_UQQaQbEtaPOI_SP"));
-        
     TString weightFlag[3] = {"w_Qu_","w_Qu^2_","w_QuQaQb_"}; 
     TH1D* pHistSumOfWeightsPtRP[3] = {NULL};                    
     TH1D* pHistSumOfWeightsEtaRP[3] = {NULL};                    
@@ -456,14 +540,41 @@ void AliFlowAnalysisWithScalarProduct::GetOutputHistograms(TList *outputListHist
      pHistSumOfWeightsEtaPOI[i] = dynamic_cast<TH1D*>(outputListHistos->FindObject(Form("Flow_SumOfWeights%sEtaPOI_SP",weightFlag[i].Data())));
     }
 
+    TH1D*     pHistQaQb     = dynamic_cast<TH1D*>(outputListHistos->FindObject("Flow_QaQb_SP"));
+    TH1D*     pHistQaQbNorm = dynamic_cast<TH1D*>(outputListHistos->FindObject("Flow_QaQbNorm_SP"));
+    TH1D*     pHistQaQbCos  = dynamic_cast<TH1D*>(outputListHistos->FindObject("Flow_QaQbCos_SP"));
+    TH1D*     pHistResolution = dynamic_cast<TH1D*>(outputListHistos->FindObject("Flow_resolution_SP"));
+    TH1D*     pHistQaNorm   = dynamic_cast<TH1D*>(outputListHistos->FindObject("Flow_QaNorm_SP"));
+    TH2D*     pHistQaNormvsMa   = dynamic_cast<TH2D*>(outputListHistos->FindObject("Flow_QaNormvsMa_SP"));
+    TH1D*     pHistQbNorm   = dynamic_cast<TH1D*>(outputListHistos->FindObject("Flow_QbNorm_SP"));
+    TH2D*     pHistQbNormvsMb   = dynamic_cast<TH2D*>(outputListHistos->FindObject("Flow_QbNormvsMb_SP"));
+    TH2D*     pHistMavsMb = dynamic_cast<TH2D*>(outputListHistos->FindObject("Flow_MavsMb_SP"));
+
     //pass the pointers to the task
-    if (pCommonHist && pCommonHistResults && pHistProQaQb &&
-       pHistProUQetaRP && pHistProUQetaPOI && pHistProUQPtRP && pHistProUQPtPOI 
-       && pHistSumOfLinearWeights && pHistSumOfQuadraticWeights && pHistProUQQaQbPtRP
-       && pHistProUQQaQbEtaRP && pHistProUQQaQbPtPOI && pHistProUQQaQbEtaPOI) {
+    if (pCommonHist && pCommonHistResults && pHistProQaQb && pHistProQaQbNorm &&
+       pHistSumOfLinearWeights && pHistSumOfQuadraticWeights && 
+       pHistQaQb && pHistQaQbNorm && pHistQaQbCos && pHistResolution &&
+       pHistQaNorm && pHistQaNormvsMa && pHistQbNorm && pHistQbNormvsMb && 
+       pHistMavsMb &&
+       pHistProUQetaRP && pHistProUQetaPOI && 
+       pHistProUQPtRP && pHistProUQPtPOI &&  
+       pHistProUQQaQbPtRP && pHistProUQQaQbEtaRP && 
+       pHistProUQQaQbPtPOI && pHistProUQQaQbEtaPOI) {
       this -> SetCommonHists(pCommonHist);
       this -> SetCommonHistsRes(pCommonHistResults);
       this -> SetHistProQaQb(pHistProQaQb);
+      this -> SetHistProQaQbNorm(pHistProQaQbNorm);
+      this -> SetHistSumOfLinearWeights(pHistSumOfLinearWeights);
+      this -> SetHistSumOfQuadraticWeights(pHistSumOfQuadraticWeights); 
+      this -> SetHistQaQb(pHistQaQb);
+      this -> SetHistQaQbNorm(pHistQaQbNorm);
+      this -> SetHistQaQbCos(pHistQaQbCos);
+      this -> SetHistResolution(pHistResolution);
+      this -> SetHistQaNorm(pHistQaNorm);
+      this -> SetHistQaNormvsMa(pHistQaNormvsMa);
+      this -> SetHistQbNorm(pHistQbNorm);
+      this -> SetHistQbNormvsMb(pHistQbNormvsMb);
+      this -> SetHistMavsMb(pHistMavsMb);
       this -> SetHistProUQetaRP(pHistProUQetaRP);
       this -> SetHistProUQetaPOI(pHistProUQetaPOI);
       this -> SetHistProUQPtRP(pHistProUQPtRP);
@@ -472,17 +583,17 @@ void AliFlowAnalysisWithScalarProduct::GetOutputHistograms(TList *outputListHist
       this -> SetHistProUQQaQbEtaRP(pHistProUQQaQbEtaRP);
       this -> SetHistProUQQaQbPtPOI(pHistProUQQaQbPtPOI);
       this -> SetHistProUQQaQbEtaPOI(pHistProUQQaQbEtaPOI);      
-      this -> SetHistSumOfLinearWeights(pHistSumOfLinearWeights);
-      this -> SetHistSumOfQuadraticWeights(pHistSumOfQuadraticWeights);      
-     }  
-   
-   for(Int_t i=0;i<3;i++)
-   {
-    if(pHistSumOfWeightsPtRP[i]) this -> SetHistSumOfWeightsPtRP(pHistSumOfWeightsPtRP[i],i);      
-    if(pHistSumOfWeightsEtaRP[i]) this -> SetHistSumOfWeightsEtaRP(pHistSumOfWeightsEtaRP[i],i);      
-    if(pHistSumOfWeightsPtPOI[i]) this -> SetHistSumOfWeightsPtPOI(pHistSumOfWeightsPtPOI[i],i);      
-    if(pHistSumOfWeightsEtaPOI[i]) this -> SetHistSumOfWeightsEtaPOI(pHistSumOfWeightsEtaPOI[i],i);      
-   }   
+         
+      for(Int_t i=0;i<3;i++) {
+       if(pHistSumOfWeightsPtRP[i]) this -> SetHistSumOfWeightsPtRP(pHistSumOfWeightsPtRP[i],i);      
+       if(pHistSumOfWeightsEtaRP[i]) this -> SetHistSumOfWeightsEtaRP(pHistSumOfWeightsEtaRP[i],i);      
+       if(pHistSumOfWeightsPtPOI[i]) this -> SetHistSumOfWeightsPtPOI(pHistSumOfWeightsPtPOI[i],i);      
+       if(pHistSumOfWeightsEtaPOI[i]) this -> SetHistSumOfWeightsEtaPOI(pHistSumOfWeightsEtaPOI[i],i);      
+      }   
+
+    } else {
+      cout<<"WARNING: Histograms needed to run Finish() in SP are not accessable!"<<endl; }
+         
   } // end of if(outputListHistos)
 }            
 
@@ -501,11 +612,20 @@ void AliFlowAnalysisWithScalarProduct::Finish() {
   Int_t iNbinsPt  = AliFlowCommonConstants::GetMaster()->GetNbinsPt();
   Int_t iNbinsEta = AliFlowCommonConstants::GetMaster()->GetNbinsEta();
    
+  //Calculate the event plane resolution
+  //----------------------------------
+  Double_t dCos2phi = fHistResolution->GetMean();
+  if (dCos2phi !=0.0){
+    Double_t dResolution = TMath::Sqrt(2*dCos2phi); 
+    cout<<"An estimate of the event plane resolution is: "<<dResolution<<endl;
+    cout<<endl;
+  }
+
   //Calculate reference flow (noname)
   //----------------------------------
   //weighted average over (QaQb/MaMb) with weight (MaMb)
-  Double_t dQaQb  = fHistProQaQb->GetBinContent(1);
-  Double_t dV = 0.;
+  Double_t dQaQb  = fHistProQaQbNorm->GetBinContent(1);
+  Double_t dV = -999.; 
   if(dQaQb>=0.)
   {
    dV = TMath::Sqrt(dQaQb); 
@@ -514,7 +634,7 @@ void AliFlowAnalysisWithScalarProduct::Finish() {
   //  statistical error = term1 * spread * term2:
   //  term1 = sqrt{sum_{i=1}^{N} w^2}/(sum_{i=1}^{N} w)
   //  term2 = 1/sqrt(1-term1^2) 
-  Double_t dSpreadQaQb = fHistProQaQb->GetBinError(1);
+  Double_t dSpreadQaQb = fHistProQaQbNorm->GetBinError(1);
   Double_t dSumOfLinearWeights = fHistSumOfLinearWeights->GetBinContent(1);
   Double_t dSumOfQuadraticWeights = fHistSumOfQuadraticWeights->GetBinContent(1);
   Double_t dTerm1 = 0.;
@@ -539,7 +659,7 @@ void AliFlowAnalysisWithScalarProduct::Finish() {
   //v as a function of eta for RP selection
   for(Int_t b=1;b<iNbinsEta+1;b++) {
     Double_t duQpro = fHistProUQetaRP->GetBinContent(b);
-    Double_t dv2pro = 0.;
+    Double_t dv2pro = -999.;
     if (dV!=0.) { dv2pro = duQpro/dV; }
     //calculate the statistical error
     Double_t dv2ProErr = CalculateStatisticalError(b, dStatErrorQaQb, fHistProUQetaRP, fHistProUQQaQbEtaRP, fHistSumOfWeightsEtaRP);
@@ -551,7 +671,7 @@ void AliFlowAnalysisWithScalarProduct::Finish() {
   //v as a function of eta for POI selection
   for(Int_t b=1;b<iNbinsEta+1;b++) {
     Double_t duQpro = fHistProUQetaPOI->GetBinContent(b);
-    Double_t dv2pro = 0.;
+    Double_t dv2pro = -999.;
     if (dV!=0.) { dv2pro = duQpro/dV; }
     //calculate the statistical error
     Double_t dv2ProErr = CalculateStatisticalError(b, dStatErrorQaQb, fHistProUQetaPOI, fHistProUQQaQbEtaPOI, fHistSumOfWeightsEtaPOI);
@@ -569,7 +689,7 @@ void AliFlowAnalysisWithScalarProduct::Finish() {
   
   for(Int_t b=1;b<iNbinsPt+1;b++) {
     Double_t duQpro = fHistProUQPtRP->GetBinContent(b);
-    Double_t dv2pro = 0.;
+    Double_t dv2pro = -999.;
     if (dV!=0.) { dv2pro = duQpro/dV; }
     //calculate the statistical error
     Double_t dv2ProErr = CalculateStatisticalError(b, dStatErrorQaQb, fHistProUQPtRP, fHistProUQQaQbPtRP, fHistSumOfWeightsPtRP);
@@ -603,7 +723,7 @@ void AliFlowAnalysisWithScalarProduct::Finish() {
   
   for(Int_t b=1;b<iNbinsPt+1;b++) {
     Double_t duQpro = fHistProUQPtPOI->GetBinContent(b);
-    Double_t dv2pro = 0.;
+    Double_t dv2pro = -999.;
     if (dV!=0.) { dv2pro = duQpro/dV; }
     //calculate the statistical error
     Double_t dv2ProErr = CalculateStatisticalError(b, dStatErrorQaQb, fHistProUQPtPOI, fHistProUQQaQbPtPOI, fHistSumOfWeightsPtPOI);
@@ -642,6 +762,7 @@ Double_t AliFlowAnalysisWithScalarProduct::CalculateStatisticalError(Int_t b, Do
   Double_t duQproSpread = pHistProUQ->GetBinError(b);
   Double_t sumOfMq = pHistSumOfWeights[0]->GetBinContent(b);
   Double_t sumOfMqSquared = pHistSumOfWeights[1]->GetBinContent(b);
+  Double_t dQaQb = fHistProQaQbNorm->GetBinContent(1);
   Double_t dTerm1 = 0.;
   Double_t dTerm2 = 0.;
   if(sumOfMq) {
@@ -660,17 +781,17 @@ Double_t AliFlowAnalysisWithScalarProduct::CalculateStatisticalError(Int_t b, Do
     Double_t dDenominator = 1.-dTerm1Cov/(dTerm2Cov*dTerm3Cov);
     Double_t dPrefactor = dTerm1Cov/(dTerm2Cov*dTerm3Cov);
     if(dDenominator!=0) {
-      Double_t dCovariance = (pHistProUQQaQb->GetBinContent(b)-fHistProQaQb->GetBinContent(1)*pHistProUQ->GetBinContent(b))/dDenominator;            
+      Double_t dCovariance = (pHistProUQQaQb->GetBinContent(b)-dQaQb*pHistProUQ->GetBinContent(b))/dDenominator;            
       dWeightedCovariance = dCovariance*dPrefactor; 
     }
   }
   
   Double_t dv2ProErr = 0.; // final statitical error: 
-  if(fHistProQaQb->GetBinContent(1)>0.) {
-    Double_t dv2ProErrorSquared = (1./4.)*pow(fHistProQaQb->GetBinContent(1),-3.)*
+  if(dQaQb>0.) {
+    Double_t dv2ProErrorSquared = (1./4.)*pow(dQaQb,-3.)*
       (pow(pHistProUQ->GetBinContent(b),2.)*pow(aStatErrorQaQb,2.)
-       + 4.*pow(fHistProQaQb->GetBinContent(1),2.)*pow(duQproErr,2.)
-       - 4.*fHistProQaQb->GetBinContent(1)*pHistProUQ->GetBinContent(b)*dWeightedCovariance);
+       + 4.*pow(dQaQb,2.)*pow(duQproErr,2.)
+       - 4.*dQaQb*pHistProUQ->GetBinContent(b)*dWeightedCovariance);
     if(dv2ProErrorSquared>0.) dv2ProErr = pow(dv2ProErrorSquared,0.5);
   } 
    
index 810090e..f952073 100644 (file)
@@ -14,6 +14,8 @@ class AliFlowEventSimple;
 class AliFlowCommonHist;
 class AliFlowCommonHistResults;
 
+class TH1D;
+class TH2D;
 class TProfile;
 class TList;
 class TFile;
@@ -40,16 +42,19 @@ class AliFlowAnalysisWithScalarProduct {
    void    WriteHistograms(TString* outputFileName);        //writes histograms locally
    void    WriteHistograms(TString outputFileName);         //writes histograms locally
    void    WriteHistograms(TDirectoryFile *outputFileName); //writes histograms locally
-   
-   void    SetDebug(Bool_t kt)   { this->fDebug = kt ; }
-   Bool_t  GetDebug() const      { return this->fDebug ; }
-
+    
    Double_t CalculateStatisticalError(Int_t bin, 
                                      Double_t aStatErrorQaQb,
                                      TProfile* aHistProUQ, 
                                      TProfile* aHistProUQQaQb, 
                                      TH1D** aHistSumOfWeights);
 
+   void    SetDebug(Bool_t kt)   { this->fDebug = kt ; }
+   Bool_t  GetDebug() const      { return this->fDebug ; }
+
+
+   void     SetRelDiffMsub(Double_t diff) { this->fRelDiffMsub = diff; }
+   Double_t GetRelDiffMsub() const        { return this->fRelDiffMsub; }
 
    //phi weights
    void    SetWeightsList(TList* const aWeightsList)  {this->fWeightsList = (TList*)aWeightsList->Clone();}
@@ -65,19 +70,32 @@ class AliFlowAnalysisWithScalarProduct {
    TProfile* GetHistProUQPtRP() const   {return this->fHistProUQPtRP;} 
    TProfile* GetHistProUQPtPOI() const  {return this->fHistProUQPtPOI;}
    TProfile* GetHistProQaQb() const     {return this->fHistProQaQb;}
+   TProfile* GetHistProQaQbNorm() const {return this->fHistProQaQbNorm;}
    TH1D*     GetHistSumOfLinearWeights() const    {return this->fHistSumOfLinearWeights;}
    TH1D*     GetHistSumOfQuadraticWeights() const {return this->fHistSumOfQuadraticWeights;}
-   TProfile* GetHistProUQQaQbPtRP() const         {return this->fHistProUQQaQbPtRP;}   
-   TProfile* GetHistProUQQaQbEtaRP() const        {return this->fHistProUQQaQbEtaRP;}   
-   TProfile* GetHistProUQQaQbPtPOI() const        {return this->fHistProUQQaQbPtPOI;}   
-   TProfile* GetHistProUQQaQbEtaPOI() const            {return this->fHistProUQQaQbEtaPOI;}   
+
+   TProfile* GetHistProUQQaQbPtRP() const   {return this->fHistProUQQaQbPtRP;}   
+   TProfile* GetHistProUQQaQbEtaRP() const  {return this->fHistProUQQaQbEtaRP;}   
+   TProfile* GetHistProUQQaQbPtPOI() const  {return this->fHistProUQQaQbPtPOI;}   
+   TProfile* GetHistProUQQaQbEtaPOI() const {return this->fHistProUQQaQbEtaPOI;} 
    TH1D*     GetHistSumOfWeightsPtRP(Int_t i) const    {return this->fHistSumOfWeightsPtRP[i];}
    TH1D*     GetHistSumOfWeightsEtaRP(Int_t i) const   {return this->fHistSumOfWeightsEtaRP[i];}
    TH1D*     GetHistSumOfWeightsPtPOI(Int_t i) const   {return this->fHistSumOfWeightsPtPOI[i];}
    TH1D*     GetHistSumOfWeightsEtaPOI(Int_t i) const  {return this->fHistSumOfWeightsEtaPOI[i];}
+
    AliFlowCommonHist*        GetCommonHists() const    {return this->fCommonHists; }
    AliFlowCommonHistResults* GetCommonHistsRes() const {return this->fCommonHistsRes; }
 
+   TH1D* GetHistQaQb() const       {return this->fHistQaQb;}
+   TH1D* GetHistQaQbNorm() const   {return this->fHistQaQbNorm;}
+   TH1D* GetHistQaQbCos() const    {return this->fHistQaQbCos;}
+   TH1D* GetHistResolution() const {return this->fHistResolution;}
+   TH1D* GetHistQaNorm() const     {return this->fHistQaNorm;}
+   TH2D* GetHistQaNormvsMa() const {return this->fHistQaNormvsMa;}
+   TH1D* GetHistQbNorm() const     {return this->fHistQbNorm;}
+   TH2D* GetHistQbNormvsMb() const {return this->fHistQbNormvsMb;}
+   TH2D* GetMavsMb() const         {return this->fHistMavsMb;}
+
    //histogram setters
    void SetHistProUQetaRP(TProfile* const aHistProUQetaRP)   
      {this->fHistProUQetaRP = aHistProUQetaRP;}
@@ -89,10 +107,13 @@ class AliFlowAnalysisWithScalarProduct {
      {this->fHistProUQPtPOI = aHistProUQPtPOI;}
    void SetHistProQaQb(TProfile* const aHistProQaQb)         
      {this->fHistProQaQb = aHistProQaQb;}
+   void SetHistProQaQbNorm(TProfile* const aHistProQaQbNorm)         
+     {this->fHistProQaQbNorm = aHistProQaQbNorm;}
    void SetHistSumOfLinearWeights(TH1D* const aHistSumOfLinearWeights) 
      {this->fHistSumOfLinearWeights = aHistSumOfLinearWeights;}
    void SetHistSumOfQuadraticWeights(TH1D* const aHistSumOfQuadraticWeights) 
      {this->fHistSumOfQuadraticWeights = aHistSumOfQuadraticWeights;}
+
    void SetHistProUQQaQbPtRP(TProfile* const aHistProUQQaQbPtRP)     
      {this->fHistProUQQaQbPtRP = aHistProUQQaQbPtRP;}   
    void SetHistProUQQaQbEtaRP(TProfile* const aHistProUQQaQbEtaRP)   
@@ -109,12 +130,31 @@ class AliFlowAnalysisWithScalarProduct {
      {this->fHistSumOfWeightsPtPOI[i] = aHistSumOfWeightsPtPOI;}  
    void SetHistSumOfWeightsEtaPOI(TH1D* const aHistSumOfWeightsEtaPOI, Int_t const i) 
      {this->fHistSumOfWeightsEtaPOI[i] = aHistSumOfWeightsEtaPOI;}   
+
    void SetCommonHists(AliFlowCommonHist* const someCommonHists)              
      {this->fCommonHists = someCommonHists; }
    void SetCommonHistsRes(AliFlowCommonHistResults* const someCommonHistsRes) 
      {this->fCommonHistsRes = someCommonHistsRes; }
-   
 
+   void SetHistQaQb(TH1D* const aHistQaQb)
+     {this->fHistQaQb = aHistQaQb;}
+   void SetHistQaQbNorm(TH1D* const aHistQaQbNorm)
+     {this->fHistQaQbNorm = aHistQaQbNorm;}
+   void SetHistQaQbCos(TH1D* const aHistQaQbCos)
+     {this->fHistQaQbCos = aHistQaQbCos;}
+   void SetHistResolution(TH1D* const aHistResolution)
+     {this->fHistResolution = aHistResolution;}
+   void SetHistQaNorm(TH1D* const aHistQaNorm)
+     {this->fHistQaNorm = aHistQaNorm;}
+   void SetHistQaNormvsMa(TH2D* const aHistQaNormvsMa)
+     {this->fHistQaNormvsMa = aHistQaNormvsMa;}
+   void SetHistQbNorm(TH1D* const aHistQbNorm)
+     {this->fHistQbNorm = aHistQbNorm;}
+   void SetHistQbNormvsMb(TH2D* const aHistQbNormvsMb)
+     {this->fHistQbNormvsMb = aHistQbNormvsMb;}
+   void SetHistMavsMb(TH2D* const aHistMavsMb)
+     {this->fHistMavsMb = aHistMavsMb;}
+   
  private:
    AliFlowAnalysisWithScalarProduct(const AliFlowAnalysisWithScalarProduct& anAnalysis);            //copy constructor
    AliFlowAnalysisWithScalarProduct& operator=(const AliFlowAnalysisWithScalarProduct& anAnalysis); //assignment operator 
@@ -122,16 +162,20 @@ class AliFlowAnalysisWithScalarProduct {
    Int_t      fEventNumber;      // event counter
    Bool_t     fDebug ;           // flag for analysis: more print statements
 
+   Double_t   fRelDiffMsub;      // the relative difference the two subevent multiplicities can have
+
    TList*     fWeightsList;      // list holding input histograms with phi weights
    Bool_t     fUsePhiWeights;    // use phi weights
    TH1F*      fPhiWeights;       // histogram holding phi weights
 
    TList*     fHistList;         //list to hold all output histograms  
+
    TProfile*  fHistProUQetaRP;   //uQ(eta) for RP
    TProfile*  fHistProUQetaPOI;  //uQ(eta) for POI
    TProfile*  fHistProUQPtRP;    //uQ(pt) for RP
    TProfile*  fHistProUQPtPOI;   //uQ(pt) for POI
-   TProfile*  fHistProQaQb;      //average of QaQb (for event plane resolution)
+   TProfile*  fHistProQaQb;      //average of QaQb 
+   TProfile*  fHistProQaQbNorm;  //average of QaQb/MaMb
    TH1D*      fHistSumOfLinearWeights;     //holds sum of Ma*Mb
    TH1D*      fHistSumOfQuadraticWeights;  //holds sum of (Ma*Mb)^2
    
@@ -146,7 +190,17 @@ class AliFlowAnalysisWithScalarProduct {
    
    AliFlowCommonHist*        fCommonHists;    //control histograms
    AliFlowCommonHistResults* fCommonHistsRes; //results histograms
-
+   
+   TH1D*      fHistQaQb;         // distribution of QaQb
+   TH1D*      fHistQaQbNorm;     // distribution of QaQb/MaMb
+   TH1D*      fHistQaQbCos;      // distribution of the angle between Qa and Qb (from Acos (va*vb))
+   TH1D*      fHistResolution;   // distribution of cos(2(phi_a - phi_b))
+   TH1D*      fHistQaNorm;       // distribution of Qa/Ma
+   TH2D*      fHistQaNormvsMa;   // distribution of Qa/Ma vs Ma
+   TH1D*      fHistQbNorm;       // distribution of Qb/Mb
+   TH2D*      fHistQbNormvsMb;   // distribution of Qb/Mb vs Mb
+   TH2D*      fHistMavsMb;       // Ma vs Mb
+      
    ClassDef(AliFlowAnalysisWithScalarProduct,0)  // macro for rootcint
      };
  
index d19c1a2..50a5c8d 100644 (file)
 * provided "as is" without express or implied warranty.                  * 
 **************************************************************************/
 
+///////////////////////////////////////////////
+// AliAnalysisTaskScalarProduct:
+//
+// analysis task for Scalar Product Method
+//
+// Author: Naomi van der Kolk (kolk@nikhef.nl)
+///////////////////////////////////////////////
+
+
 #include "Riostream.h" //needed as include
 
 class TFile;
@@ -28,14 +37,6 @@ class AliAnalysisTaskSE;
 #include "AliFlowCommonHist.h"
 #include "AliFlowCommonHistResults.h"
 
-///////////////////////////////////////////////
-// AliAnalysisTaskScalarProduct:
-//
-// analysis task for Scalar Product Method
-//
-// Author: Naomi van der Kolk (kolk@nikhef.nl)
-///////////////////////////////////////////////
-
 ClassImp(AliAnalysisTaskScalarProduct)
 
 //________________________________________________________________________
@@ -45,7 +46,8 @@ AliAnalysisTaskScalarProduct::AliAnalysisTaskScalarProduct(const char *name, Boo
   fSP(NULL),
   fListHistos(NULL),
   fUsePhiWeights(usePhiWeights),
-  fListWeights(NULL)
+  fListWeights(NULL),
+  fRelDiffMsub(1.0)
 {
   // Constructor
   cout<<"AliAnalysisTaskScalarProduct::AliAnalysisTaskScalarProduct(const char *name)"<<endl;
@@ -67,7 +69,8 @@ AliAnalysisTaskScalarProduct::AliAnalysisTaskScalarProduct() :
   fSP(NULL),
   fListHistos(NULL),
   fUsePhiWeights(kFALSE),
-  fListWeights(NULL)
+  fListWeights(NULL),
+  fRelDiffMsub(1.0)
   {
   // Constructor
   cout<<"AliAnalysisTaskScalarProduct::AliAnalysisTaskScalarProduct()"<<endl;
@@ -96,7 +99,10 @@ void AliAnalysisTaskScalarProduct::UserCreateOutputObjects()
   cout<<"AliAnalysisTaskScalarProduct::CreateOutputObjects()"<<endl;
   
   //Analyser
-  fSP  = new AliFlowAnalysisWithScalarProduct() ;
+  fSP = new AliFlowAnalysisWithScalarProduct();
+
+  //set the allowed relative difference in the subevent multiplicities
+  fSP->SetRelDiffMsub(fRelDiffMsub); 
     
   //for using phi weights:
   if(fUsePhiWeights) {
index ec336cd..ec859b0 100644 (file)
@@ -33,6 +33,9 @@ class AliAnalysisTaskScalarProduct : public AliAnalysisTaskSE {
   void   SetUsePhiWeights(Bool_t const aPhiW) {this->fUsePhiWeights = aPhiW;}
   Bool_t GetUsePhiWeights() const             {return this->fUsePhiWeights;}
 
+  void     SetRelDiffMsub(Double_t diff) { this->fRelDiffMsub = diff; }
+  Double_t GetRelDiffMsub() const        { return this->fRelDiffMsub; }
+
 
  private:
 
@@ -43,9 +46,10 @@ class AliAnalysisTaskScalarProduct : public AliAnalysisTaskSE {
   AliFlowAnalysisWithScalarProduct* fSP;            // analysis object
   TList*                            fListHistos;    // collection of output
 
-  Bool_t                            fUsePhiWeights; // use phi weights
-  TList*                            fListWeights;   // list with weights
+  Bool_t    fUsePhiWeights; // use phi weights
+  TList*    fListWeights;   // list with weights
 
+  Double_t  fRelDiffMsub;   // the relative difference the two subevent multiplicities can have
   
   ClassDef(AliAnalysisTaskScalarProduct, 0); // example of analysis
 };
index de176d8..9394188 100644 (file)
@@ -657,6 +657,7 @@ AliAnalysisTaskFlowEvent* AddTaskFlow(TString type, Bool_t* METHODS, Bool_t QA,
   //===========================================================================
   if (SP){
     AliAnalysisTaskScalarProduct *taskSP = new AliAnalysisTaskScalarProduct("TaskScalarProduct",WEIGHTS[0]);
+    taskSP->SetRelDiffMsub(0.1);
     mgr->AddTask(taskSP);
   }
   if (LYZ1SUM){