Q vector selection in MC

diff --git a/PWGLF/SPECTRA/PiKaPr/TestAOD/AliAnalysisTaskSpectraAllChAOD.cxx b/PWGLF/SPECTRA/PiKaPr/TestAOD/AliAnalysisTaskSpectraAllChAOD.cxx
index 6e3e26e7da6d846a4e8811cde6f541306100749d..be60e0032a15dcfa648028dcbd6bd217bab20b8a 100644 (file)
--- a/PWGLF/SPECTRA/PiKaPr/TestAOD/AliAnalysisTaskSpectraAllChAOD.cxx
+++ b/PWGLF/SPECTRA/PiKaPr/TestAOD/AliAnalysisTaskSpectraAllChAOD.cxx
@@ -126,11 +126,11 @@ void AliAnalysisTaskSpectraAllChAOD::UserCreateOutputObjects()
   fOutput->Add(NSparseHistTrk);
   
   //dimensions of THnSparse for stack
-  const Int_t nvarst=6;
-  //                                             pt          cent    IDgen        isph        y    etaselected
-  Int_t    binsHistRealSt[nvarst] = {      nptBins,  fnCentBins,        3,         2,        2,             1};
-  Double_t xminHistRealSt[nvarst] = {         0.,           0.,      -0.5,      -0.5,    -0.5,            0.5};
-  Double_t xmaxHistRealSt[nvarst] = {       10.,        100.,      2.5,       1.5,      0.5,           1.5};
+  const Int_t nvarst=7;
+  //                                             pt          cent    IDgen        isph        y    etaselected       Qvec gen
+  Int_t    binsHistRealSt[nvarst] = {      nptBins,  fnCentBins,        3,         2,        2,             1,      fnQvecBins};
+  Double_t xminHistRealSt[nvarst] = {         0.,           0.,      -0.5,      -0.5,    -0.5,            0.5,              0.};
+  Double_t xmaxHistRealSt[nvarst] = {       10.,        100.,      2.5,       1.5,      0.5,           1.5,      fQvecUpperLim};
   THnSparseF* NSparseHistSt = new THnSparseF("NSparseHistSt","NSparseHistSt",nvarst,binsHistRealSt,xminHistRealSt,xmaxHistRealSt);
   NSparseHistSt->GetAxis(0)->SetTitle("#it{p}_{T,gen}");
   NSparseHistSt->SetBinEdges(0,ptBins);
@@ -145,14 +145,16 @@ void AliAnalysisTaskSpectraAllChAOD::UserCreateOutputObjects()
   NSparseHistSt->GetAxis(4)->SetName("y");
   NSparseHistSt->GetAxis(5)->SetTitle("etaselected");
   NSparseHistSt->GetAxis(5)->SetName("etaselected");
+  NSparseHistSt->GetAxis(6)->SetTitle("Q vec gen");
+  NSparseHistSt->GetAxis(6)->SetName("Q_vec_gen");
   fOutput->Add(NSparseHistSt);
   
   //dimensions of THnSparse for the normalization
-  const Int_t nvarev=3;
-  //                                             cent             Q vec                Nch
-  Int_t    binsHistRealEv[nvarev] = {    fnCentBins,      fnQvecBins,           fnNchBins};
-  Double_t xminHistRealEv[nvarev] = {           0.,               0.,                   0.};
-  Double_t xmaxHistRealEv[nvarev] = {       100.,  fQvecUpperLim,               2000.};
+  const Int_t nvarev=4;
+  //                                             cent             Q vec                Nch    Qvec gen
+  Int_t    binsHistRealEv[nvarev] = {    fnCentBins,      fnQvecBins,           fnNchBins,      fnQvecBins};
+  Double_t xminHistRealEv[nvarev] = {           0.,               0.,                   0.,             0.};
+  Double_t xmaxHistRealEv[nvarev] = {       100.,  fQvecUpperLim,               2000.,       fQvecUpperLim};
   THnSparseF* NSparseHistEv = new THnSparseF("NSparseHistEv","NSparseHistEv",nvarev,binsHistRealEv,xminHistRealEv,xmaxHistRealEv);
   NSparseHistEv->GetAxis(0)->SetTitle(Form("%s cent",fEventCuts->GetCentralityMethod().Data()));
   NSparseHistEv->GetAxis(0)->SetName(Form("%s_cent",fEventCuts->GetCentralityMethod().Data()));
@@ -160,6 +162,8 @@ void AliAnalysisTaskSpectraAllChAOD::UserCreateOutputObjects()
   NSparseHistEv->GetAxis(1)->SetName("Q_vec");
   NSparseHistEv->GetAxis(2)->SetTitle("N charged");
   NSparseHistEv->GetAxis(2)->SetName("N_ch");
+  NSparseHistEv->GetAxis(3)->SetTitle("Q vec gen");
+  NSparseHistEv->GetAxis(3)->SetName("Q_vec_gen");
   fOutput->Add(NSparseHistEv);
   
   PostData(1, fOutput  );
@@ -190,13 +194,18 @@ void AliAnalysisTaskSpectraAllChAOD::UserExec(Option_t *)
   //Default TPC priors
   //if(fHelperPID->GetPIDType()==kBayes)fHelperPID->GetPIDCombined()->SetDefaultTPCPriors();//FIXME maybe this can go in the UserCreateOutputObject?
   
-  Double_t Qvec=0.;//in case of MC we save space in the memory
-  if(!fIsMC){
+  Double_t Qvec=0.;
+  if(fIsQvecCalibMode){
+    if(fVZEROside==0)Qvec=fEventCuts->GetqV0A();
+    else if (fVZEROside==1)Qvec=fEventCuts->GetqV0C();
+  }
+  else Qvec=fEventCuts->GetQvecPercentile(fVZEROside);
+  
+  Double_t QvecMC = 0.;
+  if(fIsMC){
     if(fIsQvecCalibMode){
-      if(fVZEROside==0)Qvec=fEventCuts->GetqV0A();
-      else if (fVZEROside==1)Qvec=fEventCuts->GetqV0C();
+      QvecMC = fEventCuts->CalculateQVectorMC(fVZEROside);
     }
-    else Qvec=fEventCuts->GetQvecPercentile(fVZEROside);
   }
   
   Double_t Cent=fEventCuts->GetCent();
@@ -221,13 +230,14 @@ void AliAnalysisTaskSpectraAllChAOD::UserExec(Option_t *)
          if(partMC->Eta()>=fTrackCuts->GetEtaMin() && partMC->Eta()<=fTrackCuts->GetEtaMax())etaselected=1.;
          
          //pt     cent    IDgen        isph        y
-         Double_t varSt[6];
+         Double_t varSt[7];
          varSt[0]=partMC->Pt();
          varSt[1]=Cent;
          varSt[2]=(Double_t)fHelperPID->GetParticleSpecies(partMC);
          varSt[3]=(Double_t)partMC->IsPhysicalPrimary();
          varSt[4]=partMC->Y();
          varSt[5]=etaselected;
+         varSt[6]=QvecMC;
          ((THnSparseF*)fOutput->FindObject("NSparseHistSt"))->Fill(varSt);//stack loop
          
          //Printf("a particle");
@@ -311,10 +321,11 @@ void AliAnalysisTaskSpectraAllChAOD::UserExec(Option_t *)
     Nch++;
   } // end loop on tracks
   
-  Double_t varEv[3];
+  Double_t varEv[4];
   varEv[0]=Cent;
   varEv[1]=Qvec;
   varEv[2]=Nch;
+  varEv[3]=QvecMC;
   ((THnSparseF*)fOutput->FindObject("NSparseHistEv"))->Fill(varEv);//event loop
   
   PostData(1, fOutput  );

diff --git a/PWGLF/SPECTRA/PiKaPr/TestAOD/AliSpectraAODEventCuts.cxx b/PWGLF/SPECTRA/PiKaPr/TestAOD/AliSpectraAODEventCuts.cxx
index 4d8e758e97828907d0bca62530dd2c4839d9423c..bfae432442320507147121c5e5fbb93f5c12924e 100644 (file)
--- a/PWGLF/SPECTRA/PiKaPr/TestAOD/AliSpectraAODEventCuts.cxx
+++ b/PWGLF/SPECTRA/PiKaPr/TestAOD/AliSpectraAODEventCuts.cxx
@@ -644,6 +644,45 @@ Double_t AliSpectraAODEventCuts::GetQvecPercentile(Int_t v0side){
   return percentile;
 }
 
+//______________________________________________________
+Double_t AliSpectraAODEventCuts::CalculateQVectorMC(Int_t v0side){
+  
+  if(!fIsMC) return -999.;
+  
+  // 1. get MC array
+  TClonesArray *arrayMC = (TClonesArray*) fAOD->GetList()->FindObject(AliAODMCParticle::StdBranchName());
+  
+  if (!arrayMC) AliFatal("Error: MC particles branch not found!\n");
+  
+  Double_t Qx2mc = 0., Qy2mc = 0.;
+  Int_t mult2mc = 0;
+  
+  Int_t nMC = arrayMC->GetEntries();
+      
+  // 2. loop on generated
+  for (Int_t iMC = 0; iMC < nMC; iMC++){
+    AliAODMCParticle *partMC = (AliAODMCParticle*) arrayMC->At(iMC);
+  
+    // 3. set VZERO side - FIXME Add TPC!
+    Double_t EtaVZERO[2] = {-999.,-999.};
+    if(v0side==0/*V0A*/){ EtaVZERO[0] = 2.8; EtaVZERO[1] = 5.1; } 
+    if(v0side==1/*V0C*/){ EtaVZERO[0] = -3.7; EtaVZERO[1] = -1.7; } 
+    
+    if(partMC->Eta()<EtaVZERO[0] || partMC->Eta()>EtaVZERO[1]) continue;
+    
+    // 4. Calculate Qvec components
+    
+    Qx2mc += TMath::Cos(2.*partMC->Phi());
+    Qy2mc += TMath::Sin(2.*partMC->Phi());
+    mult2mc++;
+    
+  }
+  
+  // 5. return q vector
+  return TMath::Sqrt((Qx2mc*Qx2mc + Qy2mc*Qy2mc)/mult2mc);
+  
+}
+
 //______________________________________________________
 Bool_t AliSpectraAODEventCuts::CheckSplineArray(TObjArray * splarr){
   //check TSpline array consistency

diff --git a/PWGLF/SPECTRA/PiKaPr/TestAOD/AliSpectraAODEventCuts.h b/PWGLF/SPECTRA/PiKaPr/TestAOD/AliSpectraAODEventCuts.h
index 6546c85ec1b3642a8f47ca31fc4cbc8e0621a8cb..f412f6e3b0f9bed2778b89e42afdefa30f6754dc 100644 (file)
--- a/PWGLF/SPECTRA/PiKaPr/TestAOD/AliSpectraAODEventCuts.h
+++ b/PWGLF/SPECTRA/PiKaPr/TestAOD/AliSpectraAODEventCuts.h
@@ -156,6 +156,8 @@ class AliSpectraAODEventCuts : public TNamed
   Bool_t CheckSplineArray(TObjArray * splarr);
   TObjArray *GetSplineArrayV0A() { return fSplineArrayV0A; }
   TObjArray *GetSplineArrayV0C() { return fSplineArrayV0C; }
+  
+  Double_t CalculateQVectorMC(Int_t v0side);
 
  private: