/* returns number of pions and kaons from all the annihilations */
/* together with corresponding energy spectrums                 */
int AnalyseParticles(const char* fileName = "man_output.root", int verbose = 0){
  //---------------- Open File ------->
  // output files
  // TFile* fOut = new TFile("particles.root","RECREATE");
  TFile* fIn = new TFile(fileName);
  if(!fIn->IsOpen()){
    std::cout<<"Error opening file"<<std::endl;
    return -1;
  }
  //---------------------------------<


  //------- Prepare histograms ------->
  TH1I* hParticles = new TH1I("hParticles", "Particles;Name;Quantity", 6, 0, 6);
  // change labels of values on Xaxis of the histogram
  const char* particleNames[] = {"pi0", "pi+", "pi-", "K0", "K+", "K-"};
  for(int iBin=0; iBin<6;iBin++){
    hParticles->GetXaxis()->SetBinLabel(iBin+1,particleNames[iBin]);
  }
  
  TH1D* hPi0 = new TH1D("hPi0", "Energy spectrum for #pi^{0};Energy [eV];Entries", 100, 0, 1000);
  TH1D* hPip = new TH1D("hPip", "Energy spectrum for #pi^{+};Energy [eV];Entries", 100, 0, 1000);
  TH1D* hPim = new TH1D("hPim", "Energy spectrum for #pi^{-};Energy [eV];Entries", 100, 0, 1000);
  TH1D* hK0 = new TH1D("hK0", "Energy spectrum for K^{0};Energy [eV];Entries", 100, 0, 1000);
  TH1D* hKp = new TH1D("hKp", "Energy spectrum for K^{+};Energy [eV];Entries", 100, 0, 1000);
  TH1D* hKm = new TH1D("hKm", "Energy spectrum for K^{-};Energy [eV];Entries", 100, 0, 1000);  
  //---------------------------------<
  
  //---------------- Ntuples --------->
  TBranch* row_wise_branch;
  TLeaf* leaf;
  TTree* tCounter; // pointer to the Ntuple with number of secondaries in each event
  TTree* tParticles; // pointer for the Ntuple with particles (Ntuple is a special form of a Tree)
  int nParticles;
  // values storred in the Ntuples
  char name[200];
  int pdg;
  double mass;
  double charge;
  double lifeTime;
  int leptonNumber;
  int baryonNumber;
  double energy;
  int nSecondaries;

  tParticles = (TTree*)fIn->GetEntries("fSecondaries")->Clone("fParticles"); // getting Ntuple from the file
  nParticles = tParticles->GetEntries();
  // connecting each value in the Ntuple with specific variable in the script
  // if(tParticles->SetBranchAddress("name", name)){
    // if(verbose != 0) std::cout<<"Weird stucture of the Tree, going Leaf, by leaf instead" << std::endl;
     row_wise_branch = (TBranch*)tParticles->GetBranch("row_wise_branch");
     leaf = (TLeaf*)row_wise_branch->GetLeaf("Name");
     leaf->SetAddress(name);
//  }
  // if(tParticles->SetBranchAddress("PDGEncoding", &pdg)){
    leaf = (TLeaf*)row_wise_branch->GetLeaf("PDGEncoding");
    leaf->SetAddress(&pdg);
  // }
// if(tParticles->SetBranchAddress("PDGMass", &mass)){
  leaf = (TLeaf*)row_wise_branch->GetLeaf("PDGMass");
  leaf->SetAddress(&mass);
// }
// if(tParticles->SetBranchAddress("PDGCharge", &charge)){
     leaf = (TLeaf*)row_wise_branch->GetLeaf("PDGCharge");
     leaf->SetAddress(&charge);
  // }
// if(tParticles->SetBranchAddress("PDGLifeTime", &lifeTime)){
  // leaf = (TLeaf*)row_wise_branch->GetLeaf("PDGLifeTime");
  // leaf->SetAddress(&lifeTime);
  // }
//if(tParticles->SetBranchAddress("LeptonNumber", &leptonNumber)){
  // leaf = (TLeaf*)row_wise_branch->GetLeaf("LeptonNumber");
  // leaf->SetAddress(&leptonNumber);
  // }
// if(tParticles->SetBranchAddress("BaryonNumber", &baryonNumber)){
  // leaf = (TLeaf*)row_wise_branch->GetLeaf("BaryonNumber");
  // leaf->SetAddress(&baryonNumber);
  // }
//if(tParticles->SetBranchAddress("Energy", &energy)){
  // leaf = (TLeaf*)row_wise_branch->GetLeaf("Energy");
  // leaf->SetAddress(&energy);
  // }

  tCounter = (TTree*)fIn->Get("fCounter")->Clone("tCounter"); // getting Ntuple from the file
//if(tCounter->SetBranchAddress("fNSecondaries", &nSecondaries)){
    row_wise_branch = (TBranch*)tCounter->GetBranch("row_wise_branch");
    leaf = (TLeaf*)row_wise_branch->GetLeaf("fNSecondaries");
    leaf->SetAddress(&nSecondaries);
    //}
  //---------------------------------<
       
  int nb_PimPimPip = 0;     //Set event typologies
  int nb_KmPipPim  = 0;
  int nb_KmKpPim   = 0;
  int nb_KpPipPim  = 0;
  int nb_KpKpPim   = 0;
  int nb_PipPipPim    = 0;

  int nb_KpKpKmKmPimPip =0;
  int nb_KpKpPimKmKOL =0;
  int nb_KpKpPimKOLPim=0;
  int nb_KpKpPimKOLKOL=0;
  int nb_KpKpPimKOSKm=0;

  int nb_KpKmPipPim=0; 

 
   


  
  int all_Km = 0;
  int all_Kp = 0;
  int all_Pim =0;
  int all_Pip=0; 

  int  nb_em=0;
  int  nb_gamma=0;
  int  nb_pim=0;
  int  nb_piz=0;
  int  nb_pip=0;
  int  nb_km=0;
  int  nb_kzl=0;
  int  nb_kzs=0;
  int  nb_kp=0;
  int  nb_eta=0;


  //-- loop over all annihilations -->
  int particle_offset = 0;
  for(int event=0; event<tCounter->GetEntries(); event++){
    tCounter->GetEntry(event); // getting event (annihilation)
    if(verbose != 0) printf("Annihilation #%d (%d)\n", event, nSecondaries);
         
          int nb_Pim =0;              //Set particles
          int nb_Pip =0;
          int nb_Km  =0;
          int nb_Kp  =0;
          int nb_KOS =0;
          int nb_KOL =0;
   
    for(int i=0; i<nSecondaries; i++){ // loop over all particles
      tParticles->GetEntry(particle_offset + i);
      if(verbose != 0){
	printf("------------------------------\n");
	printf("Particle #%d\n",i);
	printf("\tName: %s\n",name);
	printf("\tPDG: %d\n",pdg);
	printf("\tMass: %f\n",mass);
	printf("\tCharge: %f\n",charge);
	printf("\tLife Time: %f\n",lifeTime);
	printf("\tLepton number: %d\n",leptonNumber);
	printf("\tBaryon number: %d\n",baryonNumber);
	printf("\tEnergy [eV]: %f\n\n",energy);
      }
                        //Set particles
          
      // std::cout<<pdg<<" "<<name<<" "<<mass<<std::endl;
      //std::cout<<name<<std::endl;
          
 
    
       	   
 
       // if(strcmp(name,"pi-")==1) nb_Pim=nb_Pim+1;
       // if(strcmp(name,"pi+")==1) nb_Pip=nb_Pip+1;
       // if(strcmp(name,"kaon-")==1) nb_Km=nb_Km+1;
       // if(strcmp(name,"kaon+")==1) nb_Kp=nb_Kp+1;
          
          

       if(pdg==-211){
	nb_Pim=nb_Pim+1;
      //	std::cout<<"Pi-"<<std::endl;
       }
       if(pdg== 211) nb_Pip=nb_Pip+1;
       if(pdg== -321){
	 nb_Km=nb_Km+1;
	 // std::cout<<"K-"<<std::endl;
       }
       if(pdg==  321) nb_Kp=nb_Kp+1;
       if(pdg==  310) nb_KOS=nb_KOS+1;
       if(pdg==  130) nb_KOL=nb_KOL+1;

    
      if(pdg == 111){ // check for pi0
	hParticles->Fill(0);
	hPi0->Fill(energy);
      }else if(abs(pdg) == 211){ // check for charged pions (pi+ is numbered 211 while pi- is -211)
	if(charge > 0){
	  hParticles->Fill(1);
	  hPip->Fill(energy);
	}else{
	  hParticles->Fill(2);
	  hPim->Fill(energy);
	}
      } 
    } // end of the particles loop
    // particle_offset += nSecondaries;
      // if(nb_Pim==2&&nb_Pip==1) nb_PimPimPip=nb_PimPimPip+1;
    all_Kp += nb_Kp;
    all_Km += nb_Km;
    all_Pip+= nb_Pip;
    all_Pim+= nb_Pim;
       if(nb_Pim==1&& nb_Pip==1&&nb_Km==1) nb_KmPipPim=nb_KmPipPim+1;
       if(nb_Km==1&&nb_Kp==1&&nb_Pim==1) nb_KmKpPim=nb_KmKpPim+1;
       if(nb_Pim==2&&nb_Pip==1) nb_PimPimPip=nb_PimPimPip+1;
       if(nb_Kp==1&&nb_Pip==1&&nb_Pim==1) nb_KpPipPim=nb_KpPipPim+1;
       if(nb_Kp==1&&nb_Pip==1&&nb_Pim==1) nb_KpPipPim=nb_KpPipPim+1;

            

       if(nb_Kp==2&&nb_Pim==1){
            nb_KpKpPim=nb_KpKpPim+1;   
           // for(int i=0; i<nSecondaries; i++){
	     // tParticles->GetEntry(particle_offset + i);
             //      if(pdg==11) nb_em=nb_em+1;
             //      if(pdg==22) nb_gamma=nb_gamma+2;
             //      if(pdg==-211) nb_pim=nb_pim+1;
             //      if(pdg== 111) nb_piz=nb_piz+1;
             //      if(pdg==211) nb_pip=nb_pip+1;
             //      if(pdg==-321) nb_km=nb_km+1;
             //      if(pdg==130)  nb_kzl=nb_kzl+1;
             //      if(pdg==310) nb_kzs=nb_kzs+1;
             //      if(pdg==321) nb_kp=nb_kp+1;
             //      if(pdg==221) nb_eta=nb_eta+1;
                   
           
      // std::cout<<name<<std::endl;


             // if(pdg==-211) std::cout<<"Pi-!   " <<std::endl; 
            // }
            }
               // str Plist[1];
      
       if(nb_Pip==2&&nb_Pim==1) nb_PipPipPim=nb_PipPipPim+1;        
      
       if(nb_Kp==2&&nb_Km==2&&nb_Pim==1&&nb_Pip==1){ nb_KpKpKmKmPimPip=nb_KpKpKmKmPimPip+1;
       //  for(int i=0; i<nSecondaries; i++){
       //       tParticles->GetEntry(particle_offset + i);
        //          std::cout<<"K+K+K-K-Pi-"<<  name <<std::endl;
       //             }
                    }
       if(nb_Kp==2&&nb_Pim==1&&nb_Km==1&&nb_KOL==1){ nb_KpKpPimKmKOL=nb_KpKpPimKmKOL+1;
              for(int i=0; i<nSecondaries; i++){
              tParticles->GetEntry(particle_offset + i);
                std::cout<<"K+K+Pi-K0L"<<  name <<std::endl;
                    }
                    }

       if(nb_Kp==2&&nb_Pim==2&&nb_KOL==1) nb_KpKpPimKOLPim=nb_KpKpPimKOLPim+1;
       if(nb_Kp==2&&nb_Pim==1&&nb_KOL==2) nb_KpKpPimKOLKOL=nb_KpKpPimKOLKOL+1; 
       if(nb_Kp==2&&nb_Pim==1&&nb_KOS==1&&nb_Km==1) nb_KpKpPimKOSKm=nb_KpKpPimKOSKm+1;
       
       if(nb_Kp==1&&nb_Km==1&&nb_Pip==1&&nb_Pim==1) nb_KpKmPipPim=nb_KpKmPipPim+1; 
                particle_offset += nSecondaries;                            
                                      // }
         
  }// end of the event loop
     
      
      
          
        std::cout<<"Pi-Pi-Pi+: "<<std::endl;  std::cout<< nb_PimPimPip<<std::endl;
        std::cout<<"Pi-Pi+K-:"<<std::endl;    std::cout<<nb_KmPipPim<<std::endl;
        std::cout<<"Pi+K-K+ :"<<std::endl;     std::cout<<nb_KmKpPim<<std::endl;
        std::cout<<"K+Pi+Pi- :"<<std::endl;     std::cout<<nb_KpPipPim<<std::endl;
        std::cout<<"K+K+Pi- :"<<std::endl;     std::cout<<nb_KpKpPim<<std::endl;
       // std::cout<<"Particles for the above event typoligi "<<"e-:"<<nb_em<<" gamma:"<<nb_gamma<<" Pi-:"<<nb_pim<<" Pi0:"<<nb_piz<<" Pi+:"<<nb_pip<<" K-:"<<nb_km<<" K0L:"<<nb_kzl<<" K0S:"<<nb_kzs<<" K+:"<<nb_kp<<" eta:"<<nb_eta<< std::endl; 
        std::cout<<"Pi+Pi+Pi- :"<<std::endl;     std::cout<<nb_PipPipPim<<std::endl;
       //Below with two K+
        std::cout<<"K+K+K-K-Pi-Pi+ :"<<std::endl;     std::cout<<nb_KpKpKmKmPimPip<<std::endl;
        std::cout<<"K+K+Pi- K-K0L :"<<std::endl;     std::cout<<nb_KpKpPimKmKOL<<std::endl;
        std::cout<<"K+K+Pi-Pi-K0L :"<<std::endl;     std::cout<<nb_KpKpPimKOLPim<<std::endl;
        std::cout<<"K+K+Pi-K0LK0L :"<<std::endl;     std::cout<<nb_KpKpPimKOLKOL<<std::endl;
         std::cout<<"K+K+Pi-K0S K- :"<<std::endl;     std::cout<<nb_KpKpPimKOSKm<<std::endl;
         
         std::cout<<"K+K-Pi+Pi-"<<std::endl;          std::cout<<nb_KpKmPipPim<<std::endl;
          




	std::cout<<"K-:"<<all_Km<<"K+"<<all_Kp<<"Pi-:"<<all_Pim<<"Pi+ :"<<all_Pip<< std::endl;

  //---------------------------------<
  
  //------- save histograms --------->
//  fOut->cd();
//  hParticles->Write(hParticles->GetName(), TObject::kOverwrite);
//  hPi0->Write(hPi0->GetName(), TObject::kOverwrite);
//  hPip->Write(hPip->GetName(), TObject::kOverwrite);
//  hPim->Write(hPim->GetName(), TObject::kOverwrite);
//  hK0->Write(hK0->GetName(), TObject::kOverwrite);
//  hKp->Write(hKp->GetName(), TObject::kOverwrite);
//  hKm->Write(hKm->GetName(), TObject::kOverwrite);
  //---------------------------------<
  
 // fOut->Close();
  fIn->Close();
  return nParticles;
}