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