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