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, nFWMuons;
 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; nFWMuons = 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 pt2 = p[0]*p[0]+p[1]*p[1];
 333       Double_t momentum = TMath::Sqrt(pt2+p[2]*p[2]);
 334       Double_t fPt = TMath::Sqrt(pt2);
 335       totalP += momentum;
 336       meanPt += fPt;
 337       if(fPt > maxPt) maxPt = fPt;
 338
 339       if(esdTrack->GetSign() > 0) {
 340         nPos++;
 341         if(fPt > fLowPtCut) nCh1GeV++;
 342         if(fPt > fHighPtCut) nCh3GeV++;
 343         if(fPt > fVeryHighPtCut) nCh10GeV++;
 344       }
 345       if(esdTrack->GetSign() < 0) {
 346         nNeg++;
 347         if(fPt > fLowPtCut) nCh1GeV++;
 348         if(fPt > fHighPtCut) nCh3GeV++;
 349         if(fPt > fVeryHighPtCut) nCh10GeV++;
 350       }
 351       if(esdTrack->GetSign() == 0) nNeutr++;
 352
 353       //PID
 354       Double_t prob[5];
 355       esdTrack->GetESDpid(prob);
 356
 357       Double_t rcc = 0.0;
 358       for(Int_t i = 0; i < AliPID::kSPECIES; i++) rcc += prob[i]*partFrac[i];
 359       if(rcc == 0.0) continue;
 360       //Bayes' formula
 361       Double_t w[5];
 362       for(Int_t i = 0; i < AliPID::kSPECIES; i++) w[i] = prob[i]*partFrac[i]/rcc;
 363
 364       //protons
 365       if ((w[4]>w[3])&&(w[4]>w[2])&&(w[4]>w[1])&&(w[4]>w[0])) nProtons++;
 366       //kaons
 367       if ((w[3]>w[4])&&(w[3]>w[2])&&(w[3]>w[1])&&(w[3]>w[0])) nKaons++;
 368       //pions
 369       if ((w[2]>w[4])&&(w[2]>w[3])&&(w[2]>w[1])&&(w[2]>w[0])) nPions++;
 370       //electrons
 371       if ((w[0]>w[4])&&(w[0]>w[3])&&(w[0]>w[2])&&(w[0]>w[1])) {
 372         nElectrons++;
 373         if(fPt > fLowPtCut) nEl1GeV++;
 374         if(fPt > fHighPtCut) nEl3GeV++;
 375         if(fPt > fVeryHighPtCut) nEl10GeV++;
 376       }
 377       ntrack++;
 378     }//esd track loop
 379
 380     /////////////
 381     //muon code//
 382     ////////////
 383     Int_t nMuonTracks = esd->GetNumberOfMuonTracks();
 384     // loop over all reconstructed tracks (also first track of combination)
 385     for (Int_t iTrack = 0; iTrack <  nMuonTracks;  iTrack++) {
 386       AliESDMuonTrack* muonTrack = esd->GetMuonTrack(iTrack);
 387       if (muonTrack == 0x0) continue;
 388
 389       // Coordinates at vertex
 390       fZ = muonTrack->GetZ();
 391       fY = muonTrack->GetBendingCoor();
 392       fX = muonTrack->GetNonBendingCoor();
 393
 394       fThetaX = muonTrack->GetThetaX();
 395       fThetaY = muonTrack->GetThetaY();
 396
 397       fPyz = 1./TMath::Abs(muonTrack->GetInverseBendingMomentum());
 398       fPzRec = - fPyz / TMath::Sqrt(1.0 + TMath::Tan(fThetaY)*TMath::Tan(fThetaY));
 399       fPxRec = fPzRec * TMath::Tan(fThetaX);
 400       fPyRec = fPzRec * TMath::Tan(fThetaY);
 401       fCharge = Int_t(TMath::Sign(1.,muonTrack->GetInverseBendingMomentum()));
 402
 403       //ChiSquare of the track if needed
 404       fChisquare = muonTrack->GetChi2()/(2.0 * muonTrack->GetNHit() - 5);
 405       fEnergy = TMath::Sqrt(fMUONMASS * fMUONMASS + fPxRec * fPxRec + fPyRec * fPyRec + fPzRec * fPzRec);
 406       fEPvector.SetPxPyPzE(fPxRec, fPyRec, fPzRec, fEnergy);
 407
 408       // total number of muons inside a vertex cut
 409       if((TMath::Abs(fZ)<fZVertexCut) && (TMath::Sqrt(fY*fY+fX*fX)<fRhoVertexCut)) {
 410         nMuons++;
 411         nFWMuons++;
 412         if(fEPvector.Pt() > fLowPtCut) {
 413           nMu1GeV++;
 414           if(fEPvector.Pt() > fHighPtCut) {
 415             nMu3GeV++;
 416             if (fEPvector.Pt() > fVeryHighPtCut) {
 417               nMu10GeV++;
 418             }
 419           }
 420         }
 421       }
 422     }//muon track loop
 423
 424     // Fill the event tags
 425     if(ntrack != 0) meanPt = meanPt/ntrack;
 426
 427     //AliInfo(Form("====================================="));
 428     //AliInfo(Form("URL: %s - GUID: %s",fturl.Data(),fguid.Data()));
 429     //AliInfo(Form("====================================="));
 430
 431     evTag->SetEventId(iEventNumber+1);
 432     evTag->SetGUID(fguid);
 433     if(fSession == "grid") {
 434       evTag->SetMD5(0);
 435       evTag->SetTURL(fturl);
 436       evTag->SetSize(0);
 437     }
 438     else evTag->SetPath(fturl);
 439
 440     evTag->SetVertexX(vertexIn->GetXv());
 441     evTag->SetVertexY(vertexIn->GetYv());
 442     evTag->SetVertexZ(vertexIn->GetZv());
 443     evTag->SetVertexZError(vertexIn->GetZRes());
 444     evTag->SetVertexFlag(fVertexflag);
 445
 446     evTag->SetT0VertexZ(esd->GetT0zVertex());
 447
 448     evTag->SetTriggerMask(esd->GetTriggerMask());
 449     evTag->SetTriggerCluster(esd->GetTriggerCluster());
 450
 451     evTag->SetZDCNeutron1Energy(esd->GetZDCN1Energy());
 452     evTag->SetZDCProton1Energy(esd->GetZDCP1Energy());
 453     evTag->SetZDCEMEnergy(esd->GetZDCEMEnergy(0),esd->GetZDCEMEnergy(1));
 454     evTag->SetZDCNeutron1Energy(esd->GetZDCN2Energy());
 455     evTag->SetZDCProton1Energy(esd->GetZDCP2Energy());
 456     evTag->SetNumOfParticipants(esd->GetZDCParticipants());
 457     evTag->SetNumOfParticipants(esd->GetZDCParticipants());
 458     evTag->SetNumOfParticipants2(esd->GetZDCParticipants2());
 459
 460
 461     evTag->SetNumOfTracks(esd->GetNumberOfTracks());
 462     evTag->SetNumOfPosTracks(nPos);
 463     evTag->SetNumOfNegTracks(nNeg);
 464     evTag->SetNumOfNeutrTracks(nNeutr);
 465
 466     evTag->SetNumOfV0s(esd->GetNumberOfV0s());
 467     evTag->SetNumOfCascades(esd->GetNumberOfCascades());
 468     evTag->SetNumOfKinks(esd->GetNumberOfKinks());
 469     evTag->SetNumOfPMDTracks(esd->GetNumberOfPmdTracks());
 470
 471     evTag->SetNumOfProtons(nProtons);
 472     evTag->SetNumOfKaons(nKaons);
 473     evTag->SetNumOfPions(nPions);
 474     evTag->SetNumOfMuons(nMuons);
 475     evTag->SetNumOfFWMuons(nFWMuons);
 476     evTag->SetNumOfElectrons(nElectrons);
 477     evTag->SetNumOfPhotons(nGamas);
 478     evTag->SetNumOfPi0s(nPi0s);
 479     evTag->SetNumOfNeutrons(nNeutrons);
 480     evTag->SetNumOfKaon0s(nK0s);
 481
 482     evTag->SetNumOfChargedAbove1GeV(nCh1GeV);
 483     evTag->SetNumOfChargedAbove3GeV(nCh3GeV);
 484     evTag->SetNumOfChargedAbove10GeV(nCh10GeV);
 485     evTag->SetNumOfMuonsAbove1GeV(nMu1GeV);
 486     evTag->SetNumOfMuonsAbove3GeV(nMu3GeV);
 487     evTag->SetNumOfMuonsAbove10GeV(nMu10GeV);
 488     evTag->SetNumOfElectronsAbove1GeV(nEl1GeV);
 489     evTag->SetNumOfElectronsAbove3GeV(nEl3GeV);
 490     evTag->SetNumOfElectronsAbove10GeV(nEl10GeV);
 491
 492     evTag->SetNumOfPHOSClusters(esd->GetNumberOfPHOSClusters());
 493     evTag->SetNumOfEMCALClusters(esd->GetNumberOfEMCALClusters());
 494
 495     evTag->SetTotalMomentum(totalP);
 496     evTag->SetMeanPt(meanPt);
 497     evTag->SetMaxPt(maxPt);
 498
 499     tag->SetRunId(iInitRunNumber);
 500     if(fIsSim) tag->SetDataType(0);
 501     else tag->SetDataType(1);
 502
 503     if(fguid != fTempGuid) {
 504       fTempGuid = fguid;
 505       ttag.Fill();
 506       tag->Clear("");
 507     }
 508     tag->AddEventTag(*evTag);
 509     if(iEventNumber+1 == chain->GetEntries()) {
 510       //AliInfo(Form("File: %s",fturl.Data()));
 511       ttag.Fill();
 512       tag->Clear("");
 513     }
 514   }//event loop
 515   lastEvent = chain->GetEntries();
 516
 517   //gSystem->GetMemInfo(meminfo);
 518   //AliInfo(Form("After the event and track loop - Memory used: %d MB - Increase: %d MB",meminfo->fMemUsed,meminfo->fMemUsed - tempmem));
 519   //tempmem = meminfo->fMemUsed;
 520
 521   //chain->Delete("");
 522
 523   //gSystem->GetMemInfo(meminfo);
 524   //AliInfo(Form("After the t->Delete - Memory used: %d MB - Increase: %d MB",meminfo->fMemUsed,meminfo->fMemUsed - tempmem));
 525   //tempmem = meminfo->fMemUsed;
 526
 527   ftag->cd();
 528   tag->Clear();
 529   ttag.Write();
 530   ftag->Close();
 531
 532   //gSystem->GetMemInfo(meminfo);
 533   //AliInfo(Form("After the file closing - Memory used: %d MB - Increase: %d MB",meminfo->fMemUsed,meminfo->fMemUsed - tempmem));
 534   //tempmem = meminfo->fMemUsed;
 535
 536   delete esd;
 537   delete tag;
 538
 539   //gSystem->GetMemInfo(meminfo);
 540   //AliInfo(Form("After the delete objects - Memory used: %d MB - Increase: %d MB",meminfo->fMemUsed,meminfo->fMemUsed - tempmem));
 541 }
 542
 543 //_____________________________________________________________________________
 544 void AliESDTagCreator::CreateTag(TFile* file, const char *guid, const char *md5, const char *turl, Long64_t size, Int_t Counter) {
 545   //private method that creates tag files
 546   TString fguid = guid;
 547   TString fmd5 = md5;
 548   TString fturl = turl;
 549   /////////////
 550   //muon code//
 551   ////////////
 552   Double_t fMUONMASS = 0.105658369;
 553   //Variables
 554   Double_t fX,fY,fZ ;
 555   Double_t fThetaX, fThetaY, fPyz, fChisquare;
 556   Double_t fPxRec, fPyRec, fPzRec, fEnergy;
 557   Int_t fCharge;
 558   TLorentzVector fEPvector;
 559
 560   Float_t fZVertexCut = 10.0;
 561   Float_t fRhoVertexCut = 2.0;
 562
 563   Float_t fLowPtCut = 1.0;
 564   Float_t fHighPtCut = 3.0;
 565   Float_t fVeryHighPtCut = 10.0;
 566   ////////////
 567
 568   Double_t partFrac[5] = {0.01, 0.01, 0.85, 0.10, 0.05};
 569
 570   // Creates the tags for all the events in a given ESD file
 571   Bool_t fIsSim = kTRUE;
 572   Int_t ntrack;
 573   Int_t nProtons, nKaons, nPions, nMuons, nElectrons, nFWMuons;
 574   Int_t nPos, nNeg, nNeutr;
 575   Int_t nK0s, nNeutrons, nPi0s, nGamas;
 576   Int_t nCh1GeV, nCh3GeV, nCh10GeV;
 577   Int_t nMu1GeV, nMu3GeV, nMu10GeV;
 578   Int_t nEl1GeV, nEl3GeV, nEl10GeV;
 579   Float_t maxPt = .0, meanPt = .0, totalP = .0;
 580   Int_t fVertexflag;
 581   Int_t iRunNumber = 0;
 582   TString fVertexName;
 583
 584   AliRunTag *tag = new AliRunTag();
 585   AliEventTag *evTag = new AliEventTag();
 586   TTree ttag("T","A Tree with event tags");
 587   TBranch * btag = ttag.Branch("AliTAG", &tag);
 588   btag->SetCompressionLevel(9);
 589   gSystem->GetMemInfo(meminfo);
 590   AliInfo(Form("After the tag initialization - Memory used: %d MB",meminfo->fMemUsed));
 591   Int_t tempmem = meminfo->fMemUsed;
 592
 593   AliInfo(Form("Creating the ESD tags......."));
 594
 595   Int_t firstEvent = 0,lastEvent = 0;
 596   TTree *t = (TTree*) file->Get("esdTree");
 597   AliESDEvent *esd = new AliESDEvent();
 598   esd->ReadFromTree(t);
 599
 600   gSystem->GetMemInfo(meminfo);
 601   AliInfo(Form("After the esd initialization - Memory used: %d MB - Increase: %d MB",meminfo->fMemUsed,meminfo->fMemUsed - tempmem));
 602   tempmem = meminfo->fMemUsed;
 603
 604   t->GetEntry(0);
 605   Int_t iInitRunNumber = esd->GetRunNumber();
 606
 607   Int_t iNumberOfEvents = (Int_t)t->GetEntries();
 608   for (Int_t iEventNumber = 0; iEventNumber < iNumberOfEvents; iEventNumber++) {
 609     ntrack = 0;
 610     nPos = 0;
 611     nNeg = 0;
 612     nNeutr =0;
 613     nK0s = 0;
 614     nNeutrons = 0;
 615     nPi0s = 0;
 616     nGamas = 0;
 617     nProtons = 0;
 618     nKaons = 0;
 619     nPions = 0;
 620     nMuons = 0;
 621     nFWMuons = 0;
 622     nElectrons = 0;
 623     nCh1GeV = 0;
 624     nCh3GeV = 0;
 625     nCh10GeV = 0;
 626     nMu1GeV = 0;
 627     nMu3GeV = 0;
 628     nMu10GeV = 0;
 629     nEl1GeV = 0;
 630     nEl3GeV = 0;
 631     nEl10GeV = 0;
 632     maxPt = .0;
 633     meanPt = .0;
 634     totalP = .0;
 635     fVertexflag = 1;
 636
 637     t->GetEntry(iEventNumber);
 638     iRunNumber = esd->GetRunNumber();
 639     if(iRunNumber != iInitRunNumber) AliFatal("Inconsistency of run numbers in the AliESD!!!");
 640     const AliESDVertex * vertexIn = esd->GetVertex();
 641     fVertexName = vertexIn->GetName();
 642     if(fVertexName == "default") fVertexflag = 0;
 643
 644     for (Int_t iTrackNumber = 0; iTrackNumber < esd->GetNumberOfTracks(); iTrackNumber++) {
 645       AliESDtrack * esdTrack = esd->GetTrack(iTrackNumber);
 646       if(esdTrack->GetLabel() != 0) fIsSim = kTRUE;
 647       else if(esdTrack->GetLabel() == 0) fIsSim = kFALSE;
 648       UInt_t status = esdTrack->GetStatus();
 649
 650       //select only tracks with ITS refit
 651       if ((status&AliESDtrack::kITSrefit)==0) continue;
 652       //select only tracks with TPC refit
 653       if ((status&AliESDtrack::kTPCrefit)==0) continue;
 654
 655       //select only tracks with the "combined PID"
 656       if ((status&AliESDtrack::kESDpid)==0) continue;
 657       Double_t p[3];
 658       esdTrack->GetPxPyPz(p);
 659       Double_t pt2 = p[0]*p[0]+p[1]*p[1];
 660       Double_t momentum = TMath::Sqrt(pt2+p[2]*p[2]);
 661       Double_t fPt = TMath::Sqrt(pt2);
 662       totalP += momentum;
 663       meanPt += fPt;
 664       if(fPt > maxPt) maxPt = fPt;
 665
 666       if(esdTrack->GetSign() > 0) {
 667         nPos++;
 668         if(fPt > fLowPtCut) nCh1GeV++;
 669         if(fPt > fHighPtCut) nCh3GeV++;
 670         if(fPt > fVeryHighPtCut) nCh10GeV++;
 671       }
 672       if(esdTrack->GetSign() < 0) {
 673         nNeg++;
 674         if(fPt > fLowPtCut) nCh1GeV++;
 675         if(fPt > fHighPtCut) nCh3GeV++;
 676         if(fPt > fVeryHighPtCut) nCh10GeV++;
 677       }
 678       if(esdTrack->GetSign() == 0) nNeutr++;
 679
 680       //PID
 681       Double_t prob[5];
 682       esdTrack->GetESDpid(prob);
 683
 684       Double_t rcc = 0.0;
 685       for(Int_t i = 0; i < AliPID::kSPECIES; i++) rcc += prob[i]*partFrac[i];
 686       if(rcc == 0.0) continue;
 687       //Bayes' formula
 688       Double_t w[5];
 689       for(Int_t i = 0; i < AliPID::kSPECIES; i++) w[i] = prob[i]*partFrac[i]/rcc;
 690
 691       //protons
 692       if ((w[4]>w[3])&&(w[4]>w[2])&&(w[4]>w[1])&&(w[4]>w[0])) nProtons++;
 693       //kaons
 694       if ((w[3]>w[4])&&(w[3]>w[2])&&(w[3]>w[1])&&(w[3]>w[0])) nKaons++;
 695       //pions
 696       if ((w[2]>w[4])&&(w[2]>w[3])&&(w[2]>w[1])&&(w[2]>w[0])) nPions++;
 697       //electrons
 698       if ((w[0]>w[4])&&(w[0]>w[3])&&(w[0]>w[2])&&(w[0]>w[1])) {
 699         nElectrons++;
 700         if(fPt > fLowPtCut) nEl1GeV++;
 701         if(fPt > fHighPtCut) nEl3GeV++;
 702         if(fPt > fVeryHighPtCut) nEl10GeV++;
 703       }
 704       ntrack++;
 705     }//esd track loop
 706
 707     /////////////
 708     //muon code//
 709     ////////////
 710     Int_t nMuonTracks = esd->GetNumberOfMuonTracks();
 711     // loop over all reconstructed tracks (also first track of combination)
 712     for (Int_t iTrack = 0; iTrack <  nMuonTracks;  iTrack++) {
 713       AliESDMuonTrack* muonTrack = esd->GetMuonTrack(iTrack);
 714       if (muonTrack == 0x0) continue;
 715
 716       // Coordinates at vertex
 717       fZ = muonTrack->GetZ();
 718       fY = muonTrack->GetBendingCoor();
 719       fX = muonTrack->GetNonBendingCoor();
 720
 721       fThetaX = muonTrack->GetThetaX();
 722       fThetaY = muonTrack->GetThetaY();
 723
 724       fPyz = 1./TMath::Abs(muonTrack->GetInverseBendingMomentum());
 725       fPzRec = - fPyz / TMath::Sqrt(1.0 + TMath::Tan(fThetaY)*TMath::Tan(fThetaY));
 726       fPxRec = fPzRec * TMath::Tan(fThetaX);
 727       fPyRec = fPzRec * TMath::Tan(fThetaY);
 728       fCharge = Int_t(TMath::Sign(1.,muonTrack->GetInverseBendingMomentum()));
 729
 730       //ChiSquare of the track if needed
 731       fChisquare = muonTrack->GetChi2()/(2.0 * muonTrack->GetNHit() - 5);
 732       fEnergy = TMath::Sqrt(fMUONMASS * fMUONMASS + fPxRec * fPxRec + fPyRec * fPyRec + fPzRec * fPzRec);
 733       fEPvector.SetPxPyPzE(fPxRec, fPyRec, fPzRec, fEnergy);
 734
 735       // total number of muons inside a vertex cut
 736       if((TMath::Abs(fZ)<fZVertexCut) && (TMath::Sqrt(fY*fY+fX*fX)<fRhoVertexCut)) {
 737         nMuons++;
 738         nFWMuons++;
 739         if(fEPvector.Pt() > fLowPtCut) {
 740           nMu1GeV++;
 741           if(fEPvector.Pt() > fHighPtCut) {
 742             nMu3GeV++;
 743             if (fEPvector.Pt() > fVeryHighPtCut) {
 744               nMu10GeV++;
 745             }
 746           }
 747         }
 748       }
 749     }//muon track loop
 750
 751     // Fill the event tags
 752     if(ntrack != 0) meanPt = meanPt/ntrack;
 753
 754     evTag->SetEventId(iEventNumber+1);
 755     evTag->SetGUID(fguid);
 756     evTag->SetMD5(fmd5);
 757     evTag->SetTURL(fturl);
 758     evTag->SetSize(size);
 759     evTag->SetVertexX(vertexIn->GetXv());
 760     evTag->SetVertexY(vertexIn->GetYv());
 761     evTag->SetVertexZ(vertexIn->GetZv());
 762     evTag->SetVertexZError(vertexIn->GetZRes());
 763     evTag->SetVertexFlag(fVertexflag);
 764
 765     evTag->SetT0VertexZ(esd->GetT0zVertex());
 766
 767     evTag->SetTriggerMask(esd->GetTriggerMask());
 768     evTag->SetTriggerCluster(esd->GetTriggerCluster());
 769
 770     evTag->SetZDCNeutron1Energy(esd->GetZDCN1Energy());
 771     evTag->SetZDCProton1Energy(esd->GetZDCP1Energy());
 772     evTag->SetZDCEMEnergy(esd->GetZDCEMEnergy(0),esd->GetZDCEMEnergy(1));
 773     evTag->SetZDCNeutron1Energy(esd->GetZDCN2Energy());
 774     evTag->SetZDCProton1Energy(esd->GetZDCP2Energy());
 775     evTag->SetNumOfParticipants(esd->GetZDCParticipants());
 776
 777
 778     evTag->SetNumOfTracks(esd->GetNumberOfTracks());
 779     evTag->SetNumOfPosTracks(nPos);
 780     evTag->SetNumOfNegTracks(nNeg);
 781     evTag->SetNumOfNeutrTracks(nNeutr);
 782
 783     evTag->SetNumOfV0s(esd->GetNumberOfV0s());
 784     evTag->SetNumOfCascades(esd->GetNumberOfCascades());
 785     evTag->SetNumOfKinks(esd->GetNumberOfKinks());
 786     evTag->SetNumOfPMDTracks(esd->GetNumberOfPmdTracks());
 787
 788     evTag->SetNumOfProtons(nProtons);
 789     evTag->SetNumOfKaons(nKaons);
 790     evTag->SetNumOfPions(nPions);
 791     evTag->SetNumOfMuons(nMuons);
 792     evTag->SetNumOfFWMuons(nFWMuons);
 793     evTag->SetNumOfElectrons(nElectrons);
 794     evTag->SetNumOfPhotons(nGamas);
 795     evTag->SetNumOfPi0s(nPi0s);
 796     evTag->SetNumOfNeutrons(nNeutrons);
 797     evTag->SetNumOfKaon0s(nK0s);
 798
 799     evTag->SetNumOfChargedAbove1GeV(nCh1GeV);
 800     evTag->SetNumOfChargedAbove3GeV(nCh3GeV);
 801     evTag->SetNumOfChargedAbove10GeV(nCh10GeV);
 802     evTag->SetNumOfMuonsAbove1GeV(nMu1GeV);
 803     evTag->SetNumOfMuonsAbove3GeV(nMu3GeV);
 804     evTag->SetNumOfMuonsAbove10GeV(nMu10GeV);
 805     evTag->SetNumOfElectronsAbove1GeV(nEl1GeV);
 806     evTag->SetNumOfElectronsAbove3GeV(nEl3GeV);
 807     evTag->SetNumOfElectronsAbove10GeV(nEl10GeV);
 808
 809     evTag->SetNumOfPHOSClusters(esd->GetNumberOfPHOSClusters());
 810     evTag->SetNumOfEMCALClusters(esd->GetNumberOfEMCALClusters());
 811
 812     evTag->SetTotalMomentum(totalP);
 813     evTag->SetMeanPt(meanPt);
 814     evTag->SetMaxPt(maxPt);
 815
 816     tag->SetRunId(iInitRunNumber);
 817     if(fIsSim) tag->SetDataType(0);
 818     else tag->SetDataType(1);
 819     tag->AddEventTag(*evTag);
 820   }//event loop
 821   lastEvent = iNumberOfEvents;
 822
 823   gSystem->GetMemInfo(meminfo);
 824   AliInfo(Form("After the event and track loop - Memory used: %d MB - Increase: %d MB",meminfo->fMemUsed,meminfo->fMemUsed - tempmem));
 825   tempmem = meminfo->fMemUsed;
 826   t->Delete("");
 827
 828   gSystem->GetMemInfo(meminfo);
 829   AliInfo(Form("After the t->Delete - Memory used: %d MB - Increase: %d MB",meminfo->fMemUsed,meminfo->fMemUsed - tempmem));
 830   tempmem = meminfo->fMemUsed;
 831
 832   TString localFileName = "Run"; localFileName += tag->GetRunId();
 833   localFileName += ".Event"; localFileName += firstEvent; localFileName += "_"; localFileName += lastEvent; localFileName += "."; localFileName += Counter;
 834   localFileName += ".ESD.tag.root";
 835
 836   TString fileName;
 837
 838   if(fStorage == 0) {
 839     fileName = localFileName.Data();
 840     AliInfo(Form("Writing tags to local file: %s",fileName.Data()));
 841   }
 842   else if(fStorage == 1) {
 843     TString alienLocation = "/alien";
 844     alienLocation += gGrid->Pwd();
 845     alienLocation += fgridpath.Data();
 846     alienLocation += "/";
 847     alienLocation +=  localFileName;
 848     alienLocation += "?se=";
 849     alienLocation += fSE.Data();
 850     fileName = alienLocation.Data();
 851     AliInfo(Form("Writing tags to grid file: %s",fileName.Data()));
 852   }
 853
 854   TFile* ftag = TFile::Open(fileName, "recreate");
 855   ftag->cd();
 856   ttag.Fill();
 857   tag->Clear();
 858   ttag.Write();
 859   ftag->Close();
 860
 861   gSystem->GetMemInfo(meminfo);
 862   AliInfo(Form("After the file closing - Memory used: %d MB - Increase: %d MB",meminfo->fMemUsed,meminfo->fMemUsed - tempmem));
 863   tempmem = meminfo->fMemUsed;
 864
 865   delete ftag;
 866   delete esd;
 867
 868   delete tag;
 869   gSystem->GetMemInfo(meminfo);
 870   AliInfo(Form("After the delete objects - Memory used: %d MB - Increase: %d MB",meminfo->fMemUsed,meminfo->fMemUsed - tempmem));
 871 }
 872
 873 //_____________________________________________________________________________
 874 void AliESDTagCreator::CreateTag(TFile* file, const char *filepath, Int_t Counter) {
 875   //private method that creates tag files
 876   /////////////
 877   //muon code//
 878   ////////////
 879   Double_t fMUONMASS = 0.105658369;
 880   //Variables
 881   Double_t fX,fY,fZ ;
 882   Double_t fThetaX, fThetaY, fPyz, fChisquare;
 883   Double_t fPxRec, fPyRec, fPzRec, fEnergy;
 884   Int_t fCharge;
 885   TLorentzVector fEPvector;
 886
 887   Float_t fZVertexCut = 10.0;
 888   Float_t fRhoVertexCut = 2.0;
 889
 890   Float_t fLowPtCut = 1.0;
 891   Float_t fHighPtCut = 3.0;
 892   Float_t fVeryHighPtCut = 10.0;
 893   ////////////
 894
 895   Double_t partFrac[5] = {0.01, 0.01, 0.85, 0.10, 0.05};
 896
 897   // Creates the tags for all the events in a given ESD file
 898   Bool_t fIsSim = kTRUE;
 899   Int_t ntrack;
 900   Int_t nProtons, nKaons, nPions, nMuons, nElectrons, nFWMuons;
 901   Int_t nPos, nNeg, nNeutr;
 902   Int_t nK0s, nNeutrons, nPi0s, nGamas;
 903   Int_t nCh1GeV, nCh3GeV, nCh10GeV;
 904   Int_t nMu1GeV, nMu3GeV, nMu10GeV;
 905   Int_t nEl1GeV, nEl3GeV, nEl10GeV;
 906   Float_t maxPt = .0, meanPt = .0, totalP = .0;
 907   Int_t fVertexflag;
 908   Int_t iRunNumber = 0;
 909   TString fVertexName;
 910
 911   AliRunTag *tag = new AliRunTag();
 912   AliEventTag *evTag = new AliEventTag();
 913   TTree ttag("T","A Tree with event tags");
 914   TBranch * btag = ttag.Branch("AliTAG", &tag);
 915   btag->SetCompressionLevel(9);
 916
 917   AliInfo(Form("Creating the ESD tags......."));
 918
 919   Int_t firstEvent = 0,lastEvent = 0;
 920
 921   TTree *t = (TTree*) file->Get("esdTree");
 922   AliESDEvent *esd = new AliESDEvent();
 923   esd->ReadFromTree(t);
 924
 925   t->GetEntry(0);
 926   Int_t iInitRunNumber = esd->GetRunNumber();
 927
 928   Int_t iNumberOfEvents = (Int_t)t->GetEntries();
 929   for (Int_t iEventNumber = 0; iEventNumber < iNumberOfEvents; iEventNumber++) {
 930     ntrack = 0;
 931     nPos = 0;
 932     nNeg = 0;
 933     nNeutr =0;
 934     nK0s = 0;
 935     nNeutrons = 0;
 936     nPi0s = 0;
 937     nGamas = 0;
 938     nProtons = 0;
 939     nKaons = 0;
 940     nPions = 0;
 941     nMuons = 0;
 942     nFWMuons = 0;
 943     nElectrons = 0;
 944     nCh1GeV = 0;
 945     nCh3GeV = 0;
 946     nCh10GeV = 0;
 947     nMu1GeV = 0;
 948     nMu3GeV = 0;
 949     nMu10GeV = 0;
 950     nEl1GeV = 0;
 951     nEl3GeV = 0;
 952     nEl10GeV = 0;
 953     maxPt = .0;
 954     meanPt = .0;
 955     totalP = .0;
 956     fVertexflag = 1;
 957
 958     t->GetEntry(iEventNumber);
 959     iRunNumber = esd->GetRunNumber();
 960     if(iRunNumber != iInitRunNumber) AliFatal("Inconsistency of run numbers in the AliESD!!!");
 961     const AliESDVertex * vertexIn = esd->GetVertex();
 962     fVertexName = vertexIn->GetName();
 963     if(fVertexName == "default") fVertexflag = 0;
 964
 965     for (Int_t iTrackNumber = 0; iTrackNumber < esd->GetNumberOfTracks(); iTrackNumber++) {
 966       AliESDtrack * esdTrack = esd->GetTrack(iTrackNumber);
 967       if(esdTrack->GetLabel() != 0) fIsSim = kTRUE;
 968       else if(esdTrack->GetLabel() == 0) fIsSim = kFALSE;
 969       UInt_t status = esdTrack->GetStatus();
 970
 971       //select only tracks with ITS refit
 972       if ((status&AliESDtrack::kITSrefit)==0) continue;
 973       //select only tracks with TPC refit
 974       if ((status&AliESDtrack::kTPCrefit)==0) continue;
 975
 976       //select only tracks with the "combined PID"
 977       if ((status&AliESDtrack::kESDpid)==0) continue;
 978       Double_t p[3];
 979       esdTrack->GetPxPyPz(p);
 980       Double_t pt2 = p[0]*p[0]+p[1]*p[1];
 981       Double_t momentum = TMath::Sqrt(pt2+p[2]*p[2]);
 982       Double_t fPt = TMath::Sqrt(pt2);
 983       totalP += momentum;
 984       meanPt += fPt;
 985       if(fPt > maxPt) maxPt = fPt;
 986
 987       if(esdTrack->GetSign() > 0) {
 988         nPos++;
 989         if(fPt > fLowPtCut) nCh1GeV++;
 990         if(fPt > fHighPtCut) nCh3GeV++;
 991         if(fPt > fVeryHighPtCut) nCh10GeV++;
 992       }
 993       if(esdTrack->GetSign() < 0) {
 994         nNeg++;
 995         if(fPt > fLowPtCut) nCh1GeV++;
 996         if(fPt > fHighPtCut) nCh3GeV++;
 997         if(fPt > fVeryHighPtCut) nCh10GeV++;
 998       }
 999       if(esdTrack->GetSign() == 0) nNeutr++;
1000
1001       //PID
1002       Double_t prob[5];
1003       esdTrack->GetESDpid(prob);
1004
1005       Double_t rcc = 0.0;
1006       for(Int_t i = 0; i < AliPID::kSPECIES; i++) rcc += prob[i]*partFrac[i];
1007       if(rcc == 0.0) continue;
1008       //Bayes' formula
1009       Double_t w[5];
1010       for(Int_t i = 0; i < AliPID::kSPECIES; i++) w[i] = prob[i]*partFrac[i]/rcc;
1011
1012       //protons
1013       if ((w[4]>w[3])&&(w[4]>w[2])&&(w[4]>w[1])&&(w[4]>w[0])) nProtons++;
1014       //kaons
1015       if ((w[3]>w[4])&&(w[3]>w[2])&&(w[3]>w[1])&&(w[3]>w[0])) nKaons++;
1016       //pions
1017       if ((w[2]>w[4])&&(w[2]>w[3])&&(w[2]>w[1])&&(w[2]>w[0])) nPions++;
1018       //electrons
1019       if ((w[0]>w[4])&&(w[0]>w[3])&&(w[0]>w[2])&&(w[0]>w[1])) {
1020         nElectrons++;
1021         if(fPt > fLowPtCut) nEl1GeV++;
1022         if(fPt > fHighPtCut) nEl3GeV++;
1023         if(fPt > fVeryHighPtCut) nEl10GeV++;
1024       }
1025       ntrack++;
1026     }//esd track loop
1027
1028     /////////////
1029     //muon code//
1030     ////////////
1031     Int_t nMuonTracks = esd->GetNumberOfMuonTracks();
1032     // loop over all reconstructed tracks (also first track of combination)
1033     for (Int_t iTrack = 0; iTrack <  nMuonTracks;  iTrack++) {
1034       AliESDMuonTrack* muonTrack = esd->GetMuonTrack(iTrack);
1035       if (muonTrack == 0x0) continue;
1036
1037       // Coordinates at vertex
1038       fZ = muonTrack->GetZ();
1039       fY = muonTrack->GetBendingCoor();
1040       fX = muonTrack->GetNonBendingCoor();
1041
1042       fThetaX = muonTrack->GetThetaX();
1043       fThetaY = muonTrack->GetThetaY();
1044
1045       fPyz = 1./TMath::Abs(muonTrack->GetInverseBendingMomentum());
1046       fPzRec = - fPyz / TMath::Sqrt(1.0 + TMath::Tan(fThetaY)*TMath::Tan(fThetaY));
1047       fPxRec = fPzRec * TMath::Tan(fThetaX);
1048       fPyRec = fPzRec * TMath::Tan(fThetaY);
1049       fCharge = Int_t(TMath::Sign(1.,muonTrack->GetInverseBendingMomentum()));
1050
1051       //ChiSquare of the track if needed
1052       fChisquare = muonTrack->GetChi2()/(2.0 * muonTrack->GetNHit() - 5);
1053       fEnergy = TMath::Sqrt(fMUONMASS * fMUONMASS + fPxRec * fPxRec + fPyRec * fPyRec + fPzRec * fPzRec);
1054       fEPvector.SetPxPyPzE(fPxRec, fPyRec, fPzRec, fEnergy);
1055
1056       // total number of muons inside a vertex cut
1057       if((TMath::Abs(fZ)<fZVertexCut) && (TMath::Sqrt(fY*fY+fX*fX)<fRhoVertexCut)) {
1058         nMuons++;
1059         nFWMuons++;
1060         if(fEPvector.Pt() > fLowPtCut) {
1061           nMu1GeV++;
1062           if(fEPvector.Pt() > fHighPtCut) {
1063             nMu3GeV++;
1064             if (fEPvector.Pt() > fVeryHighPtCut) {
1065               nMu10GeV++;
1066             }
1067           }
1068         }
1069       }
1070     }//muon track loop
1071
1072     // Fill the event tags
1073     if(ntrack != 0) meanPt = meanPt/ntrack;
1074
1075     evTag->SetEventId(iEventNumber+1);
1076     evTag->SetPath(filepath);
1077
1078     evTag->SetVertexX(vertexIn->GetXv());
1079     evTag->SetVertexY(vertexIn->GetYv());
1080     evTag->SetVertexZ(vertexIn->GetZv());
1081     evTag->SetVertexZError(vertexIn->GetZRes());
1082     evTag->SetVertexFlag(fVertexflag);
1083
1084     evTag->SetT0VertexZ(esd->GetT0zVertex());
1085
1086     evTag->SetTriggerMask(esd->GetTriggerMask());
1087     evTag->SetTriggerCluster(esd->GetTriggerCluster());
1088
1089     evTag->SetZDCNeutron1Energy(esd->GetZDCN1Energy());
1090     evTag->SetZDCProton1Energy(esd->GetZDCP1Energy());
1091     evTag->SetZDCEMEnergy(esd->GetZDCEMEnergy(0),esd->GetZDCEMEnergy(1));
1092     evTag->SetZDCNeutron1Energy(esd->GetZDCN2Energy());
1093     evTag->SetZDCProton1Energy(esd->GetZDCP2Energy());
1094     evTag->SetNumOfParticipants(esd->GetZDCParticipants());
1095
1096
1097     evTag->SetNumOfTracks(esd->GetNumberOfTracks());
1098     evTag->SetNumOfPosTracks(nPos);
1099     evTag->SetNumOfNegTracks(nNeg);
1100     evTag->SetNumOfNeutrTracks(nNeutr);
1101
1102     evTag->SetNumOfV0s(esd->GetNumberOfV0s());
1103     evTag->SetNumOfCascades(esd->GetNumberOfCascades());
1104     evTag->SetNumOfKinks(esd->GetNumberOfKinks());
1105     evTag->SetNumOfPMDTracks(esd->GetNumberOfPmdTracks());
1106
1107     evTag->SetNumOfProtons(nProtons);
1108     evTag->SetNumOfKaons(nKaons);
1109     evTag->SetNumOfPions(nPions);
1110     evTag->SetNumOfMuons(nMuons);
1111     evTag->SetNumOfFWMuons(nFWMuons);
1112     evTag->SetNumOfElectrons(nElectrons);
1113     evTag->SetNumOfPhotons(nGamas);
1114     evTag->SetNumOfPi0s(nPi0s);
1115     evTag->SetNumOfNeutrons(nNeutrons);
1116     evTag->SetNumOfKaon0s(nK0s);
1117
1118     evTag->SetNumOfChargedAbove1GeV(nCh1GeV);
1119     evTag->SetNumOfChargedAbove3GeV(nCh3GeV);
1120     evTag->SetNumOfChargedAbove10GeV(nCh10GeV);
1121     evTag->SetNumOfMuonsAbove1GeV(nMu1GeV);
1122     evTag->SetNumOfMuonsAbove3GeV(nMu3GeV);
1123     evTag->SetNumOfMuonsAbove10GeV(nMu10GeV);
1124     evTag->SetNumOfElectronsAbove1GeV(nEl1GeV);
1125     evTag->SetNumOfElectronsAbove3GeV(nEl3GeV);
1126     evTag->SetNumOfElectronsAbove10GeV(nEl10GeV);
1127
1128     evTag->SetNumOfPHOSClusters(esd->GetNumberOfPHOSClusters());
1129     evTag->SetNumOfEMCALClusters(esd->GetNumberOfEMCALClusters());
1130
1131     evTag->SetTotalMomentum(totalP);
1132     evTag->SetMeanPt(meanPt);
1133     evTag->SetMaxPt(maxPt);
1134
1135     tag->SetRunId(iInitRunNumber);
1136     if(fIsSim) tag->SetDataType(0);
1137     else tag->SetDataType(1);
1138     tag->AddEventTag(*evTag);
1139   }//event loop
1140   lastEvent = iNumberOfEvents;
1141
1142   t->Delete("");
1143
1144   TString localFileName = "Run"; localFileName += tag->GetRunId();
1145   localFileName += ".Event"; localFileName += firstEvent; localFileName += "_"; localFileName += lastEvent; localFileName += "."; localFileName += Counter;
1146   localFileName += ".ESD.tag.root";
1147
1148   TString fileName;
1149
1150   if(fStorage == 0) {
1151     fileName = localFileName.Data();
1152     AliInfo(Form("Writing tags to local file: %s",fileName.Data()));
1153   }
1154   else if(fStorage == 1) {
1155     TString alienLocation = "/alien";
1156     alienLocation += gGrid->Pwd();
1157     alienLocation += fgridpath.Data();
1158     alienLocation += "/";
1159     alienLocation +=  localFileName;
1160     alienLocation += "?se=";
1161     alienLocation += fSE.Data();
1162     fileName = alienLocation.Data();
1163     AliInfo(Form("Writing tags to grid file: %s",fileName.Data()));
1164   }
1165
1166   TFile* ftag = TFile::Open(fileName, "recreate");
1167   ftag->cd();
1168   ttag.Fill();
1169   tag->Clear();
1170   ttag.Write();
1171   ftag->Close();
1172
1173   delete ftag;
1174   delete esd;
1175
1176   delete tag;
1177 }
1178
1179 //_____________________________________________________________________________
1180 void AliESDTagCreator::CreateESDTags(Int_t fFirstEvent, Int_t fLastEvent, AliGRPObject *grpData, ULong_t * qa, Bool_t * es, Int_t qalength, Int_t eslength) {
1181   //GRP
1182   Float_t lhcLuminosity = 0.0;
1183   TString lhcState = "test";
1184   //UInt_t detectorMask = 0;
1185   Int_t detectorMask = 0;
1186
1187   detectorMask = grpData->GetDetectorMask();
1188   time_t startTime = grpData->GetTimeStart();
1189   TTimeStamp *t1 = new TTimeStamp(startTime);
1190   time_t endTime = grpData->GetTimeEnd();
1191   TTimeStamp *t2 = new TTimeStamp(endTime);
1192   const char* beamtype = grpData->GetBeamType();
1193   Float_t beamenergy = grpData->GetBeamEnergy();
1194
1195
1196   /////////////
1197   //muon code//
1198   ////////////
1199   Double_t fMUONMASS = 0.105658369;
1200   //Variables
1201   Double_t fX,fY,fZ ;
1202   Double_t fThetaX, fThetaY, fPyz, fChisquare;
1203   Double_t fPxRec,fPyRec, fPzRec, fEnergy;
1204   Int_t fCharge;
1205   TLorentzVector fEPvector;
1206
1207   Float_t fZVertexCut = 10.0;
1208   Float_t fRhoVertexCut = 2.0;
1209
1210   Float_t fLowPtCut = 1.0;
1211   Float_t fHighPtCut = 3.0;
1212   Float_t fVeryHighPtCut = 10.0;
1213   ////////////
1214
1215   Double_t partFrac[5] = {0.01, 0.01, 0.85, 0.10, 0.05};
1216
1217   // Creates the tags for all the events in a given ESD file
1218   Int_t ntrack;
1219   Int_t nProtons, nKaons, nPions, nMuons, nElectrons, nFWMuons;
1220   Int_t nPos, nNeg, nNeutr;
1221   Int_t nK0s, nNeutrons, nPi0s, nGamas;
1222   Int_t nCh1GeV, nCh3GeV, nCh10GeV;
1223   Int_t nMu1GeV, nMu3GeV, nMu10GeV;
1224   Int_t nEl1GeV, nEl3GeV, nEl10GeV;
1225   Float_t maxPt = .0, meanPt = .0, totalP = .0;
1226   Int_t fVertexflag;
1227   Int_t iRunNumber = 0;
1228   TString fVertexName("default");
1229
1230   AliInfo(Form("Creating the ESD tags......."));
1231
1232   TFile *file = TFile::Open("AliESDs.root");
1233   if (!file || !file->IsOpen()) {
1234     AliError(Form("opening failed"));
1235     delete file;
1236     return ;
1237   }
1238   Int_t lastEvent = 0;
1239   TTree *b = (TTree*) file->Get("esdTree");
1240   AliESDEvent *esd = new AliESDEvent();
1241   esd->ReadFromTree(b);
1242
1243   b->GetEntry(fFirstEvent);
1244   Int_t iInitRunNumber = esd->GetRunNumber();
1245
1246   Int_t iNumberOfEvents = (Int_t)b->GetEntries();
1247   if(fLastEvent == -1) lastEvent = (Int_t)b->GetEntries();
1248   else lastEvent = fLastEvent;
1249
1250   char fileName[256];
1251   sprintf(fileName, "Run%d.Event%d_%d.ESD.tag.root",
1252           iInitRunNumber,fFirstEvent,lastEvent);
1253   AliInfo(Form("writing tags to file %s", fileName));
1254   AliDebug(1, Form("writing tags to file %s", fileName));
1255
1256   TFile* ftag = TFile::Open(fileName, "recreate");
1257
1258   AliRunTag *tag = new AliRunTag();
1259   AliEventTag *evTag = new AliEventTag();
1260   TTree ttag("T","A Tree with event tags");
1261   TBranch * btag = ttag.Branch("AliTAG", &tag);
1262   btag->SetCompressionLevel(9);
1263
1264   if(fLastEvent != -1) iNumberOfEvents = fLastEvent + 1;
1265   for (Int_t iEventNumber = fFirstEvent; iEventNumber < iNumberOfEvents; iEventNumber++) {
1266     ntrack = 0;
1267     nPos = 0;
1268     nNeg = 0;
1269     nNeutr =0;
1270     nK0s = 0;
1271     nNeutrons = 0;
1272     nPi0s = 0;
1273     nGamas = 0;
1274     nProtons = 0;
1275     nKaons = 0;
1276     nPions = 0;
1277     nMuons = 0;
1278     nFWMuons = 0;
1279     nElectrons = 0;
1280     nCh1GeV = 0;
1281     nCh3GeV = 0;
1282     nCh10GeV = 0;
1283     nMu1GeV = 0;
1284     nMu3GeV = 0;
1285     nMu10GeV = 0;
1286     nEl1GeV = 0;
1287     nEl3GeV = 0;
1288     nEl10GeV = 0;
1289     maxPt = .0;
1290     meanPt = .0;
1291     totalP = .0;
1292     fVertexflag = 0;
1293
1294     b->GetEntry(iEventNumber);
1295     iRunNumber = esd->GetRunNumber();
1296     if(iRunNumber != iInitRunNumber) AliFatal("Inconsistency of run numbers in the AliESD!!!");
1297     const AliESDVertex * vertexIn = esd->GetVertex();
1298     if (!vertexIn) AliError("ESD has not defined vertex.");
1299     if (vertexIn) fVertexName = vertexIn->GetName();
1300     if(fVertexName != "default") fVertexflag = 1;
1301     for (Int_t iTrackNumber = 0; iTrackNumber < esd->GetNumberOfTracks(); iTrackNumber++) {
1302       AliESDtrack * esdTrack = esd->GetTrack(iTrackNumber);
1303       UInt_t status = esdTrack->GetStatus();
1304
1305       //select only tracks with ITS refit
1306       if ((status&AliESDtrack::kITSrefit)==0) continue;
1307       //select only tracks with TPC refit
1308       if ((status&AliESDtrack::kTPCrefit)==0) continue;
1309
1310       //select only tracks with the "combined PID"
1311       if ((status&AliESDtrack::kESDpid)==0) continue;
1312       Double_t p[3];
1313       esdTrack->GetPxPyPz(p);
1314       Double_t pt2 = p[0]*p[0]+p[1]*p[1];
1315       Double_t momentum = TMath::Sqrt(pt2+p[2]*p[2]);
1316       Double_t fPt = TMath::Sqrt(pt2);
1317       totalP += momentum;
1318       meanPt += fPt;
1319       if(fPt > maxPt) maxPt = fPt;
1320
1321       if(esdTrack->GetSign() > 0) {
1322         nPos++;
1323         if(fPt > fLowPtCut) nCh1GeV++;
1324         if(fPt > fHighPtCut) nCh3GeV++;
1325         if(fPt > fVeryHighPtCut) nCh10GeV++;
1326       }
1327       if(esdTrack->GetSign() < 0) {
1328         nNeg++;
1329         if(fPt > fLowPtCut) nCh1GeV++;
1330         if(fPt > fHighPtCut) nCh3GeV++;
1331         if(fPt > fVeryHighPtCut) nCh10GeV++;
1332       }
1333       if(esdTrack->GetSign() == 0) nNeutr++;
1334
1335       //PID
1336       Double_t prob[5];
1337       esdTrack->GetESDpid(prob);
1338
1339       Double_t rcc = 0.0;
1340       for(Int_t i = 0; i < AliPID::kSPECIES; i++) rcc += prob[i]*partFrac[i];
1341       if(rcc == 0.0) continue;
1342       //Bayes' formula
1343       Double_t w[5];
1344       for(Int_t i = 0; i < AliPID::kSPECIES; i++) w[i] = prob[i]*partFrac[i]/rcc;
1345
1346       //protons
1347       if ((w[4]>w[3])&&(w[4]>w[2])&&(w[4]>w[1])&&(w[4]>w[0])) nProtons++;
1348       //kaons
1349       if ((w[3]>w[4])&&(w[3]>w[2])&&(w[3]>w[1])&&(w[3]>w[0])) nKaons++;
1350       //pions
1351       if ((w[2]>w[4])&&(w[2]>w[3])&&(w[2]>w[1])&&(w[2]>w[0])) nPions++;
1352       //electrons
1353       if ((w[0]>w[4])&&(w[0]>w[3])&&(w[0]>w[2])&&(w[0]>w[1])) {
1354         nElectrons++;
1355         if(fPt > fLowPtCut) nEl1GeV++;
1356         if(fPt > fHighPtCut) nEl3GeV++;
1357         if(fPt > fVeryHighPtCut) nEl10GeV++;
1358       }
1359       ntrack++;
1360     }//track loop
1361
1362     /////////////
1363     //muon code//
1364     ////////////
1365     Int_t nMuonTracks = esd->GetNumberOfMuonTracks();
1366     // loop over all reconstructed tracks (also first track of combination)
1367     for (Int_t iTrack = 0; iTrack <  nMuonTracks;  iTrack++) {
1368       AliESDMuonTrack* muonTrack = esd->GetMuonTrack(iTrack);
1369       if (muonTrack == 0x0) continue;
1370
1371       // Coordinates at vertex
1372       fZ = muonTrack->GetZ();
1373       fY = muonTrack->GetBendingCoor();
1374       fX = muonTrack->GetNonBendingCoor();
1375
1376       fThetaX = muonTrack->GetThetaX();
1377       fThetaY = muonTrack->GetThetaY();
1378
1379       fPyz = 1./TMath::Abs(muonTrack->GetInverseBendingMomentum());
1380       fPzRec = - fPyz / TMath::Sqrt(1.0 + TMath::Tan(fThetaY)*TMath::Tan(fThetaY));
1381       fPxRec = fPzRec * TMath::Tan(fThetaX);
1382       fPyRec = fPzRec * TMath::Tan(fThetaY);
1383       fCharge = Int_t(TMath::Sign(1.,muonTrack->GetInverseBendingMomentum()));
1384
1385       //ChiSquare of the track if needed
1386       fChisquare = muonTrack->GetChi2()/(2.0 * muonTrack->GetNHit() - 5);
1387       fEnergy = TMath::Sqrt(fMUONMASS * fMUONMASS + fPxRec * fPxRec + fPyRec * fPyRec + fPzRec * fPzRec);
1388       fEPvector.SetPxPyPzE(fPxRec, fPyRec, fPzRec, fEnergy);
1389
1390       // total number of muons inside a vertex cut
1391       if((TMath::Abs(fZ)<fZVertexCut) && (TMath::Sqrt(fY*fY+fX*fX)<fRhoVertexCut)) {
1392         nMuons++;
1393         nFWMuons++;
1394         if(fEPvector.Pt() > fLowPtCut) {
1395           nMu1GeV++;
1396           if(fEPvector.Pt() > fHighPtCut) {
1397             nMu3GeV++;
1398             if (fEPvector.Pt() > fVeryHighPtCut) {
1399               nMu10GeV++;
1400             }
1401           }
1402         }
1403       }
1404     }//muon track loop
1405
1406     // Fill the event tags
1407     if(ntrack != 0)
1408       meanPt = meanPt/ntrack;
1409
1410     evTag->SetEventId(iEventNumber+1);
1411     if (vertexIn) {
1412       evTag->SetVertexX(vertexIn->GetXv());
1413       evTag->SetVertexY(vertexIn->GetYv());
1414       evTag->SetVertexZ(vertexIn->GetZv());
1415       evTag->SetVertexZError(vertexIn->GetZRes());
1416     }
1417     evTag->SetVertexFlag(fVertexflag);
1418
1419     evTag->SetT0VertexZ(esd->GetT0zVertex());
1420
1421     evTag->SetTriggerMask(esd->GetTriggerMask());
1422     evTag->SetTriggerCluster(esd->GetTriggerCluster());
1423
1424     evTag->SetZDCNeutron1Energy(esd->GetZDCN1Energy());
1425     evTag->SetZDCProton1Energy(esd->GetZDCP1Energy());
1426     evTag->SetZDCNeutron2Energy(esd->GetZDCN2Energy());
1427     evTag->SetZDCProton2Energy(esd->GetZDCP2Energy());
1428     evTag->SetZDCEMEnergy(esd->GetZDCEMEnergy(0),esd->GetZDCEMEnergy(1));
1429     evTag->SetNumOfParticipants(esd->GetZDCParticipants());
1430
1431
1432     evTag->SetNumOfTracks(esd->GetNumberOfTracks());
1433     evTag->SetNumOfPosTracks(nPos);
1434     evTag->SetNumOfNegTracks(nNeg);
1435     evTag->SetNumOfNeutrTracks(nNeutr);
1436
1437     evTag->SetNumOfV0s(esd->GetNumberOfV0s());
1438     evTag->SetNumOfCascades(esd->GetNumberOfCascades());
1439     evTag->SetNumOfKinks(esd->GetNumberOfKinks());
1440     evTag->SetNumOfPMDTracks(esd->GetNumberOfPmdTracks());
1441
1442     evTag->SetNumOfProtons(nProtons);
1443     evTag->SetNumOfKaons(nKaons);
1444     evTag->SetNumOfPions(nPions);
1445     evTag->SetNumOfMuons(nMuons);
1446     evTag->SetNumOfFWMuons(nFWMuons);
1447     evTag->SetNumOfElectrons(nElectrons);
1448     evTag->SetNumOfPhotons(nGamas);
1449     evTag->SetNumOfPi0s(nPi0s);
1450     evTag->SetNumOfNeutrons(nNeutrons);
1451     evTag->SetNumOfKaon0s(nK0s);
1452
1453     evTag->SetNumOfChargedAbove1GeV(nCh1GeV);
1454     evTag->SetNumOfChargedAbove3GeV(nCh3GeV);
1455     evTag->SetNumOfChargedAbove10GeV(nCh10GeV);
1456     evTag->SetNumOfMuonsAbove1GeV(nMu1GeV);
1457     evTag->SetNumOfMuonsAbove3GeV(nMu3GeV);
1458     evTag->SetNumOfMuonsAbove10GeV(nMu10GeV);
1459     evTag->SetNumOfElectronsAbove1GeV(nEl1GeV);
1460     evTag->SetNumOfElectronsAbove3GeV(nEl3GeV);
1461     evTag->SetNumOfElectronsAbove10GeV(nEl10GeV);
1462
1463     evTag->SetNumOfPHOSClusters(esd->GetNumberOfPHOSClusters());
1464     evTag->SetNumOfEMCALClusters(esd->GetNumberOfEMCALClusters());
1465
1466     evTag->SetTotalMomentum(totalP);
1467     evTag->SetMeanPt(meanPt);
1468     evTag->SetMaxPt(maxPt);
1469
1470     tag->SetLHCTag(lhcLuminosity,lhcState);
1471     tag->SetDetectorTag(detectorMask);
1472
1473     tag->SetRunId(iInitRunNumber);
1474     tag->SetRunStartTime(t1->GetDate());
1475     tag->SetRunStopTime(t2->GetDate());
1476     tag->SetBeamEnergy(beamenergy);
1477     tag->SetBeamType(beamtype);
1478
1479     //QA setting
1480     tag->SetQA(qa, qalength) ;
1481     tag->SetEventSpecies(es, eslength) ;
1482
1483     tag->AddEventTag(*evTag);
1484   }
1485
1486   ftag->cd();
1487   ttag.Fill();
1488   tag->Clear();
1489   ttag.Write();
1490   ftag->Close();
1491   file->cd();
1492   delete tag;
1493   delete evTag;
1494 }
1495
1496 //_____________________________________________________________________________
1497 void AliESDTagCreator::SwitchOffBranches() const {
1498   //
1499   // Switch of branches on user request
1500   TObjArray * tokens = fBranches.Tokenize(" ");
1501   Int_t ntok = tokens->GetEntries();
1502   for (Int_t i = 0; i < ntok; i++)  {
1503     TString str = ((TObjString*) tokens->At(i))->GetString();
1504     fChain->SetBranchStatus(Form("%s%s%s","*", str.Data(), "*"), 0);
1505     AliInfo(Form("Branch %s switched off \n", str.Data()));
1506   }
1507 }