[u/mrichter/AliRoot.git] / STEER / AliTagCreator.cxx

/**************************************************************************
 * Author: Panos Christakoglou.                                           *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

/* $Id$ */

//-----------------------------------------------------------------
//           AliTagCreator class
//   This is the class to deal with the tag creation (post process)
//   Origin: Panos Christakoglou, UOA-CERN, Panos.Christakoglou@cern.ch
//-----------------------------------------------------------------

//ROOT
#include <TFile.h>
#include <TString.h>
#include <TTree.h>
#include <TSystem.h>
#include <TChain.h>
#include <TLorentzVector.h>

//ROOT-AliEn
#include <TGrid.h>
#include <TGridResult.h>

//AliRoot
#include "AliRunTag.h"
#include "AliEventTag.h"
#include "AliESD.h"
#include "AliESDVertex.h"
#include "AliLog.h"


#include "AliTagCreator.h"


ClassImp(AliTagCreator)


//______________________________________________________________________________
AliTagCreator::AliTagCreator() //local mode
{
  //==============Default constructor for a AliTagCreator==================
  fgridpath = "";
  fSE = "";   
  fStorage = 0; 
}

//______________________________________________________________________________
AliTagCreator::~AliTagCreator()
{
//================Default destructor for a AliTagCreator=======================
}

//______________________________________________________________________________
void AliTagCreator::SetStorage(Int_t storage)
{
  // Sets correctly the storage: 0 for local, 1 for GRID
  fStorage = storage;
  if(fStorage == 0)
    AliInfo(Form("Tags will be stored locally...."));
  if(fStorage == 1)
    AliInfo(Form("Tags will be stored in the grid...."));
  if((fStorage != 0)&&(fStorage != 1))
    {
      AliInfo(Form("Storage was not properly set!!!"));
      abort();
    }  
}


//______________________________________________________________________________
Bool_t AliTagCreator::ReadESDCollection(TGridResult *fresult)
{
  // Reads the entry of the TGridResult and creates the tags
  Int_t nEntries = fresult->GetEntries();

  TString alienUrl;
  const char *guid;
  const char *md5;
  const char *turl;
  Long64_t size = -1;

  Int_t counter = 0;
  for(Int_t i = 0; i < nEntries; i++)
    {
      alienUrl = fresult->GetKey(i,"turl");
      guid = fresult->GetKey(i,"guid");
      if(fresult->GetKey(i,"size"))
	size = atol (fresult->GetKey(i,"size"));
      md5 = fresult->GetKey(i,"md5");
      turl = fresult->GetKey(i,"turl");
      if(md5 && !strlen(guid))
	md5 = 0;
      if(guid && !strlen(guid))
	guid = 0;

      TFile *f = TFile::Open(alienUrl,"READ");
      CreateTag(f,guid,md5,turl,size,counter);
      f->Close();
      delete f;	 
      counter += 1;
    }//grid result loop

  return kTRUE;
}

//__________________________________________________________________________
Bool_t AliTagCreator::MergeTags()
{
  //Merges the tags and stores the merged tag file 
  //locally if fStorage=0 or in the grid if fStorage=1
  AliInfo(Form("Merging tags....."));
  TChain *fgChain = new TChain("T");

  if(fStorage == 0) {
    const char * tagPattern = "tag";
    // Open the working directory
    void * dirp = gSystem->OpenDirectory(gSystem->pwd());
    const char * name = 0x0;
    // Add all files matching *pattern* to the chain
    while((name = gSystem->GetDirEntry(dirp))) {
      if (strstr(name,tagPattern))       
	fgChain->Add(name);  
    }//directory loop
    AliInfo(Form("Chained tag files: %d",fgChain->GetEntries()));
  }//local mode

  else if(fStorage == 1) {
    TString alienLocation = gGrid->Pwd();
    alienLocation += fgridpath.Data();
    alienLocation += "/";

    TGridResult *tagresult = gGrid->Query(alienLocation,"*tag.root","","");
    Int_t nEntries = tagresult->GetEntries();
    for(Int_t i = 0; i < nEntries; i++) {
      TString alienUrl = tagresult->GetKey(i,"turl");
      fgChain->Add(alienUrl);
    }//grid result loop      
    AliInfo(Form("Chained tag files: %d",fgChain->GetEntries()));
  }//grid mode
 
  AliRunTag *tag = new AliRunTag;
  AliEventTag *evTag = new AliEventTag;
  fgChain->SetBranchAddress("AliTAG",&tag);
   
  //Defining new tag objects
  AliRunTag *newTag = new AliRunTag();
  TTree ttag("T","A Tree with event tags");
  TBranch * btag = ttag.Branch("AliTAG", &newTag);
  btag->SetCompressionLevel(9);
  for(Int_t iTagFiles = 0; iTagFiles < fgChain->GetEntries(); iTagFiles++) {
    fgChain->GetEntry(iTagFiles);
    newTag->SetRunId(tag->GetRunId());
    const TClonesArray *tagList = tag->GetEventTags();
    for(Int_t j = 0; j < tagList->GetEntries(); j++) {
      evTag = (AliEventTag *) tagList->At(j);
      newTag->AddEventTag(*evTag);
    }
    ttag.Fill();
    newTag->Clear();
  }//tag file loop 
  
  TString localFileName = "Run"; localFileName += tag->GetRunId(); 
  localFileName += ".Merged"; localFileName += ".ESD.tag.root";
     
  TString alienFileName = "/alien";
  alienFileName += gGrid->Pwd();
  alienFileName += fgridpath.Data();
  alienFileName += "/";
  alienFileName +=  localFileName;
  alienFileName += "?se=";
  alienFileName += fSE.Data();

  TString filename = 0x0;
  
  if(fStorage == 0) {
    filename = localFileName.Data();      
    AliInfo(Form("Writing merged tags to local file: %s",filename.Data()));
  } 
  if(fStorage == 1) {
    filename = alienFileName.Data();
    AliInfo(Form("Writing merged tags to grid file: %s",filename.Data()));     
  }
  
  TFile* ftag = TFile::Open(filename, "recreate");
  ftag->cd();
  ttag.Write();
  ftag->Close();

  delete tag;
  delete evTag;
  delete newTag;

  return kTRUE;
}

//_____________________________________________________________________________
void AliTagCreator::CreateTag(TFile* file, const char *guid, const char *md5, const char *turl, Long64_t size, Int_t Counter)
{
  /////////////
  //muon code//
  ////////////
  Double_t MUON_MASS = 0.105658369;
  //Variables
  Double_t X,Y,Z ;
  Double_t ThetaX, ThetaY, Pyz,Chisquare;
  Double_t PxRec,PyRec, PzRec, Energy;
  Int_t Charge;
  TLorentzVector EPvector;

  Float_t ZVertexCut = 10.0; 
  Float_t RhoVertexCut = 2.0; 

  Float_t LowPtCut = 1.0;
  Float_t HighPtCut = 3.0;
  Float_t VeryHighPtCut = 10.0;
  ////////////

  Double_t partFrac[5] = {0.01, 0.01, 0.85, 0.10, 0.05};

  // Creates the tags for all the events in a given ESD file
  Int_t ntrack;
  Int_t nProtons, nKaons, nPions, nMuons, nElectrons;
  Int_t nPos, nNeg, nNeutr;
  Int_t nK0s, nNeutrons, nPi0s, nGamas;
  Int_t nCh1GeV, nCh3GeV, nCh10GeV;
  Int_t nMu1GeV, nMu3GeV, nMu10GeV;
  Int_t nEl1GeV, nEl3GeV, nEl10GeV;
  Float_t maxPt = .0, meanPt = .0, totalP = .0;

  AliRunTag *tag = new AliRunTag();
  AliEventTag *evTag = new AliEventTag();
  TTree ttag("T","A Tree with event tags");
  TBranch * btag = ttag.Branch("AliTAG", &tag);
  btag->SetCompressionLevel(9);
  
  AliInfo(Form("Creating the tags......."));	
  
  Int_t firstEvent = 0,lastEvent = 0;
  TTree *t = (TTree*) file->Get("esdTree");
  TBranch * b = t->GetBranch("ESD");
  AliESD *esd = 0;
  b->SetAddress(&esd);
  
  tag->SetRunId(esd->GetRunNumber());
  
  Int_t iNumberOfEvents = b->GetEntries();
  for (Int_t iEventNumber = 0; iEventNumber < iNumberOfEvents; iEventNumber++) {
    ntrack = 0;
    nPos = 0;
    nNeg = 0;
    nNeutr =0;
    nK0s = 0;
    nNeutrons = 0;
    nPi0s = 0;
    nGamas = 0;
    nProtons = 0;
    nKaons = 0;
    nPions = 0;
    nMuons = 0;
    nElectrons = 0;	  
    nCh1GeV = 0;
    nCh3GeV = 0;
    nCh10GeV = 0;
    nMu1GeV = 0;
    nMu3GeV = 0;
    nMu10GeV = 0;
    nEl1GeV = 0;
    nEl3GeV = 0;
    nEl10GeV = 0;
    maxPt = .0;
    meanPt = .0;
    totalP = .0;
    
    b->GetEntry(iEventNumber);
    const AliESDVertex * vertexIn = esd->GetVertex();
    
    for (Int_t iTrackNumber = 0; iTrackNumber < esd->GetNumberOfTracks(); iTrackNumber++) {
      AliESDtrack * esdTrack = esd->GetTrack(iTrackNumber);
      UInt_t status = esdTrack->GetStatus();
      
      //select only tracks with ITS refit
      if ((status&AliESDtrack::kITSrefit)==0) continue;
      //select only tracks with TPC refit
      if ((status&AliESDtrack::kTPCrefit)==0) continue;
      
      //select only tracks with the "combined PID"
      if ((status&AliESDtrack::kESDpid)==0) continue;
      Double_t p[3];
      esdTrack->GetPxPyPz(p);
      Double_t momentum = sqrt(pow(p[0],2) + pow(p[1],2) + pow(p[2],2));
      Double_t fPt = sqrt(pow(p[0],2) + pow(p[1],2));
      totalP += momentum;
      meanPt += fPt;
      if(fPt > maxPt) maxPt = fPt;
      
      if(esdTrack->GetSign() > 0) {
	nPos++;
	if(fPt > LowPtCut) nCh1GeV++;
	if(fPt > HighPtCut) nCh3GeV++;
	if(fPt > VeryHighPtCut) nCh10GeV++;
      }
      if(esdTrack->GetSign() < 0) {
	nNeg++;
	if(fPt > LowPtCut) nCh1GeV++;
	if(fPt > HighPtCut) nCh3GeV++;
	if(fPt > VeryHighPtCut) nCh10GeV++;
      }
      if(esdTrack->GetSign() == 0) nNeutr++;
      
      //PID
      Double_t prob[5];
      esdTrack->GetESDpid(prob);
      
      Double_t rcc = 0.0;
      for(Int_t i = 0; i < AliPID::kSPECIES; i++) rcc += prob[i]*partFrac[i];
      if(rcc == 0.0) continue;
      //Bayes' formula
      Double_t w[5];
      for(Int_t i = 0; i < AliPID::kSPECIES; i++) w[i] = prob[i]*partFrac[i]/rcc;
      
      //protons
      if ((w[4]>w[3])&&(w[4]>w[2])&&(w[4]>w[1])&&(w[4]>w[0])) nProtons++;
      //kaons
      if ((w[3]>w[4])&&(w[3]>w[2])&&(w[3]>w[1])&&(w[3]>w[0])) nKaons++;
      //pions
      if ((w[2]>w[4])&&(w[2]>w[3])&&(w[2]>w[1])&&(w[2]>w[0])) nPions++; 
      //electrons
      if ((w[0]>w[4])&&(w[0]>w[3])&&(w[0]>w[2])&&(w[0]>w[1])) {
	nElectrons++;
	if(fPt > LowPtCut) nEl1GeV++;
	if(fPt > HighPtCut) nEl3GeV++;
	if(fPt > VeryHighPtCut) nEl10GeV++;
      }	  
      ntrack++;
    }//esd track loop
    
    /////////////
    //muon code//
    ////////////
    Int_t nMuonTracks = esd->GetNumberOfMuonTracks();
    // loop over all reconstructed tracks (also first track of combination)
    for (Int_t iTrack = 0; iTrack <  nMuonTracks;  iTrack++) {
      AliESDMuonTrack* muonTrack = esd->GetMuonTrack(iTrack);
      if (muonTrack == 0x0) continue;
      
      // Coordinates at vertex
      Z = muonTrack->GetZ(); 
      Y = muonTrack->GetBendingCoor();
      X = muonTrack->GetNonBendingCoor(); 
      
      ThetaX = muonTrack->GetThetaX();
      ThetaY = muonTrack->GetThetaY();
      
      Pyz = 1./TMath::Abs(muonTrack->GetInverseBendingMomentum());
      PzRec = - Pyz / TMath::Sqrt(1.0 + TMath::Tan(ThetaY)*TMath::Tan(ThetaY));
      PxRec = PzRec * TMath::Tan(ThetaX);
      PyRec = PzRec * TMath::Tan(ThetaY);
      Charge = Int_t(TMath::Sign(1.,muonTrack->GetInverseBendingMomentum()));
      
      //ChiSquare of the track if needed
      Chisquare = muonTrack->GetChi2()/(2.0 * muonTrack->GetNHit() - 5);
      Energy = TMath::Sqrt(MUON_MASS * MUON_MASS + PxRec * PxRec + PyRec * PyRec + PzRec * PzRec);
      EPvector.SetPxPyPzE(PxRec, PyRec, PzRec, Energy);
      
      // total number of muons inside a vertex cut 
      if((TMath::Abs(Z)<ZVertexCut) && (TMath::Sqrt(Y*Y+X*X)<RhoVertexCut)) {
	nMuons++;
	if(EPvector.Pt() > LowPtCut) {
	  nMu1GeV++; 
	  if(EPvector.Pt() > HighPtCut) {
	    nMu3GeV++; 
	    if (EPvector.Pt() > VeryHighPtCut) {
	      nMu10GeV++;
	    }
	  }
	}
      }
    }//muon track loop
    
    // Fill the event tags 
    if(ntrack != 0)
      meanPt = meanPt/ntrack;
    
    evTag->SetEventId(iEventNumber+1);
    evTag->SetGUID(guid);
    evTag->SetMD5(md5);
    evTag->SetTURL(turl);
    evTag->SetSize(size);
    evTag->SetVertexX(vertexIn->GetXv());
    evTag->SetVertexY(vertexIn->GetYv());
    evTag->SetVertexZ(vertexIn->GetZv());
    
    evTag->SetT0VertexZ(esd->GetT0zVertex());
    
    evTag->SetTrigger(esd->GetTrigger());
    
    evTag->SetZDCNeutronEnergy(esd->GetZDCNEnergy());
    evTag->SetZDCProtonEnergy(esd->GetZDCPEnergy());
    evTag->SetZDCEMEnergy(esd->GetZDCEMEnergy());
    evTag->SetNumOfParticipants(esd->GetZDCParticipants());
    
    
    evTag->SetNumOfTracks(esd->GetNumberOfTracks());
    evTag->SetNumOfPosTracks(nPos);
    evTag->SetNumOfNegTracks(nNeg);
    evTag->SetNumOfNeutrTracks(nNeutr);
    
    evTag->SetNumOfV0s(esd->GetNumberOfV0s());
    evTag->SetNumOfCascades(esd->GetNumberOfCascades());
    evTag->SetNumOfKinks(esd->GetNumberOfKinks());
    evTag->SetNumOfPMDTracks(esd->GetNumberOfPmdTracks());
    
    evTag->SetNumOfProtons(nProtons);
    evTag->SetNumOfKaons(nKaons);
    evTag->SetNumOfPions(nPions);
    evTag->SetNumOfMuons(nMuons);
    evTag->SetNumOfElectrons(nElectrons);
    evTag->SetNumOfPhotons(nGamas);
    evTag->SetNumOfPi0s(nPi0s);
    evTag->SetNumOfNeutrons(nNeutrons);
    evTag->SetNumOfKaon0s(nK0s);
    
    evTag->SetNumOfChargedAbove1GeV(nCh1GeV);
    evTag->SetNumOfChargedAbove3GeV(nCh3GeV);
    evTag->SetNumOfChargedAbove10GeV(nCh10GeV);
    evTag->SetNumOfMuonsAbove1GeV(nMu1GeV);
    evTag->SetNumOfMuonsAbove3GeV(nMu3GeV);
    evTag->SetNumOfMuonsAbove10GeV(nMu10GeV);
    evTag->SetNumOfElectronsAbove1GeV(nEl1GeV);
    evTag->SetNumOfElectronsAbove3GeV(nEl3GeV);
    evTag->SetNumOfElectronsAbove10GeV(nEl10GeV);
    
    evTag->SetNumOfPHOSTracks(esd->GetNumberOfPHOSParticles());
    evTag->SetNumOfEMCALTracks(esd->GetNumberOfEMCALParticles());
    
    evTag->SetTotalMomentum(totalP);
    evTag->SetMeanPt(meanPt);
    evTag->SetMaxPt(maxPt);
    
    tag->AddEventTag(*evTag);
  }//event loop
  lastEvent = iNumberOfEvents;
  
  t->Delete("");
  
  ttag.Fill();
  tag->Clear();

  TString localFileName = "Run"; localFileName += tag->GetRunId(); 
  localFileName += ".Event"; localFileName += firstEvent; localFileName += "_"; localFileName += lastEvent; localFileName += "."; localFileName += Counter;
  localFileName += ".ESD.tag.root";

  TString alienLocation = "/alien";
  alienLocation += gGrid->Pwd();
  alienLocation += fgridpath.Data();
  alienLocation += "/";
  alienLocation +=  localFileName;
  alienLocation += "?se=";
  alienLocation += fSE.Data();

  TString fileName;
  
  if(fStorage == 0)
    {
      fileName = localFileName.Data();      
      AliInfo(Form("Writing tags to local file: %s",fileName.Data()));
    }
  if(fStorage == 1)
    {
      fileName = alienLocation.Data();
      AliInfo(Form("Writing tags to grid file: %s",fileName.Data()));
    }

  TFile* ftag = TFile::Open(fileName, "recreate");
  ftag->cd();
  ttag.Write();
  ftag->Close();

  delete ftag;
  delete esd;

  delete tag;
  delete evTag;
}
Commit	Line	Data
b45e5084	1	/**************************************************************************
	2	* Author: Panos Christakoglou. *
	3	* Contributors are mentioned in the code where appropriate. *
	4	* *
	5	* Permission to use, copy, modify and distribute this software and its *
	6	* documentation strictly for non-commercial purposes is hereby granted *
	7	* without fee, provided that the above copyright notice appears in all *
	8	* copies and that both the copyright notice and this permission notice *
	9	* appear in the supporting documentation. The authors make no claims *
	10	* about the suitability of this software for any purpose. It is *
	11	* provided "as is" without express or implied warranty. *
	12	**************************************************************************/
	13
	14	/* $Id$ */
	15
	16	//-----------------------------------------------------------------
	17	// AliTagCreator class
	18	// This is the class to deal with the tag creation (post process)
	19	// Origin: Panos Christakoglou, UOA-CERN, Panos.Christakoglou@cern.ch
	20	//-----------------------------------------------------------------
	21
	22	//ROOT
b45e5084	23	#include <TFile.h>
b45e5084	24	#include <TString.h>
cb1645b7	25	#include <TTree.h>
4fd84456	26	#include <TSystem.h>
4fd84456	27	#include <TChain.h>
2bdb9d38	28	#include <TLorentzVector.h>
b45e5084	29
	30	//ROOT-AliEn
	31	#include <TGrid.h>
b45e5084	32	#include <TGridResult.h>
b45e5084	33
	34	//AliRoot
	35	#include "AliRunTag.h"
	36	#include "AliEventTag.h"
	37	#include "AliESD.h"
	38	#include "AliESDVertex.h"
b45e5084	39	#include "AliLog.h"
	40
	41
	42	#include "AliTagCreator.h"
	43
	44
	45	ClassImp(AliTagCreator)
	46
	47
	48	//______________________________________________________________________________
	49	AliTagCreator::AliTagCreator() //local mode
	50	{
	51	//==============Default constructor for a AliTagCreator==================
28afeb2e	52	fgridpath = "";
b45e5084	53	fSE = "";
28afeb2e	54	fStorage = 0;
b45e5084	55	}
b45e5084	56
b45e5084	57	//______________________________________________________________________________
	58	AliTagCreator::~AliTagCreator()
	59	{
	60	//================Default destructor for a AliTagCreator=======================
	61	}
	62
28afeb2e	63	//______________________________________________________________________________
	64	void AliTagCreator::SetStorage(Int_t storage)
	65	{
cb1645b7	66	// Sets correctly the storage: 0 for local, 1 for GRID
28afeb2e	67	fStorage = storage;
	68	if(fStorage == 0)
	69	AliInfo(Form("Tags will be stored locally...."));
	70	if(fStorage == 1)
	71	AliInfo(Form("Tags will be stored in the grid...."));
	72	if((fStorage != 0)&&(fStorage != 1))
	73	{
	74	AliInfo(Form("Storage was not properly set!!!"));
	75	abort();
	76	}
	77	}
	78
b45e5084	79
b45e5084	80	//______________________________________________________________________________
cb1645b7	81	Bool_t AliTagCreator::ReadESDCollection(TGridResult *fresult)
b45e5084	82	{
cb1645b7	83	// Reads the entry of the TGridResult and creates the tags
cb1645b7	84	Int_t nEntries = fresult->GetEntries();
b45e5084	85
cb1645b7	86	TString alienUrl;
b45e5084	87	const char *guid;
e16601cf	88	const char *md5;
	89	const char *turl;
	90	Long64_t size = -1;
	91
b45e5084	92	Int_t counter = 0;
cb1645b7	93	for(Int_t i = 0; i < nEntries; i++)
09c52d4a	94	{
cb1645b7	95	alienUrl = fresult->GetKey(i,"turl");
09c52d4a	96	guid = fresult->GetKey(i,"guid");
e16601cf	97	if(fresult->GetKey(i,"size"))
	98	size = atol (fresult->GetKey(i,"size"));
	99	md5 = fresult->GetKey(i,"md5");
	100	turl = fresult->GetKey(i,"turl");
	101	if(md5 && !strlen(guid))
	102	md5 = 0;
	103	if(guid && !strlen(guid))
	104	guid = 0;
	105
cb1645b7	106	TFile *f = TFile::Open(alienUrl,"READ");
e16601cf	107	CreateTag(f,guid,md5,turl,size,counter);
09c52d4a	108	f->Close();
09c52d4a	109	delete f;
b45e5084	110	counter += 1;
09c52d4a	111	}//grid result loop
b45e5084	112
	113	return kTRUE;
	114	}
	115
4fd84456	116	//__________________________________________________________________________
	117	Bool_t AliTagCreator::MergeTags()
	118	{
	119	//Merges the tags and stores the merged tag file
	120	//locally if fStorage=0 or in the grid if fStorage=1
	121	AliInfo(Form("Merging tags....."));
	122	TChain *fgChain = new TChain("T");
	123
	124	if(fStorage == 0) {
	125	const char * tagPattern = "tag";
	126	// Open the working directory
	127	void * dirp = gSystem->OpenDirectory(gSystem->pwd());
	128	const char * name = 0x0;
	129	// Add all files matching pattern to the chain
	130	while((name = gSystem->GetDirEntry(dirp))) {
	131	if (strstr(name,tagPattern))
	132	fgChain->Add(name);
	133	}//directory loop
	134	AliInfo(Form("Chained tag files: %d",fgChain->GetEntries()));
	135	}//local mode
	136
	137	else if(fStorage == 1) {
	138	TString alienLocation = gGrid->Pwd();
	139	alienLocation += fgridpath.Data();
	140	alienLocation += "/";
	141
	142	TGridResult tagresult = gGrid->Query(alienLocation,"tag.root","","");
	143	Int_t nEntries = tagresult->GetEntries();
	144	for(Int_t i = 0; i < nEntries; i++) {
	145	TString alienUrl = tagresult->GetKey(i,"turl");
	146	fgChain->Add(alienUrl);
	147	}//grid result loop
	148	AliInfo(Form("Chained tag files: %d",fgChain->GetEntries()));
	149	}//grid mode
	150
	151	AliRunTag *tag = new AliRunTag;
	152	AliEventTag *evTag = new AliEventTag;
	153	fgChain->SetBranchAddress("AliTAG",&tag);
	154
	155	//Defining new tag objects
	156	AliRunTag *newTag = new AliRunTag();
	157	TTree ttag("T","A Tree with event tags");
	158	TBranch * btag = ttag.Branch("AliTAG", &newTag);
	159	btag->SetCompressionLevel(9);
	160	for(Int_t iTagFiles = 0; iTagFiles < fgChain->GetEntries(); iTagFiles++) {
	161	fgChain->GetEntry(iTagFiles);
	162	newTag->SetRunId(tag->GetRunId());
	163	const TClonesArray *tagList = tag->GetEventTags();
	164	for(Int_t j = 0; j < tagList->GetEntries(); j++) {
	165	evTag = (AliEventTag *) tagList->At(j);
	166	newTag->AddEventTag(*evTag);
	167	}
	168	ttag.Fill();
	169	newTag->Clear();
	170	}//tag file loop
	171
	172	TString localFileName = "Run"; localFileName += tag->GetRunId();
	173	localFileName += ".Merged"; localFileName += ".ESD.tag.root";
	174
	175	TString alienFileName = "/alien";
	176	alienFileName += gGrid->Pwd();
	177	alienFileName += fgridpath.Data();
	178	alienFileName += "/";
	179	alienFileName += localFileName;
180	alienFileName += "?se=";
181	alienFileName += fSE.Data();
182
183	TString filename = 0x0;
184
185	if(fStorage == 0) {
186	filename = localFileName.Data();
187	AliInfo(Form("Writing merged tags to local file: %s",filename.Data()));
188	}
189	if(fStorage == 1) {
190	filename = alienFileName.Data();
191	AliInfo(Form("Writing merged tags to grid file: %s",filename.Data()));
192	}
193
194	TFile* ftag = TFile::Open(filename, "recreate");
195	ftag->cd();
196	ttag.Write();
197	ftag->Close();
198
199	delete tag;
200	delete evTag;
201	delete newTag;
202
203	return kTRUE;
204	}
b45e5084	205
b45e5084	206	//_____________________________________________________________________________
e16601cf	207	void AliTagCreator::CreateTag(TFile* file, const char guid, const char md5, const char *turl, Long64_t size, Int_t Counter)
b45e5084	208	{
2bdb9d38	209	/////////////
	210	//muon code//
	211	////////////
	212	Double_t MUON_MASS = 0.105658369;
	213	//Variables
	214	Double_t X,Y,Z ;
	215	Double_t ThetaX, ThetaY, Pyz,Chisquare;
	216	Double_t PxRec,PyRec, PzRec, Energy;
	217	Int_t Charge;
	218	TLorentzVector EPvector;
	219
	220	Float_t ZVertexCut = 10.0;
	221	Float_t RhoVertexCut = 2.0;
	222
	223	Float_t LowPtCut = 1.0;
	224	Float_t HighPtCut = 3.0;
	225	Float_t VeryHighPtCut = 10.0;
	226	////////////
	227
	228	Double_t partFrac[5] = {0.01, 0.01, 0.85, 0.10, 0.05};
	229
cb1645b7	230	// Creates the tags for all the events in a given ESD file
b45e5084	231	Int_t ntrack;
cb1645b7	232	Int_t nProtons, nKaons, nPions, nMuons, nElectrons;
	233	Int_t nPos, nNeg, nNeutr;
	234	Int_t nK0s, nNeutrons, nPi0s, nGamas;
	235	Int_t nCh1GeV, nCh3GeV, nCh10GeV;
	236	Int_t nMu1GeV, nMu3GeV, nMu10GeV;
	237	Int_t nEl1GeV, nEl3GeV, nEl10GeV;
	238	Float_t maxPt = .0, meanPt = .0, totalP = .0;
b45e5084	239
b45e5084	240	AliRunTag *tag = new AliRunTag();
b45e5084	241	AliEventTag *evTag = new AliEventTag();
b45e5084	242	TTree ttag("T","A Tree with event tags");
4fd84456	243	TBranch * btag = ttag.Branch("AliTAG", &tag);
b45e5084	244	btag->SetCompressionLevel(9);
	245
	246	AliInfo(Form("Creating the tags......."));
	247
	248	Int_t firstEvent = 0,lastEvent = 0;
	249	TTree t = (TTree) file->Get("esdTree");
	250	TBranch * b = t->GetBranch("ESD");
	251	AliESD *esd = 0;
	252	b->SetAddress(&esd);
	253
	254	tag->SetRunId(esd->GetRunNumber());
	255
cb1645b7	256	Int_t iNumberOfEvents = b->GetEntries();
2bdb9d38	257	for (Int_t iEventNumber = 0; iEventNumber < iNumberOfEvents; iEventNumber++) {
	258	ntrack = 0;
	259	nPos = 0;
	260	nNeg = 0;
	261	nNeutr =0;
	262	nK0s = 0;
	263	nNeutrons = 0;
	264	nPi0s = 0;
	265	nGamas = 0;
	266	nProtons = 0;
	267	nKaons = 0;
	268	nPions = 0;
	269	nMuons = 0;
	270	nElectrons = 0;
	271	nCh1GeV = 0;
	272	nCh3GeV = 0;
	273	nCh10GeV = 0;
	274	nMu1GeV = 0;
	275	nMu3GeV = 0;
	276	nMu10GeV = 0;
	277	nEl1GeV = 0;
	278	nEl3GeV = 0;
	279	nEl10GeV = 0;
	280	maxPt = .0;
	281	meanPt = .0;
	282	totalP = .0;
	283
	284	b->GetEntry(iEventNumber);
	285	const AliESDVertex * vertexIn = esd->GetVertex();
	286
	287	for (Int_t iTrackNumber = 0; iTrackNumber < esd->GetNumberOfTracks(); iTrackNumber++) {
	288	AliESDtrack * esdTrack = esd->GetTrack(iTrackNumber);
	289	UInt_t status = esdTrack->GetStatus();
b45e5084	290
2bdb9d38	291	//select only tracks with ITS refit
	292	if ((status&AliESDtrack::kITSrefit)==0) continue;
	293	//select only tracks with TPC refit
	294	if ((status&AliESDtrack::kTPCrefit)==0) continue;
b45e5084	295
2bdb9d38	296	//select only tracks with the "combined PID"
	297	if ((status&AliESDtrack::kESDpid)==0) continue;
	298	Double_t p[3];
	299	esdTrack->GetPxPyPz(p);
	300	Double_t momentum = sqrt(pow(p[0],2) + pow(p[1],2) + pow(p[2],2));
	301	Double_t fPt = sqrt(pow(p[0],2) + pow(p[1],2));
	302	totalP += momentum;
	303	meanPt += fPt;
	304	if(fPt > maxPt) maxPt = fPt;
b45e5084	305
2bdb9d38	306	if(esdTrack->GetSign() > 0) {
	307	nPos++;
	308	if(fPt > LowPtCut) nCh1GeV++;
	309	if(fPt > HighPtCut) nCh3GeV++;
	310	if(fPt > VeryHighPtCut) nCh10GeV++;
	311	}
	312	if(esdTrack->GetSign() < 0) {
	313	nNeg++;
	314	if(fPt > LowPtCut) nCh1GeV++;
	315	if(fPt > HighPtCut) nCh3GeV++;
	316	if(fPt > VeryHighPtCut) nCh10GeV++;
	317	}
	318	if(esdTrack->GetSign() == 0) nNeutr++;
b45e5084	319
2bdb9d38	320	//PID
	321	Double_t prob[5];
	322	esdTrack->GetESDpid(prob);
b45e5084	323
2bdb9d38	324	Double_t rcc = 0.0;
	325	for(Int_t i = 0; i < AliPID::kSPECIES; i++) rcc += prob[i]*partFrac[i];
	326	if(rcc == 0.0) continue;
	327	//Bayes' formula
	328	Double_t w[5];
	329	for(Int_t i = 0; i < AliPID::kSPECIES; i++) w[i] = prob[i]*partFrac[i]/rcc;
b45e5084	330
2bdb9d38	331	//protons
	332	if ((w[4]>w[3])&&(w[4]>w[2])&&(w[4]>w[1])&&(w[4]>w[0])) nProtons++;
	333	//kaons
	334	if ((w[3]>w[4])&&(w[3]>w[2])&&(w[3]>w[1])&&(w[3]>w[0])) nKaons++;
	335	//pions
	336	if ((w[2]>w[4])&&(w[2]>w[3])&&(w[2]>w[1])&&(w[2]>w[0])) nPions++;
	337	//electrons
	338	if ((w[0]>w[4])&&(w[0]>w[3])&&(w[0]>w[2])&&(w[0]>w[1])) {
	339	nElectrons++;
	340	if(fPt > LowPtCut) nEl1GeV++;
	341	if(fPt > HighPtCut) nEl3GeV++;
	342	if(fPt > VeryHighPtCut) nEl10GeV++;
	343	}
	344	ntrack++;
	345	}//esd track loop
	346
	347	/////////////
	348	//muon code//
	349	////////////
	350	Int_t nMuonTracks = esd->GetNumberOfMuonTracks();
	351	// loop over all reconstructed tracks (also first track of combination)
	352	for (Int_t iTrack = 0; iTrack < nMuonTracks; iTrack++) {
	353	AliESDMuonTrack* muonTrack = esd->GetMuonTrack(iTrack);
	354	if (muonTrack == 0x0) continue;
b45e5084	355
2bdb9d38	356	// Coordinates at vertex
	357	Z = muonTrack->GetZ();
	358	Y = muonTrack->GetBendingCoor();
	359	X = muonTrack->GetNonBendingCoor();
b45e5084	360
2bdb9d38	361	ThetaX = muonTrack->GetThetaX();
2bdb9d38	362	ThetaY = muonTrack->GetThetaY();
b45e5084	363
2bdb9d38	364	Pyz = 1./TMath::Abs(muonTrack->GetInverseBendingMomentum());
	365	PzRec = - Pyz / TMath::Sqrt(1.0 + TMath::Tan(ThetaY)*TMath::Tan(ThetaY));
	366	PxRec = PzRec * TMath::Tan(ThetaX);
	367	PyRec = PzRec * TMath::Tan(ThetaY);
	368	Charge = Int_t(TMath::Sign(1.,muonTrack->GetInverseBendingMomentum()));
b45e5084	369
2bdb9d38	370	//ChiSquare of the track if needed
	371	Chisquare = muonTrack->GetChi2()/(2.0 * muonTrack->GetNHit() - 5);
	372	Energy = TMath::Sqrt(MUON_MASS * MUON_MASS + PxRec * PxRec + PyRec * PyRec + PzRec * PzRec);
	373	EPvector.SetPxPyPzE(PxRec, PyRec, PzRec, Energy);
b45e5084	374
2bdb9d38	375	// total number of muons inside a vertex cut
	376	if((TMath::Abs(Z)<ZVertexCut) && (TMath::Sqrt(YY+XX)<RhoVertexCut)) {
	377	nMuons++;
	378	if(EPvector.Pt() > LowPtCut) {
	379	nMu1GeV++;
	380	if(EPvector.Pt() > HighPtCut) {
	381	nMu3GeV++;
	382	if (EPvector.Pt() > VeryHighPtCut) {
	383	nMu10GeV++;
	384	}
	385	}
	386	}
	387	}
	388	}//muon track loop
	389
	390	// Fill the event tags
	391	if(ntrack != 0)
	392	meanPt = meanPt/ntrack;
	393
	394	evTag->SetEventId(iEventNumber+1);
	395	evTag->SetGUID(guid);
	396	evTag->SetMD5(md5);
	397	evTag->SetTURL(turl);
	398	evTag->SetSize(size);
	399	evTag->SetVertexX(vertexIn->GetXv());
	400	evTag->SetVertexY(vertexIn->GetYv());
	401	evTag->SetVertexZ(vertexIn->GetZv());
	402
	403	evTag->SetT0VertexZ(esd->GetT0zVertex());
	404
	405	evTag->SetTrigger(esd->GetTrigger());
	406
	407	evTag->SetZDCNeutronEnergy(esd->GetZDCNEnergy());
	408	evTag->SetZDCProtonEnergy(esd->GetZDCPEnergy());
	409	evTag->SetZDCEMEnergy(esd->GetZDCEMEnergy());
	410	evTag->SetNumOfParticipants(esd->GetZDCParticipants());
	411
	412
	413	evTag->SetNumOfTracks(esd->GetNumberOfTracks());
	414	evTag->SetNumOfPosTracks(nPos);
	415	evTag->SetNumOfNegTracks(nNeg);
	416	evTag->SetNumOfNeutrTracks(nNeutr);
	417
	418	evTag->SetNumOfV0s(esd->GetNumberOfV0s());
	419	evTag->SetNumOfCascades(esd->GetNumberOfCascades());
	420	evTag->SetNumOfKinks(esd->GetNumberOfKinks());
	421	evTag->SetNumOfPMDTracks(esd->GetNumberOfPmdTracks());
	422
	423	evTag->SetNumOfProtons(nProtons);
	424	evTag->SetNumOfKaons(nKaons);
	425	evTag->SetNumOfPions(nPions);
	426	evTag->SetNumOfMuons(nMuons);
	427	evTag->SetNumOfElectrons(nElectrons);
	428	evTag->SetNumOfPhotons(nGamas);
	429	evTag->SetNumOfPi0s(nPi0s);
	430	evTag->SetNumOfNeutrons(nNeutrons);
	431	evTag->SetNumOfKaon0s(nK0s);
	432
	433	evTag->SetNumOfChargedAbove1GeV(nCh1GeV);
	434	evTag->SetNumOfChargedAbove3GeV(nCh3GeV);
	435	evTag->SetNumOfChargedAbove10GeV(nCh10GeV);
	436	evTag->SetNumOfMuonsAbove1GeV(nMu1GeV);
	437	evTag->SetNumOfMuonsAbove3GeV(nMu3GeV);
	438	evTag->SetNumOfMuonsAbove10GeV(nMu10GeV);
439	evTag->SetNumOfElectronsAbove1GeV(nEl1GeV);
440	evTag->SetNumOfElectronsAbove3GeV(nEl3GeV);
441	evTag->SetNumOfElectronsAbove10GeV(nEl10GeV);
442
443	evTag->SetNumOfPHOSTracks(esd->GetNumberOfPHOSParticles());
444	evTag->SetNumOfEMCALTracks(esd->GetNumberOfEMCALParticles());
445
446	evTag->SetTotalMomentum(totalP);
447	evTag->SetMeanPt(meanPt);
448	evTag->SetMaxPt(maxPt);
449
450	tag->AddEventTag(*evTag);
451	}//event loop
cb1645b7	452	lastEvent = iNumberOfEvents;
b45e5084	453
b45e5084	454	t->Delete("");
2bdb9d38	455
b45e5084	456	ttag.Fill();
	457	tag->Clear();
	458
cb1645b7	459	TString localFileName = "Run"; localFileName += tag->GetRunId();
	460	localFileName += ".Event"; localFileName += firstEvent; localFileName += "_"; localFileName += lastEvent; localFileName += "."; localFileName += Counter;
	461	localFileName += ".ESD.tag.root";
b45e5084	462
cb1645b7	463	TString alienLocation = "/alien";
	464	alienLocation += gGrid->Pwd();
	465	alienLocation += fgridpath.Data();
	466	alienLocation += "/";
	467	alienLocation += localFileName;
	468	alienLocation += "?se=";
	469	alienLocation += fSE.Data();
28afeb2e	470
cb1645b7	471	TString fileName;
28afeb2e	472
	473	if(fStorage == 0)
	474	{
cb1645b7	475	fileName = localFileName.Data();
cb1645b7	476	AliInfo(Form("Writing tags to local file: %s",fileName.Data()));
28afeb2e	477	}
	478	if(fStorage == 1)
	479	{
cb1645b7	480	fileName = alienLocation.Data();
cb1645b7	481	AliInfo(Form("Writing tags to grid file: %s",fileName.Data()));
28afeb2e	482	}
b45e5084	483
cb1645b7	484	TFile* ftag = TFile::Open(fileName, "recreate");
b45e5084	485	ftag->cd();
	486	ttag.Write();
	487	ftag->Close();
	488
	489	delete ftag;
	490	delete esd;
	491
	492	delete tag;
b45e5084	493	delete evTag;
	494	}
	495