[u/mrichter/AliRoot.git] / JETAN / createEvents.C

// $Id$

#ifndef __CINT__
#include <stdio.h>
#include <Riostream.h>
#include <time.h>
#include <TROOT.h>
#include <TCanvas.h>
#include <TRandom.h>
#include <TSystem.h>
#include <TFile.h>
#include <TMath.h>
#include <TChain.h>
#include <TTree.h>
#include <TStopwatch.h>
#include <TObjArray.h>
#include <TObjString.h>
#include <TString.h>
#include <TCanvas.h>
#include <TH2F.h>
#include <TPaveText.h>
#include <TStyle.h>
#include <AliStack.h>
#include <AliRunLoader.h>
#include <AliHeader.h>
#include <AliGenPythiaEventHeader.h>
#include <AliJetParticle.h>
#include <AliJetParticlesReader.h>
#include <AliJetParticlesReaderKine.h>
#include <AliJetParticlesReaderKineGoodTPC.h>
#include <AliJetParticlesReaderESD.h>
#include <AliJetParticlesReaderHLT.h>
#include <AliJetEventParticles.h>
#endif

#if 0
    kITSin=0x0001,kITSout=0x0002,kITSrefit=0x0004,kITSpid=0x0008,
    kTPCin=0x0010,kTPCout=0x0020,kTPCrefit=0x0040,kTPCpid=0x0080,
    kTRDin=0x0100,kTRDout=0x0200,kTRDrefit=0x0400,kTRDpid=0x0800,
    kTOFin=0x1000,kTOFout=0x2000,kTOFrefit=0x4000,kTOFpid=0x8000,
    kPHOSpid=0x10000, kRICHpid=0x20000, kEMCALpid=0x40000,
    kTRDbackup=0x80000,
    kTRDStop=0x20000000,
    kESDpid=0x40000000,
    kTIME=0x80000000
#endif

//#define constrained
#define useHoughITS

void createEvents(Char_t *dir, Int_t files,Char_t *sdir=".");
void createEvents(Int_t esd, TObjArray *dirs, Float_t ptcut, Char_t *savefile);
void createEvents(Int_t esd,Char_t *dir="./input",Float_t ptcut=0.05,Char_t *savefile=0);
void createEvents(Int_t esd,Char_t *dirfile,Char_t *dir,Float_t ptcut=0.05,Char_t *savefile=0);
void createMixedEvents(Char_t *fname,Char_t *hname,Char_t *savefile=0,Int_t nMaxEvents=-1);
void createDirs(Char_t *dir,Int_t files,Char_t *savefile);
void testEvents(Char_t *sdir);
void readEvents(Char_t *filename,Int_t nMaxEvents=-1);
void displayMixedEvents(Char_t *fname,Char_t *hname,Float_t ptcut=2.0,Int_t nMaxEvents=-1);

//-----------------------

void createEvents(Char_t *dir, Int_t files, Char_t *sdir)
{
  createDirs(dir,files,"./dirlist.txt");
  Char_t filename[1000];
  sprintf(filename,"%s/aliev-type0.root",sdir);
  createEvents(0,"./dirlist.txt",dir,0.4,filename);
  sprintf(filename,"%s/aliev-type1.root",sdir);
  createEvents(1,"./dirlist.txt",dir,0.4,filename);
  sprintf(filename,"%s/aliev-type2.root",sdir);
  createEvents(2,"./dirlist.txt",dir,0.4,filename);
  sprintf(filename,"%s/aliev-type3.root",sdir);
  createEvents(3,"./dirlist.txt",dir,0.4,filename);
#ifdef useHoughITS
  sprintf(filename,"%s/aliev-type4.root",sdir);
  createEvents(4,"./dirlist.txt",dir,0.4,filename);
#endif
  sprintf(filename,"%s/aliev-type10.root",sdir);
  createEvents(10,"./dirlist.txt",dir,0.4,filename);
  sprintf(filename,"%s/aliev-type11.root",sdir);
  createEvents(11,"./dirlist.txt",dir,0.4,filename);
  sprintf(filename,"%s/aliev-type12.root",sdir);
  createEvents(12,"./dirlist.txt",dir,0.4,filename);
}

void createEvents(Int_t esd,Char_t *dirfile,Char_t *dir,
                  Float_t ptcut,Char_t *savefile)
{
  TObjArray *dirs=new TObjArray(10000);

  FILE *c=fopen(dirfile,"r");
  if(!c) return;

  Int_t saveErrIgLevel=gErrorIgnoreLevel;
  gErrorIgnoreLevel=kFatal;
  while(!feof(c)){
    Char_t adddir[1024];
    fscanf(c,"%s\n",adddir);
    dirs->Add(new TObjString(adddir));
    //cout << adddir << endl;
  }
  gErrorIgnoreLevel=saveErrIgLevel;

  Char_t buffer[8096];
  if(savefile)
    sprintf(buffer,"%s",savefile);
  else
    sprintf(buffer,"%s/aliev-type%d.root",dir,esd);

  createEvents(esd,dirs,ptcut,buffer);
}

void createEvents(Int_t esd,Char_t *dir,
		  Float_t ptcut,Char_t *savefile)
{
  TObjArray *dirs=new TObjArray(1);
  dirs->Add(new TObjString(dir));

  Char_t buffer[8096];
  if(savefile)
    sprintf(buffer,"%s",savefile);
  else
    sprintf(buffer,"%s/aliev-type%d.root",dir,esd);

  createEvents(esd,dirs,ptcut,buffer);
}

void createEvents(Int_t esd, TObjArray *dirs, Float_t ptcut, Char_t *savefile)
{
  AliJetParticlesReader *reader=0;
  TString desc="";
  if (esd==0) {
    reader=new AliJetParticlesReaderKine(dirs);
    //((AliJetParticlesReaderKine*)reader)->SetUseTracks(kTRUE);
    ((AliJetParticlesReaderKine*)reader)->SetCharged(kTRUE);
    ((AliJetParticlesReaderKine*)reader)->SetEM(kFALSE);
    ((AliJetParticlesReaderKine*)reader)->SetNeutral(kFALSE);
    desc+="Kine";
  } else if(esd==1){
#ifdef constrained
    reader=new AliJetParticlesReaderESD(1,dirs);
    ((AliJetParticlesReaderESD*)reader)->SetCompareFlag(
	                   AliESDtrack::kITSrefit+AliESDtrack::kTPCrefit+AliESDtrack::kTRDrefit);
#else
    reader=new AliJetParticlesReaderESD(0,dirs);
    ((AliJetParticlesReaderESD*)reader)->SetCompareFlag(AliESDtrack::kTPCin);
#endif
    desc+="ESD";
  } else if (esd==2) {
    reader=new AliJetParticlesReaderHLT(kTRUE,dirs);
    //((AliJetParticlesReaderHLT*)reader)->SetMinHits(30);
    desc+="HLTConf";
  } else if (esd==3) {
    reader=new AliJetParticlesReaderHLT(kFALSE,dirs);
    desc+="HLTHough";
  } else if (esd==4) {
    reader=new AliJetParticlesReaderHLT(kFALSE,dirs,"AliESDits.root");
    //((AliJetParticlesReaderHLT*)reader)->SetMinHits(4);
    desc+="HLTITSHough";
  } else if (esd==10) {
    reader=new AliJetParticlesReaderKine(dirs);
    //((AliJetParticlesReaderKine*)reader)->SetUseTracks(kTRUE);
    ((AliJetParticlesReaderKine*)reader)->SetCharged(kTRUE);
    ((AliJetParticlesReaderKine*)reader)->SetEM(kTRUE);
    ((AliJetParticlesReaderKine*)reader)->SetNeutral(kTRUE);
    desc+="Kine-All";
  } else if (esd==11) {
    reader=new AliJetParticlesReaderKine(dirs);
    //((AliJetParticlesReaderKine*)reader)->SetUseTracks(kTRUE);
    ((AliJetParticlesReaderKine*)reader)->SetCharged(kTRUE);
    ((AliJetParticlesReaderKine*)reader)->SetEM(kTRUE);
    ((AliJetParticlesReaderKine*)reader)->SetNeutral(kFALSE);
    desc+="Kine-EM";
  } else if (esd==12) {
    reader=new AliJetParticlesReaderKineGoodTPC(dirs);
    desc+="GoodTPC";
  } else {
    cout << "Parameter settings: " << endl;
    cout << "0  == Kine (ch)"  << endl;
    cout << "1  == TPC"  << endl;
    cout << "2  == Conformal"  << endl;
    cout << "3  == Hough"  << endl;
    cout << "4  == ITSHough"  << endl;
    cout << "10 == Kine (all)"  << endl;
    cout << "11 == Kine (ch+em)"  << endl;
    cout << "12 == Kine (good tpc)"  << endl;
    return;
  }
  //reader->ReadEventsFromTo(1,0);
  reader->SetPtCut(ptcut,100);
  reader->SetEtaCut(-0.9,0.9);

  desc+=" Ptcut: ";
  desc+=ptcut;

  TFile* file = new TFile(savefile,"RECREATE");
  TTree* tree = new TTree("AJEPtree","AliJetEventParticles Tree");
  reader->SetTree(tree);
  Int_t counter=0;
  while(reader->Next())
    {
      AliJetEventParticles *ev=new AliJetEventParticles(*reader->GetEventParticles());
      if(gErrorIgnoreLevel<=kWarning){
	cout << "Read event: " << ++counter << endl;
	ev->Print();
      }
      delete ev;
    }
  file->Write();
  file->Close();
  delete file;
  delete reader;
}

void createDirs(Char_t *dir,Int_t files,Char_t *savefile)
{
  Char_t buffer[8096];
  if(savefile)
    sprintf(buffer,"%s",savefile);
  else
    sprintf(buffer,"%s/dirinput.txt",dir);

  FILE *c=fopen(buffer,"w");
  if(!c) return;

  Int_t saveErrIgLevel=gErrorIgnoreLevel;
  gErrorIgnoreLevel=kFatal;
  for(Int_t i=0;i<files;i++){
    Char_t adddir[1024];
    //sprintf(adddir,"%s/%05d",dir,i);
    sprintf(adddir,"%s/%d",dir,i);
    Char_t fname[1024];
    sprintf(fname,"%s/galice.root",adddir);
    TFile f(fname);
    if(!f.IsOpen()) continue;
    //cout << f.GetErrno() << endl;
    f.Close();
    AliRunLoader *r = AliRunLoader::Open(fname); 
    if(!r) continue;
    if(r->GetNumberOfEvents() <= 0){
      delete r;
      continue;
    }
    if(r->LoadKinematics()){
      delete r;
      continue;
    }
    delete r;
    sprintf(fname,"%s/AliESDs.root",adddir);
    TFile g(fname);
    if(!g.IsOpen()) continue;

    TTree *tree = dynamic_cast<TTree*>(g.Get("esdTree"));
    if(!tree || tree->GetEntries()<=0){
      g.Close();
      continue;
    }
    g.Close();
#ifdef useHoughITS
    sprintf(fname,"%s/AliESDits.root",adddir);
    TFile h(fname);
    if(!h.IsOpen()) continue;

    tree = dynamic_cast<TTree*>(h.Get("esdTree"));
    if(!tree || tree->GetEntries()<=0){
      h.Close();
      continue;
    }
    h.Close();
#endif
   
    cout << "Added " << adddir << endl;
    fprintf(c,"%s\n",adddir);
  }
  gErrorIgnoreLevel=saveErrIgLevel;
  fclose(c);
}

void createMixedEvents(Char_t *fname,Char_t *hname,Char_t *savefile,Int_t nMaxEvents)
{
  //connect to events
  TChain *theTree = new TChain("AJEPtree");
  theTree->Add(fname);
  AliJetEventParticles *ev=new AliJetEventParticles();
  theTree->SetBranchAddress("particles",&ev);

  Int_t treeentries=(Int_t)theTree->GetEntries();
  if((nMaxEvents<0) || (nMaxEvents>treeentries))
    nMaxEvents=treeentries;
  //cout << "Found " << nMaxEvents << " in " << fname << endl;

  TChain *backtheTree = new TChain("AJEPtree");
  backtheTree->Add(hname);
  AliJetEventParticles *backev=new AliJetEventParticles();
  backtheTree->SetBranchAddress("particles",&backev);

  Int_t backtreeentries=(Int_t)backtheTree->GetEntries();
  Int_t nPerBackground=nMaxEvents/backtreeentries;
  if(nPerBackground==0) nPerBackground=1;

  Char_t buffer[8096];
  if(savefile)
    sprintf(buffer,"%s",savefile);
  else
    sprintf(buffer,"./aliev-mixed.root");
  TFile* file = new TFile(buffer,"RECREATE");
  TTree* tree = new TTree("AJEPtree","AliJetEventParticles Tree");
  AliJetEventParticles *mixedev=new AliJetEventParticles(0);
  tree->Branch("particles","AliJetEventParticles",&mixedev,32000,1);

  //=========================================================================
  // start the event loop
  //=========================================================================
  Int_t nEvent = 0;
  Int_t nEventHijing = -1;
  Int_t nEventHijingCounter = nPerBackground;
  while(nEvent<nMaxEvents){

    //get the event
    theTree->GetEvent(nEvent);

    //load background if needed
    if(nEventHijingCounter==nPerBackground){    
      backev->Reset();
      nEventHijing++;
      backtheTree->GetEvent(nEventHijing);
      if(nEventHijing==backtreeentries) nEventHijing=0;
      nEventHijingCounter=0;
    }

    backev->AddSignal(*ev);
    TString dummy="Counter: ";
    dummy+=nEvent;
    dummy+=" ";
    //dummy+=ev->GetHeader();
    dummy+="(Pythia ";dummy+=nEvent;
    dummy+=" ";
    //dummy+=backev->GetHeader();
    dummy+=", Hijing ";dummy+=nEventHijing;
    dummy+=")";
    mixedev->Set(*backev);
    mixedev->SetHeader(dummy);
    if(gErrorIgnoreLevel<=kWarning){
      cout << "Read event: " << nEvent << " " << nEventHijing << endl;
      mixedev->Print();
    }

    tree->Fill();

    nEvent++;
    nEventHijingCounter++;
    ev->Reset();
    mixedev->Reset();
  } //end of nev loop

  file->Write();
  file->Close();
  delete file;
  delete ev;
  delete theTree;
  delete backtheTree;
  delete backev;
}

void testEvents(Char_t *sdir)
{
  Char_t filename[1000];
  sprintf(filename,"%s/aliev-type0.root",sdir);
  readEvents(filename,0);
  sprintf(filename,"%s/aliev-type1.root",sdir);
  readEvents(filename,0);
  sprintf(filename,"%s/aliev-type2.root",sdir);
  readEvents(filename,0);
  sprintf(filename,"%s/aliev-type3.root",sdir);
  readEvents(filename,0);
  sprintf(filename,"%s/aliev-type10.root",sdir);
  readEvents(filename,0);
  sprintf(filename,"%s/aliev-type11.root",sdir);
  readEvents(filename,0);
  sprintf(filename,"%s/aliev-type12.root",sdir);
}

void readEvents(Char_t *filename,Int_t nMaxEvents)
{
  //connect to events
  TChain *theTree = new TChain("AJEPtree");
  theTree->Add(filename);
  AliJetEventParticles *ev=new AliJetEventParticles();
  theTree->SetBranchAddress("particles",&ev);

  Int_t treeentries=(Int_t)theTree->GetEntries();
  if((nMaxEvents<0) || (nMaxEvents>treeentries))
    nMaxEvents=treeentries;

  cout << "Found " << treeentries << ", but use " << nMaxEvents << " in " << filename << endl;

  //=========================================================================
  // start the event loop
  //=========================================================================
  Int_t nEvent = 0;
  while(nEvent<nMaxEvents){
    if ((nEvent % 100) == 0) {
      cout << "Reading event " << nEvent << endl;
    }

    //connect the event
    theTree->GetEvent(nEvent);

    ev->Print();
#if 0
    const TClonesArray *p=ev->GetParticles();
    for(Int_t i=0;i<p->GetEntriesFast();i++) {
      cout << i << endl;
      ((AliJetParticle*)p->At(i))->Print();
    }
#endif
    ev->Clear();
    nEvent++;
  } //end of nev loop

  delete ev;
  delete theTree;
}

void displayMixedEvents(Char_t *fname,Char_t *hname,Float_t ptcut,Int_t nMaxEvents)
{
  //connect to events
  TChain *theTree = new TChain("AJEPtree");
  theTree->Add(fname);
  AliJetEventParticles *ev=new AliJetEventParticles();
  theTree->SetBranchAddress("particles",&ev);

  Int_t treeentries=(Int_t)theTree->GetEntries();
  if((nMaxEvents<0) || (nMaxEvents>treeentries))
    nMaxEvents=treeentries;
  //cout << "Found " << nMaxEvents << " in " << fname << endl;

  TChain *backtheTree = new TChain("AJEPtree");
  backtheTree->Add(hname);
  AliJetEventParticles *backev=new AliJetEventParticles();
  backtheTree->SetBranchAddress("particles",&backev);

  Int_t backtreeentries=(Int_t)backtheTree->GetEntries();
  Int_t nPerBackground=nMaxEvents/backtreeentries;
  if(nPerBackground==0) nPerBackground=1;

  AliJetEventParticles *mixedev=new AliJetEventParticles(0);

  Int_t neta=20;
  Int_t nphi=62;

  TCanvas *c1=new TCanvas("c1","",1000,500);
  c1->Divide(2,1);
  c1->Update();
  gStyle->SetLabelSize(0.03,"XYZ");
  gStyle->SetTitleOffset(1.0,"XYZ");
  gStyle->SetTitleOffset(1.3,"Y");
  gROOT->ForceStyle();

  TH2F *h2pp=new TH2F("h2pp","E_{T}^{ch} in PP [GeV];#phi;#eta",nphi,0,TMath::TwoPi(),neta,-1,1);
  h2pp->SetStats(0);
  TH2F *h2pb=new TH2F("h2pb","E_{T}^{ch} in PbPb [GeV];#phi;#eta;",nphi,0,TMath::TwoPi(),neta,-1,1);
  h2pb->SetStats(0);

  //=========================================================================
  // start the event loop
  //=========================================================================
  Int_t nEvent = 0;
  Int_t nEventHijing = -1;
  Int_t nEventHijingCounter = nPerBackground;
  while(nEvent<nMaxEvents){

    //get the event
    theTree->GetEvent(nEvent);

    //load background if needed
    if(nEventHijingCounter==nPerBackground){    
      backev->Reset();
      nEventHijing++;
      backtheTree->GetEvent(nEventHijing);
      if(nEventHijing==backtreeentries) nEventHijing=0;
      nEventHijingCounter=0;
    }
    backev->AddSignal(*ev);
    TString dummy="Counter: ";
    dummy+=nEvent;
    dummy+=" ";
    //dummy+=ev->GetHeader();
    dummy+="(Pythia ";dummy+=nEvent;
    dummy+=" ";
    //dummy+=backev->GetHeader();
    dummy+=", Hijing ";dummy+=nEventHijing;
    dummy+=")";
    mixedev->Set(*backev);
    mixedev->SetHeader(dummy);
    if(gErrorIgnoreLevel<=kWarning){
      cout << "Read event: " << nEvent << " " << nEventHijing << endl;
      mixedev->Print();
    }

    // loop over all particles...
    c1->cd(1);
    h2pp->Reset();
    AliJetParticle *aliparticle = NULL;
    TIterator *iter = ev->GetParticles()->MakeIterator();
    while ((aliparticle = (AliJetParticle *) iter->Next()) != NULL) {
      Float_t pt=aliparticle->Pt();
      if(pt<ptcut) continue;
      h2pp->Fill(aliparticle->Phi(),aliparticle->Eta(),pt);
    }
    delete iter;
    h2pp->Draw("Lego");
    c1->cd(2);
    h2pb->Reset();
    iter = mixedev->GetParticles()->MakeIterator();
    while ((aliparticle = (AliJetParticle *) iter->Next()) != NULL) {
      Float_t pt=aliparticle->Pt();
      if(pt<ptcut) continue;
      h2pb->Fill(aliparticle->Phi(),aliparticle->Eta(),pt);
    }
    delete iter;

    h2pb->Draw("Lego");
    c1->Update();
    Char_t c;cin >> c; 
    if(c=='s'){
      Char_t sfilenamehist[1024];
      sprintf(sfilenamehist,"./display-planehists-%d-%d.eps",nEvent,nEventHijingCounter);
      c1->SaveAs(sfilenamehist);
    }
    else if(c=='q') break;

    nEvent++;
    nEventHijingCounter++;
    ev->Reset();
    mixedev->Reset();
  } //end of nev loop


  delete ev;
  delete theTree;
  delete backtheTree;
  delete backev;
}
Commit	Line	Data
b9a6a391	1	// $Id$
	2
	3	#ifndef __CINT__
	4	#include <stdio.h>
ca10471a	5	#include <Riostream.h>
b9a6a391	6	#include <time.h>
	7	#include <TROOT.h>
	8	#include <TCanvas.h>
	9	#include <TRandom.h>
	10	#include <TSystem.h>
	11	#include <TFile.h>
	12	#include <TMath.h>
	13	#include <TChain.h>
	14	#include <TTree.h>
	15	#include <TStopwatch.h>
	16	#include <TObjArray.h>
	17	#include <TObjString.h>
	18	#include <TString.h>
	19	#include <TCanvas.h>
	20	#include <TH2F.h>
	21	#include <TPaveText.h>
	22	#include <TStyle.h>
	23	#include <AliStack.h>
	24	#include <AliRunLoader.h>
	25	#include <AliHeader.h>
	26	#include <AliGenPythiaEventHeader.h>
	27	#include <AliJetParticle.h>
	28	#include <AliJetParticlesReader.h>
	29	#include <AliJetParticlesReaderKine.h>
	30	#include <AliJetParticlesReaderKineGoodTPC.h>
	31	#include <AliJetParticlesReaderESD.h>
	32	#include <AliJetParticlesReaderHLT.h>
	33	#include <AliJetEventParticles.h>
	34	#endif
	35
	36	#if 0
	37	kITSin=0x0001,kITSout=0x0002,kITSrefit=0x0004,kITSpid=0x0008,
	38	kTPCin=0x0010,kTPCout=0x0020,kTPCrefit=0x0040,kTPCpid=0x0080,
	39	kTRDin=0x0100,kTRDout=0x0200,kTRDrefit=0x0400,kTRDpid=0x0800,
	40	kTOFin=0x1000,kTOFout=0x2000,kTOFrefit=0x4000,kTOFpid=0x8000,
	41	kPHOSpid=0x10000, kRICHpid=0x20000, kEMCALpid=0x40000,
	42	kTRDbackup=0x80000,
	43	kTRDStop=0x20000000,
	44	kESDpid=0x40000000,
	45	kTIME=0x80000000
	46	#endif
	47
	48	//#define constrained
	49	#define useHoughITS
	50
	51	void createEvents(Char_t dir, Int_t files,Char_t sdir=".");
	52	void createEvents(Int_t esd, TObjArray dirs, Float_t ptcut, Char_t savefile);
	53	void createEvents(Int_t esd,Char_t dir="./input",Float_t ptcut=0.05,Char_t savefile=0);
	54	void createEvents(Int_t esd,Char_t dirfile,Char_t dir,Float_t ptcut=0.05,Char_t *savefile=0);
	55	void createMixedEvents(Char_t fname,Char_t hname,Char_t *savefile=0,Int_t nMaxEvents=-1);
	56	void createDirs(Char_t dir,Int_t files,Char_t savefile);
	57	void testEvents(Char_t *sdir);
	58	void readEvents(Char_t *filename,Int_t nMaxEvents=-1);
	59	void displayMixedEvents(Char_t fname,Char_t hname,Float_t ptcut=2.0,Int_t nMaxEvents=-1);
	60
	61	//-----------------------
	62
	63	void createEvents(Char_t dir, Int_t files, Char_t sdir)
	64	{
	65	createDirs(dir,files,"./dirlist.txt");
	66	Char_t filename[1000];
	67	sprintf(filename,"%s/aliev-type0.root",sdir);
	68	createEvents(0,"./dirlist.txt",dir,0.4,filename);
	69	sprintf(filename,"%s/aliev-type1.root",sdir);
70	createEvents(1,"./dirlist.txt",dir,0.4,filename);
71	sprintf(filename,"%s/aliev-type2.root",sdir);
72	createEvents(2,"./dirlist.txt",dir,0.4,filename);
73	sprintf(filename,"%s/aliev-type3.root",sdir);
74	createEvents(3,"./dirlist.txt",dir,0.4,filename);
75	#ifdef useHoughITS
76	sprintf(filename,"%s/aliev-type4.root",sdir);
77	createEvents(4,"./dirlist.txt",dir,0.4,filename);
78	#endif
79	sprintf(filename,"%s/aliev-type10.root",sdir);
80	createEvents(10,"./dirlist.txt",dir,0.4,filename);
81	sprintf(filename,"%s/aliev-type11.root",sdir);
82	createEvents(11,"./dirlist.txt",dir,0.4,filename);
83	sprintf(filename,"%s/aliev-type12.root",sdir);
84	createEvents(12,"./dirlist.txt",dir,0.4,filename);
85	}
86
87	void createEvents(Int_t esd,Char_t dirfile,Char_t dir,
88	Float_t ptcut,Char_t *savefile)
89	{
90	TObjArray *dirs=new TObjArray(10000);
91
92	FILE *c=fopen(dirfile,"r");
93	if(!c) return;
94
95	Int_t saveErrIgLevel=gErrorIgnoreLevel;
96	gErrorIgnoreLevel=kFatal;
97	while(!feof(c)){
98	Char_t adddir[1024];
99	fscanf(c,"%s\n",adddir);
100	dirs->Add(new TObjString(adddir));
101	//cout << adddir << endl;
102	}
103	gErrorIgnoreLevel=saveErrIgLevel;
104
105	Char_t buffer[8096];
106	if(savefile)
107	sprintf(buffer,"%s",savefile);
108	else
109	sprintf(buffer,"%s/aliev-type%d.root",dir,esd);
110
111	createEvents(esd,dirs,ptcut,buffer);
112	}
113
114	void createEvents(Int_t esd,Char_t *dir,
115	Float_t ptcut,Char_t *savefile)
116	{
117	TObjArray *dirs=new TObjArray(1);
118	dirs->Add(new TObjString(dir));
119
120	Char_t buffer[8096];
121	if(savefile)
122	sprintf(buffer,"%s",savefile);
123	else
124	sprintf(buffer,"%s/aliev-type%d.root",dir,esd);
125
126	createEvents(esd,dirs,ptcut,buffer);
127	}
128
129	void createEvents(Int_t esd, TObjArray dirs, Float_t ptcut, Char_t savefile)
130	{
131	AliJetParticlesReader *reader=0;
132	TString desc="";
133	if (esd==0) {
134	reader=new AliJetParticlesReaderKine(dirs);
135	//((AliJetParticlesReaderKine*)reader)->SetUseTracks(kTRUE);
136	((AliJetParticlesReaderKine*)reader)->SetCharged(kTRUE);
137	((AliJetParticlesReaderKine*)reader)->SetEM(kFALSE);
138	((AliJetParticlesReaderKine*)reader)->SetNeutral(kFALSE);
139	desc+="Kine";
140	} else if(esd==1){
141	#ifdef constrained
142	reader=new AliJetParticlesReaderESD(1,dirs);
143	((AliJetParticlesReaderESD*)reader)->SetCompareFlag(
144	AliESDtrack::kITSrefit+AliESDtrack::kTPCrefit+AliESDtrack::kTRDrefit);
145	#else
146	reader=new AliJetParticlesReaderESD(0,dirs);
147	((AliJetParticlesReaderESD*)reader)->SetCompareFlag(AliESDtrack::kTPCin);
148	#endif
149	desc+="ESD";
150	} else if (esd==2) {
151	reader=new AliJetParticlesReaderHLT(kTRUE,dirs);
152	//((AliJetParticlesReaderHLT*)reader)->SetMinHits(30);
153	desc+="HLTConf";
154	} else if (esd==3) {
155	reader=new AliJetParticlesReaderHLT(kFALSE,dirs);
156	desc+="HLTHough";
157	} else if (esd==4) {
158	reader=new AliJetParticlesReaderHLT(kFALSE,dirs,"AliESDits.root");
159	//((AliJetParticlesReaderHLT*)reader)->SetMinHits(4);
160	desc+="HLTITSHough";
161	} else if (esd==10) {
162	reader=new AliJetParticlesReaderKine(dirs);
163	//((AliJetParticlesReaderKine*)reader)->SetUseTracks(kTRUE);
164	((AliJetParticlesReaderKine*)reader)->SetCharged(kTRUE);
165	((AliJetParticlesReaderKine*)reader)->SetEM(kTRUE);
166	((AliJetParticlesReaderKine*)reader)->SetNeutral(kTRUE);
167	desc+="Kine-All";
168	} else if (esd==11) {
169	reader=new AliJetParticlesReaderKine(dirs);
170	//((AliJetParticlesReaderKine*)reader)->SetUseTracks(kTRUE);
171	((AliJetParticlesReaderKine*)reader)->SetCharged(kTRUE);
172	((AliJetParticlesReaderKine*)reader)->SetEM(kTRUE);
173	((AliJetParticlesReaderKine*)reader)->SetNeutral(kFALSE);
174	desc+="Kine-EM";
175	} else if (esd==12) {
176	reader=new AliJetParticlesReaderKineGoodTPC(dirs);
177	desc+="GoodTPC";
178	} else {
179	cout << "Parameter settings: " << endl;
180	cout << "0 == Kine (ch)" << endl;
181	cout << "1 == TPC" << endl;
182	cout << "2 == Conformal" << endl;
183	cout << "3 == Hough" << endl;
184	cout << "4 == ITSHough" << endl;
185	cout << "10 == Kine (all)" << endl;
186	cout << "11 == Kine (ch+em)" << endl;
187	cout << "12 == Kine (good tpc)" << endl;
188	return;
189	}
190	//reader->ReadEventsFromTo(1,0);
191	reader->SetPtCut(ptcut,100);
192	reader->SetEtaCut(-0.9,0.9);
193
194	desc+=" Ptcut: ";
195	desc+=ptcut;
196
197	TFile* file = new TFile(savefile,"RECREATE");
198	TTree* tree = new TTree("AJEPtree","AliJetEventParticles Tree");
199	reader->SetTree(tree);
200	Int_t counter=0;
201	while(reader->Next())
202	{
203	AliJetEventParticles ev=new AliJetEventParticles(reader->GetEventParticles());
204	if(gErrorIgnoreLevel<=kWarning){
205	cout << "Read event: " << ++counter << endl;
206	ev->Print();
207	}
208	delete ev;
209	}
210	file->Write();
211	file->Close();
212	delete file;
213	delete reader;
214	}
215
216	void createDirs(Char_t dir,Int_t files,Char_t savefile)
217	{
218	Char_t buffer[8096];
219	if(savefile)
220	sprintf(buffer,"%s",savefile);
221	else
222	sprintf(buffer,"%s/dirinput.txt",dir);
223
224	FILE *c=fopen(buffer,"w");
225	if(!c) return;
226
227	Int_t saveErrIgLevel=gErrorIgnoreLevel;
228	gErrorIgnoreLevel=kFatal;
229	for(Int_t i=0;i<files;i++){
230	Char_t adddir[1024];
231	//sprintf(adddir,"%s/%05d",dir,i);
232	sprintf(adddir,"%s/%d",dir,i);
233	Char_t fname[1024];
234	sprintf(fname,"%s/galice.root",adddir);
235	TFile f(fname);
236	if(!f.IsOpen()) continue;
237	//cout << f.GetErrno() << endl;
238	f.Close();
239	AliRunLoader *r = AliRunLoader::Open(fname);
240	if(!r) continue;
241	if(r->GetNumberOfEvents() <= 0){
242	delete r;
243	continue;
244	}
245	if(r->LoadKinematics()){
246	delete r;
247	continue;
248	}
249	delete r;
250	sprintf(fname,"%s/AliESDs.root",adddir);
251	TFile g(fname);
252	if(!g.IsOpen()) continue;
253
254	TTree tree = dynamic_cast<TTree>(g.Get("esdTree"));
255	if(!tree \|\| tree->GetEntries()<=0){
256	g.Close();
257	continue;
258	}
259	g.Close();
260	#ifdef useHoughITS
261	sprintf(fname,"%s/AliESDits.root",adddir);
262	TFile h(fname);
263	if(!h.IsOpen()) continue;
264
265	tree = dynamic_cast<TTree*>(h.Get("esdTree"));
266	if(!tree \|\| tree->GetEntries()<=0){
267	h.Close();
268	continue;
269	}
270	h.Close();
271	#endif
272
273	cout << "Added " << adddir << endl;
274	fprintf(c,"%s\n",adddir);
275	}
276	gErrorIgnoreLevel=saveErrIgLevel;
277	fclose(c);
278	}
279
280	void createMixedEvents(Char_t fname,Char_t hname,Char_t *savefile,Int_t nMaxEvents)
281	{
282	//connect to events
283	TChain *theTree = new TChain("AJEPtree");
284	theTree->Add(fname);
285	AliJetEventParticles *ev=new AliJetEventParticles();
286	theTree->SetBranchAddress("particles",&ev);
287
288	Int_t treeentries=(Int_t)theTree->GetEntries();
289	if((nMaxEvents<0) \|\| (nMaxEvents>treeentries))
290	nMaxEvents=treeentries;
291	//cout << "Found " << nMaxEvents << " in " << fname << endl;
292
293	TChain *backtheTree = new TChain("AJEPtree");
294	backtheTree->Add(hname);
295	AliJetEventParticles *backev=new AliJetEventParticles();
296	backtheTree->SetBranchAddress("particles",&backev);
297
298	Int_t backtreeentries=(Int_t)backtheTree->GetEntries();
299	Int_t nPerBackground=nMaxEvents/backtreeentries;
300	if(nPerBackground==0) nPerBackground=1;
301
302	Char_t buffer[8096];
303	if(savefile)
304	sprintf(buffer,"%s",savefile);
305	else
306	sprintf(buffer,"./aliev-mixed.root");
307	TFile* file = new TFile(buffer,"RECREATE");
308	TTree* tree = new TTree("AJEPtree","AliJetEventParticles Tree");
309	AliJetEventParticles *mixedev=new AliJetEventParticles(0);
310	tree->Branch("particles","AliJetEventParticles",&mixedev,32000,1);
311
312	//=========================================================================
313	// start the event loop
314	//=========================================================================
315	Int_t nEvent = 0;
316	Int_t nEventHijing = -1;
317	Int_t nEventHijingCounter = nPerBackground;
318	while(nEvent<nMaxEvents){
319
320	//get the event
321	theTree->GetEvent(nEvent);
322
323	//load background if needed
324	if(nEventHijingCounter==nPerBackground){
325	backev->Reset();
326	nEventHijing++;
327	backtheTree->GetEvent(nEventHijing);
328	if(nEventHijing==backtreeentries) nEventHijing=0;
329	nEventHijingCounter=0;
330	}
331
332	backev->AddSignal(*ev);
333	TString dummy="Counter: ";
334	dummy+=nEvent;
335	dummy+=" ";
336	//dummy+=ev->GetHeader();
337	dummy+="(Pythia ";dummy+=nEvent;
338	dummy+=" ";
339	//dummy+=backev->GetHeader();
340	dummy+=", Hijing ";dummy+=nEventHijing;
341	dummy+=")";
342	mixedev->Set(*backev);
343	mixedev->SetHeader(dummy);
344	if(gErrorIgnoreLevel<=kWarning){
345	cout << "Read event: " << nEvent << " " << nEventHijing << endl;
346	mixedev->Print();
347	}
348
349	tree->Fill();
350
351	nEvent++;
352	nEventHijingCounter++;
353	ev->Reset();
354	mixedev->Reset();
355	} //end of nev loop
356
357	file->Write();
358	file->Close();
359	delete file;
360	delete ev;
361	delete theTree;
362	delete backtheTree;
363	delete backev;
364	}
365
366	void testEvents(Char_t *sdir)
367	{
368	Char_t filename[1000];
369	sprintf(filename,"%s/aliev-type0.root",sdir);
370	readEvents(filename,0);
371	sprintf(filename,"%s/aliev-type1.root",sdir);
372	readEvents(filename,0);
373	sprintf(filename,"%s/aliev-type2.root",sdir);
374	readEvents(filename,0);
375	sprintf(filename,"%s/aliev-type3.root",sdir);
376	readEvents(filename,0);
377	sprintf(filename,"%s/aliev-type10.root",sdir);
378	readEvents(filename,0);
379	sprintf(filename,"%s/aliev-type11.root",sdir);
380	readEvents(filename,0);
381	sprintf(filename,"%s/aliev-type12.root",sdir);
382	}
383
384	void readEvents(Char_t *filename,Int_t nMaxEvents)
385	{
386	//connect to events
387	TChain *theTree = new TChain("AJEPtree");
388	theTree->Add(filename);
389	AliJetEventParticles *ev=new AliJetEventParticles();
390	theTree->SetBranchAddress("particles",&ev);
391
392	Int_t treeentries=(Int_t)theTree->GetEntries();
393	if((nMaxEvents<0) \|\| (nMaxEvents>treeentries))
394	nMaxEvents=treeentries;
395
396	cout << "Found " << treeentries << ", but use " << nMaxEvents << " in " << filename << endl;
397
398	//=========================================================================
399	// start the event loop
400	//=========================================================================
401	Int_t nEvent = 0;
402	while(nEvent<nMaxEvents){
403	if ((nEvent % 100) == 0) {
404	cout << "Reading event " << nEvent << endl;
405	}
406
407	//connect the event
408	theTree->GetEvent(nEvent);
409
410	ev->Print();
411	#if 0
412	const TClonesArray *p=ev->GetParticles();
413	for(Int_t i=0;i<p->GetEntriesFast();i++) {
414	cout << i << endl;
415	((AliJetParticle*)p->At(i))->Print();
416	}
417	#endif
418	ev->Clear();
419	nEvent++;
420	} //end of nev loop
421
422	delete ev;
423	delete theTree;
424	}
425
426	void displayMixedEvents(Char_t fname,Char_t hname,Float_t ptcut,Int_t nMaxEvents)
427	{
428	//connect to events
429	TChain *theTree = new TChain("AJEPtree");
430	theTree->Add(fname);
431	AliJetEventParticles *ev=new AliJetEventParticles();
432	theTree->SetBranchAddress("particles",&ev);
433
434	Int_t treeentries=(Int_t)theTree->GetEntries();
435	if((nMaxEvents<0) \|\| (nMaxEvents>treeentries))
436	nMaxEvents=treeentries;
437	//cout << "Found " << nMaxEvents << " in " << fname << endl;
438
439	TChain *backtheTree = new TChain("AJEPtree");
440	backtheTree->Add(hname);
441	AliJetEventParticles *backev=new AliJetEventParticles();
442	backtheTree->SetBranchAddress("particles",&backev);
443
444	Int_t backtreeentries=(Int_t)backtheTree->GetEntries();
445	Int_t nPerBackground=nMaxEvents/backtreeentries;
446	if(nPerBackground==0) nPerBackground=1;
447
448	AliJetEventParticles *mixedev=new AliJetEventParticles(0);
449
450	Int_t neta=20;
451	Int_t nphi=62;
452
453	TCanvas *c1=new TCanvas("c1","",1000,500);
454	c1->Divide(2,1);
455	c1->Update();
456	gStyle->SetLabelSize(0.03,"XYZ");
457	gStyle->SetTitleOffset(1.0,"XYZ");
458	gStyle->SetTitleOffset(1.3,"Y");
459	gROOT->ForceStyle();
460
461	TH2F *h2pp=new TH2F("h2pp","E_{T}^{ch} in PP [GeV];#phi;#eta",nphi,0,TMath::TwoPi(),neta,-1,1);
462	h2pp->SetStats(0);
463	TH2F *h2pb=new TH2F("h2pb","E_{T}^{ch} in PbPb [GeV];#phi;#eta;",nphi,0,TMath::TwoPi(),neta,-1,1);
464	h2pb->SetStats(0);
465
466	//=========================================================================
467	// start the event loop
468	//=========================================================================
469	Int_t nEvent = 0;
470	Int_t nEventHijing = -1;
471	Int_t nEventHijingCounter = nPerBackground;
472	while(nEvent<nMaxEvents){
473
474	//get the event
475	theTree->GetEvent(nEvent);
476
477	//load background if needed
478	if(nEventHijingCounter==nPerBackground){
479	backev->Reset();
480	nEventHijing++;
481	backtheTree->GetEvent(nEventHijing);
482	if(nEventHijing==backtreeentries) nEventHijing=0;
483	nEventHijingCounter=0;
484	}
485	backev->AddSignal(*ev);
486	TString dummy="Counter: ";
487	dummy+=nEvent;
488	dummy+=" ";
489	//dummy+=ev->GetHeader();
490	dummy+="(Pythia ";dummy+=nEvent;
491	dummy+=" ";
492	//dummy+=backev->GetHeader();
493	dummy+=", Hijing ";dummy+=nEventHijing;
494	dummy+=")";
495	mixedev->Set(*backev);
496	mixedev->SetHeader(dummy);
497	if(gErrorIgnoreLevel<=kWarning){
498	cout << "Read event: " << nEvent << " " << nEventHijing << endl;
499	mixedev->Print();
500	}
501
502	// loop over all particles...
503	c1->cd(1);
504	h2pp->Reset();
505	AliJetParticle *aliparticle = NULL;
506	TIterator *iter = ev->GetParticles()->MakeIterator();
507	while ((aliparticle = (AliJetParticle *) iter->Next()) != NULL) {
508	Float_t pt=aliparticle->Pt();
509	if(pt<ptcut) continue;
510	h2pp->Fill(aliparticle->Phi(),aliparticle->Eta(),pt);
511	}
512	delete iter;
513	h2pp->Draw("Lego");
514	c1->cd(2);
515	h2pb->Reset();
516	iter = mixedev->GetParticles()->MakeIterator();
517	while ((aliparticle = (AliJetParticle *) iter->Next()) != NULL) {
518	Float_t pt=aliparticle->Pt();
519	if(pt<ptcut) continue;
520	h2pb->Fill(aliparticle->Phi(),aliparticle->Eta(),pt);
521	}
522	delete iter;
523
524	h2pb->Draw("Lego");
525	c1->Update();
526	Char_t c;cin >> c;
527	if(c=='s'){
528	Char_t sfilenamehist[1024];
529	sprintf(sfilenamehist,"./display-planehists-%d-%d.eps",nEvent,nEventHijingCounter);
530	c1->SaveAs(sfilenamehist);
531	}
532	else if(c=='q') break;
533
534	nEvent++;
535	nEventHijingCounter++;
536	ev->Reset();
537	mixedev->Reset();
538	} //end of nev loop
539
540
541	delete ev;
542	delete theTree;
543	delete backtheTree;
544	delete backev;
545	}