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