]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWG2/FLOW/AliFlowCommon/AliFlowEventSimple.cxx
modifications for subevents
[u/mrichter/AliRoot.git] / PWG2 / FLOW / AliFlowCommon / AliFlowEventSimple.cxx
index 2c03bdf9233cfff1d8b950965bb910a6b55a9c7f..1d712d1a63bfbaf344a7cdd68593f0a190533d97 100644 (file)
@@ -21,6 +21,7 @@
 #include "TH1F.h"
 #include "TH1D.h"
 #include "TProfile.h"
+#include "TBrowser.h"
 #include "AliFlowVector.h"
 #include "AliFlowTrackSimple.h"
 #include "AliFlowEventSimple.h"
@@ -39,10 +40,28 @@ ClassImp(AliFlowEventSimple)
 
 //-----------------------------------------------------------------------
 
-  AliFlowEventSimple::AliFlowEventSimple(Int_t aLenght):
-    fTrackCollection(NULL),
-    fNumberOfTracks(0),
-    fEventNSelTracksIntFlow(0)
+AliFlowEventSimple::AliFlowEventSimple():
+  fTrackCollection(NULL),
+  fNumberOfTracks(0),
+  fEventNSelTracksRP(0),
+  fMCReactionPlaneAngle(0.),
+  fNumberOfTracksWrap(NULL),
+  fEventNSelTracksRPWrap(NULL),
+  fMCReactionPlaneAngleWrap(NULL)
+{
+  cout << "AliFlowEventSimple: Default constructor to be used only by root for io" << endl;
+}
+
+//-----------------------------------------------------------------------
+
+AliFlowEventSimple::AliFlowEventSimple(Int_t aLenght):
+  fTrackCollection(NULL),
+  fNumberOfTracks(0),
+  fEventNSelTracksRP(0),
+  fMCReactionPlaneAngle(0.),
+  fNumberOfTracksWrap(NULL),
+  fEventNSelTracksRPWrap(NULL),
+  fMCReactionPlaneAngleWrap(NULL)
 {
   //constructor 
   fTrackCollection =  new TObjArray(aLenght) ;
@@ -54,7 +73,11 @@ AliFlowEventSimple::AliFlowEventSimple(const AliFlowEventSimple& anEvent):
   TObject(),
   fTrackCollection(anEvent.fTrackCollection),
   fNumberOfTracks(anEvent.fNumberOfTracks),
-  fEventNSelTracksIntFlow(anEvent.fEventNSelTracksIntFlow)
+  fEventNSelTracksRP(anEvent.fEventNSelTracksRP),
+  fMCReactionPlaneAngle(anEvent.fMCReactionPlaneAngle),
+  fNumberOfTracksWrap(anEvent.fNumberOfTracksWrap),
+  fEventNSelTracksRPWrap(anEvent.fEventNSelTracksRPWrap),
+  fMCReactionPlaneAngleWrap(anEvent.fMCReactionPlaneAngleWrap)
 {
   //copy constructor 
 }
@@ -65,8 +88,12 @@ AliFlowEventSimple& AliFlowEventSimple::operator=(const AliFlowEventSimple& anEv
 {
   *fTrackCollection = *anEvent.fTrackCollection ;
   fNumberOfTracks = anEvent.fNumberOfTracks;
-  fEventNSelTracksIntFlow = anEvent.fEventNSelTracksIntFlow;
-  
+  fEventNSelTracksRP = anEvent.fEventNSelTracksRP;
+  fMCReactionPlaneAngle = anEvent.fMCReactionPlaneAngle;
+  fNumberOfTracksWrap = anEvent.fNumberOfTracksWrap; 
+  fEventNSelTracksRPWrap = anEvent.fEventNSelTracksRPWrap;
+  fMCReactionPlaneAngleWrap=anEvent.fMCReactionPlaneAngleWrap;
+
   return *this;
 }
 
@@ -75,7 +102,10 @@ AliFlowEventSimple& AliFlowEventSimple::operator=(const AliFlowEventSimple& anEv
 AliFlowEventSimple::~AliFlowEventSimple()
 {
   //destructor
-  fTrackCollection->Delete() ; delete fTrackCollection ;
+  if (fTrackCollection) fTrackCollection->Delete(); delete fTrackCollection;
+  if (fNumberOfTracksWrap) delete fNumberOfTracksWrap;
+  if (fEventNSelTracksRPWrap) delete fEventNSelTracksRPWrap;
+  if (fMCReactionPlaneAngleWrap) delete fMCReactionPlaneAngleWrap;
 }
 
 //----------------------------------------------------------------------- 
@@ -117,15 +147,7 @@ AliFlowVector AliFlowEventSimple::GetQ(Int_t n, TList *weightsList, Bool_t usePh
   TH1F *phiWeights = NULL;
   TH1D *ptWeights  = NULL;
   TH1D *etaWeights = NULL;
-
-  Double_t dSumOfWeightsToPower2 = 0.; // sum_{i=1}^{n} pow((wPhi*wPt*wEta)_i, 2)
-  Double_t dSumOfWeightsToPower3 = 0.; // sum_{i=1}^{n} pow((wPhi*wPt*wEta)_i, 3)
-  Double_t dSumOfWeightsToPower4 = 0.; // sum_{i=1}^{n} pow((wPhi*wPt*wEta)_i, 4)
-  Double_t dSumOfWeightsToPower5 = 0.; // sum_{i=1}^{n} pow((wPhi*wPt*wEta)_i, 5)
-  Double_t dSumOfWeightsToPower6 = 0.; // sum_{i=1}^{n} pow((wPhi*wPt*wEta)_i, 6)
-  Double_t dSumOfWeightsToPower7 = 0.; // sum_{i=1}^{n} pow((wPhi*wPt*wEta)_i, 7)
-  Double_t dSumOfWeightsToPower8 = 0.; // sum_{i=1}^{n} pow((wPhi*wPt*wEta)_i, 8) 
-
+  
   if(weightsList)
   {
    if(usePhiWeights)
@@ -159,7 +181,7 @@ AliFlowVector AliFlowEventSimple::GetQ(Int_t n, TList *weightsList, Bool_t usePh
    pTrack = (AliFlowTrackSimple*)TrackCollection()->At(i); 
    if(pTrack)
    {
-    if(pTrack->UseForIntegratedFlow()) 
+    if(pTrack->InRPSelection()) 
     {
      dPhi = pTrack->Phi();
      dPt  = pTrack->Pt();
@@ -187,41 +209,167 @@ AliFlowVector AliFlowEventSimple::GetQ(Int_t n, TList *weightsList, Bool_t usePh
     
      // weighted multiplicity:
      iUsedTracks+=wPhi*wPt*wEta;
-    
-     // weights raised to various powers are summed up:
-     dSumOfWeightsToPower2+=pow(wPhi*wPt*wEta, 2); 
-     dSumOfWeightsToPower3+=pow(wPhi*wPt*wEta, 3); 
-     dSumOfWeightsToPower4+=pow(wPhi*wPt*wEta, 4); 
-     dSumOfWeightsToPower5+=pow(wPhi*wPt*wEta, 5); 
-     dSumOfWeightsToPower6+=pow(wPhi*wPt*wEta, 6); 
-     dSumOfWeightsToPower7+=pow(wPhi*wPt*wEta, 7); 
-     dSumOfWeightsToPower8+=pow(wPhi*wPt*wEta, 8); 
-     
-    } // end of if (pTrack->UseForIntegratedFlow())
+         
+    } // end of if (pTrack->InRPSelection())
    } // end of if (pTrack)
    else {cerr << "no particle!!!"<<endl;}
   } // loop over particles
     
   vQ.Set(dQX,dQY);
   vQ.SetMult(iUsedTracks);
-  vQ.SetSumOfWeightsToPower2(dSumOfWeightsToPower2);
-  vQ.SetSumOfWeightsToPower3(dSumOfWeightsToPower3);
-  vQ.SetSumOfWeightsToPower4(dSumOfWeightsToPower4);
-  vQ.SetSumOfWeightsToPower5(dSumOfWeightsToPower5);
-  vQ.SetSumOfWeightsToPower6(dSumOfWeightsToPower6);
-  vQ.SetSumOfWeightsToPower7(dSumOfWeightsToPower7);
-  vQ.SetSumOfWeightsToPower8(dSumOfWeightsToPower8);
-
+  
   return vQ;
   
 }
 
+//-----------------------------------------------------------------------   
+void AliFlowEventSimple::GetQsub(AliFlowVector* Qarray, Int_t n, TList *weightsList, Bool_t usePhiWeights, Bool_t usePtWeights, Bool_t useEtaWeights) 
+{
+  
+  // calculate Q-vector in harmonic n without weights (default harmonic n=2)  
+  Double_t dQX = 0.;
+  Double_t dQY = 0.;
+    
+  Int_t iOrder = n;
+  Double_t iUsedTracks = 0;
+  Double_t dPhi = 0.;
+  Double_t dPt  = 0.;
+  Double_t dEta = 0.;
+  
+  AliFlowTrackSimple* pTrack = NULL;
+  Int_t    nBinsPhi    = 0; 
+  Double_t dBinWidthPt = 0.;
+  Double_t dPtMin      = 0.;
+  Double_t dBinWidthEta= 0.;
+  Double_t dEtaMin     = 0.;
+  Double_t wPhi = 1.;  // weight Phi  
+  Double_t wPt  = 1.;  // weight Pt 
+  Double_t wEta = 1.;  // weight Eta 
+  
+  TH1F *phiWeights = NULL;
+  TH1D *ptWeights  = NULL;
+  TH1D *etaWeights = NULL;
+  
+  if(weightsList)
+    {
+      if(usePhiWeights)
+       {
+    phiWeights = dynamic_cast<TH1F *>(weightsList->FindObject("phi_weights"));
+    if(phiWeights) nBinsPhi = phiWeights->GetNbinsX();
+   }          
+   if(usePtWeights)
+   {
+    ptWeights = dynamic_cast<TH1D *>(weightsList->FindObject("pt_weights"));
+    if(ptWeights)
+    {
+     dBinWidthPt = ptWeights->GetBinWidth(1); // assuming that all bins have the same width
+     dPtMin = (ptWeights->GetXaxis())->GetXmin();
+    } 
+   }       
+   if(useEtaWeights)
+   {
+    etaWeights = dynamic_cast<TH1D *>(weightsList->FindObject("eta_weights"));
+    if(etaWeights)
+    {
+     dBinWidthEta = etaWeights->GetBinWidth(1); // assuming that all bins have the same width
+     dEtaMin = (etaWeights->GetXaxis())->GetXmin();
+    } 
+   }          
+  } // end of if(weightsList)
+  
+  //loop over the two subevents
+  for (Int_t s=0;s<2;s++)  
+    {
+      // loop over tracks    
+      for(Int_t i=0;i<fNumberOfTracks;i++)                               
+       {
+         pTrack = (AliFlowTrackSimple*)TrackCollection()->At(i); 
+         if(pTrack)
+           {
+             if(pTrack->InRPSelection())
+               {
+                 if (pTrack->InSubevent(s)) {
+                   dPhi = pTrack->Phi();
+                   dPt  = pTrack->Pt();
+                   dEta = pTrack->Eta();
+                 
+                   // determine Phi weight: (to be improved, I should here only access it + the treatment of gaps in the if statement)
+                   if(phiWeights && nBinsPhi)
+                     {
+                       wPhi = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(dPhi*nBinsPhi/TMath::TwoPi())));
+                     }
+                   // determine v'(pt) weight:    
+                   if(ptWeights && dBinWidthPt)
+                     {
+                       wPt=ptWeights->GetBinContent(1+(Int_t)(TMath::Floor((dPt-dPtMin)/dBinWidthPt))); 
+                     }            
+                   // determine v'(eta) weight:    
+                   if(etaWeights && dBinWidthEta)
+                     {
+                       wEta=etaWeights->GetBinContent(1+(Int_t)(TMath::Floor((dEta-dEtaMin)/dBinWidthEta))); 
+                     } 
+                   
+                   // building up the weighted Q-vector:       
+                   dQX += wPhi*wPt*wEta*TMath::Cos(iOrder*dPhi);
+                   dQY += wPhi*wPt*wEta*TMath::Sin(iOrder*dPhi); 
+                   
+                   // weighted multiplicity:
+                   iUsedTracks+=wPhi*wPt*wEta;
+                                   
+                 } // end of subevent 
+               } // end of if (pTrack->InRPSelection())
+           } // end of if (pTrack)
+         else {cerr << "no particle!!!"<<endl;}
+       } // loop over particles
+      Qarray[s].Set(dQX,dQY);
+      Qarray[s].SetMult(iUsedTracks);
+      //reset
+      iUsedTracks = 0;
+      dQX = 0.;
+      dQY = 0.;
+    }
+  
+  //return vQ;
+  
+}
 
 
+//----------------------------------------------------------------------- 
+void AliFlowEventSimple::Print(Option_t *option) const
+{
+  //   -*-*-*-*-*Print some global quantities for this histogram collection class *-*-*-*-*-*-*-*
+  //             ===============================================
+  //   printf( "TH1.Print Name  = %s, Entries= %d, Total sum= %g\n",GetName(),Int_t(fEntries),GetSumOfWeights());
+  printf( "Class.Print Name = %s, Total number of tracks= %d, Number of selected tracks= %d, MC EventPlaneAngle= %f",
+         GetName(),fNumberOfTracks, fEventNSelTracksRP, fMCReactionPlaneAngle );
 
+  if (fTrackCollection) {  
+    fTrackCollection->Print(option);
+  }
+  else {
+    printf( "Empty track collection \n");
+  }
+}
 
+//----------------------------------------------------------------------- 
+ void AliFlowEventSimple::Browse(TBrowser *b)
+{
+  if (!b) return;
+  if (!fNumberOfTracksWrap) {
+    fNumberOfTracksWrap = new TParameter<int>("fNumberOfTracks", fNumberOfTracks);
+    b->Add(fNumberOfTracksWrap);
+  }
+  if (!fEventNSelTracksRPWrap) {
+    fEventNSelTracksRPWrap = new TParameter<int>("fEventNSelTracksRP", fEventNSelTracksRP);
+    b->Add(fEventNSelTracksRPWrap);
+  }
+  if (!fMCReactionPlaneAngleWrap) {
+    fMCReactionPlaneAngleWrap = new TParameter<double>(" fMCReactionPlaneAngle",  fMCReactionPlaneAngle);
+    b->Add( fMCReactionPlaneAngleWrap);
+  }
+  if (fTrackCollection) b->Add(fTrackCollection,"AliFlowTracksSimple");
+}
 
-
-
-
-
+