JETAN/jetan2004/createEvents.C

   1 // $Id$
   2
   3 #ifndef __CINT__
   4 #include <stdio.h>
   5 #include <Riostream.h>
   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, kHMPIDpid=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 }