]> 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 3faadb8d3403a0538d2f496826dee78c4b4ae8bb..1d712d1a63bfbaf344a7cdd68593f0a190533d97 100644 (file)
@@ -43,9 +43,10 @@ ClassImp(AliFlowEventSimple)
 AliFlowEventSimple::AliFlowEventSimple():
   fTrackCollection(NULL),
   fNumberOfTracks(0),
-  fEventNSelTracksIntFlow(0),
+  fEventNSelTracksRP(0),
+  fMCReactionPlaneAngle(0.),
   fNumberOfTracksWrap(NULL),
-  fEventNSelTracksIntFlowWrap(NULL),
+  fEventNSelTracksRPWrap(NULL),
   fMCReactionPlaneAngleWrap(NULL)
 {
   cout << "AliFlowEventSimple: Default constructor to be used only by root for io" << endl;
@@ -56,9 +57,10 @@ AliFlowEventSimple::AliFlowEventSimple():
 AliFlowEventSimple::AliFlowEventSimple(Int_t aLenght):
   fTrackCollection(NULL),
   fNumberOfTracks(0),
-  fEventNSelTracksIntFlow(0),
+  fEventNSelTracksRP(0),
+  fMCReactionPlaneAngle(0.),
   fNumberOfTracksWrap(NULL),
-  fEventNSelTracksIntFlowWrap(NULL),
+  fEventNSelTracksRPWrap(NULL),
   fMCReactionPlaneAngleWrap(NULL)
 {
   //constructor 
@@ -71,10 +73,10 @@ AliFlowEventSimple::AliFlowEventSimple(const AliFlowEventSimple& anEvent):
   TObject(),
   fTrackCollection(anEvent.fTrackCollection),
   fNumberOfTracks(anEvent.fNumberOfTracks),
-  fEventNSelTracksIntFlow(anEvent.fEventNSelTracksIntFlow),
+  fEventNSelTracksRP(anEvent.fEventNSelTracksRP),
   fMCReactionPlaneAngle(anEvent.fMCReactionPlaneAngle),
   fNumberOfTracksWrap(anEvent.fNumberOfTracksWrap),
-  fEventNSelTracksIntFlowWrap(anEvent.fEventNSelTracksIntFlowWrap),
+  fEventNSelTracksRPWrap(anEvent.fEventNSelTracksRPWrap),
   fMCReactionPlaneAngleWrap(anEvent.fMCReactionPlaneAngleWrap)
 {
   //copy constructor 
@@ -86,10 +88,10 @@ AliFlowEventSimple& AliFlowEventSimple::operator=(const AliFlowEventSimple& anEv
 {
   *fTrackCollection = *anEvent.fTrackCollection ;
   fNumberOfTracks = anEvent.fNumberOfTracks;
-  fEventNSelTracksIntFlow = anEvent.fEventNSelTracksIntFlow;
+  fEventNSelTracksRP = anEvent.fEventNSelTracksRP;
   fMCReactionPlaneAngle = anEvent.fMCReactionPlaneAngle;
   fNumberOfTracksWrap = anEvent.fNumberOfTracksWrap; 
-  fEventNSelTracksIntFlowWrap = anEvent.fEventNSelTracksIntFlowWrap;
+  fEventNSelTracksRPWrap = anEvent.fEventNSelTracksRPWrap;
   fMCReactionPlaneAngleWrap=anEvent.fMCReactionPlaneAngleWrap;
 
   return *this;
@@ -102,7 +104,7 @@ AliFlowEventSimple::~AliFlowEventSimple()
   //destructor
   if (fTrackCollection) fTrackCollection->Delete(); delete fTrackCollection;
   if (fNumberOfTracksWrap) delete fNumberOfTracksWrap;
-  if (fEventNSelTracksIntFlowWrap) delete fEventNSelTracksIntFlowWrap;
+  if (fEventNSelTracksRPWrap) delete fEventNSelTracksRPWrap;
   if (fMCReactionPlaneAngleWrap) delete fMCReactionPlaneAngleWrap;
 }
 
@@ -145,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)
@@ -187,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();
@@ -215,43 +209,141 @@ 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= %d",
-         GetName(),fNumberOfTracks, fEventNSelTracksIntFlow, fMCReactionPlaneAngle );
+  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);
@@ -269,24 +361,15 @@ void AliFlowEventSimple::Print(Option_t *option) const
     fNumberOfTracksWrap = new TParameter<int>("fNumberOfTracks", fNumberOfTracks);
     b->Add(fNumberOfTracksWrap);
   }
-  if (!fEventNSelTracksIntFlowWrap) {
-    fEventNSelTracksIntFlowWrap = new TParameter<int>("fEventNSelTracksIntFlow", fEventNSelTracksIntFlow);
-    b->Add(fEventNSelTracksIntFlowWrap);
+  if (!fEventNSelTracksRPWrap) {
+    fEventNSelTracksRPWrap = new TParameter<int>("fEventNSelTracksRP", fEventNSelTracksRP);
+    b->Add(fEventNSelTracksRPWrap);
   }
   if (!fMCReactionPlaneAngleWrap) {
-     fMCReactionPlaneAngleWrap = new TParameter<double>(" fMCReactionPlaneAngle",  fMCReactionPlaneAngle);
+    fMCReactionPlaneAngleWrap = new TParameter<double>(" fMCReactionPlaneAngle",  fMCReactionPlaneAngle);
     b->Add( fMCReactionPlaneAngleWrap);
   }
   if (fTrackCollection) b->Add(fTrackCollection,"AliFlowTracksSimple");
 }
 
-
-
-
-
-
-
-
-
-
-
+