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