AliTOFanalyzeMatching.C

void AliTOFanalyzeMatching(const char* datafile)
{

  //
  // Matching efficiency and contamination
  // for different particle species
  // (pions, kaons and protons).
  // All histos are saved into a separate file
  // datafile is assumed to be the file name containing
  // the results of the matching in TNtuple format.
  //
  // Author: F. Pierella | pierella@bo.infn.it
  //
  // Use case:
  // start root
  // root[0] .L AliTOFanalyzeMatching.C
  // root[1] AliTOFanalyzeMatching("matchingNtuple.root")

  // output (histos!) filename
  char outFileName[100];
  strcpy(outFileName,"histo");
  strcat(outFileName,datafile);
  
  // dummy histos (for normalization)
  TH1F* hpitot= new TH1F("hpitot","",12,0.,3.);
  TH1F* hkatot= new TH1F("hkatot","",12,0.,3.);
  TH1F* hprtot= new TH1F("hprtot","",12,0.,3.);
  
  TH1F* hpimatched= new TH1F("hpimatched","",12,0.,3.);
  TH1F* hkamatched= new TH1F("hkamatched","",12,0.,3.);
  TH1F* hprmatched= new TH1F("hprmatched","",12,0.,3.);
  
  
  // matching efficiency histos
  TH1F* hpimatcheff= new TH1F("hpimatcheff","Matching efficiency for pions",12,0.,3.);
  TH1F* hkamatcheff= new TH1F("hkamatcheff","Matching efficiency for kaons",12,0.,3.);
  TH1F* hprmatcheff= new TH1F("hprmatcheff","Matching efficiency for protons",12,0.,3.);

  // matching contamination histos
  TH1F* hpimatchcon= new TH1F("hpimatchcon","Matching contamination for pions",12,0.,3.);
  TH1F* hkamatchcon= new TH1F("hkamatchcon","Matching contamination for kaons",12,0.,3.);
  TH1F* hprmatchcon= new TH1F("hprmatchcon","Matching contamination for protons",12,0.,3.);
  
  
  TFile *file = TFile::Open(datafile,"old");
  TNtuple* fNtuple= (TNtuple*)file->Get("Ntuple"); // get ntuple from file
  Int_t nvar = fNtuple->GetNvar(); cout <<"N of var.="<< nvar << endl;
  fNtuple->GetEvent(0);
  
  file->cd();
  Int_t nparticles = (Int_t)fNtuple->GetEntries();
  
  for (Int_t i=0; i < nparticles; i++) {
    fNtuple->GetEvent(i);
    Int_t event=fNtuple->GetLeaf("event")->GetValue();   
    Int_t pdgcode=fNtuple->GetLeaf("ipart")->GetValue(); 
    Int_t matc=fNtuple->GetLeaf("matc")->GetValue(0);
    Float_t px=fNtuple->GetLeaf("pxvtx")->GetValue(0);
    Float_t py=fNtuple->GetLeaf("pyvtx")->GetValue(0);
    Float_t pz=fNtuple->GetLeaf("pzvtx")->GetValue(0);
    
    Float_t pvtx=TMath::Sqrt(px*px+py*py+pz*pz);  
    Float_t ptvtx=TMath::Sqrt(px*px+py*py);
    Int_t abspdgcode=TMath::Abs(pdgcode);
    

    // N (1+2+3+4+(-4)) cases
    if(matc>=1 || matc==-4){
      switch(abspdgcode){
      case 211:
	hpitot->Fill(pvtx);
	break;
      case 321:
	hkatot->Fill(pvtx);
	break;
      case 2212:
	hprtot->Fill(pvtx);
	break;
      }
    }


    // N_matched (3+4) cases
    if(matc==3 || matc==4){
      switch(abspdgcode){
      case 211:
	hpimatched->Fill(pvtx);
	break;
      case 321:
	hkamatched->Fill(pvtx);
	break;
      case 2212:
	hprmatched->Fill(pvtx);
	break;
      }
    }


    // N_t (3) case
    if(matc==3){
      switch(abspdgcode){
      case 211:
	hpimatcheff->Fill(pvtx);
	break;
      case 321:
	hkamatcheff->Fill(pvtx);
	break;
      case 2212:
	hprmatcheff->Fill(pvtx);
	break;
      }
    }

    // N_w (4) case
    if(matc==4){
      switch(abspdgcode){
      case 211:
	hpimatchcon->Fill(pvtx);
	break;
      case 321:
	hkamatchcon->Fill(pvtx);
	break;
      case 2212:
	hprmatchcon->Fill(pvtx);
	break;
      }
    }

  }

  // histo normalization
  // efficiency
  hpimatcheff->Divide(hpitot);
  hkamatcheff->Divide(hkatot);
  hprmatcheff->Divide(hprtot);

  // contamination
  hpimatchcon->Divide(hpimatched);
  hkamatchcon->Divide(hkamatched);
  hprmatchcon->Divide(hprmatched);


  TFile *houtfile = new TFile(outFileName,"recreate");
  houtfile->cd();

  hpitot->Write();
  hkatot->Write();
  hprtot->Write();
    
  hpimatched->Write();
  hkamatched->Write();
  hprmatched->Write();
    
  hpimatcheff->Write();
  hkamatcheff->Write();
  hprmatcheff->Write();

  hpimatchcon->Write();
  hkamatchcon->Write();
  hprmatchcon->Write();
  houtfile->Close();   

  cout << "File " << outFileName << " with histos has been created" << endl;
}