// 15-nov-2002 -> /auto/alice/rcabrera/jets/test
//.........................................................................
// Author:  Renan Cabrera (Creighton U.)
// 8/22/2002
// 11/13/2002 19:00 
// 11/14/2002 12:00  Now we save all events. If we do not find a jet in this event
//                   then we save only partons.
// 11/14/2002 20:00  FillFromTracks in background (1,0)                    
//........................................................................

///////////////////////////////////////////////////////////////////////////
//
//  This macro "makeJMDST.C"  must be put in the directory where we
//  want to produce the output MDST files.
//  All the parameters must be defined in the following functions:
//     InitAliEMCALJetFinder  [Parameters of the JetFinder]
//     DefineInput:
//    		 InDir:    Input Directory 
//   		 OutFile:  JetMDST output filename
//               MaxEvent: Maximum number of events in each galice.root input file 
//    		 file1:    firs subdirectory input in the structure.
//    		 file2:	   last subdirectory input in the structure. 
//               EventType: Expected "Pythia" or "Hijing"
// 
// The idea of the JetMDST is to summarize all important information needed for
// direct analysis. Each event generates the information of:
//     -Partons: in the AliEMCALParton object that includes Et, Eta, Phi
//     -Number of Jets        
//     -Jets:   in the AliEMCALJet object that includes
//		Reconstructed Et, Phi, Eta of the Jet.
//		Et, Phi, Eta of the tracks in the jet cone.
//
// About the background mixing, this mixing can be carried out for the case of 
// EventType = Pythia. In this situation the JMDST will contain all Pythia events
// mixed with a Hijing background defined in the variable BackgroundFile   
// Note.- Please do not add more histograms to this macro. The right place to do that is
//        in "readJMDST.C"
//////////////////////////////////////////////////////////////////////////////

#include "TROOT.h"
#include "TFile.h" 
#include "TObjArray.h"
#include "TTree.h"
#include "TBranch.h"
#include "TParticle.h"
#include "TObject.h"

TH2F* hLego = 0   ;
TH2F* hLegoB  = 0 ;
TH2F* hLegoEMCAL=0;

TH1F* hJetEt,*hJetPhi,*hJetEta,*hPartonEt,*hPartonPhi,*hPartonEta;
TH1F* hTracksEta,*hNJets;
TH2F* hEMCALCell;
TH1F* hNOpenedEvents;

AliGenEventHeader* evHeader;
    // JetFinder Parameters
    Int_t debugJet;
    Float_t coneRadius ;
    Float_t etSeed    ;
    Float_t minJetEt  ;
    Float_t minCellEt ; 
    Float_t minPtCut  ; 
    Int_t modeSubtraction;
    Float_t minMove ;
    Float_t maxMove ;
    Float_t precBg  ;
    Bool_t hadronCorrection;

TClonesArray * Tjets = new TClonesArray("AliEMCALJet", 10);
TClonesArray * Tpartons = new TClonesArray("AliEMCALParton",4);
Int_t debug;

void makeJMDST(Int_t deb = 0,Int_t Bg = 1)
//................................................................
// deb: argument that defines the comments.
// Bg = 1 to mix background (for TypeEvent = pythia )
// Bg = 0 to not mix background.
// The "DefineInput" function defines all needed parameters.
//...............................................................
{
debug = deb;
TString InDir,nFileIn,OutFile,EventType;
Int_t file1,file2;
Int_t MaxEvent;
DefineInput(InDir,OutFile,EventType,file1,file2,MaxEvent);

TFile JetMicroDst(OutFile.Data(), "RECREATE") ; // this is the micro DST file
TTree * mdst = new TTree("mdst","jets dst") ; // this is the micro dst

InitHistograms();

AliEMCALGetter *gime = 0;

AliEMCALJetFinder* jetFinder = 0;

Int_t NJets,NPartons;      //Number of jets and partons
Int_t EventN;      //Event Number
Int_t NParticles;  //Number of particles in event 
Int_t NOpenedEvents=0; // Number of successfully opened events

TBranch * partonsbranch   = mdst->Branch("partons",&Tpartons,32000,1) ;
TBranch * jetsbranch      = mdst->Branch("jets",&Tjets,64000,1);


TBranch * NJetsbranch     = mdst->Branch("njets", &NJets, "njets/I" );
TBranch * EventNbranch    = mdst->Branch("eventn",&EventN,"eventn/I");
TBranch * NParticlesbranch= mdst->Branch("nparticles",&NParticles,"nparticles/I");


//::::::::::::::::::::::::::::::Background::::::::::::::::::::::::::
// Note.-   The file of the backaground can be changed
//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
if( Bg && EventType.Contains("Pythia") ){
  // TString BackgroundFile("/auto/alice/AliENProd/V30703/00070/00010/galice.root");

   TString BackgroundFile("/auto/alice/sahal/PPR3/3/1/galice.root");

   gime = 0 ;
   gime = AliEMCALGetter::GetInstance(BackgroundFile.Data());
   
        jetFinder =  new AliEMCALJetFinder("UA1 Jet Finder", "Pre Production");
        InitAliEMCALJetFinder(jetFinder);

   Int_t nparticles = gAlice->GetEvent(0);
   jetFinder->FillFromHits(0);
   jetFinder->FillFromTracks(1,0);
   jetFinder->SaveBackgroundEvent();

}
//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

for(Int_t ifile=file1; ifile<=file2; ifile ++){ // loop over files

   nFileIn = InDir;

 //**************************************************************************************\\ 

 if (ifile < 10 )
   {
     nFileIn += "0000"; //<- This line must be modified according to the structure of the
   } else if ( ifile > 99  )
   {
      nFileIn += "00"; //<- This line must be modified according to the structure of the
                     //   input directory
   }else
   {
      nFileIn += "000"; //<- This line must be modified according to the structure of the
   }

  //***************************************************************************************//


   nFileIn += ifile;
   nFileIn += "/galice.root";

   gime=0;
   gime = AliEMCALGetter::GetInstance( nFileIn.Data() );
   if(gime==0) 
       {printf(" Error readin file %s (skiping this file!)\n ", nFileIn.Data() ); continue;}
   printf("InputFile = %s found \n", nFileIn.Data() );

 
   if( !jetFinder ){
        if(debug>0)printf("Allocating jetFinder\n"); 
	jetFinder =  new AliEMCALJetFinder("UA1 Jet Finder", "Pre Production");
 	InitAliEMCALJetFinder(jetFinder);}


        if (debug>0) printf("===============Beginning to run over events\n");

   for (Int_t nev = 0 ; nev < MaxEvent ; nev++) { // loop over events
 	
         printf("\n Number of file=%d   Number of event = %d \n",ifile,nev);

 	 jetFinder->ResetJets();
	 Int_t nparticles = gAlice->GetEvent(nev);
         if(debug>1)printf("Number of Particles = %d\n",nparticles);

	 if (nparticles <= 0) {
              printf(" nparticles<=0, breaking event loop\n");
              break;}
              else {
              if(debug==2) printf("Event successfully opened\n ");
              NOpenedEvents++;             
		} 
         if(Bg==1) jetFinder->InitFromBackground();

         if(debug>0) printf("Beginning to fill information to find jets\n");
	
         NParticles = nparticles; //To fill Tree	
	 jetFinder->FillFromHits(1);      if(debug>0)printf("After FillFromHits\n");
	 jetFinder->FillFromTracks(1, 0); if(debug>0)printf("After FillFromTracks\n");
	 jetFinder->Find();if(debug>0)printf("After jetFinder\n");
        
	 NJets = jetFinder->Njets();

	 if(NJets==0) { FillPartonsToMDST(EventType,NPartons);hNJets->Fill(0); } 
		else{
                if(NJets>10){printf("Number of Jets=%d ,reducing to 10",NJets);NJets = 10;}
       
         	 hNJets->Fill(NJets);

        	 FillJetsToMDST(jetFinder,NJets);if(debug>0)printf("\n After FillJets\n");

        	 FillPartonsToMDST(EventType,NPartons);//printf("After FillPartons----\n");
                    }

	  EventN = nev; //Filling the event number
          mdst->Fill();
          if(NJets>0) Tjets->Delete();

   }// loop over events

}//loop over files

hNOpenedEvents->Fill(NOpenedEvents);

mdst->Print();
JetMicroDst->Write();
JetMicroDst->Close();

TCanvas *c1 = new TCanvas("c1","hLego",100,100,600,500);
c1->cd();//c1->SetLogz();
hLego->Draw("lego");        
c1->SaveAs("hLego.gif");

}

//=======================================================

void  FillPartonsToMDST(TString EventType,Int_t &NPartons )
//
//
// 
{

 Int_t npart = (gAlice->GetHeader())->GetNprimary();
     if(debug>1) printf("\n AliEMCALJetFinder::FillFromPartons : npart %i\n", npart);

 evHeader = gAlice->GetHeader()->GenEventHeader();

 TString EventType2(evHeader->ClassName());

if(EventType.Contains("Hijing")&&EventType2.Contains("Hijing")){
        if(debug>1)printf("Filling Hijing Partons\n ");
 	TLorentzVector parton[4];
        NPartons = 4;
 ((AliGenHijingEventHeader)evHeader)->GetJets(parton[0], parton[1], parton[2], parton[3]);

	for(Int_t i=0;i<4;i++){
          if(debug>1)printf(" Hijing Parton %d\n",i);
          if( parton[i].Pt()==0  )
          {printf("Error Pt=0,not filling partons and returning");  return;  }

          new((*Tpartons)[i])  AliEMCALParton( parton[i].Pt(),parton[i].Phi(),parton[i].Eta() );

	  hPartonEt->Fill(  parton[i].Pt() );
          hPartonPhi->Fill( parton[i].Phi()*180/TMath::Pi()  );
          hPartonEta->Fill( parton[i].Eta() );
         
                     } 
	 }

else if( EventType.Contains("Pythia")&& EventType2.Contains("Pythia") ){
        if(debug>1) printf("Filling Pythia Partons\n");
        NPartons = 2;
	TParticle *MPart;
        MPart = gAlice->Particle(6);
        if(debug>1) printf(" Pythia Parton 1\n");
        new((*Tpartons)[0])  AliEMCALParton( MPart->Pt(),MPart->Phi(),MPart->Eta() );

	hPartonEt->Fill(  MPart->Pt() );
        hPartonPhi->Fill( MPart->Phi()*180/TMath::Pi() );
        hPartonEta->Fill( MPart->Eta() );

	MPart = gAlice->Particle(7);
        printf(" Pythia Parton 2\n");
	new((*Tpartons)[1])  AliEMCALParton( MPart->Pt(),MPart->Phi(),MPart->Eta() );

	hPartonEt->Fill(  MPart->Pt() );
        hPartonPhi->Fill( MPart->Phi()*180/TMath::Pi()  );
        hPartonEta->Fill( MPart->Eta() );

	}

else	printf("ERROR. Incompatibility of the Header name (Pythia-Hijing)!!\n");


}
//..........................................................

void FillJetsToMDST(AliEMCALJetFinder* jetFinder,Int_t NJets  )
{
  AliEMCALJet *jet;
  Float_t PtT[50],EtaT[50],PhiT[50]; 
  Int_t PdgT[50];

  for( Int_t i=0;i<NJets;i++) {
       //Int_t i=0;
       jet = jetFinder->GetJetT(i);

       new((*Tjets)[i]) AliEMCALJet( *jet );           

	   hJetEt->Fill(  jet->Energy() );
           hJetPhi->Fill( jet->Phi()*180/TMath::Pi() );
	   hJetEta->Fill( jet->Eta() );
	
	  hLego  = new TH2F( *(jetFinder->GetLego()) );

	         Int_t NTracks = jet->NTracks();
		 jet->TrackList(PtT, EtaT, PhiT, PdgT);
		
       	 for(Int_t iT=0; iT< NTracks; iT++ )
                	{
               		 hTracksEta->Fill( EtaT[iT] );
                	}
	 //if(jet->Energy()<40. ) exit(0);
 		
    } 
}

//______________________________________________________________________

InitHistograms()
{

 hJetEt =  new TH1F("hJetEt","Transverse Energy",100,0,100);
 hJetEta = new TH1F("hJetEta","Eta",100,-0.7,0.7);
 hJetPhi = new TH1F("hJetPhi","Phi",120,0,120);

 hPartonEt =  new TH1F("hPartonEt","Transverse Energy",100,0,100);
 hPartonEta = new TH1F("hPartonEta","Eta",100,-0.7,0.7);
 hPartonPhi = new TH1F("hPartonPhi","Phi",120,0,120);

 hTracksEta = new TH1F("hTracksEta","Tracks's Eta",100,-0.7,0.7);

 hNJets =  new TH1F("hNJets","Number of Jets Distribution",10,0,10);

 hNOpenedEvents = new TH1F("hNOpenedEvents","NUmber of opened Events",100000,1,100000);

}

//______________________________________________________________________
void InitAliEMCALJetFinder(AliEMCALJetFinder* jetFinder)
{
     // Parameters for the Jet Finder	
     debugJet   = 1;
     coneRadius = 0.3;
     etSeed     = 5.;
     minJetEt   = 10.;  //
     minCellEt  = 2.;   // be carefull
     minPtCut   = 0;    //2;    // for including charge particles to file
     modeSubtraction = 1;
     minMove = 0.05;
     maxMove = 0.15;
     precBg  = 0.035;
     Bool_t hadronCorrection = kTRUE;

    jetFinder->Init();
    jetFinder->SetDebug(debugJet);
    jetFinder->SetConeRadius(coneRadius);
    jetFinder->SetEtSeed(etSeed);
    jetFinder->SetMinJetEt(minJetEt);
    jetFinder->SetMinCellEt(minCellEt);
    jetFinder->SetPtCut(minPtCut);
    jetFinder->SetHadronCorrection(0);
    if( hadronCorrection ) {
      jetFinder->SetHadronCorrection(1);
      jetFinder->SetHadronCorrector(AliEMCALHadronCorrectionv0::Instance());
    }
    jetFinder->SetSamplingFraction(12.9);
    jetFinder->SetParametersForBgSubtraction(modeSubtraction,minMove,maxMove,precBg);

    jetFinder->SetMomentumSmearing(1);
    jetFinder->SetEfficiencySim(1);
}


//_____________________________________________________________________
void DefineInput(TString& InDir, TString& OutFile,
TString& EventType,Int_t &file1,Int_t &file2,Int_t &MaxEvent )
//
// Note.- Some of the input directories are just soft links
//
{

    /*InDir = "/auto/alice/AliENProd/V30703/00054/";
    OutFile = "JetMDSTPythia50Bg.root";
    EventType = "Pythia";
    MaxEvent = 100;
    file1 = 1;
    file2 = 40;
    break;*/

    EventType = "Pythia";
    InDir = "/auto/alice/AliENProd/V30703/00062/";
    OutFile = "JetMDSTPythia50Bg.root";

    InDir     = "/auto/alice/AliENProd/V30703/00062/";
    OutFile   = "/auto/alice/pavlinov/jet/microDST/"; // my output dir
    OutFile  += "JetMDSTPythia50Bg.root";
    //
    MaxEvent = 100;
    file1 = 1;
    file2 = 10;
    break;

   /* InDir = "/auto/alice/AliENProd/V30703/00063/";
    OutFile = "JetMDSTPythia75Bg.root";
    EventType = "Pythia";
    MaxEvent = 100;
    file1 = 1;
    file2 = 40;
    break;*/

   /* InDir = "/auto/alice/AliENProd/V30703/00064/";
    OutFile = "JetMDSTPythia100Bg.root";
    EventType = "Pythia";
    MaxEvent = 100;
    file1 = 1;
    file2 = 5;
    break;*/


/*    InDir = "/auto/alice/proposal/AliENProd/2002-06/V3.08.Rev.03/00077/";
    OutFile = "JetMDSTHijing75.root";
    EventType = "Hijing";
    MaxEvent = 100;
    file1 = 91;
    file2 = 91;
    break;   */


/*  InDir = "/auto/alice/AliEn_Prod/2002-06/V3.08.Rev.03/00084";
    OutFile = "JetMDSTHijing200.root";
    EventType = "Hijing";
    MaxEvent = 100;
    file1 = 62;
    file2 = 62;
    break;  */


}