STEER/AliAODTagCreator.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 //           AliAODTagCreator 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 "AliPID.h"
  39 #include "AliAODEvent.h"
  40 #include "AliAODVertex.h"
  41 #include "AliLog.h"
  42
  43 #include "AliAODTagCreator.h"
  44
  45
  46 ClassImp(AliAODTagCreator)
  47
  48
  49 //______________________________________________________________________________
  50   AliAODTagCreator::AliAODTagCreator() :
  51     AliTagCreator() {
  52   //==============Default constructor for a AliAODTagCreator================
  53 }
  54
  55 //______________________________________________________________________________
  56 AliAODTagCreator::~AliAODTagCreator() {
  57 //================Default destructor for a AliAODTagCreator===================
  58 }
  59
  60 //______________________________________________________________________________
  61 Bool_t AliAODTagCreator::ReadGridCollection(TGridResult *fresult) {
  62   // Reads the entry of the TGridResult and creates the tags
  63   Int_t nEntries = fresult->GetEntries();
  64
  65   TString alienUrl;
  66   const char *guid;
  67   const char *md5;
  68   const char *turl;
  69   Long64_t size = -1;
  70
  71   Int_t counter = 0;
  72   for(Int_t i = 0; i < nEntries; i++) {
  73     alienUrl = fresult->GetKey(i,"turl");
  74     guid = fresult->GetKey(i,"guid");
  75     if(fresult->GetKey(i,"size")) size = atol (fresult->GetKey(i,"size"));
  76     md5 = fresult->GetKey(i,"md5");
  77     turl = fresult->GetKey(i,"turl");
  78     if(md5 && !strlen(guid)) md5 = 0;
  79     if(guid && !strlen(guid)) guid = 0;
  80
  81     TFile *f = TFile::Open(alienUrl,"READ");
  82     //CreateTag(f,guid,md5,turl,size,counter);
  83     f->Close();
  84     delete f;
  85     counter += 1;
  86   }//grid result loop
  87
  88   return kTRUE;
  89 }
  90
  91 //______________________________________________________________________________
  92 Bool_t AliAODTagCreator::ReadLocalCollection(const char *localpath) {
  93   // Checks the different subdirs of the given local path and in the
  94   // case where it finds an AliAODs.root file it creates the tags
  95
  96   void *dira =  gSystem->OpenDirectory(localpath);
  97   Char_t fPath[256];
  98   const char * dirname = 0x0;
  99   const char * filename = 0x0;
 100   const char * pattern = "AliAODs.root";
 101
 102   Int_t counter = 0;
 103   while((dirname = gSystem->GetDirEntry(dira))) {
 104     sprintf(fPath,"%s/%s",localpath,dirname);
 105     void *dirb =  gSystem->OpenDirectory(fPath);
 106     while((filename = gSystem->GetDirEntry(dirb))) {
 107       if(strstr(filename,pattern)) {
 108         TString fAODFileName;
 109         fAODFileName = fPath;
 110         fAODFileName += "/";
 111         fAODFileName += pattern;
 112         TFile *f = TFile::Open(fAODFileName,"READ");
 113         //CreateTag(f,fAODFileName,counter);
 114         f->Close();
 115         delete f;
 116
 117         counter += 1;
 118       }//pattern check
 119     }//child directory's entry loop
 120   }//parent directory's entry loop
 121
 122   return kTRUE;
 123 }
 124
 125 //______________________________________________________________________________
 126 Bool_t AliAODTagCreator::ReadCAFCollection(const char *filename) {
 127   // Temporary solution for CAF: Takes as an input the ascii file that
 128   // lists the AODs stored in the SE of the CAF and creates the tags.
 129
 130   // Open the input stream
 131   ifstream in;
 132   in.open(filename);
 133
 134   Int_t counter = 0;
 135   TString esdfile;
 136   // Read the input list of files and add them to the chain
 137   while(in.good()) {
 138     in >> esdfile;
 139     if (!esdfile.Contains("root")) continue; // protection
 140     TFile *f = TFile::Open(esdfile,"READ");
 141     //CreateTag(f,esdfile,counter);
 142     f->Close();
 143     delete f;
 144
 145     counter += 1;
 146   }
 147
 148   return kTRUE;
 149 }
 150
 151 //__________________________________________________________________________
 152 void AliAODTagCreator::CreateAODTags(Int_t fFirstEvent, Int_t fLastEvent) {
 153   //creates tag files for AODs
 154
 155   Float_t fLowPtCut = 1.0;
 156   Float_t fHighPtCut = 3.0;
 157   Float_t fVeryHighPtCut = 10.0;
 158   ////////////
 159
 160   Double_t partFrac[10] = {0.01, 0.01, 0.85, 0.10, 0.05, 0., 0., 0., 0., 0.};
 161
 162   // Creates the tags for all the events in a given AOD file
 163   Int_t ntrack;
 164   Int_t nProtons, nKaons, nPions, nMuons, nElectrons;
 165   Int_t nPos, nNeg, nNeutr;
 166   Int_t nKinks, nV0s, nCascades;
 167   Int_t nK0s, nNeutrons, nPi0s, nGamas;
 168   Int_t nCh1GeV, nCh3GeV, nCh10GeV;
 169   Int_t nMu1GeV, nMu3GeV, nMu10GeV;
 170   Int_t nEl1GeV, nEl3GeV, nEl10GeV;
 171   Float_t maxPt = .0, meanPt = .0, totalP = .0;
 172
 173   AliRunTag *tag = new AliRunTag();
 174   TTree ttag("T","A Tree with event tags");
 175   TBranch * btag = ttag.Branch("AliTAG", &tag);
 176   btag->SetCompressionLevel(9);
 177
 178   //reading the esd tag file
 179   TChain *oldTagTree = new TChain("T");
 180   const char * tagPattern = "ESD.tag";
 181   // Open the working directory
 182   void * dirp = gSystem->OpenDirectory(gSystem->pwd());
 183   const char * name = 0x0;
 184   // Add all files matching *pattern* to the chain
 185   while((name = gSystem->GetDirEntry(dirp))) {
 186     if (strstr(name,tagPattern)) oldTagTree->Add(name);
 187   }//directory loop
 188   AliInfo(Form("Chained tag files: %d",oldTagTree->GetEntries()));
 189
 190   //reading the esd tag file
 191   AliRunTag *oldtag = new AliRunTag();
 192   TString tagFilename;
 193   oldTagTree->SetBranchAddress("AliTAG",&oldtag);
 194   oldTagTree->GetEntry(0);
 195   tag->CopyStandardContent(oldtag);
 196   const TClonesArray *evTagList = oldtag->GetEventTags();
 197
 198   AliInfo(Form("Creating the AOD tags......."));
 199
 200   TFile *file = TFile::Open("AliAOD.root");
 201   if (!file || !file->IsOpen()) {
 202     AliError(Form("opening failed"));
 203     delete file;
 204     return ;
 205   }
 206   TTree *aodTree = (TTree*)file->Get("aodTree");
 207   AliAODEvent *aod = new AliAODEvent();
 208   aod->ReadFromTree(aodTree);
 209
 210   Int_t lastEvent = 0;
 211   if(fLastEvent == -1) lastEvent = (Int_t)aodTree->GetEntries();
 212   else lastEvent = fLastEvent;
 213
 214   // loop over events
 215   Int_t nEvents = aodTree->GetEntries();
 216   for (Int_t iEventNumber = 0; iEventNumber < nEvents; iEventNumber++) {
 217     AliEventTag *evTag = (AliEventTag *)evTagList->At(iEventNumber);
 218     ntrack = 0;
 219     nPos = 0; nNeg = 0; nNeutr =0;
 220     nKinks = 0; nV0s = 0; nCascades = 0;
 221     nK0s = 0; nNeutrons = 0; nPi0s = 0; nGamas = 0;
 222     nProtons = 0;  nKaons = 0; nPions = 0; nMuons = 0; nElectrons = 0;
 223     nCh1GeV = 0; nCh3GeV = 0; nCh10GeV = 0;
 224     nMu1GeV = 0; nMu3GeV = 0; nMu10GeV = 0;
 225     nEl1GeV = 0; nEl3GeV = 0; nEl10GeV = 0;
 226     maxPt = .0; meanPt = .0; totalP = .0;
 227
 228     // read events
 229     aodTree->GetEvent(iEventNumber);
 230
 231     // set pointers
 232     aod->GetStdContent();
 233
 234     Int_t nTracks = aod->GetNTracks();
 235     // loop over vertices
 236     Int_t nVtxs = aod->GetNVertices();
 237     for (Int_t nVtx = 0; nVtx < nVtxs; nVtx++) {
 238       // print track info
 239       AliAODVertex *vertex = aod->GetVertex(nVtx);
 240       if(vertex->GetType() == 1) nKinks += 1;
 241       if(vertex->GetType() == 2) nV0s += 1;
 242       if(vertex->GetType() == 3) nCascades += 1;
 243     }
 244     for (Int_t nTr = 0; nTr < nTracks; nTr++) {
 245       AliAODTrack *track = aod->GetTrack(nTr);
 246
 247       Double_t fPt = track->Pt();
 248       if(fPt > maxPt) maxPt = fPt;
 249       if(track->Charge() > 0) {
 250         nPos++;
 251         if(fPt > fLowPtCut) nCh1GeV++;
 252         if(fPt > fHighPtCut) nCh3GeV++;
 253         if(fPt > fVeryHighPtCut) nCh10GeV++;
 254       }
 255       if(track->Charge() < 0) {
 256         nNeg++;
 257         if(fPt > fLowPtCut) nCh1GeV++;
 258         if(fPt > fHighPtCut) nCh3GeV++;
 259         if(fPt > fVeryHighPtCut) nCh10GeV++;
 260       }
 261       if(track->Charge() == 0) nNeutr++;
 262
 263       //PID
 264       const Double32_t *prob = track->PID();
 265       Double_t rcc = 0.0;
 266       for(Int_t i = 0; i < AliPID::kSPECIES; i++) rcc += prob[i]*partFrac[i];
 267       if(rcc == 0.0) continue;
 268       //Bayes' formula
 269       Double_t w[10];
 270       for(Int_t i = 0; i < AliPID::kSPECIES; i++) w[i] = prob[i]*partFrac[i]/rcc;
 271
 272       //protons
 273       if ((w[4]>w[3])&&(w[4]>w[2])&&(w[4]>w[1])&&(w[4]>w[0])) nProtons++;
 274       //kaons
 275       if ((w[3]>w[4])&&(w[3]>w[2])&&(w[3]>w[1])&&(w[3]>w[0])) nKaons++;
 276       //pions
 277       if ((w[2]>w[4])&&(w[2]>w[3])&&(w[2]>w[1])&&(w[2]>w[0])) nPions++;
 278       //muons
 279       if ((w[1]>w[4])&&(w[1]>w[3])&&(w[1]>w[2])&&(w[1]>w[0])) {
 280         nMuons++;
 281         if(fPt > fLowPtCut) nMu1GeV++;
 282         if(fPt > fHighPtCut) nMu3GeV++;
 283         if(fPt > fVeryHighPtCut) nMu10GeV++;
 284         }
 285       //electrons
 286       if ((w[0]>w[4])&&(w[0]>w[3])&&(w[0]>w[2])&&(w[0]>w[1])) {
 287         nElectrons++;
 288         if(fPt > fLowPtCut) nEl1GeV++;
 289         if(fPt > fHighPtCut) nEl3GeV++;
 290         if(fPt > fVeryHighPtCut) nEl10GeV++;
 291       }
 292
 293       totalP += track->P();
 294       meanPt += fPt;
 295       ntrack++;
 296     }//track loop
 297     // Fill the event tags
 298     if(ntrack != 0)
 299       meanPt = meanPt/ntrack;
 300
 301     evTag->SetEventId(iEventNumber+1);
 302
 303     evTag->SetNumOfTracks(nTracks);
 304     evTag->SetNumOfPosTracks(nPos);
 305     evTag->SetNumOfNegTracks(nNeg);
 306     evTag->SetNumOfNeutrTracks(nNeutr);
 307
 308     evTag->SetNumOfV0s(nV0s);
 309     evTag->SetNumOfCascades(nCascades);
 310     evTag->SetNumOfKinks(nKinks);
 311
 312     evTag->SetNumOfProtons(nProtons);
 313     evTag->SetNumOfKaons(nKaons);
 314     evTag->SetNumOfPions(nPions);
 315     evTag->SetNumOfMuons(nMuons);
 316     evTag->SetNumOfElectrons(nElectrons);
 317     evTag->SetNumOfPhotons(nGamas);
 318     evTag->SetNumOfPi0s(nPi0s);
 319     evTag->SetNumOfNeutrons(nNeutrons);
 320     evTag->SetNumOfKaon0s(nK0s);
 321
 322     evTag->SetNumOfChargedAbove1GeV(nCh1GeV);
 323     evTag->SetNumOfChargedAbove3GeV(nCh3GeV);
 324     evTag->SetNumOfChargedAbove10GeV(nCh10GeV);
 325     evTag->SetNumOfMuonsAbove1GeV(nMu1GeV);
 326     evTag->SetNumOfMuonsAbove3GeV(nMu3GeV);
 327     evTag->SetNumOfMuonsAbove10GeV(nMu10GeV);
 328     evTag->SetNumOfElectronsAbove1GeV(nEl1GeV);
 329     evTag->SetNumOfElectronsAbove3GeV(nEl3GeV);
 330     evTag->SetNumOfElectronsAbove10GeV(nEl10GeV);
 331
 332     evTag->SetTotalMomentum(totalP);
 333     evTag->SetMeanPt(meanPt);
 334     evTag->SetMaxPt(maxPt);
 335     tag->AddEventTag(*evTag);
 336   }//event loop
 337   if(fLastEvent == -1) lastEvent = (Int_t)aodTree->GetEntries();
 338   else lastEvent = fLastEvent;
 339
 340   ttag.Fill();
 341   tag->Clear();
 342
 343   char fileName[256];
 344   sprintf(fileName, "Run%d.Event%d_%d.AOD.tag.root",
 345           tag->GetRunId(),fFirstEvent,lastEvent );
 346   AliInfo(Form("writing tags to file %s", fileName));
 347   AliDebug(1, Form("writing tags to file %s", fileName));
 348
 349   TFile* ftag = TFile::Open(fileName, "recreate");
 350   ftag->cd();
 351   ttag.Write();
 352   ftag->Close();
 353   file->cd();
 354   file->Close();
 355 }
 356
 357 //_____________________________________________________________________________
 358 /*void AliAODTagCreator::CreateTag(TFile* file, const char *guid, const char *md5, const char *turl, Long64_t size, Int_t Counter) {
 359   //private method that creates tag files
 360   TString fguid = guid;
 361   TString fmd5 = md5;
 362   TString fturl = turl;
 363   /////////////
 364   //muon code//
 365   ////////////
 366   Double_t fMUONMASS = 0.105658369;
 367   //Variables
 368   Double_t fX,fY,fZ ;
 369   Double_t fThetaX, fThetaY, fPyz, fChisquare;
 370   Double_t fPxRec, fPyRec, fPzRec, fEnergy;
 371   Int_t fCharge;
 372   TLorentzVector fEPvector;
 373
 374   Float_t fZVertexCut = 10.0;
 375   Float_t fRhoVertexCut = 2.0;
 376
 377   Float_t fLowPtCut = 1.0;
 378   Float_t fHighPtCut = 3.0;
 379   Float_t fVeryHighPtCut = 10.0;
 380   ////////////
 381
 382   Double_t partFrac[5] = {0.01, 0.01, 0.85, 0.10, 0.05};
 383
 384   // Creates the tags for all the events in a given AOD file
 385   Bool_t fIsSim = kTRUE;
 386   Int_t ntrack;
 387   Int_t nProtons, nKaons, nPions, nMuons, nElectrons;
 388   Int_t nPos, nNeg, nNeutr;
 389   Int_t nK0s, nNeutrons, nPi0s, nGamas;
 390   Int_t nCh1GeV, nCh3GeV, nCh10GeV;
 391   Int_t nMu1GeV, nMu3GeV, nMu10GeV;
 392   Int_t nEl1GeV, nEl3GeV, nEl10GeV;
 393   Float_t maxPt = .0, meanPt = .0, totalP = .0;
 394   Int_t fVertexflag;
 395   Int_t iRunNumber = 0;
 396   TString fVertexName;
 397
 398   AliRunTag *tag = new AliRunTag();
 399   AliEventTag *evTag = new AliEventTag();
 400   TTree ttag("T","A Tree with event tags");
 401   TBranch * btag = ttag.Branch("AliTAG", &tag);
 402   btag->SetCompressionLevel(9);
 403
 404   AliInfo(Form("Creating the tags......."));
 405
 406   Int_t firstEvent = 0,lastEvent = 0;
 407   TTree *t = (TTree*) file->Get("esdTree");
 408   AliESDEvent *esd = new AliESDEvent();
 409   esd->ReadFromTree(t);
 410
 411   t->GetEntry(0);
 412   Int_t iInitRunNumber = esd->GetRunNumber();
 413
 414   Int_t iNumberOfEvents = (Int_t)t->GetEntries();
 415   for (Int_t iEventNumber = 0; iEventNumber < iNumberOfEvents; iEventNumber++) {
 416     ntrack = 0;
 417     nPos = 0;
 418     nNeg = 0;
 419     nNeutr =0;
 420     nK0s = 0;
 421     nNeutrons = 0;
 422     nPi0s = 0;
 423     nGamas = 0;
 424     nProtons = 0;
 425     nKaons = 0;
 426     nPions = 0;
 427     nMuons = 0;
 428     nElectrons = 0;
 429     nCh1GeV = 0;
 430     nCh3GeV = 0;
 431     nCh10GeV = 0;
 432     nMu1GeV = 0;
 433     nMu3GeV = 0;
 434     nMu10GeV = 0;
 435     nEl1GeV = 0;
 436     nEl3GeV = 0;
 437     nEl10GeV = 0;
 438     maxPt = .0;
 439     meanPt = .0;
 440     totalP = .0;
 441     fVertexflag = 1;
 442
 443     t->GetEntry(iEventNumber);
 444     iRunNumber = esd->GetRunNumber();
 445     if(iRunNumber != iInitRunNumber) AliFatal("Inconsistency of run numbers in the AliESD!!!");
 446     const AliESDVertex * vertexIn = esd->GetVertex();
 447     fVertexName = vertexIn->GetName();
 448     if(fVertexName == "default") fVertexflag = 0;
 449
 450     for (Int_t iTrackNumber = 0; iTrackNumber < esd->GetNumberOfTracks(); iTrackNumber++) {
 451       AliESDtrack * esdTrack = esd->GetTrack(iTrackNumber);
 452       if(esdTrack->GetLabel() != 0) fIsSim = kTRUE;
 453       else if(esdTrack->GetLabel() == 0) fIsSim = kFALSE;
 454       UInt_t status = esdTrack->GetStatus();
 455
 456       //select only tracks with ITS refit
 457       if ((status&AliESDtrack::kITSrefit)==0) continue;
 458       //select only tracks with TPC refit
 459       if ((status&AliESDtrack::kTPCrefit)==0) continue;
 460
 461       //select only tracks with the "combined PID"
 462       if ((status&AliESDtrack::kESDpid)==0) continue;
 463       Double_t p[3];
 464       esdTrack->GetPxPyPz(p);
 465       Double_t momentum = sqrt(pow(p[0],2) + pow(p[1],2) + pow(p[2],2));
 466       Double_t fPt = sqrt(pow(p[0],2) + pow(p[1],2));
 467       totalP += momentum;
 468       meanPt += fPt;
 469       if(fPt > maxPt) maxPt = fPt;
 470
 471       if(esdTrack->GetSign() > 0) {
 472         nPos++;
 473         if(fPt > fLowPtCut) nCh1GeV++;
 474         if(fPt > fHighPtCut) nCh3GeV++;
 475         if(fPt > fVeryHighPtCut) nCh10GeV++;
 476       }
 477       if(esdTrack->GetSign() < 0) {
 478         nNeg++;
 479         if(fPt > fLowPtCut) nCh1GeV++;
 480         if(fPt > fHighPtCut) nCh3GeV++;
 481         if(fPt > fVeryHighPtCut) nCh10GeV++;
 482       }
 483       if(esdTrack->GetSign() == 0) nNeutr++;
 484
 485       //PID
 486       Double_t prob[5];
 487       esdTrack->GetESDpid(prob);
 488
 489       Double_t rcc = 0.0;
 490       for(Int_t i = 0; i < AliPID::kSPECIES; i++) rcc += prob[i]*partFrac[i];
 491       if(rcc == 0.0) continue;
 492       //Bayes' formula
 493       Double_t w[5];
 494       for(Int_t i = 0; i < AliPID::kSPECIES; i++) w[i] = prob[i]*partFrac[i]/rcc;
 495
 496       //protons
 497       if ((w[4]>w[3])&&(w[4]>w[2])&&(w[4]>w[1])&&(w[4]>w[0])) nProtons++;
 498       //kaons
 499       if ((w[3]>w[4])&&(w[3]>w[2])&&(w[3]>w[1])&&(w[3]>w[0])) nKaons++;
 500       //pions
 501       if ((w[2]>w[4])&&(w[2]>w[3])&&(w[2]>w[1])&&(w[2]>w[0])) nPions++;
 502       //electrons
 503       if ((w[0]>w[4])&&(w[0]>w[3])&&(w[0]>w[2])&&(w[0]>w[1])) {
 504         nElectrons++;
 505         if(fPt > fLowPtCut) nEl1GeV++;
 506         if(fPt > fHighPtCut) nEl3GeV++;
 507         if(fPt > fVeryHighPtCut) nEl10GeV++;
 508       }
 509       ntrack++;
 510     }//esd track loop
 511
 512     /////////////
 513     //muon code//
 514     ////////////
 515     Int_t nMuonTracks = esd->GetNumberOfMuonTracks();
 516     // loop over all reconstructed tracks (also first track of combination)
 517     for (Int_t iTrack = 0; iTrack <  nMuonTracks;  iTrack++) {
 518       AliESDMuonTrack* muonTrack = esd->GetMuonTrack(iTrack);
 519       if (muonTrack == 0x0) continue;
 520
 521       // Coordinates at vertex
 522       fZ = muonTrack->GetZ();
 523       fY = muonTrack->GetBendingCoor();
 524       fX = muonTrack->GetNonBendingCoor();
 525
 526       fThetaX = muonTrack->GetThetaX();
 527       fThetaY = muonTrack->GetThetaY();
 528
 529       fPyz = 1./TMath::Abs(muonTrack->GetInverseBendingMomentum());
 530       fPzRec = - fPyz / TMath::Sqrt(1.0 + TMath::Tan(fThetaY)*TMath::Tan(fThetaY));
 531       fPxRec = fPzRec * TMath::Tan(fThetaX);
 532       fPyRec = fPzRec * TMath::Tan(fThetaY);
 533       fCharge = Int_t(TMath::Sign(1.,muonTrack->GetInverseBendingMomentum()));
 534
 535       //ChiSquare of the track if needed
 536       fChisquare = muonTrack->GetChi2()/(2.0 * muonTrack->GetNHit() - 5);
 537       fEnergy = TMath::Sqrt(fMUONMASS * fMUONMASS + fPxRec * fPxRec + fPyRec * fPyRec + fPzRec * fPzRec);
 538       fEPvector.SetPxPyPzE(fPxRec, fPyRec, fPzRec, fEnergy);
 539
 540       // total number of muons inside a vertex cut
 541       if((TMath::Abs(fZ)<fZVertexCut) && (TMath::Sqrt(fY*fY+fX*fX)<fRhoVertexCut)) {
 542         nMuons++;
 543         if(fEPvector.Pt() > fLowPtCut) {
 544           nMu1GeV++;
 545           if(fEPvector.Pt() > fHighPtCut) {
 546             nMu3GeV++;
 547             if (fEPvector.Pt() > fVeryHighPtCut) {
 548               nMu10GeV++;
 549             }
 550           }
 551         }
 552       }
 553     }//muon track loop
 554
 555     // Fill the event tags
 556     if(ntrack != 0) meanPt = meanPt/ntrack;
 557
 558     evTag->SetEventId(iEventNumber+1);
 559     evTag->SetGUID(fguid);
 560     evTag->SetMD5(fmd5);
 561     evTag->SetTURL(fturl);
 562     evTag->SetSize(size);
 563     evTag->SetVertexX(vertexIn->GetXv());
 564     evTag->SetVertexY(vertexIn->GetYv());
 565     evTag->SetVertexZ(vertexIn->GetZv());
 566     evTag->SetVertexZError(vertexIn->GetZRes());
 567     evTag->SetVertexFlag(fVertexflag);
 568
 569     evTag->SetT0VertexZ(esd->GetT0zVertex());
 570
 571     evTag->SetTriggerMask(esd->GetTriggerMask());
 572     evTag->SetTriggerCluster(esd->GetTriggerCluster());
 573
 574     evTag->SetZDCNeutron1Energy(esd->GetZDCN1Energy());
 575     evTag->SetZDCProton1Energy(esd->GetZDCP1Energy());
 576     evTag->SetZDCEMEnergy(esd->GetZDCEMEnergy());
 577     evTag->SetZDCNeutron1Energy(esd->GetZDCN2Energy());
 578     evTag->SetZDCProton1Energy(esd->GetZDCP2Energy());
 579     evTag->SetNumOfParticipants(esd->GetZDCParticipants());
 580
 581
 582     evTag->SetNumOfTracks(esd->GetNumberOfTracks());
 583     evTag->SetNumOfPosTracks(nPos);
 584     evTag->SetNumOfNegTracks(nNeg);
 585     evTag->SetNumOfNeutrTracks(nNeutr);
 586
 587     evTag->SetNumOfV0s(esd->GetNumberOfV0s());
 588     evTag->SetNumOfCascades(esd->GetNumberOfCascades());
 589     evTag->SetNumOfKinks(esd->GetNumberOfKinks());
 590     evTag->SetNumOfPMDTracks(esd->GetNumberOfPmdTracks());
 591
 592     evTag->SetNumOfProtons(nProtons);
 593     evTag->SetNumOfKaons(nKaons);
 594     evTag->SetNumOfPions(nPions);
 595     evTag->SetNumOfMuons(nMuons);
 596     evTag->SetNumOfElectrons(nElectrons);
 597     evTag->SetNumOfPhotons(nGamas);
 598     evTag->SetNumOfPi0s(nPi0s);
 599     evTag->SetNumOfNeutrons(nNeutrons);
 600     evTag->SetNumOfKaon0s(nK0s);
 601
 602     evTag->SetNumOfChargedAbove1GeV(nCh1GeV);
 603     evTag->SetNumOfChargedAbove3GeV(nCh3GeV);
 604     evTag->SetNumOfChargedAbove10GeV(nCh10GeV);
 605     evTag->SetNumOfMuonsAbove1GeV(nMu1GeV);
 606     evTag->SetNumOfMuonsAbove3GeV(nMu3GeV);
 607     evTag->SetNumOfMuonsAbove10GeV(nMu10GeV);
 608     evTag->SetNumOfElectronsAbove1GeV(nEl1GeV);
 609     evTag->SetNumOfElectronsAbove3GeV(nEl3GeV);
 610     evTag->SetNumOfElectronsAbove10GeV(nEl10GeV);
 611
 612     evTag->SetNumOfPHOSClusters(esd->GetNumberOfPHOSClusters());
 613     evTag->SetNumOfEMCALClusters(esd->GetNumberOfEMCALClusters());
 614
 615     evTag->SetTotalMomentum(totalP);
 616     evTag->SetMeanPt(meanPt);
 617     evTag->SetMaxPt(maxPt);
 618
 619     tag->SetRunId(iInitRunNumber);
 620     if(fIsSim) tag->SetDataType(0);
 621     else tag->SetDataType(1);
 622     tag->AddEventTag(*evTag);
 623   }//event loop
 624   lastEvent = iNumberOfEvents;
 625
 626   t->Delete("");
 627
 628   ttag.Fill();
 629   tag->Clear();
 630
 631   TString localFileName = "Run"; localFileName += tag->GetRunId();
 632   localFileName += ".Event"; localFileName += firstEvent; localFileName += "_"; localFileName += lastEvent; localFileName += "."; localFileName += Counter;
 633   localFileName += ".ESD.tag.root";
 634
 635   TString fileName;
 636
 637   if(fStorage == 0) {
 638     fileName = localFileName.Data();
 639     AliInfo(Form("Writing tags to local file: %s",fileName.Data()));
 640   }
 641   else if(fStorage == 1) {
 642     TString alienLocation = "/alien";
 643     alienLocation += gGrid->Pwd();
 644     alienLocation += fgridpath.Data();
 645     alienLocation += "/";
 646     alienLocation +=  localFileName;
 647     alienLocation += "?se=";
 648     alienLocation += fSE.Data();
 649     fileName = alienLocation.Data();
 650     AliInfo(Form("Writing tags to grid file: %s",fileName.Data()));
 651   }
 652
 653   TFile* ftag = TFile::Open(fileName, "recreate");
 654   ftag->cd();
 655   ttag.Write();
 656   ftag->Close();
 657
 658   delete ftag;
 659   delete esd;
 660
 661   delete tag;
 662   }*/
 663
 664 //_____________________________________________________________________________
 665 /*void AliESDTagCreator::CreateTag(TFile* file, const char *filepath, Int_t Counter) {
 666   //private method that creates tag files
 667   /////////////
 668   //muon code//
 669   ////////////
 670   Double_t fMUONMASS = 0.105658369;
 671   //Variables
 672   Double_t fX,fY,fZ ;
 673   Double_t fThetaX, fThetaY, fPyz, fChisquare;
 674   Double_t fPxRec, fPyRec, fPzRec, fEnergy;
 675   Int_t fCharge;
 676   TLorentzVector fEPvector;
 677
 678   Float_t fZVertexCut = 10.0;
 679   Float_t fRhoVertexCut = 2.0;
 680
 681   Float_t fLowPtCut = 1.0;
 682   Float_t fHighPtCut = 3.0;
 683   Float_t fVeryHighPtCut = 10.0;
 684   ////////////
 685
 686   Double_t partFrac[5] = {0.01, 0.01, 0.85, 0.10, 0.05};
 687
 688   // Creates the tags for all the events in a given ESD file
 689   Bool_t fIsSim = kTRUE;
 690   Int_t ntrack;
 691   Int_t nProtons, nKaons, nPions, nMuons, nElectrons;
 692   Int_t nPos, nNeg, nNeutr;
 693   Int_t nK0s, nNeutrons, nPi0s, nGamas;
 694   Int_t nCh1GeV, nCh3GeV, nCh10GeV;
 695   Int_t nMu1GeV, nMu3GeV, nMu10GeV;
 696   Int_t nEl1GeV, nEl3GeV, nEl10GeV;
 697   Float_t maxPt = .0, meanPt = .0, totalP = .0;
 698   Int_t fVertexflag;
 699   Int_t iRunNumber = 0;
 700   TString fVertexName;
 701
 702   AliRunTag *tag = new AliRunTag();
 703   AliEventTag *evTag = new AliEventTag();
 704   TTree ttag("T","A Tree with event tags");
 705   TBranch * btag = ttag.Branch("AliTAG", &tag);
 706   btag->SetCompressionLevel(9);
 707
 708   AliInfo(Form("Creating the tags......."));
 709
 710   Int_t firstEvent = 0,lastEvent = 0;
 711
 712   TTree *t = (TTree*) file->Get("esdTree");
 713   AliESDEvent *esd = new AliESDEvent();
 714   esd->ReadFromTree(t);
 715
 716   t->GetEntry(0);
 717   Int_t iInitRunNumber = esd->GetRunNumber();
 718
 719   Int_t iNumberOfEvents = (Int_t)t->GetEntries();
 720   for (Int_t iEventNumber = 0; iEventNumber < iNumberOfEvents; iEventNumber++) {
 721     ntrack = 0;
 722     nPos = 0;
 723     nNeg = 0;
 724     nNeutr =0;
 725     nK0s = 0;
 726     nNeutrons = 0;
 727     nPi0s = 0;
 728     nGamas = 0;
 729     nProtons = 0;
 730     nKaons = 0;
 731     nPions = 0;
 732     nMuons = 0;
 733     nElectrons = 0;
 734     nCh1GeV = 0;
 735     nCh3GeV = 0;
 736     nCh10GeV = 0;
 737     nMu1GeV = 0;
 738     nMu3GeV = 0;
 739     nMu10GeV = 0;
 740     nEl1GeV = 0;
 741     nEl3GeV = 0;
 742     nEl10GeV = 0;
 743     maxPt = .0;
 744     meanPt = .0;
 745     totalP = .0;
 746     fVertexflag = 1;
 747
 748     t->GetEntry(iEventNumber);
 749     iRunNumber = esd->GetRunNumber();
 750     if(iRunNumber != iInitRunNumber) AliFatal("Inconsistency of run numbers in the AliESD!!!");
 751     const AliESDVertex * vertexIn = esd->GetVertex();
 752     fVertexName = vertexIn->GetName();
 753     if(fVertexName == "default") fVertexflag = 0;
 754
 755     for (Int_t iTrackNumber = 0; iTrackNumber < esd->GetNumberOfTracks(); iTrackNumber++) {
 756       AliESDtrack * esdTrack = esd->GetTrack(iTrackNumber);
 757       if(esdTrack->GetLabel() != 0) fIsSim = kTRUE;
 758       else if(esdTrack->GetLabel() == 0) fIsSim = kFALSE;
 759       UInt_t status = esdTrack->GetStatus();
 760
 761       //select only tracks with ITS refit
 762       if ((status&AliESDtrack::kITSrefit)==0) continue;
 763       //select only tracks with TPC refit
 764       if ((status&AliESDtrack::kTPCrefit)==0) continue;
 765
 766       //select only tracks with the "combined PID"
 767       if ((status&AliESDtrack::kESDpid)==0) continue;
 768       Double_t p[3];
 769       esdTrack->GetPxPyPz(p);
 770       Double_t momentum = sqrt(pow(p[0],2) + pow(p[1],2) + pow(p[2],2));
 771       Double_t fPt = sqrt(pow(p[0],2) + pow(p[1],2));
 772       totalP += momentum;
 773       meanPt += fPt;
 774       if(fPt > maxPt) maxPt = fPt;
 775
 776       if(esdTrack->GetSign() > 0) {
 777         nPos++;
 778         if(fPt > fLowPtCut) nCh1GeV++;
 779         if(fPt > fHighPtCut) nCh3GeV++;
 780         if(fPt > fVeryHighPtCut) nCh10GeV++;
 781       }
 782       if(esdTrack->GetSign() < 0) {
 783         nNeg++;
 784         if(fPt > fLowPtCut) nCh1GeV++;
 785         if(fPt > fHighPtCut) nCh3GeV++;
 786         if(fPt > fVeryHighPtCut) nCh10GeV++;
 787       }
 788       if(esdTrack->GetSign() == 0) nNeutr++;
 789
 790       //PID
 791       Double_t prob[5];
 792       esdTrack->GetESDpid(prob);
 793
 794       Double_t rcc = 0.0;
 795       for(Int_t i = 0; i < AliPID::kSPECIES; i++) rcc += prob[i]*partFrac[i];
 796       if(rcc == 0.0) continue;
 797       //Bayes' formula
 798       Double_t w[5];
 799       for(Int_t i = 0; i < AliPID::kSPECIES; i++) w[i] = prob[i]*partFrac[i]/rcc;
 800
 801       //protons
 802       if ((w[4]>w[3])&&(w[4]>w[2])&&(w[4]>w[1])&&(w[4]>w[0])) nProtons++;
 803       //kaons
 804       if ((w[3]>w[4])&&(w[3]>w[2])&&(w[3]>w[1])&&(w[3]>w[0])) nKaons++;
 805       //pions
 806       if ((w[2]>w[4])&&(w[2]>w[3])&&(w[2]>w[1])&&(w[2]>w[0])) nPions++;
 807       //electrons
 808       if ((w[0]>w[4])&&(w[0]>w[3])&&(w[0]>w[2])&&(w[0]>w[1])) {
 809         nElectrons++;
 810         if(fPt > fLowPtCut) nEl1GeV++;
 811         if(fPt > fHighPtCut) nEl3GeV++;
 812         if(fPt > fVeryHighPtCut) nEl10GeV++;
 813       }
 814       ntrack++;
 815     }//esd track loop
 816
 817     /////////////
 818     //muon code//
 819     ////////////
 820     Int_t nMuonTracks = esd->GetNumberOfMuonTracks();
 821     // loop over all reconstructed tracks (also first track of combination)
 822     for (Int_t iTrack = 0; iTrack <  nMuonTracks;  iTrack++) {
 823       AliESDMuonTrack* muonTrack = esd->GetMuonTrack(iTrack);
 824       if (muonTrack == 0x0) continue;
 825
 826       // Coordinates at vertex
 827       fZ = muonTrack->GetZ();
 828       fY = muonTrack->GetBendingCoor();
 829       fX = muonTrack->GetNonBendingCoor();
 830
 831       fThetaX = muonTrack->GetThetaX();
 832       fThetaY = muonTrack->GetThetaY();
 833
 834       fPyz = 1./TMath::Abs(muonTrack->GetInverseBendingMomentum());
 835       fPzRec = - fPyz / TMath::Sqrt(1.0 + TMath::Tan(fThetaY)*TMath::Tan(fThetaY));
 836       fPxRec = fPzRec * TMath::Tan(fThetaX);
 837       fPyRec = fPzRec * TMath::Tan(fThetaY);
 838       fCharge = Int_t(TMath::Sign(1.,muonTrack->GetInverseBendingMomentum()));
 839
 840       //ChiSquare of the track if needed
 841       fChisquare = muonTrack->GetChi2()/(2.0 * muonTrack->GetNHit() - 5);
 842       fEnergy = TMath::Sqrt(fMUONMASS * fMUONMASS + fPxRec * fPxRec + fPyRec * fPyRec + fPzRec * fPzRec);
 843       fEPvector.SetPxPyPzE(fPxRec, fPyRec, fPzRec, fEnergy);
 844
 845       // total number of muons inside a vertex cut
 846       if((TMath::Abs(fZ)<fZVertexCut) && (TMath::Sqrt(fY*fY+fX*fX)<fRhoVertexCut)) {
 847         nMuons++;
 848         if(fEPvector.Pt() > fLowPtCut) {
 849           nMu1GeV++;
 850           if(fEPvector.Pt() > fHighPtCut) {
 851             nMu3GeV++;
 852             if (fEPvector.Pt() > fVeryHighPtCut) {
 853               nMu10GeV++;
 854             }
 855           }
 856         }
 857       }
 858     }//muon track loop
 859
 860     // Fill the event tags
 861     if(ntrack != 0) meanPt = meanPt/ntrack;
 862
 863     evTag->SetEventId(iEventNumber+1);
 864     evTag->SetPath(filepath);
 865
 866     evTag->SetVertexX(vertexIn->GetXv());
 867     evTag->SetVertexY(vertexIn->GetYv());
 868     evTag->SetVertexZ(vertexIn->GetZv());
 869     evTag->SetVertexZError(vertexIn->GetZRes());
 870     evTag->SetVertexFlag(fVertexflag);
 871
 872     evTag->SetT0VertexZ(esd->GetT0zVertex());
 873
 874     evTag->SetTriggerMask(esd->GetTriggerMask());
 875     evTag->SetTriggerCluster(esd->GetTriggerCluster());
 876
 877     evTag->SetZDCNeutron1Energy(esd->GetZDCN1Energy());
 878     evTag->SetZDCProton1Energy(esd->GetZDCP1Energy());
 879     evTag->SetZDCEMEnergy(esd->GetZDCEMEnergy());
 880     evTag->SetZDCNeutron1Energy(esd->GetZDCN2Energy());
 881     evTag->SetZDCProton1Energy(esd->GetZDCP2Energy());
 882     evTag->SetNumOfParticipants(esd->GetZDCParticipants());
 883
 884
 885     evTag->SetNumOfTracks(esd->GetNumberOfTracks());
 886     evTag->SetNumOfPosTracks(nPos);
 887     evTag->SetNumOfNegTracks(nNeg);
 888     evTag->SetNumOfNeutrTracks(nNeutr);
 889
 890     evTag->SetNumOfV0s(esd->GetNumberOfV0s());
 891     evTag->SetNumOfCascades(esd->GetNumberOfCascades());
 892     evTag->SetNumOfKinks(esd->GetNumberOfKinks());
 893     evTag->SetNumOfPMDTracks(esd->GetNumberOfPmdTracks());
 894
 895     evTag->SetNumOfProtons(nProtons);
 896     evTag->SetNumOfKaons(nKaons);
 897     evTag->SetNumOfPions(nPions);
 898     evTag->SetNumOfMuons(nMuons);
 899     evTag->SetNumOfElectrons(nElectrons);
 900     evTag->SetNumOfPhotons(nGamas);
 901     evTag->SetNumOfPi0s(nPi0s);
 902     evTag->SetNumOfNeutrons(nNeutrons);
 903     evTag->SetNumOfKaon0s(nK0s);
 904
 905     evTag->SetNumOfChargedAbove1GeV(nCh1GeV);
 906     evTag->SetNumOfChargedAbove3GeV(nCh3GeV);
 907     evTag->SetNumOfChargedAbove10GeV(nCh10GeV);
 908     evTag->SetNumOfMuonsAbove1GeV(nMu1GeV);
 909     evTag->SetNumOfMuonsAbove3GeV(nMu3GeV);
 910     evTag->SetNumOfMuonsAbove10GeV(nMu10GeV);
 911     evTag->SetNumOfElectronsAbove1GeV(nEl1GeV);
 912     evTag->SetNumOfElectronsAbove3GeV(nEl3GeV);
 913     evTag->SetNumOfElectronsAbove10GeV(nEl10GeV);
 914
 915     evTag->SetNumOfPHOSClusters(esd->GetNumberOfPHOSClusters());
 916     evTag->SetNumOfEMCALClusters(esd->GetNumberOfEMCALClusters());
 917
 918     evTag->SetTotalMomentum(totalP);
 919     evTag->SetMeanPt(meanPt);
 920     evTag->SetMaxPt(maxPt);
 921
 922     tag->SetRunId(iInitRunNumber);
 923     if(fIsSim) tag->SetDataType(0);
 924     else tag->SetDataType(1);
 925     tag->AddEventTag(*evTag);
 926   }//event loop
 927   lastEvent = iNumberOfEvents;
 928
 929   t->Delete("");
 930
 931   ttag.Fill();
 932   tag->Clear();
 933
 934   TString localFileName = "Run"; localFileName += tag->GetRunId();
 935   localFileName += ".Event"; localFileName += firstEvent; localFileName += "_"; localFileName += lastEvent; localFileName += "."; localFileName += Counter;
 936   localFileName += ".ESD.tag.root";
 937
 938   TString fileName;
 939
 940   if(fStorage == 0) {
 941     fileName = localFileName.Data();
 942     AliInfo(Form("Writing tags to local file: %s",fileName.Data()));
 943   }
 944   else if(fStorage == 1) {
 945     TString alienLocation = "/alien";
 946     alienLocation += gGrid->Pwd();
 947     alienLocation += fgridpath.Data();
 948     alienLocation += "/";
 949     alienLocation +=  localFileName;
 950     alienLocation += "?se=";
 951     alienLocation += fSE.Data();
 952     fileName = alienLocation.Data();
 953     AliInfo(Form("Writing tags to grid file: %s",fileName.Data()));
 954   }
 955
 956   TFile* ftag = TFile::Open(fileName, "recreate");
 957   ftag->cd();
 958   ttag.Write();
 959   ftag->Close();
 960
 961   delete ftag;
 962   delete esd;
 963
 964   delete tag;
 965   }*/
 966