/////////////////////////////////////////////////////////////////////////
//   This macro is a small example of a ROOT macro
//   illustrating how to read the hits in ACORDE,
//   fill some histograms and store it on a file
//
//   Author:
//
//	Eleazar Cuautle <ecuautle@nucleares.unam.mx>
//
//   Also:
//
//	Mario Rodriguez Cahuantzi <mrodrigu@mail.cern.ch>
//             
/////////////////////////////////////////////////////////////////////////

void TestACORDEHits (const char *filename="galice.root",Int_t numberOfEvents=0){
   // delete existing gAlice
  // Dynamically link some shared libs
  if (gClassTable->GetID("AliRun") < 0) 
    {
      gROOT->LoadMacro("loadlibs.C");
      loadlibs();
    }
  if (gAlice)
    {
      //delete gAlice->GetRunLoader();
      delete gAlice;
      gAlice = 0x0;
    }
  
  AliRunLoader *rl = AliRunLoader::Open("galice.root",
					AliConfig::GetDefaultEventFolderName(),"read");
  if (!rl)
    {
      cerr<<"Can't load RunLoader from file! \n";      return 0x0;
    }
  
  rl->LoadgAlice();
  gAlice = rl->GetAliRun();
  if (!gAlice)
    {
      cerr << " AliRun object not found on file \n";     return 0x0; 
    }
  rl->LoadHeader();
  
  // Get the pointer to the ACORDE detector
  AliLoader *acordel = rl->GetLoader("ACORDELoader");
  AliACORDE *ACORDE = (AliACORDE *) gAlice->GetDetector("ACORDE");
  if (ACORDE == 0x0 || acordel == 0x0) {
    cerr << " Can not find ACORDE or ACORDELoader \n";    return 0x0;  
  }
  //Determine how many events would you like to analyze
  if (numberOfEvents==0) numberOfEvents=(Int_t)(rl->GetNumberOfEvents());
  
  //create histograms

  TH1F	*fACORDEEnergy = new TH1F("fACORDEEnergy","Energy distribution of muons generated by AliGenCosmicParam at z = 900 cms (Just above ACORDE); Energy (GeV)",1200,0,360);
  TH1F	*fACORDEEnergyLoss = new TH1F("fACORDEEnergyLoss","Loss energy distribution  of muons generated by AliGenCosmicParam at z = 900 cms (Just above ACORDE);Energy loss (Gev)",100,0,0.030);
  TH1F	*fACORDEAzimuth = new TH1F("fACORDEAzimuth","Azimuth angle distribution of muons generated by AliGenCosmicParam at z = 900 cms (Just above ACORDE); Azimuth angle",360,-180,180);
  TH1F	*fACORDEPolar = new TH1F("fACORDEPolar","Polar angle distribution  of muons generated by AliGenCosmicParam at z = 900 cms (Just above ACORDE); Polar angle",180,0,180);
  TH2F  *fACORDExz = new TH2F("fACORDExz" ,"Distribution XZ  of muons generated by AliGenCosmicParam at z = 900 cms (Just above ACORDE); X (cms); Z (cms) ",900,-850.,850.,1200,-500.,500.); 
  TH2F  *fACORDEAzimPol =  new TH2F("fACORDEAzimPol" ,"Azimuth vs Polar angle distribution ",360,-180.,180.,180,0.,180.);
  TH1F  *fACORDEPx = new TH1F("fACORDEPx" ,"P_{x} distribution  of muons generated by AliGenCosmicParam at z = 900 cms (Just above ACORDE);GeV/c;dN/dP_{x}",60,-120.,120.);
  TH1F  *fACORDEPy = new TH1F("fACORDEPy" ,"P_{y} distribution  of muons generated by AliGenCosmicParam at z = 900 cms (Just above ACORDE);GeV/c;dN/dP_{y}",60,-120.,120.);
  TH1F  *fACORDEPz = new TH1F("fACORDEPz" ,"P_{z} distribution  of muons generated by AliGenCosmicParam at z = 900 cms (Just above ACORDE);GeV/c;dN/dP_{z}",60,-120.,120.);
  TH1F  *fACORDEPt = new TH1F("fACORDEPt" ,"P_{t} distribution  of muons generated by AliGenCosmicParam at z = 900 cms (Just above ACORDE);GeV/c;dN/dP_{t}",1200,0.,360.);
  TH2F  *fELossAzimuth = new TH2F("fELossAzimuth","Dist. Azimuth Angle VS Energy Loss; Azimuth angle; Energy (GeV)",360,0,350,100,0,0.030); 
  TH2F  *fEnergyAzimuth = new TH2F("fEnergyAzimuth","Dist. Azimuth Angle VS Energy; Azimuth angle; Energy (GeV)",360,0,360,300,0,350); 
  TH2F  *fELossPolar = new TH2F("fELossPolar","Dist. Polar Angle VS Energy Loss; Polar angle; Energy (GeV)",180,0,180,100,0,0.030); 
  TH2F  *fEnergyPolar = new TH2F("fEnergyPolar","Dist. Polar Angle VS Energy; Polar angle; Energy (GeV)",180,0,180,300,0,350); 

 //------------------------------------------------------------------//
  
  for (Int_t ievent=0; ievent<numberOfEvents; ievent++) {
    if ((ievent%10) == 0)  printf ("Processing event %d \n", ievent);
    rl->GetEvent(ievent);
    
    // Get the pointer Hit tree
    acordel->LoadHits();
    TTree *hitTree = acordel->TreeH();
    ACORDE->SetTreeAddress();
    if (!hitTree) {
      cout << " No TreeH found" << endl;      return 0x0; //rc;
    }
    
    rl->LoadKinematics();
    Int_t nTrack = (Int_t) hitTree->GetEntries();
    //    cout <<"<AliACORDEanalyzeHits> Found "<< nTrack <<" primary particles with hits \n";
    
    // Start loop on tracks in the hits containers
    for(Int_t iTrack=0; iTrack<nTrack;iTrack++){
      ACORDE->ResetHits();
      hitTree->GetEvent(iTrack);          
      if(ACORDE) {	 
	for(acordeHit=(AliACORDEhit*)ACORDE->FirstHit(-1);acordeHit;
	    acordeHit=(AliACORDEhit*)ACORDE->NextHit()) {
		fACORDEEnergy->Fill(acordeHit->Energy());
		fACORDEEnergyLoss->Fill(acordeHit->Eloss());
		fACORDEAzimuth->Fill(acordeHit->AzimuthAngle());
		fACORDEPolar->Fill(acordeHit->PolarAngle());
		fACORDExz->Fill(acordeHit->X(),acordeHit->Z());
         	fELossAzimuth->Fill((Float_t)(acordeHit->AzimuthAngle()),(Float_t)(acordeHit->Eloss()));
 		fACORDEPx->Fill( (Float_t)(acordeHit->Px()));
	 	fACORDEPy->Fill( (Float_t)(acordeHit->Py()));
	  	fACORDEPz->Fill( (Float_t)(acordeHit->Pz()));
	 	fACORDEPt->Fill( TMath::Sqrt( (acordeHit->Px())*(acordeHit->Px())+
					     (acordeHit->Py())*(acordeHit->Py())) );
	  	fACORDEAzimPol->Fill( (Float_t)(acordeHit->AzimuthAngle()),
			        (Float_t)(acordeHit->PolarAngle()));
		fEnergyAzimuth->Fill(acordeHit->Energy(),acordeHit->AzimuthAngle());
		fEnergyPolar->Fill(acordeHit->Energy(),acordeHit->PolarAngle());
		fELossPolar->Fill(acordeHit->PolarAngle(),acordeHit->Eloss());
	//cout << "Muon triggered by ACORDE "<<endl;
	//cout << "Energy: "<<acordeHit->Energy()<<" EnergyLoss: "<<acordeHit->Eloss()<<" Polar angle: "<<acordeHit->PolarAngle()<<" Azimuth angle: "<<acordeHit->AzimuthAngle()<<endl;

	} // end for acordeHits
      } // end if ACORDE
    } // end for iTrack  
  }

 // Put the histograms in a TCanvas

	TCanvas *c1 = new TCanvas("c1","ACORDE Hits distribution of muons generated by AliGenCosmicParam at z = 900 cms (Just above ACORDE)");
	c1->Divide(2,1);

	c1->cd(1)->SetLogy();
	fACORDEEnergy->SetMarkerStyle(kFullCircle);
	fACORDEEnergy->SetMarkerSize(0.7);
	fACORDEEnergy->Draw("E");

	c1->cd(2)->SetLogy();
	fACORDEEnergyLoss->SetMarkerStyle(kFullCircle);
	fACORDEEnergyLoss->SetMarkerSize(0.7);
	fACORDEEnergyLoss->Draw("E");
	
	c1->Draw();

	TCanvas *c2 = new TCanvas("c2","ACORDE Hits distribution of muons generated by AliGenCosmicParam at z = 900 cms (Just above ACORDE)");
	c2->Divide(2,2);
	c2->cd(1);
	fACORDEAzimuth->Draw();

	c2->cd(2);
	fACORDEPolar->Draw();

	c2->cd(3);
	gStyle->SetPalette(1);
	fACORDExz->DrawCopy("colz");
	
	c2->cd(4);
	gStyle->SetPalette(1);
	fACORDEAzimPol->DrawCopy("colz");

	c2->Draw();

	TCanvas *c3 = new TCanvas("c3","ACORDE Hits distribution of muons generated by AliGenCosmicParam at z = 900 cms (Just above ACORDE)");
	c3->Divide(2,2);

	c3->cd(1)->SetLogy();
	fACORDEPx->SetMarkerStyle(kFullCircle);
	fACORDEPx->SetMarkerSize(0.7);
	fACORDEPx->Draw("E");

	c3->cd(2)->SetLogy();
	fACORDEPy->SetMarkerStyle(kFullCircle);
	fACORDEPy->SetMarkerSize(0.7);
	fACORDEPy->Draw("E");

	c3->cd(3)->SetLogy();
	fACORDEPz->SetMarkerStyle(kFullCircle);
	fACORDEPz->SetMarkerSize(0.7);
	fACORDEPz->Draw("E");

	c3->cd(4)->SetLogy();
	fACORDEPt->SetMarkerStyle(kFullCircle);
	fACORDEPt->SetMarkerSize(0.7);
	fACORDEPt->Draw("E");

	c3->Draw();

	TCanvas *c4 = new TCanvas("c4","ACORDE Hits distribution of muons generated by AliGenCosmicParam at z = 900 cms (Just above ACORDE)");
	c4->Divide(2,2);

	c4->cd(1);
	gStyle->SetPalette(1);
	fEnergyPolar->DrawCopy("colz");

	c4->cd(2);
	gStyle->SetPalette(1);
	fEnergyAzimuth->DrawCopy("colz");

	c4->cd(3);
	gStyle->SetPalette(1);
	fELossPolar->DrawCopy("colz");

	c4->cd(4);
	gStyle->SetPalette(1);
	fELossAzimuth->DrawCopy("colz");
	
	c4->Draw();


 // save histos in a root file
  TFile *fout = new TFile("ACORDE_hits.root","RECREATE");
	fACORDEEnergy->Write();
	fACORDEEnergyLoss->Write();
	fACORDEAzimuth->Write();
	fACORDEPolar->Write();
	fACORDExz->Write();
	fELossAzimuth->Write();
	fACORDEPx->Write();
	fACORDEPy->Write();
	fACORDEPz->Write();
	fACORDEPt->Write();
	fACORDEAzimPol->Write();
	c1->Write();
	c2->Write();
	c3->Write();
	c4->Write();
}