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