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