added calculation on angle
authorpluettig <philipp.luettig@cern.ch>
Thu, 3 Apr 2014 16:04:08 +0000 (18:04 +0200)
committerpluettig <philipp.luettig@cern.ch>
Thu, 3 Apr 2014 16:04:08 +0000 (18:04 +0200)
corrected error in Ep determination
added crosscheck histograms

PWGLF/SPECTRA/ChargedHadrons/dNdPt/AlidNdPtAnalysisPbPbAOD.cxx
PWGLF/SPECTRA/ChargedHadrons/dNdPt/AlidNdPtAnalysisPbPbAOD.h
PWGLF/SPECTRA/ChargedHadrons/dNdPt/macros/AddTask_dNdPt_PbPbAOD.C

index cda9b08..05e7361 100644 (file)
@@ -72,6 +72,7 @@ fCrossCheckClusterLengthAcc(0),
 fCutSettings(0),
 fEventplaneDist(0),
 fMCEventplaneDist(0),
+fCorrelEventplaneMCDATA(0),
 //global
 fIsMonteCarlo(0),
 // event cut variables
@@ -526,6 +527,11 @@ void AlidNdPtAnalysisPbPbAOD::UserCreateOutputObjects()
   fMCEventplaneDist->GetXaxis()->SetTitle("#phi (MC event plane)");
   fMCEventplaneDist->Sumw2();
   
+  fCorrelEventplaneMCDATA = new TH2F("fCorrelEventplaneMCDATA","fCorrelEventplaneMCDATA",40, -2.*TMath::Pi(), 2.*TMath::Pi(), 40, -2.*TMath::Pi(), 2.*TMath::Pi());
+  fCorrelEventplaneMCDATA->GetXaxis()->SetTitle("#phi (event plane)");
+  fCorrelEventplaneMCDATA->GetYaxis()->SetTitle("#phi (MC event plane)");
+  fCorrelEventplaneMCDATA->Sumw2();
+  
   // Add Histos, Profiles etc to List
   fOutputList->Add(fZvPtEtaCent);
   fOutputList->Add(fDeltaphiPtEtaCent);
@@ -565,6 +571,7 @@ void AlidNdPtAnalysisPbPbAOD::UserCreateOutputObjects()
   fOutputList->Add(fCutSettings);
   fOutputList->Add(fEventplaneDist);
   fOutputList->Add(fMCEventplaneDist);
+  fOutputList->Add(fCorrelEventplaneMCDATA);
   
   StoreCutSettingsToHistogram();
   
@@ -669,7 +676,7 @@ void AlidNdPtAnalysisPbPbAOD::UserExec(Option_t *option)
        
        dMCEventZv = mcHdr->GetVtxZ();
        dMCTrackZvPtEtaCent[0] = dMCEventZv;
-       dMCEventplaneAngle = genHijingHeader->ReactionPlaneAngle();
+       dMCEventplaneAngle = MoveMCEventplane(genHijingHeader->ReactionPlaneAngle());
        fEventStatistics->Fill("MC all events",1);
        fMCEventplaneDist->Fill(dMCEventplaneAngle);
   }
@@ -680,7 +687,14 @@ void AlidNdPtAnalysisPbPbAOD::UserExec(Option_t *option)
   if( dCentrality < 0 ) return;
   fEventStatistics->Fill("after centrality selection",1);
   
+  // get event plane Angle from AODHeader, default is Q
+  ep = const_cast<AliAODEvent*>(eventAOD)->GetEventplane();
+  if(ep) {
+       dEventplaneAngle = ep->GetEventplane("V0",eventAOD);
+  }
   
+  //   cout << dEventplaneAngle << endl;
+  fEventplaneDist->Fill(dEventplaneAngle);
   
   // start with MC truth analysis
   if(fIsMonteCarlo)
@@ -719,9 +733,9 @@ void AlidNdPtAnalysisPbPbAOD::UserExec(Option_t *option)
          dMCTrackZvPtEtaCent[3] = dCentrality;
          fMCGenZvPtEtaCent->Fill(dMCTrackZvPtEtaCent);
          
-         dMCTrackPhiPtEtaCent[0] = mcPart->Phi()-dMCEventplaneAngle;
-         if( dMCTrackPhiPtEtaCent[0] < 0) dMCTrackPhiPtEtaCent[0] += 2.*TMath::Pi();
-         else if( dMCTrackPhiPtEtaCent[0] > 2.*TMath::Pi()) dMCTrackPhiPtEtaCent[0] -= 2.*TMath::Pi();
+         dMCTrackPhiPtEtaCent[0] = RotatePhi(mcPart->Phi(), dEventplaneAngle); // use eventplane and not reactionplan, similar to centrality vs impact paramter
+//       if( dMCTrackPhiPtEtaCent[0] < 0) dMCTrackPhiPtEtaCent[0] += 2.*TMath::Pi();
+//       else if( dMCTrackPhiPtEtaCent[0] > 2.*TMath::Pi()) dMCTrackPhiPtEtaCent[0] -= 2.*TMath::Pi();
          dMCTrackPhiPtEtaCent[1] = mcPart->Pt();
          dMCTrackPhiPtEtaCent[2] = mcPart->Eta();
          dMCTrackPhiPtEtaCent[3] = dCentrality;
@@ -767,14 +781,7 @@ void AlidNdPtAnalysisPbPbAOD::UserExec(Option_t *option)
   
   dTrackZvPtEtaCent[0] = dEventZv;
   
-  // get event plane Angle from AODHeader, default is Q
-  ep = const_cast<AliAODEvent*>(eventAOD)->GetEventplane();
-  if(ep) {
-       dEventplaneAngle = ep->GetEventplane("V0",eventAOD);
-  }
   
-  cout << dEventplaneAngle << endl;
-  fEventplaneDist->Fill(dEventplaneAngle);
   
   if(AreRelativeCutsEnabled())
   {
@@ -810,9 +817,10 @@ void AlidNdPtAnalysisPbPbAOD::UserExec(Option_t *option)
        dTrackZvPtEtaCent[2] = track->Eta();
        dTrackZvPtEtaCent[3] = dCentrality;
        
-       dTrackPhiPtEtaCent[0] = track->Phi() - dEventplaneAngle;
-       if( dTrackPhiPtEtaCent[0] < 0) dTrackPhiPtEtaCent[0] += 2.*TMath::Pi();
-       else if( dTrackPhiPtEtaCent[0] > 2.*TMath::Pi()) dTrackPhiPtEtaCent[0] -= 2.*TMath::Pi();
+       dTrackPhiPtEtaCent[0] = RotatePhi(track->Phi(), dEventplaneAngle); 
+       
+       //      if( dTrackPhiPtEtaCent[0] < -1.0*TMath::Pi()) dTrackPhiPtEtaCent[0] += 2.*TMath::Pi();
+       //      else if( dTrackPhiPtEtaCent[0] > TMath::Pi()) dTrackPhiPtEtaCent[0] -= 2.*TMath::Pi();
        dTrackPhiPtEtaCent[1] = track->Pt();
        dTrackPhiPtEtaCent[2] = track->Eta();
        dTrackPhiPtEtaCent[3] = dCentrality;
@@ -834,9 +842,10 @@ void AlidNdPtAnalysisPbPbAOD::UserExec(Option_t *option)
          dMCTrackZvPtEtaCent[2] = mcPart->Eta();
          dMCTrackZvPtEtaCent[3] = dCentrality;
          
-         dMCTrackPhiPtEtaCent[0] = mcPart->Phi()-dMCEventplaneAngle;
-         if( dMCTrackPhiPtEtaCent[0] < 0) dMCTrackPhiPtEtaCent[0] += 2.*TMath::Pi();
-         else if( dMCTrackPhiPtEtaCent[0] > 2.*TMath::Pi()) dMCTrackPhiPtEtaCent[0] -= 2.*TMath::Pi();
+         dMCTrackPhiPtEtaCent[0] = RotatePhi(mcPart->Phi(), dEventplaneAngle); // use eventplane and not reactionplan, similar to centrality vs impact paramter
+         
+         //      if( dMCTrackPhiPtEtaCent[0] < -1.0*TMath::Pi()) dMCTrackPhiPtEtaCent[0] += 2.*TMath::Pi();
+         //      else if( dMCTrackPhiPtEtaCent[0] > TMath::Pi()) dMCTrackPhiPtEtaCent[0] -= 2.*TMath::Pi();
          dMCTrackPhiPtEtaCent[1] = mcPart->Pt();
          dMCTrackPhiPtEtaCent[2] = mcPart->Eta();
          dMCTrackPhiPtEtaCent[3] = dCentrality;
@@ -919,12 +928,75 @@ void AlidNdPtAnalysisPbPbAOD::UserExec(Option_t *option)
   Double_t dEventZvMultCent[3] = {dEventZv, iAcceptedMultiplicity, dCentrality};
   fZvMultCent->Fill(dEventZvMultCent);
   
+  // store correlation between data and MC eventplane
+  if(fIsMonteCarlo) fCorrelEventplaneMCDATA->Fill(dEventplaneAngle, dMCEventplaneAngle);
+  
   PostData(1, fOutputList);
   
   // delete pointers:
   
 }
 
+Double_t AlidNdPtAnalysisPbPbAOD::MoveMCEventplane(Double_t dMCEP)
+{
+  Double_t retval = 0;
+  retval = dMCEP;
+  
+  if( (dMCEP > 0) && (dMCEP < 1./2.*TMath::Pi()) ) 
+  {
+       return retval;
+  }
+  
+  if( (dMCEP >= 1./2.*TMath::Pi()) && (dMCEP <= 3./2.*TMath::Pi()) )
+  {
+       retval -= TMath::Pi();
+       return retval;
+  }
+  
+  if(dMCEP > 3./2.*TMath::Pi())
+  {
+       retval -= 2.*TMath::Pi();
+       return retval;
+  }
+  
+  return -9999.;
+}
+
+Double_t AlidNdPtAnalysisPbPbAOD::RotatePhi(Double_t phiTrack, Double_t phiEP)
+{
+  Double_t dPhi = 0;
+  dPhi = phiTrack - phiEP;
+  if ((dPhi >= -1./2. * TMath::Pi() ) && 
+         (dPhi <= 1./2. * TMath::Pi() ) )
+  {
+       return dPhi;
+  }
+  
+  if( (dPhi < 0) )
+  {
+       dPhi += 2.*TMath::Pi();
+  }
+  
+  if ((dPhi > 0) && 
+         (dPhi > 1./2. * TMath::Pi() ) && 
+         (dPhi <= 3./2. * TMath::Pi() ) )
+  {
+       dPhi -= TMath::Pi();
+       return dPhi;
+  }    
+  
+  if ((dPhi > 0) && 
+         (dPhi > 3./2. * TMath::Pi() )) 
+  {
+       dPhi -= 2.*TMath::Pi();
+       return dPhi;
+  }
+  
+//   Printf("[E] dphi = %.4f , phiTrack = %.4f, phiEP = %.4f", dPhi, phiTrack, phiEP);
+  
+  return -9999.;
+}
+
 Bool_t AlidNdPtAnalysisPbPbAOD::SetRelativeCuts(AliAODEvent *event)
 {
   //
@@ -1044,7 +1116,7 @@ Bool_t AlidNdPtAnalysisPbPbAOD::IsTrackAccepted(AliAODTrack *tr, Double_t dCentr
   
   Double_t dLengthInTPC = 0;
   if ( DoCutLengthInTPCPtDependent() ) { dLengthInTPC = dummy.GetLengthInActiveZone(&par,3,236, bMagZ ,0,0); }
-
+  
   Double_t dNClustersTPC = tr->GetTPCNcls();
   Double_t dCrossedRowsTPC = tr->GetTPCClusterInfo(2,1);
   Double_t dFindableClustersTPC = tr->GetTPCNclsF();
index 55ae1e3..905847c 100644 (file)
@@ -153,6 +153,8 @@ class AlidNdPtAnalysisPbPbAOD : public AliAnalysisTaskSE {
     AliGenHijingEventHeader* GetHijingEventHeader(AliAODMCHeader *header);
     AliGenPythiaEventHeader* GetPythiaEventHeader(AliAODMCHeader *header);
     
+       Double_t RotatePhi(Double_t phiTrack, Double_t phiEP);
+       Double_t MoveMCEventplane(Double_t dMCEP);
     
     Bool_t SetRelativeCuts(AliAODEvent *event);
     
@@ -206,6 +208,7 @@ class AlidNdPtAnalysisPbPbAOD : public AliAnalysisTaskSE {
     
     TH1F               *fEventplaneDist; // event plane distribution in phi
     TH1F               *fMCEventplaneDist; // MC event plane distribution in phi
+    TH2F               *fCorrelEventplaneMCDATA; // correlation between data and MC eventplane
     // global variables
     Bool_t fIsMonteCarlo;
     
index 0d6e920..8eddc00 100644 (file)
@@ -44,7 +44,10 @@ AlidNdPtAnalysisPbPbAOD *AddTask_dNdPt_PbPbAOD( UInt_t uTriggerMask = AliVEvent:
   Int_t nBinPtCheck = sizeof(binsPtCheck)/sizeof(Double_t);
   task->SetBinsPtCheck(nBinPtCheck, binsPtCheck);
   
-  Double_t binsPhi[] = {0 ,0.10472 ,0.20944 ,0.314159 ,0.418879 ,0.523599 ,0.628319 ,0.733038 ,0.837758 ,0.942478 ,1.0472 ,1.15192 ,1.25664 ,1.36136 ,1.46608 ,1.5708 ,1.67552 ,1.78024 ,1.88496 ,1.98968 ,2.0944 ,2.19911 ,2.30383 ,2.40855 ,2.51327 ,2.61799 ,2.72271 ,2.82743 ,2.93215 ,3.03687 ,3.14159 ,3.24631 ,3.35103 ,3.45575 ,3.56047 ,3.66519 ,3.76991 ,3.87463 ,3.97935 ,4.08407 ,4.18879 ,4.29351 ,4.39823 ,4.50295 ,4.60767 ,4.71239 ,4.81711 ,4.92183 ,5.02655 ,5.13127 ,5.23599 ,5.34071 ,5.44543 ,5.55015 ,5.65487 ,5.75959 ,5.86431 ,5.96903 ,6.07375 ,6.17847 ,6.28319};
+//   Double_t binsPhi[] = {0 ,0.10472 ,0.20944 ,0.314159 ,0.418879 ,0.523599 ,0.628319 ,0.733038 ,0.837758 ,0.942478 ,1.0472 ,1.15192 ,1.25664 ,1.36136 ,1.46608 ,1.5708 ,1.67552 ,1.78024 ,1.88496 ,1.98968 ,2.0944 ,2.19911 ,2.30383 ,2.40855 ,2.51327 ,2.61799 ,2.72271 ,2.82743 ,2.93215 ,3.03687 ,3.14159 ,3.24631 ,3.35103 ,3.45575 ,3.56047 ,3.66519 ,3.76991 ,3.87463 ,3.97935 ,4.08407 ,4.18879 ,4.29351 ,4.39823 ,4.50295 ,4.60767 ,4.71239 ,4.81711 ,4.92183 ,5.02655 ,5.13127 ,5.23599 ,5.34071 ,5.44543 ,5.55015 ,5.65487 ,5.75959 ,5.86431 ,5.96903 ,6.07375 ,6.17847 ,6.28319};
+  
+  Double_t binsPhi[] = {-1.*TMath::Pi(), -2.97625, -2.8109, -2.64555, -2.4802, -2.31486, -2.14951, -1.98416, -1.81882, -1.65347, -1.48812, -1.32278, -1.15743, -0.992082, -0.826735, -0.661388, -0.496041, -0.330694, -0.165347, 0, 0.165347, 0.330694, 0.496041, 0.661388, 0.826735, 0.992082, 1.15743, 1.32278, 1.48812, 1.65347, 1.81882, 1.98416, 2.14951, 2.31486, 2.4802, 2.64555, 2.8109, 2.97625, TMath::Pi()};
+  
   Int_t nBinPhi = sizeof(binsPhi)/sizeof(Double_t);
   task->SetBinsPhi(nBinPhi, binsPhi);