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