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