PWGCF/FLOW/macros/extractFlowVZERO.C

   1 extractFlowVZERO(Int_t itech=0,Int_t ibin = 3,Int_t step = 0,Int_t arm=2,Float_t pTh = 0.9,Int_t charge=0,Int_t addbins=0){
   2   TCanvas *can1 = new TCanvas("cV0Acheck","cV0Acheck");
   3
   4 /*
   5 read outVZEROv2.root or outVZEROv3.root files
   6
   7 itech=0 TPC stand alone PID and selection
   8 itech=1 TPC stand alone PID but TOF track selection
   9 itech=2 TPC&TOF PID (TOF strictly required)
  10 itech=3 TPC|TOF PID (the same of before but using TPC stand alone when TOF not available)
  11
  12 ibin = centrality bin
  13
  14 step = species (0=all charges, 1=pi, 2=K, 3=p, 4=e, 5=d, 6=t, 7=3He
  15
  16 arm = armonic 2 (v2) or 3 (v3)
  17
  18 pTh = probability threshold = 0.6 , 0.8 , 0.9 , 0.95
  19
  20 charge = 1(pos), -1(neg), 0(both)
  21
  22 addbins = to merge more centrality bin (i.e. ibin = 2, addbins = 3 merge 2,3,4,5 centrality bins) - preliminary version!!!!!!
  23
  24 */
  25
  26   // load libraries
  27   gSystem->Load("libANALYSIS");
  28   gSystem->Load("libANALYSISalice");
  29   gSystem->Load("libCORRFW.so");
  30
  31   // get objects
  32   char filein[100];
  33   sprintf(filein,"outVZEROv%i.root",arm);
  34   TFile *f=new TFile(filein);
  35   sprintf(filein,"contVZEROv%i",arm);
  36   TList *l = f->Get(filein);
  37   l->ls();
  38   AliCFContainer *c = l->At(0);
  39   AliCFContainer *c3 = l->At(1);
  40   TProfile *p1 = l->At(2);
  41   TProfile *p2 = l->At(3);
  42   TProfile *p3 = l->At(4);
  43   Int_t iVar[2] = {1,3};
  44
  45   Int_t iVarPsi[1] = {5};
  46   /* 0=centrality , 1= cos(arm*deltaphi) , 2=charge , 3=pt , 4=probability , 5=EP , 6=PIDmask*/
  47   Double_t iMin[7] = {ibin,-1,-1.5,0,pTh+0.00001,-3.14,0};
  48   Double_t iMax[7] = {ibin+addbins,1,1.5,20,1.1,3.14,1};
  49
  50   char *iTechName[4] = {"tpc","tpcInTof","tof","tpcORtof"};
  51
  52   Bool_t kOnlyTPC = kFALSE;
  53
  54   if(itech==0 && (step!=0)){
  55       kOnlyTPC = kTRUE;
  56   }
  57   else if(itech==1 && (step!=0)){
  58       iMin[6] = 2;
  59       iMax[6] = 2;
  60   }
  61   else if(itech==2 || itech==1){
  62       iMin[6] = 1;
  63       iMax[6] = 1;
  64   }
  65   else{
  66       iMin[6] = 0;
  67       iMax[6] = 1;
  68   }
  69
  70   if(charge==-1) iMax[2]=0;
  71   if(charge==1) iMin[2]=0;
  72
  73   // EP distribution needed for NUA corrections
  74   TH2F *hPhiA = l->At(5);
  75   TH2F *hPhiC = l->At(6);
  76
  77   TH1D *hPhiPA = hPhiA->ProjectionY("_py",iMin[0]+1,iMax[0]+1);
  78   TH1D *hPhiPC = hPhiC->ProjectionY("_py",iMin[0]+1,iMax[0]+1);
  79
  80   // EP resolution variables
  81   Float_t res1=0,res2=0,res3=0;
  82   Float_t eres1=0,eres2=0,eres3=0;
  83
  84   for(Int_t i=iMin[0];i<=iMax[0];i++){ // in case of more centrality bins weighting with the errors
  85     if(p1->GetBinError(i+1)){
  86       eres1 += 1./p1->GetBinError(i+1)/p1->GetBinError(i+1);
  87       res1 += p1->GetBinContent(i+1)/p1->GetBinError(i+1)/p1->GetBinError(i+1);
  88     }
  89     if(p2->GetBinError(i+1)){
  90       eres2 += 1./p2->GetBinError(i+1)/p2->GetBinError(i+1);
  91       res2 += p2->GetBinContent(i+1)/p2->GetBinError(i+1)/p2->GetBinError(i+1);
  92     }
  93     if(p3->GetBinError(i+1)){
  94       eres3 += 1./p3->GetBinError(i+1)/p3->GetBinError(i+1);
  95       res3 += p3->GetBinContent(i+1)/p3->GetBinError(i+1)/p3->GetBinError(i+1);
  96     }
  97
  98     if(eres1) res1 /= eres1;
  99     if(eres2) res2 /= eres2;
 100     if(eres3) res3 /= eres3;
 101
 102     if(eres1) eres1 = TMath::Sqrt(1./eres1);
 103     if(eres1) eres2 = TMath::Sqrt(1./eres2);
 104     if(eres1) eres3 = TMath::Sqrt(1./eres3);
 105   }
 106
 107   // NUA correction (fit to EP distribution)
 108   TF1 *fpol = new TF1("fPol","pol0");
 109   hPhiPA->Fit("fPol","","",-TMath::Pi()/arm,TMath::Pi()/arm);
 110   Float_t scalA = fPol->GetParameter(0);
 111   hPhiPC->Fit("fPol","","",-TMath::Pi()/arm,TMath::Pi()/arm);
 112   Float_t scalC = fPol->GetParameter(0);
 113
 114   AliCFContainer *c2[20];
 115   AliCFContainer *c4[20];
 116   TH2D *h[20],*h2[20];
 117
 118   AliCFContainer *c2bis[20];
 119   AliCFContainer *c4bis[20];
 120
 121   AliCFContainer *cPsi2[20];
 122   AliCFContainer *cPsi4[20];
 123   TH1D *hPsi[20],*hPsi2[20];
 124   AliCFContainer *cPsi2bis[20];
 125   AliCFContainer *cPsi4bis[20];
 126
 127   Float_t intA = 0;
 128   Float_t intC = 0;
 129
 130   for(Int_t i=5;i<15;i++){
 131     printf("%i\n",i);
 132     iMin[5] = hPhiPA->GetBinCenter(i+1);
 133     iMax[5] = hPhiPA->GetBinCenter(i+1);
 134     if(!kOnlyTPC){
 135       c2[i] = c->MakeSlice(2,iVar,iMin,iMax);
 136       c4[i]= c3->MakeSlice(2,iVar,iMin,iMax);
 137     }
 138     else{ // merge to maskPID bins (TPC standalone without TOF + TPC stand alone with TOF)
 139       iMin[6] = 0;
 140       iMax[6] = 0;
 141       c2[i] = c->MakeSlice(2,iVar,iMin,iMax);
 142       c4[i]= c3->MakeSlice(2,iVar,iMin,iMax);
 143
 144       iMin[6] = 2;
 145       iMax[6] = 2;
 146       c2bis[i] = c->MakeSlice(2,iVar,iMin,iMax);
 147       c4bis[i]= c3->MakeSlice(2,iVar,iMin,iMax);
 148    }
 149
 150
 151     if(!kOnlyTPC){
 152       cPsi2[i] = c->MakeSlice(1,iVarPsi,iMin,iMax);
 153       cPsi4[i]= c3->MakeSlice(1,iVarPsi,iMin,iMax);
 154     }
 155     else{ // merge to maskPID bins (TPC standalone without TOF + TPC stand alone with TOF)
 156       iMin[6] = 0;
 157       iMax[6] = 0;
 158       cPsi2[i] = c->MakeSlice(1,iVarPsi,iMin,iMax);
 159       cPsi4[i]= c3->MakeSlice(1,iVarPsi,iMin,iMax);
 160
 161       iMin[6] = 2;
 162       iMax[6] = 2;
 163       cPsi2bis[i] = c->MakeSlice(1,iVarPsi,iMin,iMax);
 164       cPsi4bis[i]= c3->MakeSlice(1,iVarPsi,iMin,iMax);
 165     }
 166
 167     h[i] = (TH2D *) c2[i]->Project(step,0,1);
 168     h2[i] = (TH2D *) c4[i]->Project(step,0,1);
 169     if(kOnlyTPC){
 170       TH2D *htemp = (TH2D *) c2bis[i]->Project(step,0,1);
 171       h[i]->Add(htemp);
 172       htemp = (TH2D *) c4bis[i]->Project(step,0,1);
 173       h2[i]->Add(htemp);
 174     }
 175
 176     if(hPhiPA->GetBinContent(i)) h[i]->Scale(scalA/hPhiPA->GetBinContent(i+1)); // reweighting for NUA correction
 177     if(hPhiPC->GetBinContent(i)) h2[i]->Scale(scalC/hPhiPC->GetBinContent(i+1)); // reweighting for NUA correction
 178
 179     hPsi[i] = (TH1D *) cPsi2[i]->Project(step,0);
 180     hPsi2[i] = (TH1D *) cPsi4[i]->Project(step,0);
 181
 182     if(kOnlyTPC){ // merge to maskPID bins (TPC standalone without TOF + TPC stand alone with TOF)
 183       TH1D *htemp2 = (TH1D *) cPsi2bis[i]->Project(step,0);
 184       hPsi[i]->Add(htemp2);
 185       htemp2 = (TH1D *) cPsi4bis[i]->Project(step,0);
 186       hPsi2[i]->Add(htemp2);
 187     }
 188
 189     if(hPhiPA->GetBinContent(i+1)) hPsi[i]->Scale(scalA/hPhiPA->GetBinContent(i+1)); // check of reweighting for NUA correction
 190     if(hPhiPC->GetBinContent(i+1)) hPsi2[i]->Scale(scalC/hPhiPC->GetBinContent(i+1)); // check of reweighting for NUA correction
 191
 192     intA += hPsi[i]->Integral();
 193     intC += hPsi2[i]->Integral();
 194
 195   }
 196
 197 //   hPhiPA->Scale(intA / (scalA)/10);
 198 //   hPhiPC->Scale(intC / (scalC)/10);
 199
 200   // NUA correction check V0A
 201   hPhiPA->Draw();
 202   hPsi[5]->Draw("SAME");
 203   hPsi[5]->Scale(scalA / intA * 10);
 204   for(Int_t i=6;i<15;i++){
 205     h[5]->Add(h[i]);
 206     h2[5]->Add(h2[i]);
 207     hPsi[i]->Draw("SAME");
 208     hPsi[i]->Scale(scalA / intA * 10);
 209   }
 210
 211   // NUA correction check V0C
 212   TCanvas *can2 = new TCanvas("cV0Ccheck","cV0Ccheck");
 213   hPhiPC->Draw();
 214   hPsi2[5]->Draw("SAME");
 215   for(Int_t i=6;i<15;i++){
 216     h[5]->Add(h[i]);
 217     h2[5]->Add(h2[i]);
 218     hPsi2[i]->Draw("SAME");
 219     hPsi2[i]->Scale(scalC /intC * 10);
 220   }
 221   hPsi2[5]->Scale(scalC /intC * 10);
 222
 223   // Flow for V0A and V0C separately
 224   TCanvas *can3 = new TCanvas("cFlowComp","cFlowComp");
 225   TProfile *pp = h[5]->ProfileY();
 226   pp->Draw();
 227   TProfile *pp2 = h2[5]->ProfileY();
 228   pp2->Draw("SAME");
 229
 230   printf("nev (selected within 0-80%) = %i\n",p1->GetEntries());
 231
 232
 233   // correction for resoltion
 234   Float_t scaling = sqrt(res1*res3/res2);
 235   pp->Scale(1./scaling);
 236   printf("resolution V0A = %f\n",scaling);
 237   Float_t err1_2 = eres1*eres1/res1/res1/4 +
 238                    eres2*eres2/res2/res2/4 +
 239                    eres3*eres3/res3/res3/4;
 240   Float_t err2_2 = err1_2;
 241   err1_2 /= scaling*scaling;
 242   scaling = sqrt(res2*res3/res1);
 243   err2_2 /= scaling*scaling;
 244   pp2->Scale(1./scaling);
 245   printf("resolution V0C =%f\n",scaling);
 246
 247   char title[100];
 248   sprintf(title,"VZERO EP;p_{t} (GeV/c);v_{%i}",arm);
 249   pp->SetTitle(title);
 250
 251   // Average V0A-V0C
 252   TH1D *pAll = pp->ProjectionX();
 253
 254   for(Int_t i=1;i <= pAll->GetNbinsX();i++){
 255        Float_t e1 = err1_2*pp->GetBinContent(i)*pp->GetBinContent(i) + pp->GetBinError(i)*pp->GetBinError(i);
 256        Float_t e2 = err2_2*pp2->GetBinContent(i)*pp2->GetBinContent(i) + pp2->GetBinError(i)*pp2->GetBinError(i);
 257        Float_t xval = 0,exval = 0;
 258        if(e1 >0 && e2>0){
 259          xval = (pp->GetBinContent(i)/e1 + pp2->GetBinContent(i)/e2)/(1/e1 + 1/e2);
 260          exval = 1./sqrt(1/e1 + 1/e2);
 261        }
 262        pAll->SetBinContent(i,xval);
 263        pAll->SetBinError(i,exval);
 264   }
 265   // combined measurement
 266   TCanvas *can4 = new TCanvas("cFlowCombined","cFlowCombined");
 267   pAll->Draw();
 268
 269   char name[100];
 270   sprintf(name,"out%i-%i_%i_%4.2f_%i%sv%i.root",iMin[0],iMax[0],step,pTh,charge,iTechName[itech],arm);
 271   TFile *fout = new TFile(name,"RECREATE");
 272
 273   pAll->SetName("histo");
 274   pAll->Write();
 275   can1->Write();
 276   can2->Write();
 277   can3->Write();
 278   fout->Close();
 279 }