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