PWGCF/FLOW/macros/extractFlowVZERO.C

   1 TFile *fo = NULL;
   2 TF1 *fPsi = NULL;
   3
   4 Bool_t kLib = kFALSE;
   5 Bool_t kVZEROrescorr = kTRUE; // VZERO resolution corrections
   6 Bool_t kNUAcorr = kTRUE; // NUA corrections
   7 Bool_t kNUOcorr = kFALSE; // Not Uniform Occupancy (TOF) corrections
   8 Bool_t kPIDcorr = kFALSE; // PID corrections
   9 Bool_t kFEEDcorr = kFALSE; // antiproton feed-down corrections
  10 Bool_t kCleanMemory = kTRUE; // if kFALSE take care of the large numbers of TCanvases
  11
  12 TH1D *hNUO[8]; // NUO corrections
  13
  14 void LoadLib();
  15 void extractFlowVZERO(Int_t icentr,const char *type="",Int_t arm=2,Float_t pTh=0.8,Bool_t isMC=kFALSE,Int_t addbin=0);
  16 TProfile *extractFlowVZEROsingle(Int_t icentr,Int_t spec,Int_t arm=2,Bool_t isMC=kFALSE,Float_t pTh=0.9,Int_t addbin=0,const char *nameSp="",Float_t detMin=0,Float_t detMax=1,Int_t chMin=-1,Int_t chMax=1);
  17 void compareTPCTOF(Int_t icentr,Int_t spec,Int_t arm=2,Float_t pTh=0.9,Int_t addbin=0);
  18 void plotQApid(Int_t ic,Float_t pt,Int_t addbin=0);
  19
  20 void LoadLib(){
  21   if(! kLib){
  22     gSystem->Load("libVMC.so");
  23     gSystem->Load("libPhysics.so");
  24     gSystem->Load("libTree.so");
  25     gSystem->Load("libSTEERBase.so");
  26     gSystem->Load("libANALYSIS.so");
  27     gSystem->Load("libAOD.so");
  28     gSystem->Load("libESD.so");
  29     gSystem->Load("libANALYSIS.so");
  30     gSystem->Load("libANALYSISalice.so");
  31     gSystem->Load("libCORRFW.so");
  32     gSystem->Load("libNetx.so");
  33     gSystem->Load("libPWGflowBase.so");
  34     gSystem->Load("libPWGflowTasks.so");
  35
  36     kLib = kTRUE;
  37   }
  38 }
  39
  40 void extractFlowVZERO(Int_t icentr,const char *type,Int_t arm,Float_t pTh,Bool_t isMC,Int_t addbin){
  41   LoadLib();
  42
  43   char name[100];
  44   snprintf(name,100,"AnalysisResults%s.root",type);
  45   if(!fo) fo = new TFile(name);
  46
  47   new TCanvas();
  48
  49   Int_t cMin[] = {0,5,10,20,30,40,50,60,70};
  50   Int_t cMax[] = {5,10,20,30,40,50,60,70,80};
  51
  52   if(kNUOcorr){ // Compute correction for NUO in TOF
  53     compareTPCTOF(icentr,0,arm,pTh,addbin);
  54 //     compareTPCTOF(icentr,1,arm,pTh,addbin);
  55 //     compareTPCTOF(icentr,2,arm,pTh,addbin);
  56 //     compareTPCTOF(icentr,3,arm,pTh,addbin);
  57 //     compareTPCTOF(icentr,4,arm,pTh,addbin);
  58   }
  59
  60   TProfile *pAll;
  61   pAll=extractFlowVZEROsingle(icentr,0,arm,0,pTh,addbin,"all",0,1);
  62   pAll->SetMarkerStyle(24);
  63   TProfile *pPiTOF,*pPiTPC,*pPiTPC2;
  64   pPiTOF=extractFlowVZEROsingle(icentr,1,arm,0,pTh,addbin,"piTOF",1,1);
  65   pPiTPC=extractFlowVZEROsingle(icentr,1,arm,0,pTh,addbin,"piTPC",0,0);
  66   pPiTPC2=extractFlowVZEROsingle(icentr,1,arm,0,pTh,addbin,"piTPC2",2,2);
  67   pPiTPC->Add(pPiTPC2);
  68
  69   TH1D *hPi = pPiTOF->ProjectionX("hPi");
  70   hPi->SetLineColor(4);
  71   hPi->SetMarkerColor(4);
  72   hPi->SetMarkerStyle(20);
  73   for(Int_t i=1;i <=hPi->GetNbinsX();i++){
  74     Float_t x = hPi->GetBinCenter(i);
  75     if(x < 0.2){
  76       hPi->SetBinContent(i,0);
  77       hPi->SetBinError(i,0);
  78     }
  79     else if(x < 0.5){
  80       hPi->SetBinContent(i,pPiTPC->GetBinContent(i));
  81       hPi->SetBinError(i,pPiTPC->GetBinError(i));
  82     }
  83     else{
  84       if(kNUOcorr){
  85         hPi->SetBinContent(i,pPiTOF->GetBinContent(i) + hNUO[0]->GetBinContent(i));
  86         hPi->SetBinError(i,pPiTOF->GetBinError(i));
  87       }
  88       else{
  89         hPi->SetBinContent(i,pPiTOF->GetBinContent(i));
  90         hPi->SetBinError(i,pPiTOF->GetBinError(i));
  91       }
  92     }
  93   }
  94
  95   TProfile *pElTOF,*pElTPC,*pElTPC2;
  96   pElTOF=extractFlowVZEROsingle(icentr,4,arm,0,pTh,addbin,"piTOF",1,1);
  97   pElTPC=extractFlowVZEROsingle(icentr,4,arm,0,pTh,addbin,"piTPC",0,0);
  98   pElTPC2=extractFlowVZEROsingle(icentr,4,arm,0,pTh,addbin,"piTPC2",2,2);
  99   pElTPC->Add(pElTPC2);
 100
 101   TH1D *hEl = pElTOF->ProjectionX("hEl");
 102   hEl->SetLineColor(6);
 103   hEl->SetMarkerColor(6);
 104   hEl->SetMarkerStyle(20);
 105   for(Int_t i=1;i <=hEl->GetNbinsX();i++){
 106     Float_t x = hEl->GetBinCenter(i);
 107     if(x < 0.2){
 108       hEl->SetBinContent(i,0);
 109       hEl->SetBinError(i,0);
 110     }
 111     else if(x < 0.3){
 112       hEl->SetBinContent(i,pElTPC->GetBinContent(i));
 113       hEl->SetBinError(i,pElTPC->GetBinError(i));
 114     }
 115     else{
 116       if(kNUOcorr){
 117         hEl->SetBinContent(i,pElTOF->GetBinContent(i) + hNUO[0]->GetBinContent(i));
 118         hEl->SetBinError(i,pElTOF->GetBinError(i));
 119       }
 120       else{
 121         hEl->SetBinContent(i,pElTOF->GetBinContent(i));
 122         hEl->SetBinError(i,pElTOF->GetBinError(i));
 123       }
 124     }
 125   }
 126
 127   TProfile *pKTOF,*pKTPC,*pKTPC2;
 128   pKTOF=extractFlowVZEROsingle(icentr,2,arm,0,pTh,addbin,"kaTOF",1,1);
 129   pKTPC=extractFlowVZEROsingle(icentr,2,arm,0,pTh,addbin,"kaTPC",0,0);
 130   pKTPC2=extractFlowVZEROsingle(icentr,2,arm,0,pTh,addbin,"kaTPC2",2,2);
 131   pKTPC->Add(pKTPC2);
 132
 133   TH1D *hK = pKTOF->ProjectionX("hKa");
 134   hK->SetLineColor(1);
 135   hK->SetMarkerColor(1);
 136   hK->SetMarkerStyle(21);
 137   for(Int_t i=1;i <=hK->GetNbinsX();i++){
 138     Float_t x = hK->GetBinCenter(i);
 139     if(x < 0.25){
 140       hK->SetBinContent(i,0);
 141       hK->SetBinError(i,0);
 142     }
 143     else if(x < 0.45){
 144       hK->SetBinContent(i,pKTPC->GetBinContent(i));
 145       hK->SetBinError(i,pKTPC->GetBinError(i));
 146     }
 147     else{
 148       if(kNUOcorr){
 149         hK->SetBinContent(i,pKTOF->GetBinContent(i) + hNUO[0]->GetBinContent(i));
 150         hK->SetBinError(i,pKTOF->GetBinError(i));
 151       }
 152       else{
 153         hK->SetBinContent(i,pKTOF->GetBinContent(i));
 154         hK->SetBinError(i,pKTOF->GetBinError(i));
 155       }
 156     }
 157   }
 158
 159   TProfile *pPrTOF,*pPrTOF2,*pPrTPC,*pPrTPC2;
 160   pPrTOF=extractFlowVZEROsingle(icentr,3,arm,0,pTh,addbin,"prTOF",1,1,-1,-1);
 161   pPrTOF2=extractFlowVZEROsingle(icentr,3,arm,0,pTh,addbin,"prTOF2",1,1,-1,1);
 162   pPrTPC=extractFlowVZEROsingle(icentr,3,arm,0,pTh,addbin,"prTPC",0,0,-1,-1);
 163   pPrTPC2=extractFlowVZEROsingle(icentr,3,arm,0,pTh,addbin,"prTPC2",2,2,-1,-1);
 164   pPrTPC->Add(pPrTPC2);
 165
 166   TH1D *hPr = pPrTOF->ProjectionX("hPr");
 167   hPr->SetLineColor(2);
 168   hPr->SetMarkerColor(2);
 169   hPr->SetMarkerStyle(22);
 170   for(Int_t i=1;i <=hPr->GetNbinsX();i++){
 171     Float_t x = hPr->GetBinCenter(i);
 172     if(x < 0.3){
 173       hPr->SetBinContent(i,0);
 174       hPr->SetBinError(i,0);
 175     }
 176     else if(x < 1.0){
 177       hPr->SetBinContent(i,pPrTPC->GetBinContent(i));
 178       hPr->SetBinError(i,pPrTPC->GetBinError(i));
 179     }
 180     else{
 181       if(x < 3){
 182         if(kNUOcorr){
 183           hPr->SetBinContent(i,pPrTOF->GetBinContent(i) + hNUO[0]->GetBinContent(i));
 184           hPr->SetBinError(i,pPrTOF->GetBinError(i));
 185         }
 186         else{
 187           hPr->SetBinContent(i,pPrTOF->GetBinContent(i));
 188           hPr->SetBinError(i,pPrTOF->GetBinError(i));
 189         }
 190       }
 191       else{
 192         if(kNUOcorr){
 193           hPr->SetBinContent(i,pPrTOF2->GetBinContent(i) + hNUO[0]->GetBinContent(i));
 194           hPr->SetBinError(i,pPrTOF2->GetBinError(i));
 195         }
 196         else{
 197           hPr->SetBinContent(i,pPrTOF2->GetBinContent(i));
 198           hPr->SetBinError(i,pPrTOF2->GetBinError(i));
 199         }
 200       }
 201     }
 202   }
 203
 204   pAll->Draw();
 205   hPi->Draw("SAME");
 206   hK->Draw("SAME");
 207   hPr->Draw("SAME");
 208
 209
 210   char name[100];
 211
 212
 213   // PID correction
 214   if(kPIDcorr){
 215     TFile *fPidTOF = new TFile("$ALICE_ROOT/PWGCF/FLOW/other/BayesianPIDcontTPCTOF.root");
 216     TFile *fPidTPC = new TFile("$ALICE_ROOT/PWGCF/FLOW/other/BayesianPIDcontTPC.root");
 217     // pi histos
 218     sprintf(name,"Pi_IDas_Picentr%i_pth%4.2f",icentr,pTh);
 219     TH1D *hPidPiPi=(TH1D *) fPidTOF->Get(name);
 220     sprintf(name,"Pi_IDas_Elcentr%i_pth%4.2f",icentr,pTh);
 221     TH1D *hPidPiEl=(TH1D *) fPidTOF->Get(name);
 222     sprintf(name,"Pi_IDas_Kacentr%i_pth%4.2f",icentr,pTh);
 223     TH1D *hPidPiKa=(TH1D *) fPidTOF->Get(name);
 224     sprintf(name,"Pi_IDas_Prcentr%i_pth%4.2f",icentr,pTh);
 225     TH1D *hPidPiPr=(TH1D *) fPidTOF->Get(name);
 226     TH1D *hPidAll = new TH1D(*hPidPiPi);
 227     hPidAll->Add(hPidPiKa);
 228     hPidAll->Add(hPidPiPr);
 229     hPidAll->Add(hPidPiEl);
 230     hPidPiPi->Divide(hPidAll);
 231     hPidPiKa->Divide(hPidAll);
 232     hPidPiPr->Divide(hPidAll);
 233     hPidPiEl->Divide(hPidAll);
 234
 235     sprintf(name,"Pi_IDas_Picentr%i_pth%4.2f",icentr,pTh);
 236     TH1D *hPidPiPiTPC=(TH1D *) fPidTPC->Get(name);
 237     sprintf(name,"Pi_IDas_Elcentr%i_pth%4.2f",icentr,pTh);
 238     TH1D *hPidPiElTPC=(TH1D *) fPidTPC->Get(name);
 239     sprintf(name,"Pi_IDas_Kacentr%i_pth%4.2f",icentr,pTh);
 240     TH1D *hPidPiKaTPC=(TH1D *) fPidTPC->Get(name);
 241     sprintf(name,"Pi_IDas_Prcentr%i_pth%4.2f",icentr,pTh);
 242     TH1D *hPidPiPrTPC=(TH1D *) fPidTPC->Get(name);
 243     hPidAll->Reset();
 244     hPidAll->Add(hPidPiPiTPC);
 245     hPidAll->Add(hPidPiKaTPC);
 246     hPidAll->Add(hPidPiPrTPC);
 247     hPidAll->Add(hPidPiElTPC);
 248     hPidPiPiTPC->Divide(hPidAll);
 249     hPidPiKaTPC->Divide(hPidAll);
 250     hPidPiPrTPC->Divide(hPidAll);
 251     hPidPiElTPC->Divide(hPidAll);
 252
 253     // K histos
 254     sprintf(name,"Ka_IDas_Picentr%i_pth%4.2f",icentr,pTh);
 255     TH1D *hPidKaPi=(TH1D *) fPidTOF->Get(name);
 256     sprintf(name,"Ka_IDas_Elcentr%i_pth%4.2f",icentr,pTh);
 257     TH1D *hPidKaEl=(TH1D *) fPidTOF->Get(name);
 258     sprintf(name,"Ka_IDas_Kacentr%i_pth%4.2f",icentr,pTh);
 259     TH1D *hPidKaKa=(TH1D *) fPidTOF->Get(name);
 260     sprintf(name,"Ka_IDas_Prcentr%i_pth%4.2f",icentr,pTh);
 261     TH1D *hPidKaPr=(TH1D *) fPidTOF->Get(name);
 262     hPidAll->Reset();
 263     hPidAll->Add(hPidKaPi);
 264     hPidAll->Add(hPidKaKa);
 265     hPidAll->Add(hPidKaPr);
 266     hPidAll->Add(hPidKaEl);
 267     hPidKaPi->Divide(hPidAll);
 268     hPidKaKa->Divide(hPidAll);
 269     hPidKaPr->Divide(hPidAll);
 270     hPidKaEl->Divide(hPidAll);
 271
 272     sprintf(name,"Ka_IDas_Picentr%i_pth%4.2f",icentr,pTh);
 273     TH1D *hPidKaPiTPC=(TH1D *) fPidTPC->Get(name);
 274     sprintf(name,"Ka_IDas_Elcentr%i_pth%4.2f",icentr,pTh);
 275     TH1D *hPidKaElTPC=(TH1D *) fPidTPC->Get(name);
 276     sprintf(name,"Ka_IDas_Kacentr%i_pth%4.2f",icentr,pTh);
 277     TH1D *hPidKaKaTPC=(TH1D *) fPidTPC->Get(name);
 278     sprintf(name,"Ka_IDas_Prcentr%i_pth%4.2f",icentr,pTh);
 279     TH1D *hPidKaPrTPC=(TH1D *) fPidTPC->Get(name);
 280     hPidAll->Reset();
 281     hPidAll->Add(hPidKaPiTPC);
 282     hPidAll->Add(hPidKaKaTPC);
 283     hPidAll->Add(hPidKaPrTPC);
 284     hPidAll->Add(hPidKaElTPC);
 285     hPidKaPiTPC->Divide(hPidAll);
 286     hPidKaKaTPC->Divide(hPidAll);
 287     hPidKaPrTPC->Divide(hPidAll);
 288     hPidKaElTPC->Divide(hPidAll);
 289
 290     // pr histos
 291     sprintf(name,"Pr_IDas_Picentr%i_pth%4.2f",icentr,pTh);
 292     TH1D *hPidPrPi=(TH1D *) fPidTOF->Get(name);
 293     sprintf(name,"Pr_IDas_Elcentr%i_pth%4.2f",icentr,pTh);
 294     TH1D *hPidPrEl=(TH1D *) fPidTOF->Get(name);
 295     sprintf(name,"Pr_IDas_Kacentr%i_pth%4.2f",icentr,pTh);
 296     TH1D *hPidPrKa=(TH1D *) fPidTOF->Get(name);
 297     sprintf(name,"Pr_IDas_Prcentr%i_pth%4.2f",icentr,pTh);
 298     TH1D *hPidPrPr=(TH1D *) fPidTOF->Get(name);
 299     hPidAll->Reset();
 300     hPidAll->Add(hPidPrPi);
 301     hPidAll->Add(hPidPrKa);
 302     hPidAll->Add(hPidPrPr);
 303     hPidAll->Add(hPidPrEl);
 304     hPidPrPi->Divide(hPidAll);
 305     hPidPrKa->Divide(hPidAll);
 306     hPidPrPr->Divide(hPidAll);
 307     hPidPrEl->Divide(hPidAll);
 308
 309     sprintf(name,"Pr_IDas_Picentr%i_pth%4.2f",icentr,pTh);
 310     TH1D *hPidPrPiTPC=(TH1D *) fPidTPC->Get(name);
 311     sprintf(name,"Pr_IDas_Elcentr%i_pth%4.2f",icentr,pTh);
 312     TH1D *hPidPrElTPC=(TH1D *) fPidTPC->Get(name);
 313     sprintf(name,"Pr_IDas_Kacentr%i_pth%4.2f",icentr,pTh);
 314     TH1D *hPidPrKaTPC=(TH1D *) fPidTPC->Get(name);
 315     sprintf(name,"Pr_IDas_Prcentr%i_pth%4.2f",icentr,pTh);
 316     TH1D *hPidPrPrTPC=(TH1D *) fPidTPC->Get(name);
 317     hPidAll->Reset();
 318     hPidAll->Add(hPidPrPiTPC);
 319     hPidAll->Add(hPidPrKaTPC);
 320     hPidAll->Add(hPidPrPrTPC);
 321     hPidAll->Add(hPidPrElTPC);
 322     hPidPrPiTPC->Divide(hPidAll);
 323     hPidPrKaTPC->Divide(hPidAll);
 324     hPidPrPrTPC->Divide(hPidAll);
 325     hPidPrElTPC->Divide(hPidAll);
 326
 327     for(Int_t k=1;k <=hPi->GetNbinsX();k++){
 328       Float_t pt = hPi->GetBinCenter(k);
 329       Float_t xPi=hPi->GetBinContent(k)*hPidPiPi->Interpolate(pt) + hK->GetBinContent(k)*hPidPiKa->Interpolate(pt) + hPr->GetBinContent(k)*hPidPiPr->Interpolate(pt) + hEl->GetBinContent(k)*hPidPiEl->Interpolate(pt);
 330       if(pt < 0.5)
 331         xPi=hPi->GetBinContent(k)*hPidPiPiTPC->Interpolate(pt) + hK->GetBinContent(k)*hPidPiKaTPC->Interpolate(pt) + hPr->GetBinContent(k)*hPidPiPrTPC->Interpolate(pt) + hEl->GetBinContent(k)*hPidPiElTPC->Interpolate(pt);
 332       Float_t xKa=hPi->GetBinContent(k)*hPidKaPi->Interpolate(pt) + hK->GetBinContent(k)*hPidKaKa->Interpolate(pt) + hPr->GetBinContent(k)*hPidKaPr->Interpolate(pt) + hEl->GetBinContent(k)*hPidKaEl->Interpolate(pt);
 333       if(pt < 0.45)
 334         xKa=hPi->GetBinContent(k)*hPidKaPiTPC->Interpolate(pt) + hK->GetBinContent(k)*hPidKaKaTPC->Interpolate(pt) + hPr->GetBinContent(k)*hPidKaPrTPC->Interpolate(pt) + hEl->GetBinContent(k)*hPidKaElTPC->Interpolate(pt);
 335       Float_t xPr=hPi->GetBinContent(k)*hPidPrPi->Interpolate(pt) + hK->GetBinContent(k)*hPidPrKa->Interpolate(pt) + hPr->GetBinContent(k)*hPidPrPr->Interpolate(pt) + hEl->GetBinContent(k)*hPidPrEl->Interpolate(pt);
 336       if(pt < 1)
 337         xPr=hPi->GetBinContent(k)*hPidPrPiTPC->Interpolate(pt) + hK->GetBinContent(k)*hPidPrKaTPC->Interpolate(pt) + hPr->GetBinContent(k)*hPidPrPrTPC->Interpolate(pt) + hEl->GetBinContent(k)*hPidPrElTPC->Interpolate(pt);
 338
 339       hPi->SetBinContent(k,hPi->GetBinContent(k)*2 - xPi);
 340       hK->SetBinContent(k,hK->GetBinContent(k)*2 - xKa);
 341       hPr->SetBinContent(k,hPr->GetBinContent(k)*2 - xPr);
 342     }
 343   }
 344
 345
 346   // antiproton Feed down
 347   TProfile *pFromLambda = extractFlowVZEROsingle(icentr,11,arm,0,pTh,addbin,"pFromLambda",1,1);
 348   TProfile *piFromK = extractFlowVZEROsingle(icentr,12,arm,0,pTh,addbin,"piFromK",1,1,1,1);
 349   TProfile *pFromLambda2 = extractFlowVZEROsingle(icentr,11,arm,0,0.6,addbin,"pFromLambdanoPID",0,1);
 350   TProfile *piFromK2 = extractFlowVZEROsingle(icentr,12,arm,0,0.6,addbin,"piFromKnoPID",0,1);
 351   TProfile *piFromK3 = extractFlowVZEROsingle(icentr,12,arm,0,0.6,addbin,"piFromKnoPIDtof",1,1);
 352
 353   TH1D *hFeedSyst = NULL;
 354
 355   if(kFEEDcorr){
 356     hFeedSyst = new TH1D(*hPr);
 357     hFeedSyst->SetName("hFeedSyst");
 358     hFeedSyst->Reset();
 359     for(Int_t k=1;k <=hPr->GetNbinsX();k++){
 360       Float_t contam = 3.23174e-01 * TMath::Exp(- 9.46743e-01 *  hPr->GetBinCenter(k));
 361       Float_t corr = contam * pFromLambda->GetBinContent(k)/(1-contam);
 362       Float_t corrErr = contam * pFromLambda->GetBinError(k)/(1-contam);
 363       Float_t value = hPr->GetBinContent(k)/(1-contam) - corr;
 364       Float_t valueErr = hPr->GetBinError(k)/(1-contam);
 365
 366       hFeedSyst->SetBinContent(k,hPr->GetBinContent(k) - value);
 367       hFeedSyst->SetBinContent(k,sqrt(corrErr*corrErr + valueErr*valueErr - hPr->GetBinError(k)*hPr->GetBinError(k)));
 368
 369       hPr->SetBinContent(k,value);
 370       hPr->SetBinError(k,sqrt(corrErr*corrErr + valueErr*valueErr));
 371
 372     }
 373     hFeedSyst->Divide(hPr);
 374   }
 375
 376   // write output
 377   snprintf(name,100,"results%03i-%03iv%i_pTh%3.1f%s.root",cMin[icentr],cMax[icentr+addbin],arm,pTh,type);
 378   TFile *fout = new TFile(name,"RECREATE");
 379   pAll->ProjectionX()->Write();
 380   hPi->Write();
 381   hK->Write();
 382   hPr->Write();
 383   if(isMC){
 384     TH1D *pTmp = extractFlowVZEROsingle(icentr,0,arm,1,pTh,addbin,"allMC",1,1,-1,1)->ProjectionX();
 385     pTmp->SetLineColor(6);
 386     pTmp->SetMarkerColor(6);
 387     pTmp->SetMarkerStyle(24);
 388     pTmp->Write();
 389     pTmp = extractFlowVZEROsingle(icentr,1,arm,1,pTh,addbin,"piMC",1,1,-1,1)->ProjectionX();
 390     pTmp->SetLineColor(4);
 391     pTmp->SetMarkerColor(4);
 392     pTmp->SetMarkerStyle(24);
 393     pTmp->Write();
 394     pTmp = extractFlowVZEROsingle(icentr,2,arm,1,pTh,addbin,"kaMC",1,1,-1,1)->ProjectionX();
 395     pTmp->SetLineColor(1);
 396     pTmp->SetMarkerColor(1);
 397     pTmp->SetMarkerStyle(25);
 398     pTmp->Write();
 399     pTmp = extractFlowVZEROsingle(icentr,3,arm,1,pTh,addbin,"prMC",1,1,-1,-1)->ProjectionX();
 400     pTmp->SetLineColor(2);
 401     pTmp->SetMarkerColor(2);
 402     pTmp->SetMarkerStyle(26);
 403     pTmp->Write();
 404   }
 405   extractFlowVZEROsingle(icentr,2,arm,0,pTh,addbin,"kProf")->Write();
 406   extractFlowVZEROsingle(icentr,9,arm,0,pTh,addbin,"ks",0,1,1,1)->Write();
 407   extractFlowVZEROsingle(icentr,9,arm,0,pTh,addbin,"ksMy",0,1,-1,-1)->Write();
 408   extractFlowVZEROsingle(icentr,10,arm,0,pTh,addbin,"lambda")->Write();
 409   pFromLambda->Write();
 410   piFromK->Write();
 411   pFromLambda2->Write();
 412   piFromK2->Write();
 413   piFromK3->Write();
 414   if(hFeedSyst) hFeedSyst->Write();
 415   fout->Close();
 416 }
 417
 418 TProfile *extractFlowVZEROsingle(Int_t icentr,Int_t spec,Int_t arm,Bool_t isMC,Float_t pTh,Int_t addbin,const char *nameSp,Float_t detMin,Float_t detMax,Int_t chMin,Int_t chMax){
 419   LoadLib();
 420
 421   pTh += 0.00001;
 422
 423   // NUA correction currently are missing
 424   char name[100];
 425   char stringa[200];
 426
 427   snprintf(name,100,"AnalysisResults.root");
 428   if(!fo) fo = new TFile(name);
 429   snprintf(name,100,"contVZEROv%i",arm);
 430   TList *cont = (TList *) fo->Get(name);
 431
 432   cont->ls();
 433
 434   Float_t xMin[5] = {icentr/*centrality bin*/,chMin/*charge*/,pTh/*prob*/,-TMath::Pi()/arm/*Psi*/,detMin/*PID mask*/};
 435   Float_t xMax[5] = {icentr+addbin,chMax,1.0,TMath::Pi()/arm,detMax};
 436
 437   cont->ls();
 438
 439   TProfile *p1 = cont->At(2);
 440   TProfile *p2 = cont->At(3);
 441   TProfile *p3 = cont->At(4);
 442
 443   TH2F *hPsi2DA = cont->At(5);
 444   TH2F *hPsi2DC = cont->At(6);
 445
 446   TH1D *hPsiA = hPsi2DA->ProjectionY("PsiA",icentr+1,icentr+addbin+1);
 447   TH1D *hPsiC = hPsi2DC->ProjectionY("PsiC",icentr+1,icentr+addbin+1);
 448
 449   if(!fPsi) fPsi = new TF1("fPsi","pol0",-TMath::Pi()/arm,TMath::Pi()/arm);
 450   hPsiA->Fit(fPsi,"0");
 451   Float_t offsetA = fPsi->GetParameter(0);
 452   hPsiC->Fit(fPsi,"0");
 453   Float_t offsetC = fPsi->GetParameter(0);
 454
 455   Int_t nbinPsi = hPsiA->GetNbinsX();
 456
 457   Float_t *NUAcorrA = new Float_t[nbinPsi];
 458   Float_t *NUAcorrC = new Float_t[nbinPsi];
 459
 460   for(Int_t i=0;i < nbinPsi;i++){
 461     NUAcorrA[i] = offsetA/(hPsiA->GetBinContent(i+1));
 462     NUAcorrC[i] = offsetC/(hPsiC->GetBinContent(i+1));
 463   }
 464
 465   Float_t res1=0,res2=0,res3=0;
 466   Float_t eres1=0,eres2=0,eres3=0;
 467
 468   for(Int_t i = icentr; i <= icentr+addbin;i++){
 469     if(p1->GetBinError(i+1)){
 470       eres1 += 1./p1->GetBinError(i+1)/p1->GetBinError(i+1);
 471       res1 += p1->GetBinContent(i+1)/p1->GetBinError(i+1)/p1->GetBinError(i+1);
 472     }
 473     if(p2->GetBinError(i+1)){
 474       eres2 += 1./p2->GetBinError(i+1)/p2->GetBinError(i+1);
 475       res2 += p2->GetBinContent(i+1)/p2->GetBinError(i+1)/p2->GetBinError(i+1);
 476     }
 477     if(p3->GetBinError(i+1)){
 478       eres3 += 1./p3->GetBinError(i+1)/p3->GetBinError(i+1);
 479       res3 += p3->GetBinContent(i+1)/p3->GetBinError(i+1)/p3->GetBinError(i+1);
 480     }
 481   }
 482
 483   res1 /= eres1;
 484   res2 /= eres2;
 485   res3 /= eres3;
 486
 487   eres1 = sqrt(1./eres1);
 488   eres2 = sqrt(1./eres2);
 489   eres3 = sqrt(1./eres3);
 490
 491   AliFlowVZEROResults *a = (AliFlowVZEROResults *) cont->At(0);
 492   AliFlowVZEROResults *b = (AliFlowVZEROResults *) cont->At(1);
 493   TProfile *pp,*pp2;
 494   if(kNUAcorr){ // with NUA corrections
 495     pp = a->GetV2reweight(spec,xMin,xMax,3,NUAcorrA);
 496     pp2 = b->GetV2reweight(spec,xMin,xMax,3,NUAcorrC);
 497   }
 498   else{
 499     pp = a->GetV2(spec,xMin,xMax);
 500     pp2 = b->GetV2(spec,xMin,xMax);
 501   }
 502
 503   Float_t scaling = sqrt(res1*res3/res2);
 504   if(kVZEROrescorr){
 505     pp->Scale(1./scaling);
 506   }
 507
 508   Float_t err1_2 = eres1*eres1/res1/res1/4 +
 509     eres2*eres2/res2/res2/4 +
 510     eres3*eres3/res3/res3/4;
 511   Float_t err2_2 = err1_2;
 512   err1_2 /= scaling*scaling;
 513   printf("resolution V0A = %f +/- %f\n",scaling,err1_2);
 514   scaling = sqrt(res2*res3/res1);
 515   err2_2 /= scaling*scaling;
 516   if(kVZEROrescorr){
 517     pp2->Scale(1./scaling);
 518   }
 519   printf("resolution V0C =%f +/- %f\n",scaling,err2_2);
 520
 521   pp->SetName("V0A");
 522   pp2->SetName("V0C");
 523
 524   if(! kCleanMemory){
 525     new TCanvas();
 526     pp->Draw();
 527     pp2->Draw("SAME");
 528   }
 529
 530   TProfile *pData = new TProfile(*pp);
 531   pData->Add(pp2);
 532   snprintf(stringa,100,"%sData",nameSp);
 533   pData->SetName(stringa);
 534
 535   TProfile *pMc = NULL;
 536
 537   TProfile *ppMC;
 538   TProfile *ppMC2;
 539
 540   if(isMC){
 541     snprintf(name,100,"contVZEROmc");
 542     cont = (TList *) fo->Get(name);
 543     cont->ls();
 544     if(arm == 2){
 545       AliFlowVZEROResults *c = (AliFlowVZEROResults *) cont->At(0);
 546       if(! kCleanMemory) c->GetV2(spec,xMin,xMax)->Draw("SAME");
 547     }
 548     AliFlowVZEROResults *cA;
 549     if(fo->Get("contVZEROv2")) cA = (AliFlowVZEROResults *) cont->At(1+2*(arm==3));
 550     else cA = (AliFlowVZEROResults *) cont->At(0);
 551     AliFlowVZEROResults *cC;
 552     if(fo->Get("contVZEROv2")) cC = (AliFlowVZEROResults *) cont->At(2+2*(arm==3));
 553     else cC = (AliFlowVZEROResults *) cont->At(1);
 554
 555     TProfile *p1mc = cont->At(5+3*(arm==3));
 556     TProfile *p2mc = cont->At(6+3*(arm==3));
 557     TProfile *p3mc = cont->At(7+3*(arm==3));
 558
 559     Float_t resMC1=0,resMC2=0,resMC3=0;
 560     Float_t eresMC1=0,eresMC2=0,eresMC3=0;
 561
 562     for(Int_t i = icentr; i <= icentr+addbin;i++){
 563       if(p1mc->GetBinError(i+1)){
 564         eresMC1 += 1./p1mc->GetBinError(i+1)/p1mc->GetBinError(i+1);
 565         resMC1 += p1mc->GetBinContent(i+1)/p1mc->GetBinError(i+1)/p1mc->GetBinError(i+1);
 566       }
 567       if(p2mc->GetBinError(i+1)){
 568         eresMC2 += 1./p2mc->GetBinError(i+1)/p2mc->GetBinError(i+1);
 569         resMC2 += p2mc->GetBinContent(i+1)/p2mc->GetBinError(i+1)/p2mc->GetBinError(i+1);
 570       }
 571       if(p3mc->GetBinError(i+1)){
 572         eresMC3 += 1./p3mc->GetBinError(i+1)/p3mc->GetBinError(i+1);
 573         resMC3 += p3mc->GetBinContent(i+1)/p3mc->GetBinError(i+1)/p3mc->GetBinError(i+1);
 574       }
 575     }
 576
 577     resMC1 /= eresMC1;
 578     resMC2 /= eresMC2;
 579     resMC3 /= eresMC3;
 580
 581     eresMC1 = sqrt(1./eresMC1);
 582     eresMC2 = sqrt(1./eresMC2);
 583     eresMC3 = sqrt(1./eresMC3);
 584
 585     ppMC = cA->GetV2(spec,xMin,xMax);
 586     ppMC2 = cC->GetV2(spec,xMin,xMax);
 587
 588     scaling = sqrt(resMC1*resMC3/resMC2);
 589     ppMC->Scale(1./scaling);
 590
 591     err1_2 = eresMC1*eresMC1/resMC1/resMC1/4 +
 592       eresMC2*eresMC2/resMC2/resMC2/4 +
 593       eresMC3*eresMC3/resMC3/resMC3/4;
 594     err2_2 = err1_2;
 595     err1_2 /= scaling*scaling;
 596     printf("resolution V0A (MC) = %f +/- %f\n",scaling,err1_2);
 597     scaling = sqrt(resMC2*resMC3/resMC1);
 598     err2_2 /= scaling*scaling;
 599     ppMC2->Scale(1./scaling);
 600     printf("resolution V0C (MC) =%f +/- %f\n",scaling,err2_2);
 601
 602     ppMC->SetName("V0Amc");
 603     ppMC2->SetName("V0Cmc");
 604
 605     if(! kCleanMemory){
 606       ppMC->Draw("SAME");
 607       ppMC2->Draw("SAME");
 608     }
 609
 610     pMc = new TProfile(*ppMC);
 611     pMc->Add(ppMC2);
 612     snprintf(stringa,100,"%sMC",nameSp);
 613     pMc->SetName(stringa);
 614     pMc->SetLineColor(2);
 615   }
 616
 617   if(! kCleanMemory){
 618     new TCanvas();
 619     pData->Draw();
 620   }
 621   if(pMc && !kCleanMemory){
 622     pMc->Draw("SAME");
 623     TH1D *hData = pData->ProjectionX();
 624     TH1D *hMc = pMc->ProjectionX();
 625     hData->Divide(hMc);
 626     new TCanvas();
 627     hData->Draw();
 628   }
 629
 630   delete[] NUAcorrA;
 631   delete[] NUAcorrC;
 632
 633   if(kCleanMemory){
 634     if(pp) delete pp;
 635     if(pp2) delete pp2;
 636     if(ppMC) delete ppMC;
 637     if(ppMC2) delete ppMC2;
 638   }
 639
 640   delete cont;
 641
 642   if(isMC) return pMc;
 643   return pData;
 644 }
 645 void compareTPCTOF(Int_t icentr,Int_t spec,Int_t arm,Float_t pTh,Int_t addbin){
 646   LoadLib();
 647   Float_t ptMaxFit[8] = {20,0.7,0.55,1.2,2,2,2,2};
 648
 649   TProfile *pTPC = extractFlowVZEROsingle(icentr,spec,arm,0,pTh,addbin,"TPC",0,0,-1,-1+2*(spec!=3)); // TPC && !TOF (TPC PID)
 650   TProfile *pTPC2 = extractFlowVZEROsingle(icentr,spec,arm,0,pTh,addbin,"TPCextra",2,2,-1,-1+2*(spec!=3)); //TPC && TOF (TPC PID)
 651
 652   TProfile *pTOF = extractFlowVZEROsingle(icentr,spec,arm,0,pTh,addbin,"TOF",1,1,-1,-1+2*(spec!=3)); //TPC && TOF (TPC&TOF PID)
 653
 654   if(spec > 0) pTPC->Add(pTPC2);
 655   else pTPC->Add(pTOF);
 656
 657   char name[100];
 658   snprintf(name,100,"NUO_%i",spec);
 659
 660   hNUO[spec] = pTPC->ProjectionX(name);
 661   hNUO[spec]->Add(pTOF,-1);
 662
 663   if(spec && hNUO[0]) hNUO[spec]->Add(hNUO[0],-1);
 664
 665   if(! kCleanMemory){
 666     new TCanvas();
 667     pTPC->Draw();
 668     pTOF->Draw("SAME");
 669
 670     new TCanvas();
 671     pTPC->SetLineColor(1);
 672     pTOF->SetLineColor(2);
 673   }
 674
 675   hNUO[spec]->Draw();
 676
 677   if(kCleanMemory){
 678     if(pTPC) delete pTPC;
 679     if(pTPC2) delete pTPC2;
 680     if(pTOF) delete pTOF;
 681   }
 682 }
 683
 684 void plotQApid(Int_t ic,Float_t pt,Int_t addbin){
 685   char name[100];
 686   char stringa[200];
 687   LoadLib();
 688
 689   snprintf(name,100,"AnalysisResults.root");
 690   if(!fo) fo = new TFile(name);
 691   snprintf(name,100,"contVZEROv2");
 692   TList *cont = (TList *) fo->Get(name);
 693   AliFlowVZEROResults *pidqa = cont->FindObject("qaPID");
 694   Float_t xval[2] = {2.,pt+0.00001};
 695   Float_t xval2[2] = {2.+addbin,pt+0.00001};
 696
 697   TProfile *proTPCpi = pidqa->GetV2(0,xval,xval2);
 698   TProfile *proTOFpi = pidqa->GetV2(1,xval,xval2);
 699   TProfile *proTPCka = pidqa->GetV2(2,xval,xval2);
 700   TProfile *proTOFka = pidqa->GetV2(3,xval,xval2);
 701   TProfile *proTPCpr = pidqa->GetV2(4,xval,xval2);
 702   TProfile *proTOFpr = pidqa->GetV2(5,xval,xval2);
 703
 704   proTPCpi->SetName("hPiTPC");
 705   proTOFpi->SetName("hPiTOF");
 706   proTPCka->SetName("hKaTPC");
 707   proTOFka->SetName("hKaTOF");
 708   proTPCpr->SetName("hPrTPC");
 709   proTOFpr->SetName("hPrTOF");
 710
 711   proTPCpi->GetXaxis()->SetTitle("dE/dx - dE/dx_{calc}^{#pi} (a.u)");
 712   proTOFpi->GetXaxis()->SetTitle("t_{tof} - t_{calc}^{#pi} (ps)");
 713   proTPCka->GetXaxis()->SetTitle("dE/dx - dE/dx_{calc}^{K} (a.u)");
 714   proTOFka->GetXaxis()->SetTitle("t_{tof} - t_{calc}^{K} (ps)");
 715   proTPCpr->GetXaxis()->SetTitle("dE/dx - dE/dx_{calc}^{p} (a.u)");
 716   proTOFpr->GetXaxis()->SetTitle("t_{tof} - t_{calc}^{p} (ps)");
 717
 718   TCanvas *c = new TCanvas("cVproj","cVproj");
 719   c->Divide(2,3);
 720   c->cd(1);
 721   proTPCpi->Draw();
 722   c->cd(2);
 723   proTOFpi->Draw();
 724   c->cd(3);
 725   proTPCka->Draw();
 726   c->cd(4);
 727   proTOFka->Draw();
 728   c->cd(5);
 729   proTPCpr->Draw();
 730   c->cd(6);
 731   proTOFpr->Draw();
 732 }