From: prino Date: Mon, 4 Jun 2012 13:40:52 +0000 (+0000) Subject: Possibility to get the TPC event plane resolution from 3 subevents (Giacomo) X-Git-Url: http://git.uio.no/git/?a=commitdiff_plain;h=b356f0ee6ed70422c83c68a0e58cfcd0470a9650;p=u%2Fmrichter%2FAliRoot.git Possibility to get the TPC event plane resolution from 3 subevents (Giacomo) --- diff --git a/PWGHF/vertexingHF/charmFlow/AliAnalysisTaskSEHFv2.cxx b/PWGHF/vertexingHF/charmFlow/AliAnalysisTaskSEHFv2.cxx index 635f9c9cac7..9e80b0c1950 100644 --- a/PWGHF/vertexingHF/charmFlow/AliAnalysisTaskSEHFv2.cxx +++ b/PWGHF/vertexingHF/charmFlow/AliAnalysisTaskSEHFv2.cxx @@ -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 diff --git a/PWGHF/vertexingHF/charmFlow/AliAnalysisTaskSEHFv2.h b/PWGHF/vertexingHF/charmFlow/AliAnalysisTaskSEHFv2.h index 9851dc3bcdd..2f43b0a04af 100644 --- a/PWGHF/vertexingHF/charmFlow/AliAnalysisTaskSEHFv2.h +++ b/PWGHF/vertexingHF/charmFlow/AliAnalysisTaskSEHFv2.h @@ -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 }; diff --git a/PWGHF/vertexingHF/charmFlow/DmesonsFlowAnalysis.C b/PWGHF/vertexingHF/charmFlow/DmesonsFlowAnalysis.C index df8f30d8c85..85313ddade2 100644 --- a/PWGHF/vertexingHF/charmFlow/DmesonsFlowAnalysis.C +++ b/PWGHF/vertexingHF/charmFlow/DmesonsFlowAnalysis.C @@ -32,11 +32,14 @@ #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 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;icentAdd((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;iresGetMean(); - resolFull=TMath::Sqrt(resolSub[2]*resolSub[0]/resolSub[1]); + Double_t errors[3]; + for(Int_t ires=0;iresGetMean(); + 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;iSetPoint(iPoint,i+(Float_t)step/2.,resolFull); + gresovscentr->SetPointError(iPoint,5./2.,errory); iPoint++; }