[u/mrichter/AliRoot.git] / PWG2 / EBYE / macros / readBalanceFunction.C

const TString gBFAnalysisType[7] = {"y","eta","qlong","qout","qside","qinv","phi"};

const Double_t cent[9]  = {382.8,329.7,260.5,186.4,128.9,85.,52.8,30.,15.8};   // hard coded at the moment for centrality percentiles 
const Double_t centE[9] = {3.1,4.6,4.4,3.9,3.3,2.6,2.0,1.3,0.6};               // (0-5,5-10,10-20,20-30,...,70-80)

void readBalanceFunction(Bool_t bHistos = kTRUE, TString inFile = "AnalysisResults.root",Int_t fStartBinBFWidth = 1, Int_t fRebin = 1) {
  // Macro to read the output of the BF analysis:  MW: CHANGE THIS!!!!
  //i) Prints and draws the final BF output
  //ii) Plots the QA part of the analysis
  //iii) store BF in output file
  //Author: Panos.Christakoglou@cern.ch, m.weber@cern.ch
  //Loading the needed libraries
  gSystem->Load("libProofPlayer.so");
  gSystem->Load("libANALYSIS.so");
  gSystem->Load("libANALYSISalice.so");
  gSystem->Load("libPWG2ebye.so");

  //Draw BF              
  drawBF(bHistos,inFile, fStartBinBFWidth, fRebin);

  //Merge the output
  //mergeOutput("/alice/cern.ch/user/p/pchrist/Balance/pp/7TeV/LHC10b/output/");
}

//___________________________________________________________//
void drawBF(Bool_t bHistos = kTRUE, TString inFile = "AnalysisResults.root", Int_t fStartBinBFWidth = 1, Int_t fRebin = 1) {
  //Function to draw the BF objects and write them into the output file

  Int_t maximumCanvases = 10;
  Int_t iCanvas         = 0;
  TCanvas *cQA[10];
  TCanvas *cQAV0M = new TCanvas("cQAV0M","V0M multiplicities");
  cQAV0M->Divide(4,3);
  TCanvas *cQARef = new TCanvas("cQARef","reference track multiplicities");
  cQARef->Divide(4,3);
  TCanvas *cBF[10];
  TCanvas *cBFS[10];

  // get the file
  TFile *f = TFile::Open(inFile.Data());
  if(!f) {
    Printf("File not found!!!");
    break;
  }

  // get the BF output directory
  TDirectoryFile *dir = dynamic_cast<TDirectoryFile *>(f->Get("PWG2EbyE.outputBalanceFunctionAnalysis"));
  if(!dir) {
    Printf("Output directory not found!!!");
    break;
  }

  // loop over all lists and plot the BF and QA
  TList *list = NULL;
  TString listName;
  TIter nextkey( dir->GetListOfKeys() );
  TKey *key;

  AliBalance *bf[10][7];
  AliBalance *bfs[10][7];
  TH1D *gbf[10][7];
  TH1D *gbfs[10][7];

  for(Int_t i = 0; i < 10; i++){
    for(Int_t j = 0; j < 7; j++){
      gbf[i][j] = NULL;
      gbfs[i][j] = NULL;
    }
  }

  TH1D *fHistP[7]; //N+
  TH1D *fHistN[7]; //N-
  TH1D *fHistPN[7]; //N+-
  TH1D *fHistNP[7]; //N-+
  TH1D *fHistPP[7]; //N++
  TH1D *fHistNN[7]; //N--

  TH1D *fHistPS[7]; //N+
  TH1D *fHistNS[7]; //N-
  TH1D *fHistPNS[7]; //N+-
  TH1D *fHistNPS[7]; //N-+
  TH1D *fHistPPS[7]; //N++
  TH1D *fHistNNS[7]; //N--

  Double_t WM[10];     // weighted mean for eta (recalculated from fStartBin)
  Double_t WME[10];    // error
  Double_t WMS[10];     // weighted mean for eta (recalculated from fStartBin) (shuffled)
  Double_t WMSE[10];    // error (shuffled)

  while ( (key = (TKey*)nextkey())) {

    list = (TList*)key->ReadObj();
    listName = TString(list->GetName());

    cout<<"Processing list "<<listName<<endl;

    // ----------------------------------------------------
    // plot QA histograms
    if(listName.Contains("QA")){     

      iCanvas ++;
      if(iCanvas == 10) {cout<<"TOO MANY LISTS --> increase MAXIMUM"<<endl; return;}
      cQA[iCanvas] = new TCanvas(listName,listName);
      cQA[iCanvas]->Divide(4,3);

     cQA[iCanvas]->cd(1);
      TH1F* histEventStats = (TH1F*)list->FindObject("fHistEventStats");
      if(histEventStats){
	histEventStats->SetFillColor(9);
	histEventStats->Draw();
      }

      cQA[iCanvas]->cd(2);
      TH1F* histTriggerStats = (TH1F*)list->FindObject("fHistTriggerStats");
      if(histTriggerStats){
	histTriggerStats->SetFillColor(9);
	histTriggerStats->Draw();
      }

     cQA[iCanvas]->cd(3);
      TH1F* histTrackStats = (TH1F*)list->FindObject("fHistTrackStats");
      if(histTrackStats){
	histTrackStats->SetFillColor(9);
	histTrackStats->Draw();
      }

      cQA[iCanvas]->cd(4);
      TH1F* histCentStats = (TH1F*)list->FindObject("fHistCentStats");
      if(histCentStats){
	histCentStats->SetFillColor(9);
	histCentStats->Draw("colz");
      }


      cQA[iCanvas]->cd(5);
      TH1F* histVx = (TH1F*)list->FindObject("fHistVx");
      if(histVx){
	histVx->SetFillColor(9);
	histVx->Draw();
      }
      cQA[iCanvas]->cd(6);
      TH1F* histVy = (TH1F*)list->FindObject("fHistVy");
      if(histVy){
	histVy->SetFillColor(9);
	histVy->Draw();
      }
      
      cQA[iCanvas]->cd(7);
      TH1F* histVz = (TH1F*)list->FindObject("fHistVz");
      if(histVz){
	histVz->SetFillColor(9);
	histVz->Draw();
      }

      cQA[iCanvas]->cd(8);
      cQA[iCanvas]->cd(8)->SetLogz();
      TH2F* histDCA = (TH2F*)list->FindObject("fHistDCA");  
      if(histDCA) histDCA->Draw("colz");
      

      cQA[iCanvas]->cd(9);
      cQA[iCanvas]->cd(9)->SetLogz();
      TH2F* histClus = (TH2F*)list->FindObject("fHistClus");
      if(histClus) histClus->Draw("colz");
      
      cQA[iCanvas]->cd(10);
      TH1F* histPt = (TH1F*)list->FindObject("fHistPt");
      if(histPt){
	histPt->SetFillColor(9);
	histPt->Draw();
      }
      
      cQA[iCanvas]->cd(11);
      TH1F* histEta = (TH1F*)list->FindObject("fHistEta");
      if(histEta){
	histEta->SetFillColor(9);
	histEta->Draw();
      }
      
      cQA[iCanvas]->cd(12);
      TH1F* histPhi = (TH1F*)list->FindObject("fHistPhi");
      if(histPhi){
	histPhi->SetFillColor(9);
	histPhi->Draw();
      }

      // centrality estimator QA
      cQAV0M->cd(iCanvas);
      cQAV0M->cd(iCanvas)->SetLogz();
      TH1F* histV0M = (TH1F*)list->FindObject("fHistV0M");
      if(histV0M){
	histV0M->Draw("colz");
      }

      cQARef->cd(iCanvas);
      cQARef->cd(iCanvas)->SetLogz();
      TH1F* histRef = (TH1F*)list->FindObject("fHistRefTracks");
      if(histRef){
	histRef->Draw("colz");
      }
    }
    // ----------------------------------------------------

    // ----------------------------------------------------
    // calculate and plot BF 
    if(listName.Contains("BF_")){
      cBF[iCanvas] = new TCanvas(listName,listName);
      cBF[iCanvas]->Divide(3,3);

      for(Int_t a = 0; a < 7; a++){

	cout<<"ANALYSE "<<gBFAnalysisType[a]<<endl;

	// create the BF object
	bf[iCanvas][a]  = new AliBalance();

	fHistP[a] = (TH1D*)list->FindObject(Form("fHistP%s",gBFAnalysisType[a].Data()));
	fHistN[a] = (TH1D*)list->FindObject(Form("fHistN%s",gBFAnalysisType[a].Data()));
	fHistPP[a] = (TH1D*)list->FindObject(Form("fHistPP%s",gBFAnalysisType[a].Data()));
	fHistPN[a] = (TH1D*)list->FindObject(Form("fHistPN%s",gBFAnalysisType[a].Data()));
	fHistNP[a] = (TH1D*)list->FindObject(Form("fHistNP%s",gBFAnalysisType[a].Data()));
	fHistNN[a] = (TH1D*)list->FindObject(Form("fHistNN%s",gBFAnalysisType[a].Data()));

	// rebin histograms (be careful with divider!)
	fHistP[a]->Rebin(fRebin);
	fHistN[a]->Rebin(fRebin);
	fHistPP[a]->Rebin(fRebin);
	fHistPN[a]->Rebin(fRebin);
	fHistNP[a]->Rebin(fRebin);
	fHistNN[a]->Rebin(fRebin);

	// set histograms in AliBalance object
	bf[iCanvas][a]->SetHistNp(a, fHistP[a]);
	bf[iCanvas][a]->SetHistNn(a, fHistN[a]);
	bf[iCanvas][a]->SetHistNpp(a, fHistPP[a]);
	bf[iCanvas][a]->SetHistNpn(a, fHistPN[a]);
	bf[iCanvas][a]->SetHistNnp(a, fHistNP[a]);
	bf[iCanvas][a]->SetHistNnn(a, fHistNN[a]);

	gbf[iCanvas][a] = bf[iCanvas][a]->GetBalanceFunctionHistogram(a);
	gbf[iCanvas][a]->SetName(Form("%s_BF_%s",listName.Data(),gBFAnalysisType[a].Data()));

	cBF[iCanvas]->cd(a+1);
	gbf[iCanvas][a]->SetMarkerStyle(20);

	if(!bHistos){
	  gbf[iCanvas][a]->DrawCopy("AP");
	  if(a==1) GetWeightedMean(gbf[iCanvas][a],fStartBinBFWidth,WM[iCanvas-1],WME[iCanvas-1]); // for eta recalculate width (from 0.1 only!)
	}
	else{
	  fHistPN[a]->SetLineColor(2);
	  fHistPN[a]->Draw();
	  fHistPP[a]->SetLineColor(1);
	  fHistPP[a]->Draw("same");
	  fHistNP[a]->SetLineColor(4);
	  fHistNP[a]->Draw("same");
	  fHistNN[a]->SetLineColor(8);
	  fHistNN[a]->Draw("same");
	}
      }
    }
    // ----------------------------------------------------

    // ----------------------------------------------------
    // calculate and plot BF (shuffled)
    if(listName.Contains("BFShuffled")){
      cBFS[iCanvas] = new TCanvas(listName,listName);
      cBFS[iCanvas]->Divide(3,3);

      for(Int_t a = 0; a < 7; a++){

	// create the BF object
	bfs[iCanvas][a]  = new AliBalance();

	fHistPS[a] = (TH1D*)list->FindObject(Form("fHistP%s_shuffle",gBFAnalysisType[a].Data()));
	fHistNS[a] = (TH1D*)list->FindObject(Form("fHistP%s_shuffle",gBFAnalysisType[a].Data()));
	fHistPPS[a] = (TH1D*)list->FindObject(Form("fHistPP%s_shuffle",gBFAnalysisType[a].Data()));
	fHistPNS[a] = (TH1D*)list->FindObject(Form("fHistPN%s_shuffle",gBFAnalysisType[a].Data()));
	fHistNPS[a] = (TH1D*)list->FindObject(Form("fHistNP%s_shuffle",gBFAnalysisType[a].Data()));
	fHistNNS[a] = (TH1D*)list->FindObject(Form("fHistNN%s_shuffle",gBFAnalysisType[a].Data()));

	// rebin histograms (be careful with divider!)
	fHistPS[a]->Rebin(fRebin);
	fHistNS[a]->Rebin(fRebin);
	fHistPPS[a]->Rebin(fRebin);
	fHistPNS[a]->Rebin(fRebin);
	fHistNPS[a]->Rebin(fRebin);
	fHistNNS[a]->Rebin(fRebin);

	// set histograms in AliBalance object
	bfs[iCanvas][a]->SetHistNp(a, fHistPS[a]);
	bfs[iCanvas][a]->SetHistNn(a, fHistNS[a]);
	bfs[iCanvas][a]->SetHistNpp(a, fHistPPS[a]);
	bfs[iCanvas][a]->SetHistNpn(a, fHistPNS[a]);
	bfs[iCanvas][a]->SetHistNnp(a, fHistNPS[a]);
	bfs[iCanvas][a]->SetHistNnn(a, fHistNNS[a]);

	gbfs[iCanvas][a] = bfs[iCanvas][a]->GetBalanceFunctionHistogram(a);
	gbfs[iCanvas][a]->SetName(Form("%s_BF_%s",listName.Data(),gBFAnalysisType[a].Data()));

	cBFS[iCanvas]->cd(a+1);
	gbfs[iCanvas][a]->SetMarkerStyle(20);
	if(!bHistos){
	  gbfs[iCanvas][a]->DrawCopy("AP");
	  if(a==1) GetWeightedMean(gbfs[iCanvas][a],fStartBinBFWidth,WMS[iCanvas-1],WMSE[iCanvas-1]); // for eta recalculate width (from 0.1 only!)
	}
	else{
	  fHistPNS[a]->SetLineColor(2);
	  fHistPNS[a]->Draw();
	  fHistPPS[a]->SetLineColor(1);
	  fHistPPS[a]->Draw("same");
	  fHistNPS[a]->SetLineColor(4);
	  fHistNPS[a]->Draw("same");
	  fHistNNS[a]->SetLineColor(8);
	  fHistNNS[a]->Draw("same");
	}
      }
    }
    // ----------------------------------------------------
  }

  // for BF calculation create also graphs with weighted mean for eta
  TGraphErrors *gWM  = NULL;
  TGraphErrors *gWMS = NULL;
  if(!bHistos && gbf[1][1]){
    gWM = new TGraphErrors(iCanvas,cent,WM,centE,WME);
    gWM->SetName("gCentrality");    
  }
  if(!bHistos && gbfs[1][1]){
    gWMS = new TGraphErrors(iCanvas,cent,WMS,centE,WMSE);  
    gWMS->SetName("gCentralityS");
  }

  TFile *fOut = TFile::Open(Form("Histograms_WMstart%d_rebin%d_%s", fStartBinBFWidth, fRebin,inFile.Data()),"RECREATE");
  fOut->cd();
  for(Int_t i = 0; i <= iCanvas; i++){
    for(Int_t a = 0; a < 7; a++){
      if(gbf[i][a]){
	gbf[i][a]->Write();
	gbf[i][a]->Delete();
      }
      if(gbfs[i][a]){
	gbfs[i][a]->Write();
	gbfs[i][a]->Delete();	
      }
    }
  }

  if(gWM) gWM->Write();
  if(gWMS) gWMS->Write();

  fOut->Close();
  f->Close();
  gROOT->Reset();
}

//____________________________________________________________________//
void GetWeightedMean(TH1D *gHistBalance, Int_t fStartBin = 1, Double_t &WM, Double_t &WME) {

  //Prints the calculated width of the BF and its error
  Double_t gSumXi = 0.0, gSumBi = 0.0, gSumBiXi = 0.0;
  Double_t gSumBiXi2 = 0.0, gSumBi2Xi2 = 0.0;
  Double_t gSumDeltaBi2 = 0.0, gSumXi2DeltaBi2 = 0.0;
  Double_t deltaBalP2 = 0.0, integral = 0.0;
  Double_t deltaErrorNew = 0.0;

  //Retrieve this variables from Histogram
  Int_t fNumberOfBins = gHistBalance->GetNbinsX();
  Double_t fP2Step    = gHistBalance->GetBinWidth(1); // assume equal binning!
  
  cout<<"=================================================="<<endl;
  cout<<"RECALCULATION OF BF WIDTH (StartBin = "<<fStartBin<<")"<<endl;
  cout<<"HISTOGRAM has "<<fNumberOfBins<<" bins with bin size of "<<fP2Step<<endl;
  for(Int_t i = fStartBin; i <= fNumberOfBins; i++) { 
    cout<<"B: "<<gHistBalance->GetBinContent(i)<<"\t Error: "<<gHistBalance->GetBinError(i)<<"\t bin: "<<gHistBalance->GetBinCenter(i)<<endl;
  } 
  cout<<"=================================================="<<endl;
  for(Int_t i = fStartBin; i <= fNumberOfBins; i++) {
    gSumXi += gHistBalance->GetBinCenter(i);
    gSumBi += gHistBalance->GetBinContent(i);
    gSumBiXi += gHistBalance->GetBinContent(i)*gHistBalance->GetBinCenter(i);
    gSumBiXi2 += gHistBalance->GetBinContent(i)*TMath::Power(gHistBalance->GetBinCenter(i),2);
    gSumBi2Xi2 += TMath::Power(gHistBalance->GetBinContent(i),2)*TMath::Power(gHistBalance->GetBinCenter(i),2);
    gSumDeltaBi2 +=  TMath::Power(gHistBalance->GetBinError(i),2);
    gSumXi2DeltaBi2 += TMath::Power(gHistBalance->GetBinCenter(i),2) * TMath::Power(gHistBalance->GetBinError(i),2);
    
    deltaBalP2 += fP2Step*TMath::Power(gHistBalance->GetBinError(i),2);
    integral += fP2Step*gHistBalance->GetBinContent(i);
  }
  for(Int_t i = 1; i < fNumberOfBins; i++)
    deltaErrorNew += gHistBalance->GetBinError(i)*(gHistBalance->GetBinCenter(i)*gSumBi - gSumBiXi)/TMath::Power(gSumBi,2);
  
  Double_t integralError = TMath::Sqrt(deltaBalP2);
  
  Double_t delta = gSumBiXi / gSumBi;
  Double_t deltaError = (gSumBiXi / gSumBi) * TMath::Sqrt(TMath::Power((TMath::Sqrt(gSumXi2DeltaBi2)/gSumBiXi),2) + TMath::Power((gSumDeltaBi2/gSumBi),2) );
  
  cout<<"Width: "<<delta<<"\t Error: "<<deltaError<<endl;
  cout<<"New error: "<<deltaErrorNew<<endl;
  cout<<"Integral: "<<integral<<"\t Error: "<<integralError<<endl;
  cout<<"=================================================="<<endl;

  WM  = delta;
  WME = deltaError;
}

//___________________________________________________________//
// NOT USED any more
//___________________________________________________________//
void mergeOutput(const char* outputDir) {
  //Function to merge the output of the sub-jobs
  //Create a BF object
  AliBalance *bf = new AliBalance();

  //connect to AliEn's API services
  TGrid::Connect("alien://"); 

  //Getting the output dir from the env. variable 
  //(JDL field: JDLVariables={"OutputDir"};)
  TGridResult* result = gGrid->Query(outputDir,"*/root_archive.zip","","-l 1000");
  
  Int_t nEntries = result->GetEntries();

  TString alienUrl;
  TDirectoryFile *dirSubJob;

  TString gCutName[4] = {"Total","Offline trigger",
                         "Vertex","Analyzed"};
  TH1F *fHistEventStats = new TH1F("fHistEventStats",
				   "Event statistics;;N_{events}",
				   4,0.5,4.5);
  for(Int_t i = 1; i <= 4; i++)
    fHistEventStats->GetXaxis()->SetBinLabel(i,gCutName[i-1].Data());

  AliESDtrackCuts *bfTrackCuts = new AliESDtrackCuts("bfTrackCuts");
  for(Int_t i = 0; i < nEntries; i++) {
    alienUrl = result->GetKey(i,"turl");
    alienUrl += "#AnalysisResults.root";
    Printf("Opening file: %s",alienUrl.Data());
    TFile *file = TFile::Open(alienUrl.Data());
    dirSubJob = dynamic_cast<TDirectoryFile *>(file->Get("PWG2EbyE.outputBalanceFunctionAnalysis.root"));

    //merge BF
    AliBalance *bfSubJob = dynamic_cast<AliBalance *>(dirSubJob->Get("AliBalance"));
    //bfSubJob->PrintResults();
    bf->Merge(bfSubJob);
    //delete bfSubJob;

    //merge event stats
    TList *listSubJob = dynamic_cast<TList *>(dirSubJob->Get("listQA"));
    fHistEventStats->Add(dynamic_cast<TH1F *>(listSubJob->At(0)));

    bfTrackCuts = dynamic_cast<AliESDtrackCuts *>(listSubJob->At(1));
    delete listSubJob;
  }

  //Create the output file
  TString outputFile = "AnalysisResults.Merged.root";
  TFile *foutput = TFile::Open(outputFile.Data(),"recreate");
  TDirectoryFile *dirOutput = new TDirectoryFile();
  dirOutput->SetName("PWG2EbyE.outputBalanceFunctionAnalysis.root");
  //dirOutput->cd();
  dirOutput->Add(bf);
  TList *list = new TList();
  list->SetName("listQA");
  list->Add(fHistEventStats);
  list->Add(bfTrackCuts);
  dirOutput->Add(list);
  dirOutput->Write();
  bf->Write();
  list->Write();
  foutput->Close();

    //cout<<alienUrl.Data()<<endl;
}
Commit	Line	Data
6432ac6a	1	const TString gBFAnalysisType[7] = {"y","eta","qlong","qout","qside","qinv","phi"};
6432ac6a	2
a0ce3c86	3	const Double_t cent[9] = {382.8,329.7,260.5,186.4,128.9,85.,52.8,30.,15.8}; // hard coded at the moment for centrality percentiles
	4	const Double_t centE[9] = {3.1,4.6,4.4,3.9,3.3,2.6,2.0,1.3,0.6}; // (0-5,5-10,10-20,20-30,...,70-80)
	5
	6	void readBalanceFunction(Bool_t bHistos = kTRUE, TString inFile = "AnalysisResults.root",Int_t fStartBinBFWidth = 1, Int_t fRebin = 1) {
6432ac6a	7	// Macro to read the output of the BF analysis: MW: CHANGE THIS!!!!
	8	//i) Prints and draws the final BF output
	9	//ii) Plots the QA part of the analysis
	10	//iii) store BF in output file
	11	//Author: Panos.Christakoglou@cern.ch, m.weber@cern.ch
a114d234	12	//Loading the needed libraries
	13	gSystem->Load("libProofPlayer.so");
	14	gSystem->Load("libANALYSIS.so");
	15	gSystem->Load("libANALYSISalice.so");
	16	gSystem->Load("libPWG2ebye.so");
	17
a114d234	18	//Draw BF
a0ce3c86	19	drawBF(bHistos,inFile, fStartBinBFWidth, fRebin);
61483909	20
61483909	21	//Merge the output
a114d234	22	//mergeOutput("/alice/cern.ch/user/p/pchrist/Balance/pp/7TeV/LHC10b/output/");
	23	}
	24
	25	//___________________________________________________________//
a0ce3c86	26	void drawBF(Bool_t bHistos = kTRUE, TString inFile = "AnalysisResults.root", Int_t fStartBinBFWidth = 1, Int_t fRebin = 1) {
6432ac6a	27	//Function to draw the BF objects and write them into the output file
	28
	29	Int_t maximumCanvases = 10;
	30	Int_t iCanvas = 0;
	31	TCanvas *cQA[10];
a0ce3c86	32	TCanvas *cQAV0M = new TCanvas("cQAV0M","V0M multiplicities");
	33	cQAV0M->Divide(4,3);
	34	TCanvas *cQARef = new TCanvas("cQARef","reference track multiplicities");
	35	cQARef->Divide(4,3);
6432ac6a	36	TCanvas *cBF[10];
	37	TCanvas *cBFS[10];
	38
	39	// get the file
	40	TFile *f = TFile::Open(inFile.Data());
a114d234	41	if(!f) {
	42	Printf("File not found!!!");
	43	break;
	44	}
	45
6432ac6a	46	// get the BF output directory
6432ac6a	47	TDirectoryFile dir = dynamic_cast<TDirectoryFile >(f->Get("PWG2EbyE.outputBalanceFunctionAnalysis"));
a114d234	48	if(!dir) {
	49	Printf("Output directory not found!!!");
	50	break;
	51	}
	52
6432ac6a	53	// loop over all lists and plot the BF and QA
	54	TList *list = NULL;
	55	TString listName;
	56	TIter nextkey( dir->GetListOfKeys() );
	57	TKey *key;
	58
	59	AliBalance *bf[10][7];
	60	AliBalance *bfs[10][7];
32a4f715	61	TH1D *gbf[10][7];
32a4f715	62	TH1D *gbfs[10][7];
6432ac6a	63
a0ce3c86	64	for(Int_t i = 0; i < 10; i++){
	65	for(Int_t j = 0; j < 7; j++){
	66	gbf[i][j] = NULL;
	67	gbfs[i][j] = NULL;
	68	}
	69	}
	70
6432ac6a	71	TH1D *fHistP[7]; //N+
	72	TH1D *fHistN[7]; //N-
	73	TH1D *fHistPN[7]; //N+-
	74	TH1D *fHistNP[7]; //N-+
	75	TH1D *fHistPP[7]; //N++
	76	TH1D *fHistNN[7]; //N--
	77
32a4f715	78	TH1D *fHistPS[7]; //N+
	79	TH1D *fHistNS[7]; //N-
	80	TH1D *fHistPNS[7]; //N+-
	81	TH1D *fHistNPS[7]; //N-+
	82	TH1D *fHistPPS[7]; //N++
	83	TH1D *fHistNNS[7]; //N--
	84
a0ce3c86	85	Double_t WM[10]; // weighted mean for eta (recalculated from fStartBin)
	86	Double_t WME[10]; // error
	87	Double_t WMS[10]; // weighted mean for eta (recalculated from fStartBin) (shuffled)
	88	Double_t WMSE[10]; // error (shuffled)
	89
6432ac6a	90	while ( (key = (TKey*)nextkey())) {
	91
	92	list = (TList*)key->ReadObj();
	93	listName = TString(list->GetName());
	94
	95	cout<<"Processing list "<<listName<<endl;
	96
	97	// ----------------------------------------------------
	98	// plot QA histograms
	99	if(listName.Contains("QA")){
	100
	101	iCanvas ++;
	102	if(iCanvas == 10) {cout<<"TOO MANY LISTS --> increase MAXIMUM"<<endl; return;}
	103	cQA[iCanvas] = new TCanvas(listName,listName);
a0ce3c86	104	cQA[iCanvas]->Divide(4,3);
	105
	106	cQA[iCanvas]->cd(1);
	107	TH1F* histEventStats = (TH1F*)list->FindObject("fHistEventStats");
	108	if(histEventStats){
	109	histEventStats->SetFillColor(9);
	110	histEventStats->Draw();
	111	}
	112
	113	cQA[iCanvas]->cd(2);
	114	TH1F* histTriggerStats = (TH1F*)list->FindObject("fHistTriggerStats");
	115	if(histTriggerStats){
	116	histTriggerStats->SetFillColor(9);
	117	histTriggerStats->Draw();
	118	}
6432ac6a	119
a0ce3c86	120	cQA[iCanvas]->cd(3);
	121	TH1F* histTrackStats = (TH1F*)list->FindObject("fHistTrackStats");
	122	if(histTrackStats){
	123	histTrackStats->SetFillColor(9);
	124	histTrackStats->Draw();
	125	}
	126
	127	cQA[iCanvas]->cd(4);
	128	TH1F* histCentStats = (TH1F*)list->FindObject("fHistCentStats");
	129	if(histCentStats){
	130	histCentStats->SetFillColor(9);
	131	histCentStats->Draw("colz");
	132	}
	133
	134
	135
	136	cQA[iCanvas]->cd(5);
6432ac6a	137	TH1F* histVx = (TH1F*)list->FindObject("fHistVx");
	138	if(histVx){
	139	histVx->SetFillColor(9);
	140	histVx->Draw();
	141	}
a0ce3c86	142	cQA[iCanvas]->cd(6);
6432ac6a	143	TH1F* histVy = (TH1F*)list->FindObject("fHistVy");
	144	if(histVy){
	145	histVy->SetFillColor(9);
	146	histVy->Draw();
	147	}
	148
a0ce3c86	149	cQA[iCanvas]->cd(7);
6432ac6a	150	TH1F* histVz = (TH1F*)list->FindObject("fHistVz");
	151	if(histVz){
	152	histVz->SetFillColor(9);
	153	histVz->Draw();
	154	}
a0ce3c86	155
	156	cQA[iCanvas]->cd(8);
	157	cQA[iCanvas]->cd(8)->SetLogz();
	158	TH2F* histDCA = (TH2F*)list->FindObject("fHistDCA");
	159	if(histDCA) histDCA->Draw("colz");
6432ac6a	160
a0ce3c86	161
	162	cQA[iCanvas]->cd(9);
	163	cQA[iCanvas]->cd(9)->SetLogz();
6432ac6a	164	TH2F* histClus = (TH2F*)list->FindObject("fHistClus");
	165	if(histClus) histClus->Draw("colz");
	166
a0ce3c86	167	cQA[iCanvas]->cd(10);
6432ac6a	168	TH1F* histPt = (TH1F*)list->FindObject("fHistPt");
	169	if(histPt){
	170	histPt->SetFillColor(9);
	171	histPt->Draw();
	172	}
	173
a0ce3c86	174	cQA[iCanvas]->cd(11);
6432ac6a	175	TH1F* histEta = (TH1F*)list->FindObject("fHistEta");
	176	if(histEta){
	177	histEta->SetFillColor(9);
	178	histEta->Draw();
	179	}
	180
a0ce3c86	181	cQA[iCanvas]->cd(12);
6432ac6a	182	TH1F* histPhi = (TH1F*)list->FindObject("fHistPhi");
	183	if(histPhi){
	184	histPhi->SetFillColor(9);
	185	histPhi->Draw();
	186	}
a0ce3c86	187
	188	// centrality estimator QA
	189	cQAV0M->cd(iCanvas);
	190	cQAV0M->cd(iCanvas)->SetLogz();
	191	TH1F* histV0M = (TH1F*)list->FindObject("fHistV0M");
	192	if(histV0M){
	193	histV0M->Draw("colz");
	194	}
	195
	196	cQARef->cd(iCanvas);
	197	cQARef->cd(iCanvas)->SetLogz();
	198	TH1F* histRef = (TH1F*)list->FindObject("fHistRefTracks");
	199	if(histRef){
	200	histRef->Draw("colz");
	201	}
6432ac6a	202	}
	203	// ----------------------------------------------------
	204
	205	// ----------------------------------------------------
	206	// calculate and plot BF
	207	if(listName.Contains("BF_")){
	208	cBF[iCanvas] = new TCanvas(listName,listName);
	209	cBF[iCanvas]->Divide(3,3);
	210
	211	for(Int_t a = 0; a < 7; a++){
	212
	213	cout<<"ANALYSE "<<gBFAnalysisType[a]<<endl;
	214
	215	// create the BF object
	216	bf[iCanvas][a] = new AliBalance();
	217
	218	fHistP[a] = (TH1D*)list->FindObject(Form("fHistP%s",gBFAnalysisType[a].Data()));
a5ab9ca3	219	fHistN[a] = (TH1D*)list->FindObject(Form("fHistN%s",gBFAnalysisType[a].Data()));
6432ac6a	220	fHistPP[a] = (TH1D*)list->FindObject(Form("fHistPP%s",gBFAnalysisType[a].Data()));
	221	fHistPN[a] = (TH1D*)list->FindObject(Form("fHistPN%s",gBFAnalysisType[a].Data()));
	222	fHistNP[a] = (TH1D*)list->FindObject(Form("fHistNP%s",gBFAnalysisType[a].Data()));
	223	fHistNN[a] = (TH1D*)list->FindObject(Form("fHistNN%s",gBFAnalysisType[a].Data()));
	224
a0ce3c86	225	// rebin histograms (be careful with divider!)
	226	fHistP[a]->Rebin(fRebin);
	227	fHistN[a]->Rebin(fRebin);
	228	fHistPP[a]->Rebin(fRebin);
	229	fHistPN[a]->Rebin(fRebin);
	230	fHistNP[a]->Rebin(fRebin);
	231	fHistNN[a]->Rebin(fRebin);
	232
33a86282	233	// set histograms in AliBalance object
	234	bf[iCanvas][a]->SetHistNp(a, fHistP[a]);
	235	bf[iCanvas][a]->SetHistNn(a, fHistN[a]);
	236	bf[iCanvas][a]->SetHistNpp(a, fHistPP[a]);
	237	bf[iCanvas][a]->SetHistNpn(a, fHistPN[a]);
	238	bf[iCanvas][a]->SetHistNnp(a, fHistNP[a]);
	239	bf[iCanvas][a]->SetHistNnn(a, fHistNN[a]);
6432ac6a	240
33a86282	241	gbf[iCanvas][a] = bf[iCanvas][a]->GetBalanceFunctionHistogram(a);
6432ac6a	242	gbf[iCanvas][a]->SetName(Form("%s_BF_%s",listName.Data(),gBFAnalysisType[a].Data()));
	243
	244	cBF[iCanvas]->cd(a+1);
	245	gbf[iCanvas][a]->SetMarkerStyle(20);
	246
	247	if(!bHistos){
a0ce3c86	248	gbf[iCanvas][a]->DrawCopy("AP");
a0ce3c86	249	if(a==1) GetWeightedMean(gbf[iCanvas][a],fStartBinBFWidth,WM[iCanvas-1],WME[iCanvas-1]); // for eta recalculate width (from 0.1 only!)
6432ac6a	250	}
6432ac6a	251	else{
6432ac6a	252	fHistPN[a]->SetLineColor(2);
32a4f715	253	fHistPN[a]->Draw();
	254	fHistPP[a]->SetLineColor(1);
	255	fHistPP[a]->Draw("same");
6432ac6a	256	fHistNP[a]->SetLineColor(4);
	257	fHistNP[a]->Draw("same");
	258	fHistNN[a]->SetLineColor(8);
	259	fHistNN[a]->Draw("same");
	260	}
	261	}
	262	}
	263	// ----------------------------------------------------
	264
	265	// ----------------------------------------------------
	266	// calculate and plot BF (shuffled)
	267	if(listName.Contains("BFShuffled")){
6432ac6a	268	cBFS[iCanvas] = new TCanvas(listName,listName);
	269	cBFS[iCanvas]->Divide(3,3);
	270
	271	for(Int_t a = 0; a < 7; a++){
	272
	273	// create the BF object
	274	bfs[iCanvas][a] = new AliBalance();
	275
32a4f715	276	fHistPS[a] = (TH1D*)list->FindObject(Form("fHistP%s_shuffle",gBFAnalysisType[a].Data()));
	277	fHistNS[a] = (TH1D*)list->FindObject(Form("fHistP%s_shuffle",gBFAnalysisType[a].Data()));
	278	fHistPPS[a] = (TH1D*)list->FindObject(Form("fHistPP%s_shuffle",gBFAnalysisType[a].Data()));
	279	fHistPNS[a] = (TH1D*)list->FindObject(Form("fHistPN%s_shuffle",gBFAnalysisType[a].Data()));
	280	fHistNPS[a] = (TH1D*)list->FindObject(Form("fHistNP%s_shuffle",gBFAnalysisType[a].Data()));
	281	fHistNNS[a] = (TH1D*)list->FindObject(Form("fHistNN%s_shuffle",gBFAnalysisType[a].Data()));
6432ac6a	282
a0ce3c86	283	// rebin histograms (be careful with divider!)
	284	fHistPS[a]->Rebin(fRebin);
	285	fHistNS[a]->Rebin(fRebin);
	286	fHistPPS[a]->Rebin(fRebin);
	287	fHistPNS[a]->Rebin(fRebin);
	288	fHistNPS[a]->Rebin(fRebin);
	289	fHistNNS[a]->Rebin(fRebin);
	290
33a86282	291	// set histograms in AliBalance object
32a4f715	292	bfs[iCanvas][a]->SetHistNp(a, fHistPS[a]);
	293	bfs[iCanvas][a]->SetHistNn(a, fHistNS[a]);
	294	bfs[iCanvas][a]->SetHistNpp(a, fHistPPS[a]);
	295	bfs[iCanvas][a]->SetHistNpn(a, fHistPNS[a]);
	296	bfs[iCanvas][a]->SetHistNnp(a, fHistNPS[a]);
	297	bfs[iCanvas][a]->SetHistNnn(a, fHistNNS[a]);
	298
	299	gbfs[iCanvas][a] = bfs[iCanvas][a]->GetBalanceFunctionHistogram(a);
6432ac6a	300	gbfs[iCanvas][a]->SetName(Form("%s_BF_%s",listName.Data(),gBFAnalysisType[a].Data()));
	301
	302	cBFS[iCanvas]->cd(a+1);
	303	gbfs[iCanvas][a]->SetMarkerStyle(20);
	304	if(!bHistos){
a0ce3c86	305	gbfs[iCanvas][a]->DrawCopy("AP");
a0ce3c86	306	if(a==1) GetWeightedMean(gbfs[iCanvas][a],fStartBinBFWidth,WMS[iCanvas-1],WMSE[iCanvas-1]); // for eta recalculate width (from 0.1 only!)
6432ac6a	307	}
6432ac6a	308	else{
32a4f715	309	fHistPNS[a]->SetLineColor(2);
	310	fHistPNS[a]->Draw();
	311	fHistPPS[a]->SetLineColor(1);
	312	fHistPPS[a]->Draw("same");
	313	fHistNPS[a]->SetLineColor(4);
	314	fHistNPS[a]->Draw("same");
	315	fHistNNS[a]->SetLineColor(8);
	316	fHistNNS[a]->Draw("same");
6432ac6a	317	}
	318	}
	319	}
	320	// ----------------------------------------------------
a114d234	321	}
a114d234	322
a0ce3c86	323	// for BF calculation create also graphs with weighted mean for eta
	324	TGraphErrors *gWM = NULL;
	325	TGraphErrors *gWMS = NULL;
	326	if(!bHistos && gbf[1][1]){
	327	gWM = new TGraphErrors(iCanvas,cent,WM,centE,WME);
	328	gWM->SetName("gCentrality");
	329	}
	330	if(!bHistos && gbfs[1][1]){
	331	gWMS = new TGraphErrors(iCanvas,cent,WMS,centE,WMSE);
	332	gWMS->SetName("gCentralityS");
	333	}
	334
	335	TFile *fOut = TFile::Open(Form("Histograms_WMstart%d_rebin%d_%s", fStartBinBFWidth, fRebin,inFile.Data()),"RECREATE");
6432ac6a	336	fOut->cd();
a0ce3c86	337	for(Int_t i = 0; i <= iCanvas; i++){
6432ac6a	338	for(Int_t a = 0; a < 7; a++){
a0ce3c86	339	if(gbf[i][a]){
	340	gbf[i][a]->Write();
	341	gbf[i][a]->Delete();
	342	}
	343	if(gbfs[i][a]){
	344	gbfs[i][a]->Write();
	345	gbfs[i][a]->Delete();
	346	}
6432ac6a	347	}
6432ac6a	348	}
a0ce3c86	349
	350	if(gWM) gWM->Write();
	351	if(gWMS) gWMS->Write();
	352
	353	fOut->Close();
	354	f->Close();
	355	gROOT->Reset();
	356	}
	357
	358	//____________________________________________________________________//
	359	void GetWeightedMean(TH1D *gHistBalance, Int_t fStartBin = 1, Double_t &WM, Double_t &WME) {
	360
	361	//Prints the calculated width of the BF and its error
	362	Double_t gSumXi = 0.0, gSumBi = 0.0, gSumBiXi = 0.0;
	363	Double_t gSumBiXi2 = 0.0, gSumBi2Xi2 = 0.0;
	364	Double_t gSumDeltaBi2 = 0.0, gSumXi2DeltaBi2 = 0.0;
	365	Double_t deltaBalP2 = 0.0, integral = 0.0;
	366	Double_t deltaErrorNew = 0.0;
	367
	368	//Retrieve this variables from Histogram
	369	Int_t fNumberOfBins = gHistBalance->GetNbinsX();
	370	Double_t fP2Step = gHistBalance->GetBinWidth(1); // assume equal binning!
	371
	372	cout<<"=================================================="<<endl;
	373	cout<<"RECALCULATION OF BF WIDTH (StartBin = "<<fStartBin<<")"<<endl;
	374	cout<<"HISTOGRAM has "<<fNumberOfBins<<" bins with bin size of "<<fP2Step<<endl;
	375	for(Int_t i = fStartBin; i <= fNumberOfBins; i++) {
	376	cout<<"B: "<<gHistBalance->GetBinContent(i)<<"\t Error: "<<gHistBalance->GetBinError(i)<<"\t bin: "<<gHistBalance->GetBinCenter(i)<<endl;
	377	}
	378	cout<<"=================================================="<<endl;
	379	for(Int_t i = fStartBin; i <= fNumberOfBins; i++) {
	380	gSumXi += gHistBalance->GetBinCenter(i);
	381	gSumBi += gHistBalance->GetBinContent(i);
	382	gSumBiXi += gHistBalance->GetBinContent(i)*gHistBalance->GetBinCenter(i);
	383	gSumBiXi2 += gHistBalance->GetBinContent(i)*TMath::Power(gHistBalance->GetBinCenter(i),2);
	384	gSumBi2Xi2 += TMath::Power(gHistBalance->GetBinContent(i),2)*TMath::Power(gHistBalance->GetBinCenter(i),2);
	385	gSumDeltaBi2 += TMath::Power(gHistBalance->GetBinError(i),2);
	386	gSumXi2DeltaBi2 += TMath::Power(gHistBalance->GetBinCenter(i),2) * TMath::Power(gHistBalance->GetBinError(i),2);
	387
	388	deltaBalP2 += fP2Step*TMath::Power(gHistBalance->GetBinError(i),2);
	389	integral += fP2Step*gHistBalance->GetBinContent(i);
	390	}
	391	for(Int_t i = 1; i < fNumberOfBins; i++)
	392	deltaErrorNew += gHistBalance->GetBinError(i)(gHistBalance->GetBinCenter(i)gSumBi - gSumBiXi)/TMath::Power(gSumBi,2);
	393
	394	Double_t integralError = TMath::Sqrt(deltaBalP2);
	395
	396	Double_t delta = gSumBiXi / gSumBi;
	397	Double_t deltaError = (gSumBiXi / gSumBi) * TMath::Sqrt(TMath::Power((TMath::Sqrt(gSumXi2DeltaBi2)/gSumBiXi),2) + TMath::Power((gSumDeltaBi2/gSumBi),2) );
	398
	399	cout<<"Width: "<<delta<<"\t Error: "<<deltaError<<endl;
	400	cout<<"New error: "<<deltaErrorNew<<endl;
	401	cout<<"Integral: "<<integral<<"\t Error: "<<integralError<<endl;
	402	cout<<"=================================================="<<endl;
	403
	404	WM = delta;
	405	WME = deltaError;
61483909	406	}
61483909	407
33a86282	408	//___________________________________________________________//
33a86282	409	// NOT USED any more
61483909	410	//___________________________________________________________//
	411	void mergeOutput(const char* outputDir) {
	412	//Function to merge the output of the sub-jobs
61483909	413	//Create a BF object
	414	AliBalance *bf = new AliBalance();
	415
	416	//connect to AliEn's API services
	417	TGrid::Connect("alien://");
	418
	419	//Getting the output dir from the env. variable
	420	//(JDL field: JDLVariables={"OutputDir"};)
	421	TGridResult* result = gGrid->Query(outputDir,"*/root_archive.zip","","-l 1000");
	422
	423	Int_t nEntries = result->GetEntries();
	424
	425	TString alienUrl;
	426	TDirectoryFile *dirSubJob;
	427
	428	TString gCutName[4] = {"Total","Offline trigger",
	429	"Vertex","Analyzed"};
	430	TH1F *fHistEventStats = new TH1F("fHistEventStats",
	431	"Event statistics;;N_{events}",
	432	4,0.5,4.5);
	433	for(Int_t i = 1; i <= 4; i++)
	434	fHistEventStats->GetXaxis()->SetBinLabel(i,gCutName[i-1].Data());
	435
	436	AliESDtrackCuts *bfTrackCuts = new AliESDtrackCuts("bfTrackCuts");
	437	for(Int_t i = 0; i < nEntries; i++) {
	438	alienUrl = result->GetKey(i,"turl");
	439	alienUrl += "#AnalysisResults.root";
	440	Printf("Opening file: %s",alienUrl.Data());
	441	TFile *file = TFile::Open(alienUrl.Data());
	442	dirSubJob = dynamic_cast<TDirectoryFile *>(file->Get("PWG2EbyE.outputBalanceFunctionAnalysis.root"));
	443
	444	//merge BF
	445	AliBalance bfSubJob = dynamic_cast<AliBalance >(dirSubJob->Get("AliBalance"));
a114d234	446	//bfSubJob->PrintResults();
61483909	447	bf->Merge(bfSubJob);
	448	//delete bfSubJob;
	449
	450	//merge event stats
	451	TList listSubJob = dynamic_cast<TList >(dirSubJob->Get("listQA"));
	452	fHistEventStats->Add(dynamic_cast<TH1F *>(listSubJob->At(0)));
	453
	454	bfTrackCuts = dynamic_cast<AliESDtrackCuts *>(listSubJob->At(1));
	455	delete listSubJob;
	456	}
	457
	458	//Create the output file
	459	TString outputFile = "AnalysisResults.Merged.root";
	460	TFile *foutput = TFile::Open(outputFile.Data(),"recreate");
	461	TDirectoryFile *dirOutput = new TDirectoryFile();
	462	dirOutput->SetName("PWG2EbyE.outputBalanceFunctionAnalysis.root");
	463	//dirOutput->cd();
	464	dirOutput->Add(bf);
	465	TList *list = new TList();
	466	list->SetName("listQA");
	467	list->Add(fHistEventStats);
	468	list->Add(bfTrackCuts);
	469	dirOutput->Add(list);
	470	dirOutput->Write();
	471	bf->Write();
	472	list->Write();
	473	foutput->Close();
	474
	475	//cout<<alienUrl.Data()<<endl;
	476	}