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