Possibility to get the TPC event plane resolution from 3 subevents (Giacomo)
authorprino <prino@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 4 Jun 2012 13:40:52 +0000 (13:40 +0000)
committerprino <prino@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 4 Jun 2012 13:40:52 +0000 (13:40 +0000)
PWGHF/vertexingHF/charmFlow/AliAnalysisTaskSEHFv2.cxx
PWGHF/vertexingHF/charmFlow/AliAnalysisTaskSEHFv2.h
PWGHF/vertexingHF/charmFlow/DmesonsFlowAnalysis.C

index 635f9c9..9e80b0c 100644 (file)
@@ -97,7 +97,7 @@ AliAnalysisTaskSE(),
   fMinCentr(20),
   fMaxCentr(80),
   fEtaGap(kFALSE),
-  fSubEvents(0)
+  fSubEvents(2)
 {
   // Default constructor
 }
@@ -122,7 +122,7 @@ AliAnalysisTaskSEHFv2::AliAnalysisTaskSEHFv2(const char *name,AliRDHFCuts *rdCut
   fMinCentr(20),
   fMaxCentr(80),
   fEtaGap(kFALSE),
-  fSubEvents(0)
+  fSubEvents(2)
 {
   // standard constructor
   Int_t pdg=421;
@@ -429,6 +429,7 @@ void AliAnalysisTaskSEHFv2::UserExec(Option_t */*option*/)
   AliEventplane *pl=aod->GetEventplane();
   if(!pl){
     AliError("AliAnalysisTaskSEHFv2::UserExec:no eventplane! v2 analysis without eventplane not possible!\n");
+    fhEventsInfo->Fill(6);
     return;
   }
   
@@ -438,6 +439,7 @@ void AliAnalysisTaskSEHFv2::UserExec(Option_t */*option*/)
   Double_t rpangleVZERO=0;
   Double_t planereso=0;
   Double_t deltaPsi=0;
+  Double_t rpangleeventC=0;
   Double_t rpangleeventB=0;
   Double_t rpangleeventA=0;
 
@@ -463,54 +465,44 @@ void AliAnalysisTaskSEHFv2::UserExec(Option_t */*option*/)
     delete g;g=0x0;
     if(fUseAfterBurner)fAfterBurner->SetEventPlane((Double_t)eventplane);
   }else{
-    if(fEventPlaneMeth!=kTPC){
-
-      //VZERO EP and resolution
-      rpangleVZERO=GetPhi0Pi(pl->GetEventplane("V0",aod,fV0EPorder));
-      if(fEventPlaneMeth>kTPCVZERO){
-               //Using V0A/C for VZERO resolution
-               if(fEventPlaneMeth==kVZEROA){
-                 rpangleVZERO=GetPhi0Pi(pl->GetEventplane("V0A",aod,fV0EPorder));
-         if(fSubEvents!=kSingleV0Side)rpangleeventB=GetPhi0Pi(pl->GetEventplane("V0C",aod,fV0EPorder));
-               }
-               else if(fEventPlaneMeth==kVZEROC){
-         rpangleVZERO=GetPhi0Pi(pl->GetEventplane("V0C",aod,fV0EPorder));
-         if(fSubEvents!=kSingleV0Side)rpangleeventB=GetPhi0Pi(pl->GetEventplane("V0A",aod,fV0EPorder));
-       }
-               // deltaPsi =rpangleeventA-rpangleeventB;
-               eventplane=rpangleVZERO;
-      }
-    }
-
-    // TPC event plane
+   // TPC event plane
     rpangleTPC = pl->GetEventplane("Q");
     if(rpangleTPC<0){
       fhEventsInfo->Fill(6);
       return;
     }
     rpangleeventA = rpangleTPC;
-    if(fSubEvents>kFullTPC||fEtaGap||fEventPlaneMeth==kVZERO){
+    if(fSubEvents==2||fEventPlaneMeth==kVZERO){
       qsub1 = pl->GetQsub1();
       qsub2 = pl->GetQsub2();
       if(!qsub1 || !qsub2){
        fhEventsInfo->Fill(6);
        return;
       }
-      if(fSubEvents==kPosTPC||fSubEvents==kSingleV0Side||fEventPlaneMeth==kVZERO){
-       rpangleeventA = qsub1->Phi()/2.;
-       if(fSubEvents==kSingleV0Side||fEventPlaneMeth==kVZERO)rpangleeventB = qsub2->Phi()/2.;
+      rpangleeventA = qsub1->Phi()/2.;
+      rpangleeventB = qsub2->Phi()/2.;
+    }
+    if(fEventPlaneMeth!=kTPC){
+      //VZERO EP and resolution
+      rpangleVZERO=GetPhi0Pi(pl->GetEventplane("V0",aod,fV0EPorder));
+      rpangleeventC=rpangleVZERO;
+      if(fEventPlaneMeth==kVZEROA||fEventPlaneMeth==kVZEROC||(fEventPlaneMeth==kTPCVZERO&&fSubEvents==3)){
+       rpangleeventB=GetPhi0Pi(pl->GetEventplane("V0A",aod,fV0EPorder));
+       rpangleeventC=GetPhi0Pi(pl->GetEventplane("V0C",aod,fV0EPorder));
+       if(fEventPlaneMeth==kVZEROA)rpangleVZERO=rpangleeventB;
+       else if(fEventPlaneMeth==kVZEROC)rpangleVZERO=rpangleeventC;
       }
-      if(fSubEvents==kNegTPC)rpangleeventA = qsub2->Phi()/2.;
     }
-    
-    if(fEventPlaneMeth<=kTPCVZERO){
+
+
+    if(fEventPlaneMeth>kTPCVZERO)eventplane=rpangleVZERO;
+    else{
       q = pl->GetQVector();
-      deltaPsi = pl->GetQsubRes();
-    }else{
-      deltaPsi=eventplane-rpangleeventA;
+      eventplane=rpangleTPC;
     }
+    deltaPsi =rpangleeventA-rpangleeventB;
   }
-  
+
   if(TMath::Abs(deltaPsi)>TMath::Pi()/2.){
     if(deltaPsi>0.) deltaPsi-=TMath::Pi();
     else deltaPsi +=TMath::Pi();
@@ -526,14 +518,14 @@ void AliAnalysisTaskSEHFv2::UserExec(Option_t */*option*/)
     ((TH1F*)fOutput->FindObject(Form("hEvPlaneA%s",centrbinname.Data())))->Fill(rpangleeventA); //Angle of first subevent
     ((TH1F*)fOutput->FindObject(Form("hEvPlaneB%s",centrbinname.Data())))->Fill(rpangleeventB); //Angle of second subevent
   }
-  if(fEventPlaneMeth>kTPCVZERO){
-    Double_t deltaSub=rpangleeventA-rpangleeventB;
+  if(fEventPlaneMeth>kTPCVZERO||fSubEvents==3){
+    Double_t deltaSub=rpangleeventA-rpangleeventC;
     if(TMath::Abs(deltaSub)>TMath::Pi()/2.){// difference of subevents reaction plane angle cannot be bigger than phi/2
       if(deltaSub>0.) deltaSub-=TMath::Pi();
       else deltaSub +=TMath::Pi();
     } 
     ((TH1F*)fOutput->FindObject(Form("hEvPlaneReso2%s",centrbinname.Data())))->Fill(TMath::Cos(2.*deltaSub)); //RP resolution   
-    deltaSub =eventplane-rpangleeventB;
+    deltaSub =rpangleeventB-rpangleeventC;
     if(TMath::Abs(deltaSub)>TMath::Pi()/2.){// difference of subevents reaction plane angle cannot be bigger than phi/2
       if(deltaSub>0.) deltaSub-=TMath::Pi();
       else deltaSub +=TMath::Pi();
@@ -779,7 +771,18 @@ void AliAnalysisTaskSEHFv2::SetEventPlaneMethod(Int_t method){
   }
   fEventPlaneMeth=method;
 }
-
+//________________________________________________________________________
+void AliAnalysisTaskSEHFv2::SetNTPCSubEvents(Int_t nsub){
+  if(nsub>3||nsub<2){
+    AliWarning("Only 2 or 3 subevents implemented. Setting 2 subevents\n");
+    nsub=2;
+  }
+  if(fEventPlaneMeth==kTPC&&nsub==3){
+    AliWarning("V0 must be enabled to run 3 TPC subevents. Enabling...\n");
+    fEventPlaneMeth=kTPCVZERO;
+  }
+  fSubEvents=nsub;
+}
 //________________________________________________________________________
 Float_t AliAnalysisTaskSEHFv2::GetPhi0Pi(Float_t phi){
   // Sets the phi angle in the range 0-pi
index 9851dc3..2f43b0a 100644 (file)
@@ -35,8 +35,8 @@ class AliAnalysisTaskSEHFv2 : public AliAnalysisTaskSE
 
   enum DecChannel{kDplustoKpipi,kD0toKpi,kDstartoKpipi}; //more particles can be added
   enum EventPlaneMeth{kTPC,kTPCVZERO,kVZERO,kVZEROA,kVZEROC}; //Event plane to be calculated in the task
-  enum SubEvents{kFullTPC,kPosTPC,kNegTPC,kSingleV0Side}; //Sub-events for V0 EP
-
+  //  enum SubEvents{kFullTPC,kPosTPC,kNegTPC,kSingleV0Side}; //Sub-events for V0 EP
+  
   AliAnalysisTaskSEHFv2();
   AliAnalysisTaskSEHFv2(const char *name, AliRDHFCuts *rdCuts, Int_t decaychannel);
  
@@ -52,6 +52,10 @@ class AliAnalysisTaskSEHFv2 : public AliAnalysisTaskSE
   void SetUseAfterBurner(Bool_t ab){fUseAfterBurner=ab;}
   void SetAfterBurner(AliHFAfterBurner *ab){fAfterBurner=ab;}
   void SetEtaGapFeatureForEventplaneFromTracks (Bool_t etaGap) {fEtaGap = etaGap;}
+
+  void SetNTPCSubEvents(Int_t nsub);
+  void Set2TPCEPSubEvents(){SetNTPCSubEvents(2);}
+  void Set3TPCEPSubEvents(){SetNTPCSubEvents(3);}
   void SetEventPlaneMethod(Int_t epmethod);
   void SetTPCEPOnly(){SetEventPlaneMethod(kTPC);}
   void SetVZEROEP(){SetEventPlaneMethod(kVZERO);}
@@ -59,9 +63,10 @@ class AliAnalysisTaskSEHFv2 : public AliAnalysisTaskSE
   void SetVZEROCEP(){SetEventPlaneMethod(kVZEROC);}
   void SetTPCEP(){SetEventPlaneMethod(kTPCVZERO);}
   void SetEventPlanesCompatibility(Float_t comp) {fEventPlanesComp=comp;}
-  void SetSubEvents(Int_t subev){if(subev>=kFullTPC&&subev<=kSingleV0Side)fSubEvents=subev;}
+  //void SetSubEvents(Int_t subev){if(subev>=kFullTPC&&subev<=kSingleV0Side)fSubEvents=subev;}
 
   Int_t GetEventPlaneMethod()const {return fEventPlaneMeth;}
+  Int_t GetNTPCSubEvents()const {return fSubEvents;}
   Float_t GetEventPlanesCompatibility()const {return fEventPlanesComp;}
   Float_t GetUpperMassLimit()const {return fUpmasslimit;}
   Float_t GetLowerMassLimit()const {return fLowmasslimit;}
@@ -107,7 +112,7 @@ class AliAnalysisTaskSEHFv2 : public AliAnalysisTaskSE
   Int_t fMinCentr;              //minimum centrality
   Int_t fMaxCentr;              //maximum centrality
   Bool_t fEtaGap;               // Eta gap feature for Eventplane from tracks; be careful that you do the correct settings in AddTaskEventPlane.C !!!!
-  Int_t fSubEvents;             //Sub-events definition for V0 EP
+  Int_t fSubEvents;             //Sub-events definition for TPC EP
 
   ClassDef(AliAnalysisTaskSEHFv2,2); // AliAnalysisTaskSE for the HF v2 analysis
 };
index df8f30d..85313dd 100644 (file)
 
 #endif
 
+/* $Id$ */
+
 //methods for the analysis of AliAnalysisTaskSEHFv2 output
 //Authors: Chiara Bianchin cbianchi@pd.infn.it
 //         Giacomo Ortona  ortona@to.infn.it
 //         Francesco Prino prino@to.infn.it
 
+//evPlane flag: -1=V0C,0=V0,1=V0A,2=TPC2subevs,3=TPC3subevs
 //global variables to be set
 Bool_t gnopng=kTRUE; //don't save in png format (only root and eps)
 const Int_t nptbinsnew=3;
@@ -49,16 +52,15 @@ Float_t phibinslim[nphibins+1]={0,TMath::Pi()/4,TMath::Pi()/2,3*TMath::Pi()/4,TM
 Int_t minPtBin[nptbinsnew]={-1,-1,-1};
 Int_t maxPtBin[nptbinsnew]={-1,-1,-1};
 Double_t mass;
+
 //methods
-//Bool_t ReadFile(TList* &list,TH1F* &hstat,AliRDHFCuts* &cutobj,TString listname,TString partname,TString path="./",TString filename="AnalysisResults.root");
 TList *LoadMassHistos(TList *inputlist,Int_t minCent,Int_t maxCent,Bool_t inoutanis);
 Int_t FindPtBin(Int_t nbins, Float_t* array,Float_t value);
 void FillSignalGraph(TList *histlist,TGraphAsymmErrors **gSignal,TGraphAsymmErrors **gSignalfs,Bool_t inoutanis);
 void DmesonsFlowAnalysis(Bool_t inoutanis=kTRUE,Int_t minC=30,Int_t maxC=50,TString partname="Dplus",Int_t evPlane=2);
 void DrawEventPlane(Int_t mincentr=0,Int_t maxcentr=0,TString filename="AnalysisResults.root",TString dirname="PWG3_D2H_HFv2",TString listname="coutputv2",Int_t evPlane=2);
-Double_t DrawEventPlane(TList *list,Int_t mincentr=0,Int_t maxcentr=0,Int_t evPlane=2);
+Double_t DrawEventPlane(TList *list,Int_t mincentr=0,Int_t maxcentr=0,Int_t evPlane=2,Double_t &error);
 void DrawEventPlaneAllCentralities(TList *list,Int_t startCentr=20,Int_t endCentr=80,Int_t evPlane=2);
-//Aggiungere <pt> method
 
 //methods implementation
 //________________________________________________________________________________
@@ -216,6 +218,11 @@ void DmesonsFlowAnalysis(Bool_t inoutanis,Int_t minC,Int_t maxC,TString partname
   TString dirname="PWG3_D2H_HFv2";
   TString listname="coutputv2";
 
+  if(evPlane>3||evPlane<-1){
+    printf("Wrong EP flag %d \n",evPlane);
+    return;
+  }
+
   AliRDHFCuts *cutsobj=0x0;
   //Load input data from AliAnalysisTaskSEHFv2
   TFile *f = TFile::Open(filename.Data());
@@ -259,30 +266,30 @@ void DmesonsFlowAnalysis(Bool_t inoutanis,Int_t minC,Int_t maxC,TString partname
 
   Int_t nphi=nphibins;
   if(inoutanis)nphi=2;
-  Int_t nSubRes=1;
-  if(evPlane<2)nSubRes=3;//3 sub-events method for VZERO EP
 
   //EP resolution
   TH1F* hevplresos=(TH1F*)list->FindObject(Form("hEvPlaneResocentr%d_%d",minC,minC+5));
   TH1F* hevplresos2=0;
   TH1F* hevplresos3=0;
-  if(evPlane<2){
+  if(evPlane!=2){
     hevplresos2=(TH1F*)list->FindObject(Form("hEvPlaneReso2centr%d_%d",minC,minC+5));
     hevplresos3=(TH1F*)list->FindObject(Form("hEvPlaneReso3centr%d_%d",minC,minC+5));
   }
   Double_t resol=0;
   for(Int_t icent=minC+5;icent<maxC;icent=icent+5){
     hevplresos->Add((TH1F*)list->FindObject(Form("hEvPlaneResocentr%d_%d",icent,icent+5)));
-    if(evPlane<2){
+    if(evPlane!=2){
       hevplresos2->Add((TH1F*)list->FindObject(Form("hEvPlaneReso2centr%d_%d",icent,icent+5)));
       hevplresos3->Add((TH1F*)list->FindObject(Form("hEvPlaneReso3centr%d_%d",icent,icent+5)));
     }
   }
-  if(evPlane>=2){
+  if(evPlane==2){
     resol=AliVertexingHFUtils::GetFullEvResol(hevplresos);
   }else{
     Double_t resolSub[3]={hevplresos->GetMean(),hevplresos2->GetMean(),hevplresos3->GetMean()};
-    resol=TMath::Sqrt(resolSub[0]*resolSub[2]/resolSub[1]);
+    if(evPlane<=0)resol=TMath::Sqrt(resolSub[1]*resolSub[2]/resolSub[0]);
+    else if(evPlane==1) resol=TMath::Sqrt(resolSub[0]*resolSub[2]/resolSub[1]);
+    else if(evPlane==3) resol=TMath::Sqrt(resolSub[0]*resolSub[1]/resolSub[2]);
   }
 
   printf("average pt for pt bin \n");
@@ -422,10 +429,10 @@ void DrawEventPlane(Int_t mincentr,Int_t maxcentr,TString filename,TString dirna
   if(!list){
     printf("list %s not found in file, please check list name\n",listname.Data());return;
   }
-  DrawEventPlane(list,mincentr,maxcentr,evPlane);
+  DrawEventPlane(list,mincentr,maxcentr,evPlane,0);
 }
 //_______________________________________________________________________________________________________________________________
-Double_t DrawEventPlane(TList *list,Int_t mincentr,Int_t maxcentr,Int_t evPlane){
+Double_t DrawEventPlane(TList *list,Int_t mincentr,Int_t maxcentr,Int_t evPlane,Double_t &error){
   
   //draw the histograms correlated to the event plane, returns event plane resolution
   
@@ -437,7 +444,7 @@ Double_t DrawEventPlane(TList *list,Int_t mincentr,Int_t maxcentr,Int_t evPlane)
 
   Double_t resolFull=0; 
   Int_t nSubRes=1;
-  if(evPlane<2)nSubRes=3;//3 sub-events method for VZERO EP
+  if(evPlane!=2)nSubRes=3;//3 sub-events method for VZERO EP
   TString namereso[3]={"Reso","Reso2","Reso3"};
   TString suffixcentr=Form("centr%d_%d",mincentr,maxcentr);
   TH2F* hevpls=(TH2F*)list->FindObject(Form("hEvPlanecentr%d_%d",mincentr,mincentr+5));
@@ -483,9 +490,9 @@ Double_t DrawEventPlane(TList *list,Int_t mincentr,Int_t maxcentr,Int_t evPlane)
   if(nSubRes>1){
     hevplresos[1]->SetLineColor(2);
     hevplresos[2]->SetLineColor(3);
-    hevplresos[0]->SetTitle("cos2(#Psi_{V0A}-#Psi_{V0B})");
-    hevplresos[1]->SetTitle("cos2(#Psi_{TPC}-#Psi_{V0A})");
-    hevplresos[2]->SetTitle("cos2(#Psi_{TPC}-#Psi_{V0B})");
+    hevplresos[0]->SetTitle("cos2(#Psi_{TPC}-#Psi_{V0A})");
+    hevplresos[1]->SetTitle("cos2(#Psi_{TPC}-#Psi_{V0C})");
+    hevplresos[2]->SetTitle("cos2(#Psi_{V0A}-#Psi_{V0C})");
     hevplresos[1]->Draw("SAME");
     hevplresos[2]->Draw("SAME");
   }
@@ -493,11 +500,30 @@ Double_t DrawEventPlane(TList *list,Int_t mincentr,Int_t maxcentr,Int_t evPlane)
   leg->SetLineColor(0);
   leg->SetFillColor(0);
   
-  if(nSubRes<=1) resolFull=AliVertexingHFUtils::GetFullEvResol(hevplresos[0]);
+  if(nSubRes<=1){
+    resolFull=AliVertexingHFUtils::GetFullEvResol(hevplresos[0]);
+    error = TMath::Abs(resolFull-AliVertexingHFUtils::GetSubEvResolLowLim(hevplresos[0]));
+  }
   else{
     Double_t resolSub[3];
-    for(Int_t ires=0;ires<nSubRes;ires++)resolSub[ires]=hevplresos[ires]->GetMean();
-    resolFull=TMath::Sqrt(resolSub[2]*resolSub[0]/resolSub[1]);
+    Double_t errors[3];
+    for(Int_t ires=0;ires<nSubRes;ires++){
+      resolSub[ires]=hevplresos[ires]->GetMean();
+      errors[ires]=hevplresos[ires]->GetMeanError();
+    }
+    Double_t lowlim[3];for(Int_t ie=0;ie<3;ie++)lowlim[ie]=TMath::Abs(resolSub[ie]-errors[ie]);
+    if(evPlane<=0){
+      resolFull=TMath::Sqrt(resolSub[1]*resolSub[2]/resolSub[0]);
+      error=resolFull-TMath::Sqrt(lowlim[2]*lowlim[1]/lowlim[0]);
+    }
+    else if(evPlane==1){
+      resolFull=TMath::Sqrt(resolSub[0]*resolSub[2]/resolSub[1]);
+      error=resolFull-TMath::Sqrt(lowlim[2]*lowlim[0]/lowlim[1]);
+    }
+    else if(evPlane==3){
+      resolFull=TMath::Sqrt(resolSub[0]*resolSub[1]/resolSub[2]);
+      error=resolFull-TMath::Sqrt(lowlim[0]*lowlim[1]/lowlim[2]);
+    }
   }
   
   TPaveText* pvreso=new TPaveText(0.1,0.1,0.6,0.2,"NDC");
@@ -516,17 +542,18 @@ Double_t DrawEventPlane(TList *list,Int_t mincentr,Int_t maxcentr,Int_t evPlane)
 
 //_______________________________________________________________________________________________________________________________
 void DrawEventPlaneAllCentralities(TList *list,Int_t startCentr,Int_t endCentr,Int_t evPlane){
-  TGraph* gresovscentr=new TGraph(0);
+  TGraphErrors* gresovscentr=new TGraphErrors(0);
   gresovscentr->SetName("gresovscentr");
   gresovscentr->SetTitle(Form("Resolution vs Centrality;centrality (%s);Resolution","%"));
 
   Int_t iPoint=0;
   Int_t step=5;//5% centrality step
   //  Int_t nCC=(endCentr-startCentr)/step;
-
+  Double_t errory=0;
   for(Int_t i=startCentr;i<endCentr;i=i+step){
-    Double_t resolFull=DrawEventPlane(list,i,i+step,evPlane);
+    Double_t resolFull=DrawEventPlane(list,i,i+step,evPlane,errory);
     gresovscentr->SetPoint(iPoint,i+(Float_t)step/2.,resolFull);
+    gresovscentr->SetPointError(iPoint,5./2.,errory);
     iPoint++;
   }