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