coverity + some patch for pp
[u/mrichter/AliRoot.git] / PWG2 / FLOW / AliFlowCommon / AliFlowEventSimple.cxx
index c54eb1a..285a583 100644 (file)
@@ -53,6 +53,7 @@ AliFlowEventSimple::AliFlowEventSimple():
   fMCReactionPlaneAngle(0.),
   fMCReactionPlaneAngleIsSet(kFALSE),
   fAfterBurnerPrecision(0.001),
+  fUserModified(kFALSE),
   fNumberOfTracksWrap(NULL),
   fNumberOfRPsWrap(NULL),
   fMCReactionPlaneAngleWrap(NULL)
@@ -61,42 +62,32 @@ AliFlowEventSimple::AliFlowEventSimple():
 }
 
 //-----------------------------------------------------------------------
-AliFlowEventSimple::AliFlowEventSimple(Int_t aLength):
-  fTrackCollection(new TObjArray(aLength)),
-  fReferenceMultiplicity(0),
-  fNumberOfTracks(0),
-  fNumberOfRPs(0),
-  fMCReactionPlaneAngle(0.),
-  fMCReactionPlaneAngleIsSet(kFALSE),
-  fAfterBurnerPrecision(0.001),
-  fNumberOfTracksWrap(NULL),
-  fNumberOfRPsWrap(NULL),
-  fMCReactionPlaneAngleWrap(NULL)
-{
-  //constructor
-}
-
-//-----------------------------------------------------------------------
-AliFlowEventSimple::AliFlowEventSimple( Int_t nParticles,
+AliFlowEventSimple::AliFlowEventSimple( Int_t n,
+                                        ConstructionMethod method,
                                         TF1* ptDist,
                                         Double_t phiMin,
                                         Double_t phiMax,
                                         Double_t etaMin,
                                         Double_t etaMax):
-  fTrackCollection(new TObjArray(nParticles)),
-  fReferenceMultiplicity(nParticles),
+  fTrackCollection(new TObjArray(n)),
+  fReferenceMultiplicity(0),
   fNumberOfTracks(0),
   fNumberOfRPs(0),
   fMCReactionPlaneAngle(0.),
   fMCReactionPlaneAngleIsSet(kFALSE),
   fAfterBurnerPrecision(0.001),
+  fUserModified(kFALSE),
   fNumberOfTracksWrap(NULL),
   fNumberOfRPsWrap(NULL),
   fMCReactionPlaneAngleWrap(NULL)
 {
-  //ctor. generates nParticles random tracks with given Pt distribution
-  //phi and eta are uniform
-  Generate(nParticles,ptDist,phiMin,phiMax,etaMin,etaMax);
+  //ctor
+  // if second argument is set to AliFlowEventSimple::kGenerate
+  // it generates n random tracks with given Pt distribution
+  // (a sane default is provided), phi and eta are uniform
+
+  if (method==kGenerate)
+    Generate(n,ptDist,phiMin,phiMax,etaMin,etaMax);
 }
 
 //-----------------------------------------------------------------------
@@ -109,6 +100,7 @@ AliFlowEventSimple::AliFlowEventSimple(const AliFlowEventSimple& anEvent):
   fMCReactionPlaneAngle(anEvent.fMCReactionPlaneAngle),
   fMCReactionPlaneAngleIsSet(anEvent.fMCReactionPlaneAngleIsSet),
   fAfterBurnerPrecision(anEvent.fAfterBurnerPrecision),
+  fUserModified(anEvent.fUserModified),
   fNumberOfTracksWrap(anEvent.fNumberOfTracksWrap),
   fNumberOfRPsWrap(anEvent.fNumberOfRPsWrap),
   fMCReactionPlaneAngleWrap(anEvent.fMCReactionPlaneAngleWrap)
@@ -120,6 +112,7 @@ AliFlowEventSimple::AliFlowEventSimple(const AliFlowEventSimple& anEvent):
 AliFlowEventSimple& AliFlowEventSimple::operator=(const AliFlowEventSimple& anEvent)
 {
   //assignment operator
+  if (&anEvent==this) return *this; //check self-assignment
   if (fTrackCollection) fTrackCollection->Delete();
   delete fTrackCollection;
   fTrackCollection = (TObjArray*)(anEvent.fTrackCollection)->Clone(); //deep copy
@@ -129,6 +122,7 @@ AliFlowEventSimple& AliFlowEventSimple::operator=(const AliFlowEventSimple& anEv
   fMCReactionPlaneAngle = anEvent.fMCReactionPlaneAngle;
   fMCReactionPlaneAngleIsSet = anEvent.fMCReactionPlaneAngleIsSet;
   fAfterBurnerPrecision = anEvent.fAfterBurnerPrecision;
+  fUserModified=anEvent.fUserModified;
   fNumberOfTracksWrap = anEvent.fNumberOfTracksWrap;
   fNumberOfRPsWrap = anEvent.fNumberOfRPsWrap;
   fMCReactionPlaneAngleWrap=anEvent.fMCReactionPlaneAngleWrap;
@@ -141,9 +135,9 @@ AliFlowEventSimple::~AliFlowEventSimple()
   //destructor
   if (fTrackCollection) fTrackCollection->Delete();
   delete fTrackCollection;
-  if (fNumberOfTracksWrap) delete fNumberOfTracksWrap;
-  if (fNumberOfRPsWrap) delete fNumberOfRPsWrap;
-  if (fMCReactionPlaneAngleWrap) delete fMCReactionPlaneAngleWrap;
+  delete fNumberOfTracksWrap;
+  delete fNumberOfRPsWrap;
+  delete fMCReactionPlaneAngleWrap;
 }
 
 //-----------------------------------------------------------------------
@@ -161,6 +155,7 @@ void AliFlowEventSimple::Generate(Int_t nParticles,
     static TF1 ptdistribution("ptSpectra","x*TMath::Exp(-pow(0.13957*0.13957+x*x,0.5)/0.4)",0.1,10.);
     ptDist=&ptdistribution;
   }
+
   for (Int_t i=0; i<nParticles; i++)
   {
     AliFlowTrackSimple* track = new AliFlowTrackSimple();
@@ -171,6 +166,8 @@ void AliFlowEventSimple::Generate(Int_t nParticles,
     AddTrack(track);
   }
   fMCReactionPlaneAngle=gRandom->Uniform(0.0,TMath::TwoPi());
+  fMCReactionPlaneAngleIsSet=kTRUE;
+  SetUserModified();
 }
 
 //-----------------------------------------------------------------------
@@ -185,8 +182,13 @@ AliFlowTrackSimple* AliFlowEventSimple::GetTrack(Int_t i)
 //-----------------------------------------------------------------------
 void AliFlowEventSimple::AddTrack( AliFlowTrackSimple* track )
 {
-  //add a track
-  fTrackCollection->AddLast(track);
+  //add a track, delete the old one if necessary
+  if (fNumberOfTracks < fTrackCollection->GetEntriesFast())
+  {
+    TObject* o = fTrackCollection->At(fNumberOfTracks);
+    if (o) delete o;
+  }
+  fTrackCollection->AddAtAndExpand(track,fNumberOfTracks);
   fNumberOfTracks++;
 }
 
@@ -200,7 +202,7 @@ AliFlowVector AliFlowEventSimple::GetQ(Int_t n, TList *weightsList, Bool_t usePh
   vQ.Set(0.,0.);
 
   Int_t iOrder = n;
-  Double_t iUsedTracks = 0;
+  Double_t sumOfWeights = 0.;
   Double_t dPhi = 0.;
   Double_t dPt = 0.;
   Double_t dEta = 0.;
@@ -283,7 +285,7 @@ AliFlowVector AliFlowEventSimple::GetQ(Int_t n, TList *weightsList, Bool_t usePh
         dQY += dWeight*wPhi*wPt*wEta*TMath::Sin(iOrder*dPhi);
 
         // weighted multiplicity:
-        iUsedTracks += dWeight*wPhi*wPt*wEta;
+        sumOfWeights += dWeight*wPhi*wPt*wEta;
 
       } // end of if (pTrack->InRPSelection())
     } // end of if (pTrack)
@@ -294,7 +296,7 @@ AliFlowVector AliFlowEventSimple::GetQ(Int_t n, TList *weightsList, Bool_t usePh
   } // loop over particles
 
   vQ.Set(dQX,dQY);
-  vQ.SetMult(iUsedTracks);
+  vQ.SetMult(sumOfWeights);
 
   return vQ;
 
@@ -309,7 +311,7 @@ void AliFlowEventSimple::Get2Qsub(AliFlowVector* Qarray, Int_t n, TList *weights
   Double_t dQY = 0.;
 
   Int_t iOrder = n;
-  Double_t iUsedTracks = 0;
+  Double_t sumOfWeights = 0.;
   Double_t dPhi = 0.;
   Double_t dPt  = 0.;
   Double_t dEta = 0.;
@@ -338,12 +340,12 @@ void AliFlowEventSimple::Get2Qsub(AliFlowVector* Qarray, Int_t n, TList *weights
     if(usePhiWeights)
     {
       phiWeightsSub0 = dynamic_cast<TH1F *>(weightsList->FindObject("phi_weights_sub0"));
-      phiWeightsSub1 = dynamic_cast<TH1F *>(weightsList->FindObject("phi_weights_sub1"));
       if(phiWeightsSub0) {
-        iNbinsPhiSub0 = phiWeightsSub0->GetNbinsX();
+       iNbinsPhiSub0 = phiWeightsSub0->GetNbinsX();
       }
+      phiWeightsSub1 = dynamic_cast<TH1F *>(weightsList->FindObject("phi_weights_sub1"));
       if(phiWeightsSub1) {
-        iNbinsPhiSub1 = phiWeightsSub1->GetNbinsX();
+       iNbinsPhiSub1 = phiWeightsSub1->GetNbinsX();
       }
     }
     if(usePtWeights)
@@ -385,16 +387,25 @@ void AliFlowEventSimple::Get2Qsub(AliFlowVector* Qarray, Int_t n, TList *weights
            dWeight = pTrack->Weight();
 
             // determine Phi weight: (to be improved, I should here only access it + the treatment of gaps in the if statement)
-            if(s == 0) { //subevent 0
-             if(phiWeightsSub0 && iNbinsPhiSub0){
-               dWphi = phiWeightsSub0->GetBinContent(1+(Int_t)(TMath::Floor(dPhi*iNbinsPhiSub0/TMath::TwoPi())));
-             } 
-           } else if (s == 1) {
-             if(phiWeightsSub1 && iNbinsPhiSub1){
-               dWphi = phiWeightsSub1->GetBinContent(1+(Int_t)(TMath::Floor(dPhi*iNbinsPhiSub1/TMath::TwoPi())));
+           //subevent 0
+           if(s == 0)  { 
+             if(phiWeightsSub0 && iNbinsPhiSub0)  {
+               Int_t phiBin = 1+(Int_t)(TMath::Floor(dPhi*iNbinsPhiSub0/TMath::TwoPi()));
+               //use the phi value at the center of the bin
+               dPhi  = phiWeightsSub0->GetBinCenter(phiBin);
+               dWphi = phiWeightsSub0->GetBinContent(phiBin);
+             }
+           } 
+           //subevent 1
+           else if (s == 1) { 
+             if(phiWeightsSub1 && iNbinsPhiSub1) {
+               Int_t phiBin = 1+(Int_t)(TMath::Floor(dPhi*iNbinsPhiSub1/TMath::TwoPi()));
+               //use the phi value at the center of the bin
+               dPhi  = phiWeightsSub1->GetBinCenter(phiBin);
+               dWphi = phiWeightsSub1->GetBinContent(phiBin);
              } 
            }
-
+           
             // determine v'(pt) weight:
             if(ptWeights && dBinWidthPt)
             {
@@ -412,7 +423,7 @@ void AliFlowEventSimple::Get2Qsub(AliFlowVector* Qarray, Int_t n, TList *weights
             dQY += dWeight*dWphi*dWpt*dWeta*TMath::Sin(iOrder*dPhi);
 
             // weighted multiplicity:
-            iUsedTracks+=dWeight*dWphi*dWpt*dWeta;
+            sumOfWeights+=dWeight*dWphi*dWpt*dWeta;
 
           } // end of subevent
         } // end of if (pTrack->InRPSelection())
@@ -423,9 +434,9 @@ void AliFlowEventSimple::Get2Qsub(AliFlowVector* Qarray, Int_t n, TList *weights
       }
     } // loop over particles
     Qarray[s].Set(dQX,dQY);
-    Qarray[s].SetMult(iUsedTracks);
+    Qarray[s].SetMult(sumOfWeights);
     //reset
-    iUsedTracks = 0;
+    sumOfWeights = 0.;
     dQX = 0.;
     dQY = 0.;
   }
@@ -457,6 +468,7 @@ void AliFlowEventSimple::Print(Option_t *option) const
 //-----------------------------------------------------------------------
 void AliFlowEventSimple::Browse(TBrowser *b)
 {
+  //browse in TBrowser
   if (!b) return;
   if (!fNumberOfTracksWrap)
   {
@@ -487,6 +499,7 @@ AliFlowEventSimple::AliFlowEventSimple( TTree* inputTree,
   fMCReactionPlaneAngle(0.),
   fMCReactionPlaneAngleIsSet(kFALSE),
   fAfterBurnerPrecision(0.001),
+  fUserModified(kFALSE),
   fNumberOfTracksWrap(NULL),
   fNumberOfRPsWrap(NULL),
   fMCReactionPlaneAngleWrap(NULL)
@@ -551,6 +564,7 @@ void AliFlowEventSimple::CloneTracks(Int_t n)
       AddTrack(static_cast<AliFlowTrackSimple*>(track->Clone()));
     }
   }
+  SetUserModified();
 }
 
 //_____________________________________________________________________________
@@ -562,6 +576,7 @@ void AliFlowEventSimple::ResolutionPt(Double_t res)
     AliFlowTrackSimple* track = static_cast<AliFlowTrackSimple*>(fTrackCollection->At(i));
     if (track) track->ResolutionPt(res);
   }
+  SetUserModified();
 }
 
 //_____________________________________________________________________________
@@ -574,6 +589,8 @@ void AliFlowEventSimple::TagSubeventsInEta( Double_t etaMinA,
   for (Int_t i=0; i<fNumberOfTracks; i++)
   {
     AliFlowTrackSimple* track = static_cast<AliFlowTrackSimple*>(fTrackCollection->At(i));
+    if (!track) continue;
+    track->ResetSubEventTags();
     Double_t eta=track->Eta();
     if (eta >= etaMinA && eta <= etaMaxA) track->SetForSubevent(0);
     if (eta >= etaMinB && eta <= etaMaxB) track->SetForSubevent(1);
@@ -581,6 +598,21 @@ void AliFlowEventSimple::TagSubeventsInEta( Double_t etaMinA,
 }
 
 //_____________________________________________________________________________
+void AliFlowEventSimple::TagSubeventsByCharge()
+{
+  //Flag two subevents in given eta ranges
+  for (Int_t i=0; i<fNumberOfTracks; i++)
+  {
+    AliFlowTrackSimple* track = static_cast<AliFlowTrackSimple*>(fTrackCollection->At(i));
+    if (!track) continue;
+    track->ResetSubEventTags();
+    Int_t charge=track->Charge();
+    if (charge<0) track->SetForSubevent(0);
+    if (charge>0) track->SetForSubevent(1);
+  }
+}
+
+//_____________________________________________________________________________
 void AliFlowEventSimple::AddV1( Double_t v1 )
 {
   //add v2 to all tracks wrt the reaction plane angle
@@ -589,6 +621,7 @@ void AliFlowEventSimple::AddV1( Double_t v1 )
     AliFlowTrackSimple* track = static_cast<AliFlowTrackSimple*>(fTrackCollection->At(i));
     if (track) track->AddV1(v1, fMCReactionPlaneAngle, fAfterBurnerPrecision);
   }
+  SetUserModified();
 }
 
 //_____________________________________________________________________________
@@ -600,6 +633,19 @@ void AliFlowEventSimple::AddV2( Double_t v2 )
     AliFlowTrackSimple* track = static_cast<AliFlowTrackSimple*>(fTrackCollection->At(i));
     if (track) track->AddV2(v2, fMCReactionPlaneAngle, fAfterBurnerPrecision);
   }
+  SetUserModified();
+}
+
+//_____________________________________________________________________________
+void AliFlowEventSimple::AddV3( Double_t v3 )
+{
+  //add v3 to all tracks wrt the reaction plane angle
+  for (Int_t i=0; i<fNumberOfTracks; i++)
+  {
+    AliFlowTrackSimple* track = static_cast<AliFlowTrackSimple*>(fTrackCollection->At(i));
+    if (track) track->AddV3(v3, fMCReactionPlaneAngle, fAfterBurnerPrecision);
+  }
+  SetUserModified();
 }
 
 //_____________________________________________________________________________
@@ -611,21 +657,62 @@ void AliFlowEventSimple::AddV4( Double_t v4 )
     AliFlowTrackSimple* track = static_cast<AliFlowTrackSimple*>(fTrackCollection->At(i));
     if (track) track->AddV4(v4, fMCReactionPlaneAngle, fAfterBurnerPrecision);
   }
+  SetUserModified();
 }
 
 //_____________________________________________________________________________
-void AliFlowEventSimple::AddFlow( Double_t v1, Double_t v2, Double_t v3, Double_t v4 )
+void AliFlowEventSimple::AddV5( Double_t v5 )
+{
+  //add v4 to all tracks wrt the reaction plane angle
+  for (Int_t i=0; i<fNumberOfTracks; i++)
+  {
+    AliFlowTrackSimple* track = static_cast<AliFlowTrackSimple*>(fTrackCollection->At(i));
+    if (track) track->AddV5(v5, fMCReactionPlaneAngle, fAfterBurnerPrecision);
+  }
+  SetUserModified();
+}
+
+//_____________________________________________________________________________
+void AliFlowEventSimple::AddFlow( Double_t v1, Double_t v2, Double_t v3, Double_t v4, Double_t v5,
+                                  Double_t rp1, Double_t rp2, Double_t rp3, Double_t rp4, Double_t rp5 )
+{
+  //add flow to all tracks wrt the reaction plane angle, for all harmonic separate angle
+  for (Int_t i=0; i<fNumberOfTracks; i++)
+  {
+    AliFlowTrackSimple* track = static_cast<AliFlowTrackSimple*>(fTrackCollection->At(i));
+    if (track) track->AddFlow(v1,v2,v3,v4,v5,rp1,rp2,rp3,rp4,rp5,fAfterBurnerPrecision);
+  }
+  SetUserModified();
+}
+
+//_____________________________________________________________________________
+void AliFlowEventSimple::AddFlow( Double_t v1, Double_t v2, Double_t v3, Double_t v4, Double_t v5 )
 {
   //add flow to all tracks wrt the reaction plane angle
   for (Int_t i=0; i<fNumberOfTracks; i++)
   {
     AliFlowTrackSimple* track = static_cast<AliFlowTrackSimple*>(fTrackCollection->At(i));
-    if (track) track->AddFlow(v1,v2,v3,v4,fMCReactionPlaneAngle, fAfterBurnerPrecision);
+    if (track) track->AddFlow(v1,v2,v3,v4,v5,fMCReactionPlaneAngle, fAfterBurnerPrecision);
   }
+  SetUserModified();
 }
 
 //_____________________________________________________________________________
-void AliFlowEventSimple::TagRP( AliFlowTrackSimpleCuts* cuts )
+void AliFlowEventSimple::AddV2( TF1* ptDepV2 )
+{
+  //add v2 to all tracks wrt the reaction plane angle
+  for (Int_t i=0; i<fNumberOfTracks; i++)
+  {
+    AliFlowTrackSimple* track = static_cast<AliFlowTrackSimple*>(fTrackCollection->At(i));
+    if (!track) continue;
+    Double_t v2 = ptDepV2->Eval(track->Pt());
+    track->AddV2(v2, fMCReactionPlaneAngle, fAfterBurnerPrecision);
+  }
+  SetUserModified();
+}
+
+//_____________________________________________________________________________
+void AliFlowEventSimple::TagRP( const AliFlowTrackSimpleCuts* cuts )
 {
   //tag tracks as reference particles (RPs)
   for (Int_t i=0; i<fNumberOfTracks; i++)
@@ -633,16 +720,21 @@ void AliFlowEventSimple::TagRP( AliFlowTrackSimpleCuts* cuts )
     AliFlowTrackSimple* track = static_cast<AliFlowTrackSimple*>(fTrackCollection->At(i));
     if (!track) continue;
     Bool_t pass=cuts->PassesCuts(track);
-    track->SetForRPSelection(pass);
+    Bool_t rpTrack=track->InRPSelection();
     if (pass) 
-      fNumberOfRPs++;
+    {
+      if (!rpTrack) fNumberOfRPs++; //only increase if not already tagged
+    }
     else
-      fNumberOfRPs--;
+    {
+      if (rpTrack) fNumberOfRPs--; //only decrease if detagging
+    }
+    track->SetForRPSelection(pass);
   }
 }
 
 //_____________________________________________________________________________
-void AliFlowEventSimple::TagPOI( AliFlowTrackSimpleCuts* cuts )
+void AliFlowEventSimple::TagPOI( const AliFlowTrackSimpleCuts* cuts )
 {
   //tag tracks as particles of interest (POIs)
   for (Int_t i=0; i<fNumberOfTracks; i++)
@@ -684,8 +776,37 @@ Int_t AliFlowEventSimple::CleanUpDeadTracks()
   for (Int_t i=0; i<fNumberOfTracks; i++)
   {
     AliFlowTrackSimple* track = static_cast<AliFlowTrackSimple*>(fTrackCollection->At(i));
-    if (track->IsDead()) {delete track;ncleaned++;}
+    if (!track) continue;
+    if (track->IsDead()) {delete track;track=NULL;ncleaned++;}
   }
   fTrackCollection->Compress(); //clean up empty slots
+  fNumberOfTracks-=ncleaned; //update number of tracks
   return ncleaned;
 }
+
+//_____________________________________________________________________________
+TF1* AliFlowEventSimple::SimplePtDepV2()
+{
+  //return a standard pt dependent v2 formula, user has to clean up!
+  return new TF1("StandardPtDepV2","((x<1.0)*(0.05/1.0)*x+(x>=1.0)*0.05)");
+}
+
+//_____________________________________________________________________________
+TF1* AliFlowEventSimple::SimplePtSpectrum()
+{
+  //return a standard pt spectrum, user has to clean up!
+  return new TF1("StandardPtSpectrum","x*TMath::Exp(-pow(0.13957*0.13957+x*x,0.5)/0.4)",0.1,10.);
+}
+
+//_____________________________________________________________________________
+void AliFlowEventSimple::ClearFast()
+{
+  //clear the counter without deleting allocated objects so they can be reused
+  fReferenceMultiplicity = 0;
+  fNumberOfTracks = 0;
+  fNumberOfRPs = 0;
+  fMCReactionPlaneAngle = 0.0;
+  fMCReactionPlaneAngleIsSet = kFALSE;
+  fAfterBurnerPrecision = 0.001;
+  fUserModified = kFALSE;
+}