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