STEER/AliESDTagCreator.cxx

   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 //           AliESDTagCreator 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
  23 #include <Riostream.h>
  24 #include <TFile.h>
  25 #include <TString.h>
  26 #include <TTree.h>
  27 #include <TSystem.h>
  28 #include <TChain.h>
  29 #include <TLorentzVector.h>
  30
  31 //ROOT-AliEn
  32 #include <TGrid.h>
  33 #include <TGridResult.h>
  34
  35 //AliRoot
  36 #include "AliRunTag.h"
  37 #include "AliEventTag.h"
  38 #include "AliESD.h"
  39 #include "AliESDEvent.h"
  40 #include "AliESDVertex.h"
  41 #include "AliLog.h"
  42
  43 #include "AliESDTagCreator.h"
  44
  45
  46 ClassImp(AliESDTagCreator)
  47
  48
  49 //______________________________________________________________________________
  50   AliESDTagCreator::AliESDTagCreator() :
  51     AliTagCreator() {
  52   //==============Default constructor for a AliESDTagCreator================
  53 }
  54
  55 //______________________________________________________________________________
  56 AliESDTagCreator::~AliESDTagCreator() {
  57 //================Default destructor for a AliESDTagCreator===================
  58 }
  59
  60 //______________________________________________________________________________
  61 Bool_t AliESDTagCreator::ReadGridCollection(TGridResult *fresult) {
  62   // Reads the entry of the TGridResult and creates the tags
  63   Int_t nEntries = fresult->GetEntries();
  64
  65   TString alienUrl;
  66   const char *guid;
  67   const char *md5;
  68   const char *turl;
  69   Long64_t size = -1;
  70
  71   Int_t counter = 0;
  72   for(Int_t i = 0; i < nEntries; i++) {
  73     alienUrl = fresult->GetKey(i,"turl");
  74     guid = fresult->GetKey(i,"guid");
  75     if(fresult->GetKey(i,"size")) size = atol (fresult->GetKey(i,"size"));
  76     md5 = fresult->GetKey(i,"md5");
  77     turl = fresult->GetKey(i,"turl");
  78     if(md5 && !strlen(guid)) md5 = 0;
  79     if(guid && !strlen(guid)) guid = 0;
  80
  81     TFile *f = TFile::Open(alienUrl,"READ");
  82     CreateTag(f,guid,md5,turl,size,counter);
  83     f->Close();
  84     delete f;
  85     counter += 1;
  86   }//grid result loop
  87
  88   return kTRUE;
  89 }
  90
  91 //______________________________________________________________________________
  92 Bool_t AliESDTagCreator::ReadLocalCollection(const char *localpath) {
  93   // Checks the different subdirs of the given local path and in the
  94   // case where it finds an AliESDs.root file it creates the tags
  95
  96   void *dira =  gSystem->OpenDirectory(localpath);
  97   Char_t fPath[256];
  98   const char * dirname = 0x0;
  99   const char * filename = 0x0;
 100   const char * pattern = "AliESDs.root";
 101
 102   Int_t counter = 0;
 103   while((dirname = gSystem->GetDirEntry(dira))) {
 104     sprintf(fPath,"%s/%s",localpath,dirname);
 105     void *dirb =  gSystem->OpenDirectory(fPath);
 106     while((filename = gSystem->GetDirEntry(dirb))) {
 107       if(strstr(filename,pattern)) {
 108         TString fESDFileName;
 109         fESDFileName = fPath;
 110         fESDFileName += "/";
 111         fESDFileName += pattern;
 112         TFile *f = TFile::Open(fESDFileName,"READ");
 113         CreateTag(f,fESDFileName,counter);
 114         f->Close();
 115         delete f;
 116
 117         counter += 1;
 118       }//pattern check
 119     }//child directory's entry loop
 120   }//parent directory's entry loop
 121
 122   return kTRUE;
 123 }
 124
 125 //______________________________________________________________________________
 126 Bool_t AliESDTagCreator::ReadCAFCollection(const char *filename) {
 127   // Temporary solution for CAF: Takes as an input the ascii file that
 128   // lists the ESDs stored in the SE of the CAF and creates the tags.
 129
 130   // Open the input stream
 131   ifstream in;
 132   in.open(filename);
 133
 134   Int_t counter = 0;
 135   TString esdfile;
 136   // Read the input list of files and add them to the chain
 137   while(in.good()) {
 138     in >> esdfile;
 139     if (!esdfile.Contains("root")) continue; // protection
 140     TFile *f = TFile::Open(esdfile,"READ");
 141     CreateTag(f,esdfile,counter);
 142     f->Close();
 143     delete f;
 144
 145     counter += 1;
 146   }
 147
 148   return kTRUE;
 149 }
 150
 151 //_____________________________________________________________________________
 152 void AliESDTagCreator::CreateTag(TFile* file, const char *guid, const char *md5, const char *turl, Long64_t size, Int_t Counter) {
 153   //private method that creates tag files
 154   TString fguid = guid;
 155   TString fmd5 = md5;
 156   TString fturl = turl;
 157   /////////////
 158   //muon code//
 159   ////////////
 160   Double_t fMUONMASS = 0.105658369;
 161   //Variables
 162   Double_t fX,fY,fZ ;
 163   Double_t fThetaX, fThetaY, fPyz, fChisquare;
 164   Double_t fPxRec, fPyRec, fPzRec, fEnergy;
 165   Int_t fCharge;
 166   TLorentzVector fEPvector;
 167
 168   Float_t fZVertexCut = 10.0;
 169   Float_t fRhoVertexCut = 2.0;
 170
 171   Float_t fLowPtCut = 1.0;
 172   Float_t fHighPtCut = 3.0;
 173   Float_t fVeryHighPtCut = 10.0;
 174   ////////////
 175
 176   Double_t partFrac[5] = {0.01, 0.01, 0.85, 0.10, 0.05};
 177
 178   // Creates the tags for all the events in a given ESD file
 179   Bool_t fIsSim = kTRUE;
 180   Int_t ntrack;
 181   Int_t nProtons, nKaons, nPions, nMuons, nElectrons;
 182   Int_t nPos, nNeg, nNeutr;
 183   Int_t nK0s, nNeutrons, nPi0s, nGamas;
 184   Int_t nCh1GeV, nCh3GeV, nCh10GeV;
 185   Int_t nMu1GeV, nMu3GeV, nMu10GeV;
 186   Int_t nEl1GeV, nEl3GeV, nEl10GeV;
 187   Float_t maxPt = .0, meanPt = .0, totalP = .0;
 188   Int_t fVertexflag;
 189   Int_t iRunNumber = 0;
 190   TString fVertexName;
 191
 192   AliRunTag *tag = new AliRunTag();
 193   AliEventTag *evTag = new AliEventTag();
 194   TTree ttag("T","A Tree with event tags");
 195   TBranch * btag = ttag.Branch("AliTAG", &tag);
 196   btag->SetCompressionLevel(9);
 197
 198   AliInfo(Form("Creating the tags......."));
 199
 200   Int_t firstEvent = 0,lastEvent = 0;
 201   TTree *t = (TTree*) file->Get("esdTree");
 202   AliESDEvent *esd = new AliESDEvent();
 203   esd->ReadFromTree(t);
 204
 205   t->GetEntry(0);
 206   Int_t iInitRunNumber = esd->GetRunNumber();
 207
 208   Int_t iNumberOfEvents = (Int_t)t->GetEntries();
 209   for (Int_t iEventNumber = 0; iEventNumber < iNumberOfEvents; iEventNumber++) {
 210     ntrack = 0;
 211     nPos = 0;
 212     nNeg = 0;
 213     nNeutr =0;
 214     nK0s = 0;
 215     nNeutrons = 0;
 216     nPi0s = 0;
 217     nGamas = 0;
 218     nProtons = 0;
 219     nKaons = 0;
 220     nPions = 0;
 221     nMuons = 0;
 222     nElectrons = 0;
 223     nCh1GeV = 0;
 224     nCh3GeV = 0;
 225     nCh10GeV = 0;
 226     nMu1GeV = 0;
 227     nMu3GeV = 0;
 228     nMu10GeV = 0;
 229     nEl1GeV = 0;
 230     nEl3GeV = 0;
 231     nEl10GeV = 0;
 232     maxPt = .0;
 233     meanPt = .0;
 234     totalP = .0;
 235     fVertexflag = 1;
 236
 237     t->GetEntry(iEventNumber);
 238     iRunNumber = esd->GetRunNumber();
 239     if(iRunNumber != iInitRunNumber) AliFatal("Inconsistency of run numbers in the AliESD!!!");
 240     const AliESDVertex * vertexIn = esd->GetVertex();
 241     fVertexName = vertexIn->GetName();
 242     if(fVertexName == "default") fVertexflag = 0;
 243
 244     for (Int_t iTrackNumber = 0; iTrackNumber < esd->GetNumberOfTracks(); iTrackNumber++) {
 245       AliESDtrack * esdTrack = esd->GetTrack(iTrackNumber);
 246       if(esdTrack->GetLabel() != 0) fIsSim = kTRUE;
 247       else if(esdTrack->GetLabel() == 0) fIsSim = kFALSE;
 248       UInt_t status = esdTrack->GetStatus();
 249
 250       //select only tracks with ITS refit
 251       if ((status&AliESDtrack::kITSrefit)==0) continue;
 252       //select only tracks with TPC refit
 253       if ((status&AliESDtrack::kTPCrefit)==0) continue;
 254
 255       //select only tracks with the "combined PID"
 256       if ((status&AliESDtrack::kESDpid)==0) continue;
 257       Double_t p[3];
 258       esdTrack->GetPxPyPz(p);
 259       Double_t momentum = sqrt(pow(p[0],2) + pow(p[1],2) + pow(p[2],2));
 260       Double_t fPt = sqrt(pow(p[0],2) + pow(p[1],2));
 261       totalP += momentum;
 262       meanPt += fPt;
 263       if(fPt > maxPt) maxPt = fPt;
 264
 265       if(esdTrack->GetSign() > 0) {
 266         nPos++;
 267         if(fPt > fLowPtCut) nCh1GeV++;
 268         if(fPt > fHighPtCut) nCh3GeV++;
 269         if(fPt > fVeryHighPtCut) nCh10GeV++;
 270       }
 271       if(esdTrack->GetSign() < 0) {
 272         nNeg++;
 273         if(fPt > fLowPtCut) nCh1GeV++;
 274         if(fPt > fHighPtCut) nCh3GeV++;
 275         if(fPt > fVeryHighPtCut) nCh10GeV++;
 276       }
 277       if(esdTrack->GetSign() == 0) nNeutr++;
 278
 279       //PID
 280       Double_t prob[5];
 281       esdTrack->GetESDpid(prob);
 282
 283       Double_t rcc = 0.0;
 284       for(Int_t i = 0; i < AliPID::kSPECIES; i++) rcc += prob[i]*partFrac[i];
 285       if(rcc == 0.0) continue;
 286       //Bayes' formula
 287       Double_t w[5];
 288       for(Int_t i = 0; i < AliPID::kSPECIES; i++) w[i] = prob[i]*partFrac[i]/rcc;
 289
 290       //protons
 291       if ((w[4]>w[3])&&(w[4]>w[2])&&(w[4]>w[1])&&(w[4]>w[0])) nProtons++;
 292       //kaons
 293       if ((w[3]>w[4])&&(w[3]>w[2])&&(w[3]>w[1])&&(w[3]>w[0])) nKaons++;
 294       //pions
 295       if ((w[2]>w[4])&&(w[2]>w[3])&&(w[2]>w[1])&&(w[2]>w[0])) nPions++;
 296       //electrons
 297       if ((w[0]>w[4])&&(w[0]>w[3])&&(w[0]>w[2])&&(w[0]>w[1])) {
 298         nElectrons++;
 299         if(fPt > fLowPtCut) nEl1GeV++;
 300         if(fPt > fHighPtCut) nEl3GeV++;
 301         if(fPt > fVeryHighPtCut) nEl10GeV++;
 302       }
 303       ntrack++;
 304     }//esd track loop
 305
 306     /////////////
 307     //muon code//
 308     ////////////
 309     Int_t nMuonTracks = esd->GetNumberOfMuonTracks();
 310     // loop over all reconstructed tracks (also first track of combination)
 311     for (Int_t iTrack = 0; iTrack <  nMuonTracks;  iTrack++) {
 312       AliESDMuonTrack* muonTrack = esd->GetMuonTrack(iTrack);
 313       if (muonTrack == 0x0) continue;
 314
 315       // Coordinates at vertex
 316       fZ = muonTrack->GetZ();
 317       fY = muonTrack->GetBendingCoor();
 318       fX = muonTrack->GetNonBendingCoor();
 319
 320       fThetaX = muonTrack->GetThetaX();
 321       fThetaY = muonTrack->GetThetaY();
 322
 323       fPyz = 1./TMath::Abs(muonTrack->GetInverseBendingMomentum());
 324       fPzRec = - fPyz / TMath::Sqrt(1.0 + TMath::Tan(fThetaY)*TMath::Tan(fThetaY));
 325       fPxRec = fPzRec * TMath::Tan(fThetaX);
 326       fPyRec = fPzRec * TMath::Tan(fThetaY);
 327       fCharge = Int_t(TMath::Sign(1.,muonTrack->GetInverseBendingMomentum()));
 328
 329       //ChiSquare of the track if needed
 330       fChisquare = muonTrack->GetChi2()/(2.0 * muonTrack->GetNHit() - 5);
 331       fEnergy = TMath::Sqrt(fMUONMASS * fMUONMASS + fPxRec * fPxRec + fPyRec * fPyRec + fPzRec * fPzRec);
 332       fEPvector.SetPxPyPzE(fPxRec, fPyRec, fPzRec, fEnergy);
 333
 334       // total number of muons inside a vertex cut
 335       if((TMath::Abs(fZ)<fZVertexCut) && (TMath::Sqrt(fY*fY+fX*fX)<fRhoVertexCut)) {
 336         nMuons++;
 337         if(fEPvector.Pt() > fLowPtCut) {
 338           nMu1GeV++;
 339           if(fEPvector.Pt() > fHighPtCut) {
 340             nMu3GeV++;
 341             if (fEPvector.Pt() > fVeryHighPtCut) {
 342               nMu10GeV++;
 343             }
 344           }
 345         }
 346       }
 347     }//muon track loop
 348
 349     // Fill the event tags
 350     if(ntrack != 0) meanPt = meanPt/ntrack;
 351
 352     evTag->SetEventId(iEventNumber+1);
 353     evTag->SetGUID(fguid);
 354     evTag->SetMD5(fmd5);
 355     evTag->SetTURL(fturl);
 356     evTag->SetSize(size);
 357     evTag->SetVertexX(vertexIn->GetXv());
 358     evTag->SetVertexY(vertexIn->GetYv());
 359     evTag->SetVertexZ(vertexIn->GetZv());
 360     evTag->SetVertexZError(vertexIn->GetZRes());
 361     evTag->SetVertexFlag(fVertexflag);
 362
 363     evTag->SetT0VertexZ(esd->GetT0zVertex());
 364
 365     evTag->SetTriggerMask(esd->GetTriggerMask());
 366     evTag->SetTriggerCluster(esd->GetTriggerCluster());
 367
 368     evTag->SetZDCNeutron1Energy(esd->GetZDCN1Energy());
 369     evTag->SetZDCProton1Energy(esd->GetZDCP1Energy());
 370     evTag->SetZDCEMEnergy(esd->GetZDCEMEnergy());
 371     evTag->SetZDCNeutron1Energy(esd->GetZDCN2Energy());
 372     evTag->SetZDCProton1Energy(esd->GetZDCP2Energy());
 373     evTag->SetNumOfParticipants(esd->GetZDCParticipants());
 374
 375
 376     evTag->SetNumOfTracks(esd->GetNumberOfTracks());
 377     evTag->SetNumOfPosTracks(nPos);
 378     evTag->SetNumOfNegTracks(nNeg);
 379     evTag->SetNumOfNeutrTracks(nNeutr);
 380
 381     evTag->SetNumOfV0s(esd->GetNumberOfV0s());
 382     evTag->SetNumOfCascades(esd->GetNumberOfCascades());
 383     evTag->SetNumOfKinks(esd->GetNumberOfKinks());
 384     evTag->SetNumOfPMDTracks(esd->GetNumberOfPmdTracks());
 385
 386     evTag->SetNumOfProtons(nProtons);
 387     evTag->SetNumOfKaons(nKaons);
 388     evTag->SetNumOfPions(nPions);
 389     evTag->SetNumOfMuons(nMuons);
 390     evTag->SetNumOfElectrons(nElectrons);
 391     evTag->SetNumOfPhotons(nGamas);
 392     evTag->SetNumOfPi0s(nPi0s);
 393     evTag->SetNumOfNeutrons(nNeutrons);
 394     evTag->SetNumOfKaon0s(nK0s);
 395
 396     evTag->SetNumOfChargedAbove1GeV(nCh1GeV);
 397     evTag->SetNumOfChargedAbove3GeV(nCh3GeV);
 398     evTag->SetNumOfChargedAbove10GeV(nCh10GeV);
 399     evTag->SetNumOfMuonsAbove1GeV(nMu1GeV);
 400     evTag->SetNumOfMuonsAbove3GeV(nMu3GeV);
 401     evTag->SetNumOfMuonsAbove10GeV(nMu10GeV);
 402     evTag->SetNumOfElectronsAbove1GeV(nEl1GeV);
 403     evTag->SetNumOfElectronsAbove3GeV(nEl3GeV);
 404     evTag->SetNumOfElectronsAbove10GeV(nEl10GeV);
 405
 406     evTag->SetNumOfPHOSClusters(esd->GetNumberOfPHOSClusters());
 407     evTag->SetNumOfEMCALClusters(esd->GetNumberOfEMCALClusters());
 408
 409     evTag->SetTotalMomentum(totalP);
 410     evTag->SetMeanPt(meanPt);
 411     evTag->SetMaxPt(maxPt);
 412
 413     tag->SetRunId(iInitRunNumber);
 414     if(fIsSim) tag->SetDataType(0);
 415     else tag->SetDataType(1);
 416     tag->AddEventTag(*evTag);
 417   }//event loop
 418   lastEvent = iNumberOfEvents;
 419
 420   t->Delete("");
 421
 422   ttag.Fill();
 423   tag->Clear();
 424
 425   TString localFileName = "Run"; localFileName += tag->GetRunId();
 426   localFileName += ".Event"; localFileName += firstEvent; localFileName += "_"; localFileName += lastEvent; localFileName += "."; localFileName += Counter;
 427   localFileName += ".ESD.tag.root";
 428
 429   TString fileName;
 430
 431   if(fStorage == 0) {
 432     fileName = localFileName.Data();
 433     AliInfo(Form("Writing tags to local file: %s",fileName.Data()));
 434   }
 435   else if(fStorage == 1) {
 436     TString alienLocation = "/alien";
 437     alienLocation += gGrid->Pwd();
 438     alienLocation += fgridpath.Data();
 439     alienLocation += "/";
 440     alienLocation +=  localFileName;
 441     alienLocation += "?se=";
 442     alienLocation += fSE.Data();
 443     fileName = alienLocation.Data();
 444     AliInfo(Form("Writing tags to grid file: %s",fileName.Data()));
 445   }
 446
 447   TFile* ftag = TFile::Open(fileName, "recreate");
 448   ftag->cd();
 449   ttag.Write();
 450   ftag->Close();
 451
 452   delete ftag;
 453   delete esd;
 454
 455   delete tag;
 456 }
 457
 458 //_____________________________________________________________________________
 459 void AliESDTagCreator::CreateTag(TFile* file, const char *filepath, Int_t Counter) {
 460   //private method that creates tag files
 461   /////////////
 462   //muon code//
 463   ////////////
 464   Double_t fMUONMASS = 0.105658369;
 465   //Variables
 466   Double_t fX,fY,fZ ;
 467   Double_t fThetaX, fThetaY, fPyz, fChisquare;
 468   Double_t fPxRec, fPyRec, fPzRec, fEnergy;
 469   Int_t fCharge;
 470   TLorentzVector fEPvector;
 471
 472   Float_t fZVertexCut = 10.0;
 473   Float_t fRhoVertexCut = 2.0;
 474
 475   Float_t fLowPtCut = 1.0;
 476   Float_t fHighPtCut = 3.0;
 477   Float_t fVeryHighPtCut = 10.0;
 478   ////////////
 479
 480   Double_t partFrac[5] = {0.01, 0.01, 0.85, 0.10, 0.05};
 481
 482   // Creates the tags for all the events in a given ESD file
 483   Bool_t fIsSim = kTRUE;
 484   Int_t ntrack;
 485   Int_t nProtons, nKaons, nPions, nMuons, nElectrons;
 486   Int_t nPos, nNeg, nNeutr;
 487   Int_t nK0s, nNeutrons, nPi0s, nGamas;
 488   Int_t nCh1GeV, nCh3GeV, nCh10GeV;
 489   Int_t nMu1GeV, nMu3GeV, nMu10GeV;
 490   Int_t nEl1GeV, nEl3GeV, nEl10GeV;
 491   Float_t maxPt = .0, meanPt = .0, totalP = .0;
 492   Int_t fVertexflag;
 493   Int_t iRunNumber = 0;
 494   TString fVertexName;
 495
 496   AliRunTag *tag = new AliRunTag();
 497   AliEventTag *evTag = new AliEventTag();
 498   TTree ttag("T","A Tree with event tags");
 499   TBranch * btag = ttag.Branch("AliTAG", &tag);
 500   btag->SetCompressionLevel(9);
 501
 502   AliInfo(Form("Creating the tags......."));
 503
 504   Int_t firstEvent = 0,lastEvent = 0;
 505
 506   TTree *t = (TTree*) file->Get("esdTree");
 507   AliESDEvent *esd = new AliESDEvent();
 508   esd->ReadFromTree(t);
 509
 510   t->GetEntry(0);
 511   Int_t iInitRunNumber = esd->GetRunNumber();
 512
 513   Int_t iNumberOfEvents = (Int_t)t->GetEntries();
 514   for (Int_t iEventNumber = 0; iEventNumber < iNumberOfEvents; iEventNumber++) {
 515     ntrack = 0;
 516     nPos = 0;
 517     nNeg = 0;
 518     nNeutr =0;
 519     nK0s = 0;
 520     nNeutrons = 0;
 521     nPi0s = 0;
 522     nGamas = 0;
 523     nProtons = 0;
 524     nKaons = 0;
 525     nPions = 0;
 526     nMuons = 0;
 527     nElectrons = 0;
 528     nCh1GeV = 0;
 529     nCh3GeV = 0;
 530     nCh10GeV = 0;
 531     nMu1GeV = 0;
 532     nMu3GeV = 0;
 533     nMu10GeV = 0;
 534     nEl1GeV = 0;
 535     nEl3GeV = 0;
 536     nEl10GeV = 0;
 537     maxPt = .0;
 538     meanPt = .0;
 539     totalP = .0;
 540     fVertexflag = 1;
 541
 542     t->GetEntry(iEventNumber);
 543     iRunNumber = esd->GetRunNumber();
 544     if(iRunNumber != iInitRunNumber) AliFatal("Inconsistency of run numbers in the AliESD!!!");
 545     const AliESDVertex * vertexIn = esd->GetVertex();
 546     fVertexName = vertexIn->GetName();
 547     if(fVertexName == "default") fVertexflag = 0;
 548
 549     for (Int_t iTrackNumber = 0; iTrackNumber < esd->GetNumberOfTracks(); iTrackNumber++) {
 550       AliESDtrack * esdTrack = esd->GetTrack(iTrackNumber);
 551       if(esdTrack->GetLabel() != 0) fIsSim = kTRUE;
 552       else if(esdTrack->GetLabel() == 0) fIsSim = kFALSE;
 553       UInt_t status = esdTrack->GetStatus();
 554
 555       //select only tracks with ITS refit
 556       if ((status&AliESDtrack::kITSrefit)==0) continue;
 557       //select only tracks with TPC refit
 558       if ((status&AliESDtrack::kTPCrefit)==0) continue;
 559
 560       //select only tracks with the "combined PID"
 561       if ((status&AliESDtrack::kESDpid)==0) continue;
 562       Double_t p[3];
 563       esdTrack->GetPxPyPz(p);
 564       Double_t momentum = sqrt(pow(p[0],2) + pow(p[1],2) + pow(p[2],2));
 565       Double_t fPt = sqrt(pow(p[0],2) + pow(p[1],2));
 566       totalP += momentum;
 567       meanPt += fPt;
 568       if(fPt > maxPt) maxPt = fPt;
 569
 570       if(esdTrack->GetSign() > 0) {
 571         nPos++;
 572         if(fPt > fLowPtCut) nCh1GeV++;
 573         if(fPt > fHighPtCut) nCh3GeV++;
 574         if(fPt > fVeryHighPtCut) nCh10GeV++;
 575       }
 576       if(esdTrack->GetSign() < 0) {
 577         nNeg++;
 578         if(fPt > fLowPtCut) nCh1GeV++;
 579         if(fPt > fHighPtCut) nCh3GeV++;
 580         if(fPt > fVeryHighPtCut) nCh10GeV++;
 581       }
 582       if(esdTrack->GetSign() == 0) nNeutr++;
 583
 584       //PID
 585       Double_t prob[5];
 586       esdTrack->GetESDpid(prob);
 587
 588       Double_t rcc = 0.0;
 589       for(Int_t i = 0; i < AliPID::kSPECIES; i++) rcc += prob[i]*partFrac[i];
 590       if(rcc == 0.0) continue;
 591       //Bayes' formula
 592       Double_t w[5];
 593       for(Int_t i = 0; i < AliPID::kSPECIES; i++) w[i] = prob[i]*partFrac[i]/rcc;
 594
 595       //protons
 596       if ((w[4]>w[3])&&(w[4]>w[2])&&(w[4]>w[1])&&(w[4]>w[0])) nProtons++;
 597       //kaons
 598       if ((w[3]>w[4])&&(w[3]>w[2])&&(w[3]>w[1])&&(w[3]>w[0])) nKaons++;
 599       //pions
 600       if ((w[2]>w[4])&&(w[2]>w[3])&&(w[2]>w[1])&&(w[2]>w[0])) nPions++;
 601       //electrons
 602       if ((w[0]>w[4])&&(w[0]>w[3])&&(w[0]>w[2])&&(w[0]>w[1])) {
 603         nElectrons++;
 604         if(fPt > fLowPtCut) nEl1GeV++;
 605         if(fPt > fHighPtCut) nEl3GeV++;
 606         if(fPt > fVeryHighPtCut) nEl10GeV++;
 607       }
 608       ntrack++;
 609     }//esd track loop
 610
 611     /////////////
 612     //muon code//
 613     ////////////
 614     Int_t nMuonTracks = esd->GetNumberOfMuonTracks();
 615     // loop over all reconstructed tracks (also first track of combination)
 616     for (Int_t iTrack = 0; iTrack <  nMuonTracks;  iTrack++) {
 617       AliESDMuonTrack* muonTrack = esd->GetMuonTrack(iTrack);
 618       if (muonTrack == 0x0) continue;
 619
 620       // Coordinates at vertex
 621       fZ = muonTrack->GetZ();
 622       fY = muonTrack->GetBendingCoor();
 623       fX = muonTrack->GetNonBendingCoor();
 624
 625       fThetaX = muonTrack->GetThetaX();
 626       fThetaY = muonTrack->GetThetaY();
 627
 628       fPyz = 1./TMath::Abs(muonTrack->GetInverseBendingMomentum());
 629       fPzRec = - fPyz / TMath::Sqrt(1.0 + TMath::Tan(fThetaY)*TMath::Tan(fThetaY));
 630       fPxRec = fPzRec * TMath::Tan(fThetaX);
 631       fPyRec = fPzRec * TMath::Tan(fThetaY);
 632       fCharge = Int_t(TMath::Sign(1.,muonTrack->GetInverseBendingMomentum()));
 633
 634       //ChiSquare of the track if needed
 635       fChisquare = muonTrack->GetChi2()/(2.0 * muonTrack->GetNHit() - 5);
 636       fEnergy = TMath::Sqrt(fMUONMASS * fMUONMASS + fPxRec * fPxRec + fPyRec * fPyRec + fPzRec * fPzRec);
 637       fEPvector.SetPxPyPzE(fPxRec, fPyRec, fPzRec, fEnergy);
 638
 639       // total number of muons inside a vertex cut
 640       if((TMath::Abs(fZ)<fZVertexCut) && (TMath::Sqrt(fY*fY+fX*fX)<fRhoVertexCut)) {
 641         nMuons++;
 642         if(fEPvector.Pt() > fLowPtCut) {
 643           nMu1GeV++;
 644           if(fEPvector.Pt() > fHighPtCut) {
 645             nMu3GeV++;
 646             if (fEPvector.Pt() > fVeryHighPtCut) {
 647               nMu10GeV++;
 648             }
 649           }
 650         }
 651       }
 652     }//muon track loop
 653
 654     // Fill the event tags
 655     if(ntrack != 0) meanPt = meanPt/ntrack;
 656
 657     evTag->SetEventId(iEventNumber+1);
 658     evTag->SetPath(filepath);
 659
 660     evTag->SetVertexX(vertexIn->GetXv());
 661     evTag->SetVertexY(vertexIn->GetYv());
 662     evTag->SetVertexZ(vertexIn->GetZv());
 663     evTag->SetVertexZError(vertexIn->GetZRes());
 664     evTag->SetVertexFlag(fVertexflag);
 665
 666     evTag->SetT0VertexZ(esd->GetT0zVertex());
 667
 668     evTag->SetTriggerMask(esd->GetTriggerMask());
 669     evTag->SetTriggerCluster(esd->GetTriggerCluster());
 670
 671     evTag->SetZDCNeutron1Energy(esd->GetZDCN1Energy());
 672     evTag->SetZDCProton1Energy(esd->GetZDCP1Energy());
 673     evTag->SetZDCEMEnergy(esd->GetZDCEMEnergy());
 674     evTag->SetZDCNeutron1Energy(esd->GetZDCN2Energy());
 675     evTag->SetZDCProton1Energy(esd->GetZDCP2Energy());
 676     evTag->SetNumOfParticipants(esd->GetZDCParticipants());
 677
 678
 679     evTag->SetNumOfTracks(esd->GetNumberOfTracks());
 680     evTag->SetNumOfPosTracks(nPos);
 681     evTag->SetNumOfNegTracks(nNeg);
 682     evTag->SetNumOfNeutrTracks(nNeutr);
 683
 684     evTag->SetNumOfV0s(esd->GetNumberOfV0s());
 685     evTag->SetNumOfCascades(esd->GetNumberOfCascades());
 686     evTag->SetNumOfKinks(esd->GetNumberOfKinks());
 687     evTag->SetNumOfPMDTracks(esd->GetNumberOfPmdTracks());
 688
 689     evTag->SetNumOfProtons(nProtons);
 690     evTag->SetNumOfKaons(nKaons);
 691     evTag->SetNumOfPions(nPions);
 692     evTag->SetNumOfMuons(nMuons);
 693     evTag->SetNumOfElectrons(nElectrons);
 694     evTag->SetNumOfPhotons(nGamas);
 695     evTag->SetNumOfPi0s(nPi0s);
 696     evTag->SetNumOfNeutrons(nNeutrons);
 697     evTag->SetNumOfKaon0s(nK0s);
 698
 699     evTag->SetNumOfChargedAbove1GeV(nCh1GeV);
 700     evTag->SetNumOfChargedAbove3GeV(nCh3GeV);
 701     evTag->SetNumOfChargedAbove10GeV(nCh10GeV);
 702     evTag->SetNumOfMuonsAbove1GeV(nMu1GeV);
 703     evTag->SetNumOfMuonsAbove3GeV(nMu3GeV);
 704     evTag->SetNumOfMuonsAbove10GeV(nMu10GeV);
 705     evTag->SetNumOfElectronsAbove1GeV(nEl1GeV);
 706     evTag->SetNumOfElectronsAbove3GeV(nEl3GeV);
 707     evTag->SetNumOfElectronsAbove10GeV(nEl10GeV);
 708
 709     evTag->SetNumOfPHOSClusters(esd->GetNumberOfPHOSClusters());
 710     evTag->SetNumOfEMCALClusters(esd->GetNumberOfEMCALClusters());
 711
 712     evTag->SetTotalMomentum(totalP);
 713     evTag->SetMeanPt(meanPt);
 714     evTag->SetMaxPt(maxPt);
 715
 716     tag->SetRunId(iInitRunNumber);
 717     if(fIsSim) tag->SetDataType(0);
 718     else tag->SetDataType(1);
 719     tag->AddEventTag(*evTag);
 720   }//event loop
 721   lastEvent = iNumberOfEvents;
 722
 723   t->Delete("");
 724
 725   ttag.Fill();
 726   tag->Clear();
 727
 728   TString localFileName = "Run"; localFileName += tag->GetRunId();
 729   localFileName += ".Event"; localFileName += firstEvent; localFileName += "_"; localFileName += lastEvent; localFileName += "."; localFileName += Counter;
 730   localFileName += ".ESD.tag.root";
 731
 732   TString fileName;
 733
 734   if(fStorage == 0) {
 735     fileName = localFileName.Data();
 736     AliInfo(Form("Writing tags to local file: %s",fileName.Data()));
 737   }
 738   else if(fStorage == 1) {
 739     TString alienLocation = "/alien";
 740     alienLocation += gGrid->Pwd();
 741     alienLocation += fgridpath.Data();
 742     alienLocation += "/";
 743     alienLocation +=  localFileName;
 744     alienLocation += "?se=";
 745     alienLocation += fSE.Data();
 746     fileName = alienLocation.Data();
 747     AliInfo(Form("Writing tags to grid file: %s",fileName.Data()));
 748   }
 749
 750   TFile* ftag = TFile::Open(fileName, "recreate");
 751   ftag->cd();
 752   ttag.Write();
 753   ftag->Close();
 754
 755   delete ftag;
 756   delete esd;
 757
 758   delete tag;
 759 }
 760
 761 //_____________________________________________________________________________
 762 void AliESDTagCreator::CreateESDTags(Int_t fFirstEvent, Int_t fLastEvent) {
 763   //GRP
 764   Float_t lhcLuminosity = 0.0;
 765   TString lhcState = "test";
 766   UInt_t detectorMask = 0;
 767
 768   /////////////
 769   //muon code//
 770   ////////////
 771   Double_t fMUONMASS = 0.105658369;
 772   //Variables
 773   Double_t fX,fY,fZ ;
 774   Double_t fThetaX, fThetaY, fPyz, fChisquare;
 775   Double_t fPxRec,fPyRec, fPzRec, fEnergy;
 776   Int_t fCharge;
 777   TLorentzVector fEPvector;
 778
 779   Float_t fZVertexCut = 10.0;
 780   Float_t fRhoVertexCut = 2.0;
 781
 782   Float_t fLowPtCut = 1.0;
 783   Float_t fHighPtCut = 3.0;
 784   Float_t fVeryHighPtCut = 10.0;
 785   ////////////
 786
 787   Double_t partFrac[5] = {0.01, 0.01, 0.85, 0.10, 0.05};
 788
 789   // Creates the tags for all the events in a given ESD file
 790   Int_t ntrack;
 791   Int_t nProtons, nKaons, nPions, nMuons, nElectrons;
 792   Int_t nPos, nNeg, nNeutr;
 793   Int_t nK0s, nNeutrons, nPi0s, nGamas;
 794   Int_t nCh1GeV, nCh3GeV, nCh10GeV;
 795   Int_t nMu1GeV, nMu3GeV, nMu10GeV;
 796   Int_t nEl1GeV, nEl3GeV, nEl10GeV;
 797   Float_t maxPt = .0, meanPt = .0, totalP = .0;
 798   Int_t fVertexflag;
 799   Int_t iRunNumber = 0;
 800   TString fVertexName("default");
 801
 802   AliRunTag *tag = new AliRunTag();
 803   AliEventTag *evTag = new AliEventTag();
 804   TTree ttag("T","A Tree with event tags");
 805   TBranch * btag = ttag.Branch("AliTAG", &tag);
 806   btag->SetCompressionLevel(9);
 807
 808   AliInfo(Form("Creating the tags......."));
 809
 810   TFile *file = TFile::Open("AliESDs.root");
 811   if (!file || !file->IsOpen()) {
 812     AliError(Form("opening failed"));
 813     delete file;
 814     return ;
 815   }
 816   Int_t lastEvent = 0;
 817   TTree *b = (TTree*) file->Get("esdTree");
 818   AliESDEvent *esd = new AliESDEvent();
 819   esd->ReadFromTree(b);
 820
 821   b->GetEntry(fFirstEvent);
 822   Int_t iInitRunNumber = esd->GetRunNumber();
 823
 824   Int_t iNumberOfEvents = (Int_t)b->GetEntries();
 825   if(fLastEvent != -1) iNumberOfEvents = fLastEvent + 1;
 826   for (Int_t iEventNumber = fFirstEvent; iEventNumber < iNumberOfEvents; iEventNumber++) {
 827     ntrack = 0;
 828     nPos = 0;
 829     nNeg = 0;
 830     nNeutr =0;
 831     nK0s = 0;
 832     nNeutrons = 0;
 833     nPi0s = 0;
 834     nGamas = 0;
 835     nProtons = 0;
 836     nKaons = 0;
 837     nPions = 0;
 838     nMuons = 0;
 839     nElectrons = 0;
 840     nCh1GeV = 0;
 841     nCh3GeV = 0;
 842     nCh10GeV = 0;
 843     nMu1GeV = 0;
 844     nMu3GeV = 0;
 845     nMu10GeV = 0;
 846     nEl1GeV = 0;
 847     nEl3GeV = 0;
 848     nEl10GeV = 0;
 849     maxPt = .0;
 850     meanPt = .0;
 851     totalP = .0;
 852     fVertexflag = 0;
 853
 854     b->GetEntry(iEventNumber);
 855     iRunNumber = esd->GetRunNumber();
 856     if(iRunNumber != iInitRunNumber) AliFatal("Inconsistency of run numbers in the AliESD!!!");
 857     const AliESDVertex * vertexIn = esd->GetVertex();
 858     if (!vertexIn) AliError("ESD has not defined vertex.");
 859     if (vertexIn) fVertexName = vertexIn->GetName();
 860     if(fVertexName != "default") fVertexflag = 1;
 861     for (Int_t iTrackNumber = 0; iTrackNumber < esd->GetNumberOfTracks(); iTrackNumber++) {
 862       AliESDtrack * esdTrack = esd->GetTrack(iTrackNumber);
 863       UInt_t status = esdTrack->GetStatus();
 864
 865       //select only tracks with ITS refit
 866       if ((status&AliESDtrack::kITSrefit)==0) continue;
 867       //select only tracks with TPC refit
 868       if ((status&AliESDtrack::kTPCrefit)==0) continue;
 869
 870       //select only tracks with the "combined PID"
 871       if ((status&AliESDtrack::kESDpid)==0) continue;
 872       Double_t p[3];
 873       esdTrack->GetPxPyPz(p);
 874       Double_t momentum = sqrt(pow(p[0],2) + pow(p[1],2) + pow(p[2],2));
 875       Double_t fPt = sqrt(pow(p[0],2) + pow(p[1],2));
 876       totalP += momentum;
 877       meanPt += fPt;
 878       if(fPt > maxPt) maxPt = fPt;
 879
 880       if(esdTrack->GetSign() > 0) {
 881         nPos++;
 882         if(fPt > fLowPtCut) nCh1GeV++;
 883         if(fPt > fHighPtCut) nCh3GeV++;
 884         if(fPt > fVeryHighPtCut) nCh10GeV++;
 885       }
 886       if(esdTrack->GetSign() < 0) {
 887         nNeg++;
 888         if(fPt > fLowPtCut) nCh1GeV++;
 889         if(fPt > fHighPtCut) nCh3GeV++;
 890         if(fPt > fVeryHighPtCut) nCh10GeV++;
 891       }
 892       if(esdTrack->GetSign() == 0) nNeutr++;
 893
 894       //PID
 895       Double_t prob[5];
 896       esdTrack->GetESDpid(prob);
 897
 898       Double_t rcc = 0.0;
 899       for(Int_t i = 0; i < AliPID::kSPECIES; i++) rcc += prob[i]*partFrac[i];
 900       if(rcc == 0.0) continue;
 901       //Bayes' formula
 902       Double_t w[5];
 903       for(Int_t i = 0; i < AliPID::kSPECIES; i++) w[i] = prob[i]*partFrac[i]/rcc;
 904
 905       //protons
 906       if ((w[4]>w[3])&&(w[4]>w[2])&&(w[4]>w[1])&&(w[4]>w[0])) nProtons++;
 907       //kaons
 908       if ((w[3]>w[4])&&(w[3]>w[2])&&(w[3]>w[1])&&(w[3]>w[0])) nKaons++;
 909       //pions
 910       if ((w[2]>w[4])&&(w[2]>w[3])&&(w[2]>w[1])&&(w[2]>w[0])) nPions++;
 911       //electrons
 912       if ((w[0]>w[4])&&(w[0]>w[3])&&(w[0]>w[2])&&(w[0]>w[1])) {
 913         nElectrons++;
 914         if(fPt > fLowPtCut) nEl1GeV++;
 915         if(fPt > fHighPtCut) nEl3GeV++;
 916         if(fPt > fVeryHighPtCut) nEl10GeV++;
 917       }
 918       ntrack++;
 919     }//track loop
 920
 921     /////////////
 922     //muon code//
 923     ////////////
 924     Int_t nMuonTracks = esd->GetNumberOfMuonTracks();
 925     // loop over all reconstructed tracks (also first track of combination)
 926     for (Int_t iTrack = 0; iTrack <  nMuonTracks;  iTrack++) {
 927       AliESDMuonTrack* muonTrack = esd->GetMuonTrack(iTrack);
 928       if (muonTrack == 0x0) continue;
 929
 930       // Coordinates at vertex
 931       fZ = muonTrack->GetZ();
 932       fY = muonTrack->GetBendingCoor();
 933       fX = muonTrack->GetNonBendingCoor();
 934
 935       fThetaX = muonTrack->GetThetaX();
 936       fThetaY = muonTrack->GetThetaY();
 937
 938       fPyz = 1./TMath::Abs(muonTrack->GetInverseBendingMomentum());
 939       fPzRec = - fPyz / TMath::Sqrt(1.0 + TMath::Tan(fThetaY)*TMath::Tan(fThetaY));
 940       fPxRec = fPzRec * TMath::Tan(fThetaX);
 941       fPyRec = fPzRec * TMath::Tan(fThetaY);
 942       fCharge = Int_t(TMath::Sign(1.,muonTrack->GetInverseBendingMomentum()));
 943
 944       //ChiSquare of the track if needed
 945       fChisquare = muonTrack->GetChi2()/(2.0 * muonTrack->GetNHit() - 5);
 946       fEnergy = TMath::Sqrt(fMUONMASS * fMUONMASS + fPxRec * fPxRec + fPyRec * fPyRec + fPzRec * fPzRec);
 947       fEPvector.SetPxPyPzE(fPxRec, fPyRec, fPzRec, fEnergy);
 948
 949       // total number of muons inside a vertex cut
 950       if((TMath::Abs(fZ)<fZVertexCut) && (TMath::Sqrt(fY*fY+fX*fX)<fRhoVertexCut)) {
 951         nMuons++;
 952         if(fEPvector.Pt() > fLowPtCut) {
 953           nMu1GeV++;
 954           if(fEPvector.Pt() > fHighPtCut) {
 955             nMu3GeV++;
 956             if (fEPvector.Pt() > fVeryHighPtCut) {
 957               nMu10GeV++;
 958             }
 959           }
 960         }
 961       }
 962     }//muon track loop
 963
 964     // Fill the event tags
 965     if(ntrack != 0)
 966       meanPt = meanPt/ntrack;
 967
 968     evTag->SetEventId(iEventNumber+1);
 969     if (vertexIn) {
 970       evTag->SetVertexX(vertexIn->GetXv());
 971       evTag->SetVertexY(vertexIn->GetYv());
 972       evTag->SetVertexZ(vertexIn->GetZv());
 973       evTag->SetVertexZError(vertexIn->GetZRes());
 974     }
 975     evTag->SetVertexFlag(fVertexflag);
 976
 977     evTag->SetT0VertexZ(esd->GetT0zVertex());
 978
 979     evTag->SetTriggerMask(esd->GetTriggerMask());
 980     evTag->SetTriggerCluster(esd->GetTriggerCluster());
 981
 982     evTag->SetZDCNeutron1Energy(esd->GetZDCN1Energy());
 983     evTag->SetZDCProton1Energy(esd->GetZDCP1Energy());
 984     evTag->SetZDCNeutron2Energy(esd->GetZDCN2Energy());
 985     evTag->SetZDCProton2Energy(esd->GetZDCP2Energy());
 986     evTag->SetZDCEMEnergy(esd->GetZDCEMEnergy());
 987     evTag->SetNumOfParticipants(esd->GetZDCParticipants());
 988
 989
 990     evTag->SetNumOfTracks(esd->GetNumberOfTracks());
 991     evTag->SetNumOfPosTracks(nPos);
 992     evTag->SetNumOfNegTracks(nNeg);
 993     evTag->SetNumOfNeutrTracks(nNeutr);
 994
 995     evTag->SetNumOfV0s(esd->GetNumberOfV0s());
 996     evTag->SetNumOfCascades(esd->GetNumberOfCascades());
 997     evTag->SetNumOfKinks(esd->GetNumberOfKinks());
 998     evTag->SetNumOfPMDTracks(esd->GetNumberOfPmdTracks());
 999
1000     evTag->SetNumOfProtons(nProtons);
1001     evTag->SetNumOfKaons(nKaons);
1002     evTag->SetNumOfPions(nPions);
1003     evTag->SetNumOfMuons(nMuons);
1004     evTag->SetNumOfElectrons(nElectrons);
1005     evTag->SetNumOfPhotons(nGamas);
1006     evTag->SetNumOfPi0s(nPi0s);
1007     evTag->SetNumOfNeutrons(nNeutrons);
1008     evTag->SetNumOfKaon0s(nK0s);
1009
1010     evTag->SetNumOfChargedAbove1GeV(nCh1GeV);
1011     evTag->SetNumOfChargedAbove3GeV(nCh3GeV);
1012     evTag->SetNumOfChargedAbove10GeV(nCh10GeV);
1013     evTag->SetNumOfMuonsAbove1GeV(nMu1GeV);
1014     evTag->SetNumOfMuonsAbove3GeV(nMu3GeV);
1015     evTag->SetNumOfMuonsAbove10GeV(nMu10GeV);
1016     evTag->SetNumOfElectronsAbove1GeV(nEl1GeV);
1017     evTag->SetNumOfElectronsAbove3GeV(nEl3GeV);
1018     evTag->SetNumOfElectronsAbove10GeV(nEl10GeV);
1019
1020     evTag->SetNumOfPHOSClusters(esd->GetNumberOfPHOSClusters());
1021     evTag->SetNumOfEMCALClusters(esd->GetNumberOfEMCALClusters());
1022
1023     evTag->SetTotalMomentum(totalP);
1024     evTag->SetMeanPt(meanPt);
1025     evTag->SetMaxPt(maxPt);
1026
1027     tag->SetLHCTag(lhcLuminosity,lhcState);
1028     tag->SetDetectorTag(detectorMask);
1029
1030     tag->SetRunId(iInitRunNumber);
1031     tag->AddEventTag(*evTag);
1032   }
1033   if(fLastEvent == -1) lastEvent = (Int_t)b->GetEntries();
1034   else lastEvent = fLastEvent;
1035
1036   ttag.Fill();
1037   tag->Clear();
1038
1039   char fileName[256];
1040   sprintf(fileName, "Run%d.Event%d_%d.ESD.tag.root",
1041           tag->GetRunId(),fFirstEvent,lastEvent );
1042   AliInfo(Form("writing tags to file %s", fileName));
1043   AliDebug(1, Form("writing tags to file %s", fileName));
1044
1045   TFile* ftag = TFile::Open(fileName, "recreate");
1046   ftag->cd();
1047   ttag.Write();
1048   ftag->Close();
1049   file->cd();
1050   delete tag;
1051   delete evTag;
1052 }
1053
1054 //__________________________________________________________________________
1055 /*void AliTagCreator::CreateAODTags(Int_t fFirstEvent, Int_t fLastEvent) {
1056   //creates tag files for AODs
1057
1058   Float_t fLowPtCut = 1.0;
1059   Float_t fHighPtCut = 3.0;
1060   Float_t fVeryHighPtCut = 10.0;
1061   ////////////
1062
1063   Double_t partFrac[10] = {0.01, 0.01, 0.85, 0.10, 0.05, 0., 0., 0., 0., 0.};
1064
1065   // Creates the tags for all the events in a given ESD file
1066   Int_t ntrack;
1067   Int_t nProtons, nKaons, nPions, nMuons, nElectrons;
1068   Int_t nPos, nNeg, nNeutr;
1069   Int_t nKinks, nV0s, nCascades;
1070   Int_t nK0s, nNeutrons, nPi0s, nGamas;
1071   Int_t nCh1GeV, nCh3GeV, nCh10GeV;
1072   Int_t nMu1GeV, nMu3GeV, nMu10GeV;
1073   Int_t nEl1GeV, nEl3GeV, nEl10GeV;
1074   Float_t maxPt = .0, meanPt = .0, totalP = .0;
1075
1076   AliRunTag *tag = new AliRunTag();
1077   TTree ttag("T","A Tree with event tags");
1078   TBranch * btag = ttag.Branch("AliTAG", &tag);
1079   btag->SetCompressionLevel(9);
1080
1081   //reading the esd tag file
1082   TChain *oldTagTree = new TChain("T");
1083   const char * tagPattern = "ESD.tag";
1084   // Open the working directory
1085   void * dirp = gSystem->OpenDirectory(gSystem->pwd());
1086   const char * name = 0x0;
1087   // Add all files matching *pattern* to the chain
1088   while((name = gSystem->GetDirEntry(dirp))) {
1089     if (strstr(name,tagPattern)) oldTagTree->Add(name);
1090   }//directory loop
1091   AliInfo(Form("Chained tag files: %d",oldTagTree->GetEntries()));
1092
1093   //reading the esd tag file
1094   AliRunTag *oldtag = new AliRunTag();
1095   TString tagFilename;
1096   oldTagTree->SetBranchAddress("AliTAG",&oldtag);
1097   oldTagTree->GetEntry(0);
1098   tag->CopyStandardContent(oldtag);
1099   const TClonesArray *evTagList = oldtag->GetEventTags();
1100
1101   AliInfo(Form("Creating the tags......."));
1102
1103   TFile *file = TFile::Open("AliAOD.root");
1104   if (!file || !file->IsOpen()) {
1105     AliError(Form("opening failed"));
1106     delete file;
1107     return ;
1108   }
1109   TTree *aodTree = (TTree*)file->Get("aodTree");
1110   AliAODEvent *aod = (AliAODEvent*)aodTree->GetUserInfo()->FindObject("AliAODEvent");
1111   TIter next(aod->GetList());
1112   TObject *el;
1113
1114   Int_t lastEvent = 0;
1115   if(fLastEvent == -1) lastEvent = (Int_t)aodTree->GetEntries();
1116   else lastEvent = fLastEvent;
1117
1118   while((el=(TNamed*)next()))
1119     aodTree->SetBranchAddress(el->GetName(),aod->GetList()->GetObjectRef(el));
1120
1121   // loop over events
1122   Int_t nEvents = aodTree->GetEntries();
1123   for (Int_t iEventNumber = 0; iEventNumber < nEvents; iEventNumber++) {
1124     AliEventTag *evTag = (AliEventTag *)evTagList->At(iEventNumber);
1125     ntrack = 0;
1126     nPos = 0; nNeg = 0; nNeutr =0;
1127     nKinks = 0; nV0s = 0; nCascades = 0;
1128     nK0s = 0; nNeutrons = 0; nPi0s = 0; nGamas = 0;
1129     nProtons = 0;  nKaons = 0; nPions = 0; nMuons = 0; nElectrons = 0;
1130     nCh1GeV = 0; nCh3GeV = 0; nCh10GeV = 0;
1131     nMu1GeV = 0; nMu3GeV = 0; nMu10GeV = 0;
1132     nEl1GeV = 0; nEl3GeV = 0; nEl10GeV = 0;
1133     maxPt = .0; meanPt = .0; totalP = .0;
1134
1135     // read events
1136     aodTree->GetEvent(iEventNumber);
1137
1138     // set pointers
1139     aod->GetStdContent();
1140
1141     Int_t nTracks = aod->GetNTracks();
1142     // loop over vertices
1143     Int_t nVtxs = aod->GetNVertices();
1144     for (Int_t nVtx = 0; nVtx < nVtxs; nVtx++) {
1145       // print track info
1146       AliAODVertex *vertex = aod->GetVertex(nVtx);
1147       if(vertex->GetType() == 1) nKinks += 1;
1148       if(vertex->GetType() == 2) nV0s += 1;
1149       if(vertex->GetType() == 3) nCascades += 1;
1150     }
1151     for (Int_t nTr = 0; nTr < nTracks; nTr++) {
1152       AliAODTrack *track = aod->GetTrack(nTr);
1153
1154       Double_t fPt = track->Pt();
1155       if(fPt > maxPt) maxPt = fPt;
1156       if(track->Charge() > 0) {
1157         nPos++;
1158         if(fPt > fLowPtCut) nCh1GeV++;
1159         if(fPt > fHighPtCut) nCh3GeV++;
1160         if(fPt > fVeryHighPtCut) nCh10GeV++;
1161       }
1162       if(track->Charge() < 0) {
1163         nNeg++;
1164         if(fPt > fLowPtCut) nCh1GeV++;
1165         if(fPt > fHighPtCut) nCh3GeV++;
1166         if(fPt > fVeryHighPtCut) nCh10GeV++;
1167       }
1168       if(track->Charge() == 0) nNeutr++;
1169
1170       //PID
1171       const Double32_t *prob = track->PID();
1172       Double_t rcc = 0.0;
1173       for(Int_t i = 0; i < AliPID::kSPECIES; i++) rcc += prob[i]*partFrac[i];
1174       if(rcc == 0.0) continue;
1175       //Bayes' formula
1176       Double_t w[10];
1177       for(Int_t i = 0; i < AliPID::kSPECIES; i++) w[i] = prob[i]*partFrac[i]/rcc;
1178
1179       //protons
1180       if ((w[4]>w[3])&&(w[4]>w[2])&&(w[4]>w[1])&&(w[4]>w[0])) nProtons++;
1181       //kaons
1182       if ((w[3]>w[4])&&(w[3]>w[2])&&(w[3]>w[1])&&(w[3]>w[0])) nKaons++;
1183       //pions
1184       if ((w[2]>w[4])&&(w[2]>w[3])&&(w[2]>w[1])&&(w[2]>w[0])) nPions++;
1185       //muons
1186       if ((w[1]>w[4])&&(w[1]>w[3])&&(w[1]>w[2])&&(w[1]>w[0])) {
1187         nMuons++;
1188         if(fPt > fLowPtCut) nMu1GeV++;
1189         if(fPt > fHighPtCut) nMu3GeV++;
1190         if(fPt > fVeryHighPtCut) nMu10GeV++;
1191         }
1192       //electrons
1193       if ((w[0]>w[4])&&(w[0]>w[3])&&(w[0]>w[2])&&(w[0]>w[1])) {
1194         nElectrons++;
1195         if(fPt > fLowPtCut) nEl1GeV++;
1196         if(fPt > fHighPtCut) nEl3GeV++;
1197         if(fPt > fVeryHighPtCut) nEl10GeV++;
1198       }
1199
1200       totalP += track->P();
1201       meanPt += fPt;
1202       ntrack++;
1203     }//track loop
1204     // Fill the event tags
1205     if(ntrack != 0)
1206       meanPt = meanPt/ntrack;
1207
1208     evTag->SetEventId(iEventNumber+1);
1209
1210     evTag->SetNumOfTracks(nTracks);
1211     evTag->SetNumOfPosTracks(nPos);
1212     evTag->SetNumOfNegTracks(nNeg);
1213     evTag->SetNumOfNeutrTracks(nNeutr);
1214
1215     evTag->SetNumOfV0s(nV0s);
1216     evTag->SetNumOfCascades(nCascades);
1217     evTag->SetNumOfKinks(nKinks);
1218
1219     evTag->SetNumOfProtons(nProtons);
1220     evTag->SetNumOfKaons(nKaons);
1221     evTag->SetNumOfPions(nPions);
1222     evTag->SetNumOfMuons(nMuons);
1223     evTag->SetNumOfElectrons(nElectrons);
1224     evTag->SetNumOfPhotons(nGamas);
1225     evTag->SetNumOfPi0s(nPi0s);
1226     evTag->SetNumOfNeutrons(nNeutrons);
1227     evTag->SetNumOfKaon0s(nK0s);
1228
1229     evTag->SetNumOfChargedAbove1GeV(nCh1GeV);
1230     evTag->SetNumOfChargedAbove3GeV(nCh3GeV);
1231     evTag->SetNumOfChargedAbove10GeV(nCh10GeV);
1232     evTag->SetNumOfMuonsAbove1GeV(nMu1GeV);
1233     evTag->SetNumOfMuonsAbove3GeV(nMu3GeV);
1234     evTag->SetNumOfMuonsAbove10GeV(nMu10GeV);
1235     evTag->SetNumOfElectronsAbove1GeV(nEl1GeV);
1236     evTag->SetNumOfElectronsAbove3GeV(nEl3GeV);
1237     evTag->SetNumOfElectronsAbove10GeV(nEl10GeV);
1238
1239     evTag->SetTotalMomentum(totalP);
1240     evTag->SetMeanPt(meanPt);
1241     evTag->SetMaxPt(maxPt);
1242     tag->AddEventTag(*evTag);
1243   }//event loop
1244   if(fLastEvent == -1) lastEvent = (Int_t)aodTree->GetEntries();
1245   else lastEvent = fLastEvent;
1246
1247   ttag.Fill();
1248   tag->Clear();
1249
1250   char fileName[256];
1251   sprintf(fileName, "Run%d.Event%d_%d.AOD.tag.root",
1252           tag->GetRunId(),fFirstEvent,lastEvent );
1253   AliInfo(Form("writing tags to file %s", fileName));
1254   AliDebug(1, Form("writing tags to file %s", fileName));
1255
1256   TFile* ftag = TFile::Open(fileName, "recreate");
1257   ftag->cd();
1258   ttag.Write();
1259   ftag->Close();
1260   file->cd();
1261   file->Close();
1262   }*/