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 <TList.h>
  30 #include <TObjString.h>
  31 #include <TLorentzVector.h>
  32 #include <TMap.h>
  33 #include <TTimeStamp.h>
  34
  35 //ROOT-AliEn
  36 #include <TGrid.h>
  37 #include <TGridResult.h>
  38
  39 //AliRoot
  40 #include "AliRunTag.h"
  41 #include "AliEventTag.h"
  42 #include "AliESD.h"
  43 #include "AliESDEvent.h"
  44 #include "AliESDVertex.h"
  45 #include "AliLog.h"
  46
  47 #include "AliESDTagCreator.h"
  48
  49
  50 ClassImp(AliESDTagCreator)
  51
  52
  53 //______________________________________________________________________________
  54   AliESDTagCreator::AliESDTagCreator() :
  55     AliTagCreator(),
  56     fChain(new TChain("esdTree")), fGUIDList(new TList()),
  57     fMD5List(new TList()), fTURLList(new TList()), fBranches(""),
  58     meminfo(new MemInfo_t) {
  59   //==============Default constructor for a AliESDTagCreator================
  60 }
  61
  62 //______________________________________________________________________________
  63 AliESDTagCreator::~AliESDTagCreator() {
  64 //================Default destructor for a AliESDTagCreator===================
  65   delete fChain;
  66   delete fGUIDList;
  67   delete fMD5List;
  68   delete fTURLList;
  69   delete meminfo;
  70 }
  71
  72 //______________________________________________________________________________
  73 Bool_t AliESDTagCreator::ReadGridCollection(TGridResult *fresult) {
  74   // Reads the entry of the TGridResult and creates the tags
  75   Int_t nEntries = fresult->GetEntries();
  76
  77   TString alienUrl;
  78   const char* guid;
  79   const char* md5;
  80   const char* turl;
  81   Long64_t size = -1;
  82
  83   Int_t counter = 0;
  84   for(Int_t i = 0; i < nEntries; i++) {
  85     alienUrl = fresult->GetKey(i,"turl");
  86     guid = fresult->GetKey(i,"guid");
  87     if(fresult->GetKey(i,"size")) size = atol (fresult->GetKey(i,"size"));
  88     md5 = fresult->GetKey(i,"md5");
  89     turl = fresult->GetKey(i,"turl");
  90     if(md5 && !strlen(guid)) md5 = 0;
  91     if(guid && !strlen(guid)) guid = 0;
  92
  93     fChain->Add(alienUrl);
  94     //fGUIDList->Add(new TObjString(guid));
  95     //fMD5List->Add(new TObjString(md5));
  96     //fTURLList->Add(new TObjString(turl));
  97
  98     //TFile *f = TFile::Open(alienUrl,"READ");
  99     //CreateTag(f,guid,md5,turl,size,counter);
 100     //f->Close();
 101     //delete f;
 102     counter += 1;
 103   }//grid result loop
 104
 105   AliInfo(Form("ESD chain created......."));
 106   AliInfo(Form("Chain entries: %d",fChain->GetEntries()));
 107   // Switch of branches on user request
 108   SwitchOffBranches();
 109   CreateTag(fChain,"grid");
 110
 111   return kTRUE;
 112 }
 113
 114 //______________________________________________________________________________
 115 Bool_t AliESDTagCreator::ReadLocalCollection(const char *localpath) {
 116   // Checks the different subdirs of the given local path and in the
 117   // case where it finds an AliESDs.root file it creates the tags
 118
 119   void *dira =  gSystem->OpenDirectory(localpath);
 120   Char_t fPath[256];
 121   const char * dirname = 0x0;
 122   const char * filename = 0x0;
 123   const char * pattern = "AliESDs.root";
 124
 125   Int_t counter = 0;
 126   while((dirname = gSystem->GetDirEntry(dira))) {
 127     sprintf(fPath,"%s/%s",localpath,dirname);
 128     void *dirb =  gSystem->OpenDirectory(fPath);
 129     while((filename = gSystem->GetDirEntry(dirb))) {
 130       if(strstr(filename,pattern)) {
 131         TString fESDFileName;
 132         fESDFileName = fPath;
 133         fESDFileName += "/";
 134         fESDFileName += pattern;
 135
 136         fChain->Add(fESDFileName);
 137
 138         //TFile *f = TFile::Open(fESDFileName,"READ");
 139         //CreateTag(f,fESDFileName,counter);
 140         //f->Close();
 141         //delete f;
 142
 143         counter += 1;
 144       }//pattern check
 145     }//child directory's entry loop
 146   }//parent directory's entry loop
 147
 148   AliInfo(Form("ESD chain created......."));
 149   AliInfo(Form("Chain entries: %d",fChain->GetEntries()));
 150   // Switch of branches on user request
 151   SwitchOffBranches();
 152   CreateTag(fChain,"local");
 153
 154   return kTRUE;
 155 }
 156
 157 //______________________________________________________________________________
 158 Bool_t AliESDTagCreator::ReadCAFCollection(const char *filename) {
 159   // Temporary solution for CAF: Takes as an input the ascii file that
 160   // lists the ESDs stored in the SE of the CAF and creates the tags.
 161
 162   // Open the input stream
 163   ifstream in;
 164   in.open(filename);
 165
 166   Int_t counter = 0;
 167   TString esdfile;
 168   // Read the input list of files and add them to the chain
 169   while(in.good()) {
 170     in >> esdfile;
 171     if (!esdfile.Contains("root")) continue; // protection
 172
 173     fChain->Add(esdfile);
 174
 175     //TFile *f = TFile::Open(esdfile,"READ");
 176     //CreateTag(f,esdfile,counter);
 177     //f->Close();
 178     //delete f;
 179
 180     counter += 1;
 181   }
 182
 183   AliInfo(Form("ESD chain created......."));
 184   AliInfo(Form("Chain entries: %d",fChain->GetEntries()));
 185   // Switch of branches on user request
 186   SwitchOffBranches();
 187   CreateTag(fChain,"proof");
 188
 189   return kTRUE;
 190 }
 191
 192 //_____________________________________________________________________________
 193 void AliESDTagCreator::CreateTag(TChain* chain, const char *type) {
 194   //private method that creates tag files
 195   TString fSession = type;
 196   TString fguid, fmd5, fturl;
 197   TString fTempGuid = 0;
 198
 199   /////////////
 200   //muon code//
 201   ////////////
 202   Double_t fMUONMASS = 0.105658369;
 203   //Variables
 204   Double_t fX,fY,fZ ;
 205   Double_t fThetaX, fThetaY, fPyz, fChisquare;
 206   Double_t fPxRec, fPyRec, fPzRec, fEnergy;
 207   Int_t fCharge;
 208   TLorentzVector fEPvector;
 209
 210   Float_t fZVertexCut = 10.0;
 211   Float_t fRhoVertexCut = 2.0;
 212
 213   Float_t fLowPtCut = 1.0;
 214   Float_t fHighPtCut = 3.0;
 215   Float_t fVeryHighPtCut = 10.0;
 216   ////////////
 217
 218   Double_t partFrac[5] = {0.01, 0.01, 0.85, 0.10, 0.05};
 219
 220   // Creates the tags for all the events in a given ESD file
 221   Bool_t fIsSim = kTRUE;
 222   Int_t ntrack;
 223   Int_t nProtons, nKaons, nPions, nMuons, nElectrons;
 224   Int_t nPos, nNeg, nNeutr;
 225   Int_t nK0s, nNeutrons, nPi0s, nGamas;
 226   Int_t nCh1GeV, nCh3GeV, nCh10GeV;
 227   Int_t nMu1GeV, nMu3GeV, nMu10GeV;
 228   Int_t nEl1GeV, nEl3GeV, nEl10GeV;
 229   Float_t maxPt = .0, meanPt = .0, totalP = .0;
 230   Int_t fVertexflag;
 231   Int_t iRunNumber = 0;
 232   TString fVertexName;
 233
 234   //gSystem->GetMemInfo(meminfo);
 235   //AliInfo(Form("After the tag initialization - Memory used: %d MB",meminfo->fMemUsed));
 236   //Int_t tempmem = meminfo->fMemUsed;
 237
 238   AliInfo(Form("Creating the ESD tags......."));
 239
 240   Int_t firstEvent = 0,lastEvent = 0;
 241   AliESDEvent *esd = new AliESDEvent();
 242   esd->ReadFromTree(chain);
 243   AliESD *esdold = 0x0;
 244
 245   //gSystem->GetMemInfo(meminfo);
 246   //AliInfo(Form("After the esd initialization - Memory used: %d MB - Increase: %d MB",meminfo->fMemUsed,meminfo->fMemUsed - tempmem));
 247   //tempmem = meminfo->fMemUsed;
 248
 249   Int_t iInitRunNumber = -1;
 250   chain->GetEntry(0);
 251   TFile *f = chain->GetFile();
 252   fTempGuid = f->GetUUID().AsString();
 253
 254   TString localFileName = "Run"; localFileName += esd->GetRunNumber();
 255   localFileName += ".Event"; localFileName += firstEvent; localFileName += "_"; localFileName += chain->GetEntries(); //localFileName += "."; localFileName += Counter;
 256   localFileName += ".ESD.tag.root";
 257
 258   TString fileName;
 259
 260   if(fStorage == 0) {
 261     fileName = localFileName.Data();
 262     AliInfo(Form("Writing tags to local file: %s",fileName.Data()));
 263   }
 264   else if(fStorage == 1) {
 265     TString alienLocation = "/alien";
 266     alienLocation += gGrid->Pwd();
 267     alienLocation += fgridpath.Data();
 268     alienLocation += "/";
 269     alienLocation +=  localFileName;
 270     alienLocation += "?se=";
 271     alienLocation += fSE.Data();
 272     fileName = alienLocation.Data();
 273     AliInfo(Form("Writing tags to grid file: %s",fileName.Data()));
 274   }
 275
 276   TFile* ftag = TFile::Open(fileName, "recreate");
 277
 278   AliRunTag *tag = new AliRunTag();
 279   AliEventTag *evTag = new AliEventTag();
 280   TTree ttag("T","A Tree with event tags");
 281   TBranch * btag = ttag.Branch("AliTAG", &tag);
 282   btag->SetCompressionLevel(9);
 283
 284   for(Int_t iEventNumber = 0; iEventNumber < chain->GetEntries(); iEventNumber++) {
 285     ntrack = 0; nPos = 0; nNeg = 0; nNeutr =0;
 286     nK0s = 0; nNeutrons = 0; nPi0s = 0;
 287     nGamas = 0; nProtons = 0; nKaons = 0;
 288     nPions = 0; nMuons = 0; nElectrons = 0;
 289     nCh1GeV = 0; nCh3GeV = 0; nCh10GeV = 0;
 290     nMu1GeV = 0; nMu3GeV = 0; nMu10GeV = 0;
 291     nEl1GeV = 0; nEl3GeV = 0; nEl10GeV = 0;
 292     maxPt = .0; meanPt = .0; totalP = .0;
 293     fVertexflag = 1;
 294
 295     chain->GetEntry(iEventNumber);
 296     esdold = esd->GetAliESDOld();
 297     if(esdold) esd->CopyFromOldESD();
 298
 299     TFile *file = chain->GetFile();
 300     const TUrl *url = file->GetEndpointUrl();
 301     fguid = file->GetUUID().AsString();
 302     if(fSession == "grid") {
 303       TString fturltemp = "alien://"; fturltemp += url->GetFile();
 304       fturl = fturltemp(0,fturltemp.Index(".root",5,0,TString::kExact)+5);
 305     }
 306     else fturl = url->GetFile();
 307
 308     if(iEventNumber == 0) iInitRunNumber = esd->GetRunNumber();
 309     iRunNumber = esd->GetRunNumber();
 310     if(iRunNumber != iInitRunNumber) AliFatal("Inconsistency of run numbers in the AliESD - You are trying to merge different runs!!!");
 311
 312     const AliESDVertex * vertexIn = esd->GetVertex();
 313     fVertexName = vertexIn->GetName();
 314     if(fVertexName == "default") fVertexflag = 0;
 315
 316     for (Int_t iTrackNumber = 0; iTrackNumber < esd->GetNumberOfTracks(); iTrackNumber++) {
 317       AliESDtrack * esdTrack = esd->GetTrack(iTrackNumber);
 318       if(esdTrack->GetLabel() != 0) fIsSim = kTRUE;
 319       else if(esdTrack->GetLabel() == 0) fIsSim = kFALSE;
 320       UInt_t status = esdTrack->GetStatus();
 321
 322       //select only tracks with ITS refit
 323       if ((status&AliESDtrack::kITSrefit)==0) continue;
 324       //select only tracks with TPC refit
 325       if ((status&AliESDtrack::kTPCrefit)==0) continue;
 326
 327       //select only tracks with the "combined PID"
 328       if ((status&AliESDtrack::kESDpid)==0) continue;
 329       Double_t p[3];
 330       esdTrack->GetPxPyPz(p);
 331       Double_t momentum = sqrt(pow(p[0],2) + pow(p[1],2) + pow(p[2],2));
 332       Double_t fPt = sqrt(pow(p[0],2) + pow(p[1],2));
 333       totalP += momentum;
 334       meanPt += fPt;
 335       if(fPt > maxPt) maxPt = fPt;
 336
 337       if(esdTrack->GetSign() > 0) {
 338         nPos++;
 339         if(fPt > fLowPtCut) nCh1GeV++;
 340         if(fPt > fHighPtCut) nCh3GeV++;
 341         if(fPt > fVeryHighPtCut) nCh10GeV++;
 342       }
 343       if(esdTrack->GetSign() < 0) {
 344         nNeg++;
 345         if(fPt > fLowPtCut) nCh1GeV++;
 346         if(fPt > fHighPtCut) nCh3GeV++;
 347         if(fPt > fVeryHighPtCut) nCh10GeV++;
 348       }
 349       if(esdTrack->GetSign() == 0) nNeutr++;
 350
 351       //PID
 352       Double_t prob[5];
 353       esdTrack->GetESDpid(prob);
 354
 355       Double_t rcc = 0.0;
 356       for(Int_t i = 0; i < AliPID::kSPECIES; i++) rcc += prob[i]*partFrac[i];
 357       if(rcc == 0.0) continue;
 358       //Bayes' formula
 359       Double_t w[5];
 360       for(Int_t i = 0; i < AliPID::kSPECIES; i++) w[i] = prob[i]*partFrac[i]/rcc;
 361
 362       //protons
 363       if ((w[4]>w[3])&&(w[4]>w[2])&&(w[4]>w[1])&&(w[4]>w[0])) nProtons++;
 364       //kaons
 365       if ((w[3]>w[4])&&(w[3]>w[2])&&(w[3]>w[1])&&(w[3]>w[0])) nKaons++;
 366       //pions
 367       if ((w[2]>w[4])&&(w[2]>w[3])&&(w[2]>w[1])&&(w[2]>w[0])) nPions++;
 368       //electrons
 369       if ((w[0]>w[4])&&(w[0]>w[3])&&(w[0]>w[2])&&(w[0]>w[1])) {
 370         nElectrons++;
 371         if(fPt > fLowPtCut) nEl1GeV++;
 372         if(fPt > fHighPtCut) nEl3GeV++;
 373         if(fPt > fVeryHighPtCut) nEl10GeV++;
 374       }
 375       ntrack++;
 376     }//esd track loop
 377
 378     /////////////
 379     //muon code//
 380     ////////////
 381     Int_t nMuonTracks = esd->GetNumberOfMuonTracks();
 382     // loop over all reconstructed tracks (also first track of combination)
 383     for (Int_t iTrack = 0; iTrack <  nMuonTracks;  iTrack++) {
 384       AliESDMuonTrack* muonTrack = esd->GetMuonTrack(iTrack);
 385       if (muonTrack == 0x0) continue;
 386
 387       // Coordinates at vertex
 388       fZ = muonTrack->GetZ();
 389       fY = muonTrack->GetBendingCoor();
 390       fX = muonTrack->GetNonBendingCoor();
 391
 392       fThetaX = muonTrack->GetThetaX();
 393       fThetaY = muonTrack->GetThetaY();
 394
 395       fPyz = 1./TMath::Abs(muonTrack->GetInverseBendingMomentum());
 396       fPzRec = - fPyz / TMath::Sqrt(1.0 + TMath::Tan(fThetaY)*TMath::Tan(fThetaY));
 397       fPxRec = fPzRec * TMath::Tan(fThetaX);
 398       fPyRec = fPzRec * TMath::Tan(fThetaY);
 399       fCharge = Int_t(TMath::Sign(1.,muonTrack->GetInverseBendingMomentum()));
 400
 401       //ChiSquare of the track if needed
 402       fChisquare = muonTrack->GetChi2()/(2.0 * muonTrack->GetNHit() - 5);
 403       fEnergy = TMath::Sqrt(fMUONMASS * fMUONMASS + fPxRec * fPxRec + fPyRec * fPyRec + fPzRec * fPzRec);
 404       fEPvector.SetPxPyPzE(fPxRec, fPyRec, fPzRec, fEnergy);
 405
 406       // total number of muons inside a vertex cut
 407       if((TMath::Abs(fZ)<fZVertexCut) && (TMath::Sqrt(fY*fY+fX*fX)<fRhoVertexCut)) {
 408         nMuons++;
 409         if(fEPvector.Pt() > fLowPtCut) {
 410           nMu1GeV++;
 411           if(fEPvector.Pt() > fHighPtCut) {
 412             nMu3GeV++;
 413             if (fEPvector.Pt() > fVeryHighPtCut) {
 414               nMu10GeV++;
 415             }
 416           }
 417         }
 418       }
 419     }//muon track loop
 420
 421     // Fill the event tags
 422     if(ntrack != 0) meanPt = meanPt/ntrack;
 423
 424     //AliInfo(Form("====================================="));
 425     //AliInfo(Form("URL: %s - GUID: %s",fturl.Data(),fguid.Data()));
 426     //AliInfo(Form("====================================="));
 427
 428     evTag->SetEventId(iEventNumber+1);
 429     evTag->SetGUID(fguid);
 430     if(fSession == "grid") {
 431       evTag->SetMD5(0);
 432       evTag->SetTURL(fturl);
 433       evTag->SetSize(0);
 434     }
 435     else evTag->SetPath(fturl);
 436
 437     evTag->SetVertexX(vertexIn->GetXv());
 438     evTag->SetVertexY(vertexIn->GetYv());
 439     evTag->SetVertexZ(vertexIn->GetZv());
 440     evTag->SetVertexZError(vertexIn->GetZRes());
 441     evTag->SetVertexFlag(fVertexflag);
 442
 443     evTag->SetT0VertexZ(esd->GetT0zVertex());
 444
 445     evTag->SetTriggerMask(esd->GetTriggerMask());
 446     evTag->SetTriggerCluster(esd->GetTriggerCluster());
 447
 448     evTag->SetZDCNeutron1Energy(esd->GetZDCN1Energy());
 449     evTag->SetZDCProton1Energy(esd->GetZDCP1Energy());
 450     evTag->SetZDCEMEnergy(esd->GetZDCEMEnergy(0),esd->GetZDCEMEnergy(1));
 451     evTag->SetZDCNeutron1Energy(esd->GetZDCN2Energy());
 452     evTag->SetZDCProton1Energy(esd->GetZDCP2Energy());
 453     evTag->SetNumOfParticipants(esd->GetZDCParticipants());
 454
 455
 456     evTag->SetNumOfTracks(esd->GetNumberOfTracks());
 457     evTag->SetNumOfPosTracks(nPos);
 458     evTag->SetNumOfNegTracks(nNeg);
 459     evTag->SetNumOfNeutrTracks(nNeutr);
 460
 461     evTag->SetNumOfV0s(esd->GetNumberOfV0s());
 462     evTag->SetNumOfCascades(esd->GetNumberOfCascades());
 463     evTag->SetNumOfKinks(esd->GetNumberOfKinks());
 464     evTag->SetNumOfPMDTracks(esd->GetNumberOfPmdTracks());
 465
 466     evTag->SetNumOfProtons(nProtons);
 467     evTag->SetNumOfKaons(nKaons);
 468     evTag->SetNumOfPions(nPions);
 469     evTag->SetNumOfMuons(nMuons);
 470     evTag->SetNumOfElectrons(nElectrons);
 471     evTag->SetNumOfPhotons(nGamas);
 472     evTag->SetNumOfPi0s(nPi0s);
 473     evTag->SetNumOfNeutrons(nNeutrons);
 474     evTag->SetNumOfKaon0s(nK0s);
 475
 476     evTag->SetNumOfChargedAbove1GeV(nCh1GeV);
 477     evTag->SetNumOfChargedAbove3GeV(nCh3GeV);
 478     evTag->SetNumOfChargedAbove10GeV(nCh10GeV);
 479     evTag->SetNumOfMuonsAbove1GeV(nMu1GeV);
 480     evTag->SetNumOfMuonsAbove3GeV(nMu3GeV);
 481     evTag->SetNumOfMuonsAbove10GeV(nMu10GeV);
 482     evTag->SetNumOfElectronsAbove1GeV(nEl1GeV);
 483     evTag->SetNumOfElectronsAbove3GeV(nEl3GeV);
 484     evTag->SetNumOfElectronsAbove10GeV(nEl10GeV);
 485
 486     evTag->SetNumOfPHOSClusters(esd->GetNumberOfPHOSClusters());
 487     evTag->SetNumOfEMCALClusters(esd->GetNumberOfEMCALClusters());
 488
 489     evTag->SetTotalMomentum(totalP);
 490     evTag->SetMeanPt(meanPt);
 491     evTag->SetMaxPt(maxPt);
 492
 493     tag->SetRunId(iInitRunNumber);
 494     if(fIsSim) tag->SetDataType(0);
 495     else tag->SetDataType(1);
 496
 497     if(fguid != fTempGuid) {
 498       fTempGuid = fguid;
 499       ttag.Fill();
 500       tag->Clear("");
 501     }
 502     tag->AddEventTag(*evTag);
 503     if(iEventNumber+1 == chain->GetEntries()) {
 504       //AliInfo(Form("File: %s",fturl.Data()));
 505       ttag.Fill();
 506       tag->Clear("");
 507     }
 508   }//event loop
 509   lastEvent = chain->GetEntries();
 510
 511   //gSystem->GetMemInfo(meminfo);
 512   //AliInfo(Form("After the event and track loop - Memory used: %d MB - Increase: %d MB",meminfo->fMemUsed,meminfo->fMemUsed - tempmem));
 513   //tempmem = meminfo->fMemUsed;
 514
 515   //chain->Delete("");
 516
 517   //gSystem->GetMemInfo(meminfo);
 518   //AliInfo(Form("After the t->Delete - Memory used: %d MB - Increase: %d MB",meminfo->fMemUsed,meminfo->fMemUsed - tempmem));
 519   //tempmem = meminfo->fMemUsed;
 520
 521   ftag->cd();
 522   tag->Clear();
 523   ttag.Write();
 524   ftag->Close();
 525
 526   //gSystem->GetMemInfo(meminfo);
 527   //AliInfo(Form("After the file closing - Memory used: %d MB - Increase: %d MB",meminfo->fMemUsed,meminfo->fMemUsed - tempmem));
 528   //tempmem = meminfo->fMemUsed;
 529
 530   delete esd;
 531   delete tag;
 532
 533   //gSystem->GetMemInfo(meminfo);
 534   //AliInfo(Form("After the delete objects - Memory used: %d MB - Increase: %d MB",meminfo->fMemUsed,meminfo->fMemUsed - tempmem));
 535 }
 536
 537 //_____________________________________________________________________________
 538 void AliESDTagCreator::CreateTag(TFile* file, const char *guid, const char *md5, const char *turl, Long64_t size, Int_t Counter) {
 539   //private method that creates tag files
 540   TString fguid = guid;
 541   TString fmd5 = md5;
 542   TString fturl = turl;
 543   /////////////
 544   //muon code//
 545   ////////////
 546   Double_t fMUONMASS = 0.105658369;
 547   //Variables
 548   Double_t fX,fY,fZ ;
 549   Double_t fThetaX, fThetaY, fPyz, fChisquare;
 550   Double_t fPxRec, fPyRec, fPzRec, fEnergy;
 551   Int_t fCharge;
 552   TLorentzVector fEPvector;
 553
 554   Float_t fZVertexCut = 10.0;
 555   Float_t fRhoVertexCut = 2.0;
 556
 557   Float_t fLowPtCut = 1.0;
 558   Float_t fHighPtCut = 3.0;
 559   Float_t fVeryHighPtCut = 10.0;
 560   ////////////
 561
 562   Double_t partFrac[5] = {0.01, 0.01, 0.85, 0.10, 0.05};
 563
 564   // Creates the tags for all the events in a given ESD file
 565   Bool_t fIsSim = kTRUE;
 566   Int_t ntrack;
 567   Int_t nProtons, nKaons, nPions, nMuons, nElectrons;
 568   Int_t nPos, nNeg, nNeutr;
 569   Int_t nK0s, nNeutrons, nPi0s, nGamas;
 570   Int_t nCh1GeV, nCh3GeV, nCh10GeV;
 571   Int_t nMu1GeV, nMu3GeV, nMu10GeV;
 572   Int_t nEl1GeV, nEl3GeV, nEl10GeV;
 573   Float_t maxPt = .0, meanPt = .0, totalP = .0;
 574   Int_t fVertexflag;
 575   Int_t iRunNumber = 0;
 576   TString fVertexName;
 577
 578   AliRunTag *tag = new AliRunTag();
 579   AliEventTag *evTag = new AliEventTag();
 580   TTree ttag("T","A Tree with event tags");
 581   TBranch * btag = ttag.Branch("AliTAG", &tag);
 582   btag->SetCompressionLevel(9);
 583   gSystem->GetMemInfo(meminfo);
 584   AliInfo(Form("After the tag initialization - Memory used: %d MB",meminfo->fMemUsed));
 585   Int_t tempmem = meminfo->fMemUsed;
 586
 587   AliInfo(Form("Creating the ESD tags......."));
 588
 589   Int_t firstEvent = 0,lastEvent = 0;
 590   TTree *t = (TTree*) file->Get("esdTree");
 591   AliESDEvent *esd = new AliESDEvent();
 592   esd->ReadFromTree(t);
 593
 594   gSystem->GetMemInfo(meminfo);
 595   AliInfo(Form("After the esd initialization - Memory used: %d MB - Increase: %d MB",meminfo->fMemUsed,meminfo->fMemUsed - tempmem));
 596   tempmem = meminfo->fMemUsed;
 597
 598   t->GetEntry(0);
 599   Int_t iInitRunNumber = esd->GetRunNumber();
 600
 601   Int_t iNumberOfEvents = (Int_t)t->GetEntries();
 602   for (Int_t iEventNumber = 0; iEventNumber < iNumberOfEvents; iEventNumber++) {
 603     ntrack = 0;
 604     nPos = 0;
 605     nNeg = 0;
 606     nNeutr =0;
 607     nK0s = 0;
 608     nNeutrons = 0;
 609     nPi0s = 0;
 610     nGamas = 0;
 611     nProtons = 0;
 612     nKaons = 0;
 613     nPions = 0;
 614     nMuons = 0;
 615     nElectrons = 0;
 616     nCh1GeV = 0;
 617     nCh3GeV = 0;
 618     nCh10GeV = 0;
 619     nMu1GeV = 0;
 620     nMu3GeV = 0;
 621     nMu10GeV = 0;
 622     nEl1GeV = 0;
 623     nEl3GeV = 0;
 624     nEl10GeV = 0;
 625     maxPt = .0;
 626     meanPt = .0;
 627     totalP = .0;
 628     fVertexflag = 1;
 629
 630     t->GetEntry(iEventNumber);
 631     iRunNumber = esd->GetRunNumber();
 632     if(iRunNumber != iInitRunNumber) AliFatal("Inconsistency of run numbers in the AliESD!!!");
 633     const AliESDVertex * vertexIn = esd->GetVertex();
 634     fVertexName = vertexIn->GetName();
 635     if(fVertexName == "default") fVertexflag = 0;
 636
 637     for (Int_t iTrackNumber = 0; iTrackNumber < esd->GetNumberOfTracks(); iTrackNumber++) {
 638       AliESDtrack * esdTrack = esd->GetTrack(iTrackNumber);
 639       if(esdTrack->GetLabel() != 0) fIsSim = kTRUE;
 640       else if(esdTrack->GetLabel() == 0) fIsSim = kFALSE;
 641       UInt_t status = esdTrack->GetStatus();
 642
 643       //select only tracks with ITS refit
 644       if ((status&AliESDtrack::kITSrefit)==0) continue;
 645       //select only tracks with TPC refit
 646       if ((status&AliESDtrack::kTPCrefit)==0) continue;
 647
 648       //select only tracks with the "combined PID"
 649       if ((status&AliESDtrack::kESDpid)==0) continue;
 650       Double_t p[3];
 651       esdTrack->GetPxPyPz(p);
 652       Double_t momentum = sqrt(pow(p[0],2) + pow(p[1],2) + pow(p[2],2));
 653       Double_t fPt = sqrt(pow(p[0],2) + pow(p[1],2));
 654       totalP += momentum;
 655       meanPt += fPt;
 656       if(fPt > maxPt) maxPt = fPt;
 657
 658       if(esdTrack->GetSign() > 0) {
 659         nPos++;
 660         if(fPt > fLowPtCut) nCh1GeV++;
 661         if(fPt > fHighPtCut) nCh3GeV++;
 662         if(fPt > fVeryHighPtCut) nCh10GeV++;
 663       }
 664       if(esdTrack->GetSign() < 0) {
 665         nNeg++;
 666         if(fPt > fLowPtCut) nCh1GeV++;
 667         if(fPt > fHighPtCut) nCh3GeV++;
 668         if(fPt > fVeryHighPtCut) nCh10GeV++;
 669       }
 670       if(esdTrack->GetSign() == 0) nNeutr++;
 671
 672       //PID
 673       Double_t prob[5];
 674       esdTrack->GetESDpid(prob);
 675
 676       Double_t rcc = 0.0;
 677       for(Int_t i = 0; i < AliPID::kSPECIES; i++) rcc += prob[i]*partFrac[i];
 678       if(rcc == 0.0) continue;
 679       //Bayes' formula
 680       Double_t w[5];
 681       for(Int_t i = 0; i < AliPID::kSPECIES; i++) w[i] = prob[i]*partFrac[i]/rcc;
 682
 683       //protons
 684       if ((w[4]>w[3])&&(w[4]>w[2])&&(w[4]>w[1])&&(w[4]>w[0])) nProtons++;
 685       //kaons
 686       if ((w[3]>w[4])&&(w[3]>w[2])&&(w[3]>w[1])&&(w[3]>w[0])) nKaons++;
 687       //pions
 688       if ((w[2]>w[4])&&(w[2]>w[3])&&(w[2]>w[1])&&(w[2]>w[0])) nPions++;
 689       //electrons
 690       if ((w[0]>w[4])&&(w[0]>w[3])&&(w[0]>w[2])&&(w[0]>w[1])) {
 691         nElectrons++;
 692         if(fPt > fLowPtCut) nEl1GeV++;
 693         if(fPt > fHighPtCut) nEl3GeV++;
 694         if(fPt > fVeryHighPtCut) nEl10GeV++;
 695       }
 696       ntrack++;
 697     }//esd track loop
 698
 699     /////////////
 700     //muon code//
 701     ////////////
 702     Int_t nMuonTracks = esd->GetNumberOfMuonTracks();
 703     // loop over all reconstructed tracks (also first track of combination)
 704     for (Int_t iTrack = 0; iTrack <  nMuonTracks;  iTrack++) {
 705       AliESDMuonTrack* muonTrack = esd->GetMuonTrack(iTrack);
 706       if (muonTrack == 0x0) continue;
 707
 708       // Coordinates at vertex
 709       fZ = muonTrack->GetZ();
 710       fY = muonTrack->GetBendingCoor();
 711       fX = muonTrack->GetNonBendingCoor();
 712
 713       fThetaX = muonTrack->GetThetaX();
 714       fThetaY = muonTrack->GetThetaY();
 715
 716       fPyz = 1./TMath::Abs(muonTrack->GetInverseBendingMomentum());
 717       fPzRec = - fPyz / TMath::Sqrt(1.0 + TMath::Tan(fThetaY)*TMath::Tan(fThetaY));
 718       fPxRec = fPzRec * TMath::Tan(fThetaX);
 719       fPyRec = fPzRec * TMath::Tan(fThetaY);
 720       fCharge = Int_t(TMath::Sign(1.,muonTrack->GetInverseBendingMomentum()));
 721
 722       //ChiSquare of the track if needed
 723       fChisquare = muonTrack->GetChi2()/(2.0 * muonTrack->GetNHit() - 5);
 724       fEnergy = TMath::Sqrt(fMUONMASS * fMUONMASS + fPxRec * fPxRec + fPyRec * fPyRec + fPzRec * fPzRec);
 725       fEPvector.SetPxPyPzE(fPxRec, fPyRec, fPzRec, fEnergy);
 726
 727       // total number of muons inside a vertex cut
 728       if((TMath::Abs(fZ)<fZVertexCut) && (TMath::Sqrt(fY*fY+fX*fX)<fRhoVertexCut)) {
 729         nMuons++;
 730         if(fEPvector.Pt() > fLowPtCut) {
 731           nMu1GeV++;
 732           if(fEPvector.Pt() > fHighPtCut) {
 733             nMu3GeV++;
 734             if (fEPvector.Pt() > fVeryHighPtCut) {
 735               nMu10GeV++;
 736             }
 737           }
 738         }
 739       }
 740     }//muon track loop
 741
 742     // Fill the event tags
 743     if(ntrack != 0) meanPt = meanPt/ntrack;
 744
 745     evTag->SetEventId(iEventNumber+1);
 746     evTag->SetGUID(fguid);
 747     evTag->SetMD5(fmd5);
 748     evTag->SetTURL(fturl);
 749     evTag->SetSize(size);
 750     evTag->SetVertexX(vertexIn->GetXv());
 751     evTag->SetVertexY(vertexIn->GetYv());
 752     evTag->SetVertexZ(vertexIn->GetZv());
 753     evTag->SetVertexZError(vertexIn->GetZRes());
 754     evTag->SetVertexFlag(fVertexflag);
 755
 756     evTag->SetT0VertexZ(esd->GetT0zVertex());
 757
 758     evTag->SetTriggerMask(esd->GetTriggerMask());
 759     evTag->SetTriggerCluster(esd->GetTriggerCluster());
 760
 761     evTag->SetZDCNeutron1Energy(esd->GetZDCN1Energy());
 762     evTag->SetZDCProton1Energy(esd->GetZDCP1Energy());
 763     evTag->SetZDCEMEnergy(esd->GetZDCEMEnergy(0),esd->GetZDCEMEnergy(1));
 764     evTag->SetZDCNeutron1Energy(esd->GetZDCN2Energy());
 765     evTag->SetZDCProton1Energy(esd->GetZDCP2Energy());
 766     evTag->SetNumOfParticipants(esd->GetZDCParticipants());
 767
 768
 769     evTag->SetNumOfTracks(esd->GetNumberOfTracks());
 770     evTag->SetNumOfPosTracks(nPos);
 771     evTag->SetNumOfNegTracks(nNeg);
 772     evTag->SetNumOfNeutrTracks(nNeutr);
 773
 774     evTag->SetNumOfV0s(esd->GetNumberOfV0s());
 775     evTag->SetNumOfCascades(esd->GetNumberOfCascades());
 776     evTag->SetNumOfKinks(esd->GetNumberOfKinks());
 777     evTag->SetNumOfPMDTracks(esd->GetNumberOfPmdTracks());
 778
 779     evTag->SetNumOfProtons(nProtons);
 780     evTag->SetNumOfKaons(nKaons);
 781     evTag->SetNumOfPions(nPions);
 782     evTag->SetNumOfMuons(nMuons);
 783     evTag->SetNumOfElectrons(nElectrons);
 784     evTag->SetNumOfPhotons(nGamas);
 785     evTag->SetNumOfPi0s(nPi0s);
 786     evTag->SetNumOfNeutrons(nNeutrons);
 787     evTag->SetNumOfKaon0s(nK0s);
 788
 789     evTag->SetNumOfChargedAbove1GeV(nCh1GeV);
 790     evTag->SetNumOfChargedAbove3GeV(nCh3GeV);
 791     evTag->SetNumOfChargedAbove10GeV(nCh10GeV);
 792     evTag->SetNumOfMuonsAbove1GeV(nMu1GeV);
 793     evTag->SetNumOfMuonsAbove3GeV(nMu3GeV);
 794     evTag->SetNumOfMuonsAbove10GeV(nMu10GeV);
 795     evTag->SetNumOfElectronsAbove1GeV(nEl1GeV);
 796     evTag->SetNumOfElectronsAbove3GeV(nEl3GeV);
 797     evTag->SetNumOfElectronsAbove10GeV(nEl10GeV);
 798
 799     evTag->SetNumOfPHOSClusters(esd->GetNumberOfPHOSClusters());
 800     evTag->SetNumOfEMCALClusters(esd->GetNumberOfEMCALClusters());
 801
 802     evTag->SetTotalMomentum(totalP);
 803     evTag->SetMeanPt(meanPt);
 804     evTag->SetMaxPt(maxPt);
 805
 806     tag->SetRunId(iInitRunNumber);
 807     if(fIsSim) tag->SetDataType(0);
 808     else tag->SetDataType(1);
 809     tag->AddEventTag(*evTag);
 810   }//event loop
 811   lastEvent = iNumberOfEvents;
 812
 813   gSystem->GetMemInfo(meminfo);
 814   AliInfo(Form("After the event and track loop - Memory used: %d MB - Increase: %d MB",meminfo->fMemUsed,meminfo->fMemUsed - tempmem));
 815   tempmem = meminfo->fMemUsed;
 816   t->Delete("");
 817
 818   gSystem->GetMemInfo(meminfo);
 819   AliInfo(Form("After the t->Delete - Memory used: %d MB - Increase: %d MB",meminfo->fMemUsed,meminfo->fMemUsed - tempmem));
 820   tempmem = meminfo->fMemUsed;
 821
 822   TString localFileName = "Run"; localFileName += tag->GetRunId();
 823   localFileName += ".Event"; localFileName += firstEvent; localFileName += "_"; localFileName += lastEvent; localFileName += "."; localFileName += Counter;
 824   localFileName += ".ESD.tag.root";
 825
 826   TString fileName;
 827
 828   if(fStorage == 0) {
 829     fileName = localFileName.Data();
 830     AliInfo(Form("Writing tags to local file: %s",fileName.Data()));
 831   }
 832   else if(fStorage == 1) {
 833     TString alienLocation = "/alien";
 834     alienLocation += gGrid->Pwd();
 835     alienLocation += fgridpath.Data();
 836     alienLocation += "/";
 837     alienLocation +=  localFileName;
 838     alienLocation += "?se=";
 839     alienLocation += fSE.Data();
 840     fileName = alienLocation.Data();
 841     AliInfo(Form("Writing tags to grid file: %s",fileName.Data()));
 842   }
 843
 844   TFile* ftag = TFile::Open(fileName, "recreate");
 845   ftag->cd();
 846   ttag.Fill();
 847   tag->Clear();
 848   ttag.Write();
 849   ftag->Close();
 850
 851   gSystem->GetMemInfo(meminfo);
 852   AliInfo(Form("After the file closing - Memory used: %d MB - Increase: %d MB",meminfo->fMemUsed,meminfo->fMemUsed - tempmem));
 853   tempmem = meminfo->fMemUsed;
 854
 855   delete ftag;
 856   delete esd;
 857
 858   delete tag;
 859   gSystem->GetMemInfo(meminfo);
 860   AliInfo(Form("After the delete objects - Memory used: %d MB - Increase: %d MB",meminfo->fMemUsed,meminfo->fMemUsed - tempmem));
 861 }
 862
 863 //_____________________________________________________________________________
 864 void AliESDTagCreator::CreateTag(TFile* file, const char *filepath, Int_t Counter) {
 865   //private method that creates tag files
 866   /////////////
 867   //muon code//
 868   ////////////
 869   Double_t fMUONMASS = 0.105658369;
 870   //Variables
 871   Double_t fX,fY,fZ ;
 872   Double_t fThetaX, fThetaY, fPyz, fChisquare;
 873   Double_t fPxRec, fPyRec, fPzRec, fEnergy;
 874   Int_t fCharge;
 875   TLorentzVector fEPvector;
 876
 877   Float_t fZVertexCut = 10.0;
 878   Float_t fRhoVertexCut = 2.0;
 879
 880   Float_t fLowPtCut = 1.0;
 881   Float_t fHighPtCut = 3.0;
 882   Float_t fVeryHighPtCut = 10.0;
 883   ////////////
 884
 885   Double_t partFrac[5] = {0.01, 0.01, 0.85, 0.10, 0.05};
 886
 887   // Creates the tags for all the events in a given ESD file
 888   Bool_t fIsSim = kTRUE;
 889   Int_t ntrack;
 890   Int_t nProtons, nKaons, nPions, nMuons, nElectrons;
 891   Int_t nPos, nNeg, nNeutr;
 892   Int_t nK0s, nNeutrons, nPi0s, nGamas;
 893   Int_t nCh1GeV, nCh3GeV, nCh10GeV;
 894   Int_t nMu1GeV, nMu3GeV, nMu10GeV;
 895   Int_t nEl1GeV, nEl3GeV, nEl10GeV;
 896   Float_t maxPt = .0, meanPt = .0, totalP = .0;
 897   Int_t fVertexflag;
 898   Int_t iRunNumber = 0;
 899   TString fVertexName;
 900
 901   AliRunTag *tag = new AliRunTag();
 902   AliEventTag *evTag = new AliEventTag();
 903   TTree ttag("T","A Tree with event tags");
 904   TBranch * btag = ttag.Branch("AliTAG", &tag);
 905   btag->SetCompressionLevel(9);
 906
 907   AliInfo(Form("Creating the ESD tags......."));
 908
 909   Int_t firstEvent = 0,lastEvent = 0;
 910
 911   TTree *t = (TTree*) file->Get("esdTree");
 912   AliESDEvent *esd = new AliESDEvent();
 913   esd->ReadFromTree(t);
 914
 915   t->GetEntry(0);
 916   Int_t iInitRunNumber = esd->GetRunNumber();
 917
 918   Int_t iNumberOfEvents = (Int_t)t->GetEntries();
 919   for (Int_t iEventNumber = 0; iEventNumber < iNumberOfEvents; iEventNumber++) {
 920     ntrack = 0;
 921     nPos = 0;
 922     nNeg = 0;
 923     nNeutr =0;
 924     nK0s = 0;
 925     nNeutrons = 0;
 926     nPi0s = 0;
 927     nGamas = 0;
 928     nProtons = 0;
 929     nKaons = 0;
 930     nPions = 0;
 931     nMuons = 0;
 932     nElectrons = 0;
 933     nCh1GeV = 0;
 934     nCh3GeV = 0;
 935     nCh10GeV = 0;
 936     nMu1GeV = 0;
 937     nMu3GeV = 0;
 938     nMu10GeV = 0;
 939     nEl1GeV = 0;
 940     nEl3GeV = 0;
 941     nEl10GeV = 0;
 942     maxPt = .0;
 943     meanPt = .0;
 944     totalP = .0;
 945     fVertexflag = 1;
 946
 947     t->GetEntry(iEventNumber);
 948     iRunNumber = esd->GetRunNumber();
 949     if(iRunNumber != iInitRunNumber) AliFatal("Inconsistency of run numbers in the AliESD!!!");
 950     const AliESDVertex * vertexIn = esd->GetVertex();
 951     fVertexName = vertexIn->GetName();
 952     if(fVertexName == "default") fVertexflag = 0;
 953
 954     for (Int_t iTrackNumber = 0; iTrackNumber < esd->GetNumberOfTracks(); iTrackNumber++) {
 955       AliESDtrack * esdTrack = esd->GetTrack(iTrackNumber);
 956       if(esdTrack->GetLabel() != 0) fIsSim = kTRUE;
 957       else if(esdTrack->GetLabel() == 0) fIsSim = kFALSE;
 958       UInt_t status = esdTrack->GetStatus();
 959
 960       //select only tracks with ITS refit
 961       if ((status&AliESDtrack::kITSrefit)==0) continue;
 962       //select only tracks with TPC refit
 963       if ((status&AliESDtrack::kTPCrefit)==0) continue;
 964
 965       //select only tracks with the "combined PID"
 966       if ((status&AliESDtrack::kESDpid)==0) continue;
 967       Double_t p[3];
 968       esdTrack->GetPxPyPz(p);
 969       Double_t momentum = sqrt(pow(p[0],2) + pow(p[1],2) + pow(p[2],2));
 970       Double_t fPt = sqrt(pow(p[0],2) + pow(p[1],2));
 971       totalP += momentum;
 972       meanPt += fPt;
 973       if(fPt > maxPt) maxPt = fPt;
 974
 975       if(esdTrack->GetSign() > 0) {
 976         nPos++;
 977         if(fPt > fLowPtCut) nCh1GeV++;
 978         if(fPt > fHighPtCut) nCh3GeV++;
 979         if(fPt > fVeryHighPtCut) nCh10GeV++;
 980       }
 981       if(esdTrack->GetSign() < 0) {
 982         nNeg++;
 983         if(fPt > fLowPtCut) nCh1GeV++;
 984         if(fPt > fHighPtCut) nCh3GeV++;
 985         if(fPt > fVeryHighPtCut) nCh10GeV++;
 986       }
 987       if(esdTrack->GetSign() == 0) nNeutr++;
 988
 989       //PID
 990       Double_t prob[5];
 991       esdTrack->GetESDpid(prob);
 992
 993       Double_t rcc = 0.0;
 994       for(Int_t i = 0; i < AliPID::kSPECIES; i++) rcc += prob[i]*partFrac[i];
 995       if(rcc == 0.0) continue;
 996       //Bayes' formula
 997       Double_t w[5];
 998       for(Int_t i = 0; i < AliPID::kSPECIES; i++) w[i] = prob[i]*partFrac[i]/rcc;
 999
1000       //protons
1001       if ((w[4]>w[3])&&(w[4]>w[2])&&(w[4]>w[1])&&(w[4]>w[0])) nProtons++;
1002       //kaons
1003       if ((w[3]>w[4])&&(w[3]>w[2])&&(w[3]>w[1])&&(w[3]>w[0])) nKaons++;
1004       //pions
1005       if ((w[2]>w[4])&&(w[2]>w[3])&&(w[2]>w[1])&&(w[2]>w[0])) nPions++;
1006       //electrons
1007       if ((w[0]>w[4])&&(w[0]>w[3])&&(w[0]>w[2])&&(w[0]>w[1])) {
1008         nElectrons++;
1009         if(fPt > fLowPtCut) nEl1GeV++;
1010         if(fPt > fHighPtCut) nEl3GeV++;
1011         if(fPt > fVeryHighPtCut) nEl10GeV++;
1012       }
1013       ntrack++;
1014     }//esd track loop
1015
1016     /////////////
1017     //muon code//
1018     ////////////
1019     Int_t nMuonTracks = esd->GetNumberOfMuonTracks();
1020     // loop over all reconstructed tracks (also first track of combination)
1021     for (Int_t iTrack = 0; iTrack <  nMuonTracks;  iTrack++) {
1022       AliESDMuonTrack* muonTrack = esd->GetMuonTrack(iTrack);
1023       if (muonTrack == 0x0) continue;
1024
1025       // Coordinates at vertex
1026       fZ = muonTrack->GetZ();
1027       fY = muonTrack->GetBendingCoor();
1028       fX = muonTrack->GetNonBendingCoor();
1029
1030       fThetaX = muonTrack->GetThetaX();
1031       fThetaY = muonTrack->GetThetaY();
1032
1033       fPyz = 1./TMath::Abs(muonTrack->GetInverseBendingMomentum());
1034       fPzRec = - fPyz / TMath::Sqrt(1.0 + TMath::Tan(fThetaY)*TMath::Tan(fThetaY));
1035       fPxRec = fPzRec * TMath::Tan(fThetaX);
1036       fPyRec = fPzRec * TMath::Tan(fThetaY);
1037       fCharge = Int_t(TMath::Sign(1.,muonTrack->GetInverseBendingMomentum()));
1038
1039       //ChiSquare of the track if needed
1040       fChisquare = muonTrack->GetChi2()/(2.0 * muonTrack->GetNHit() - 5);
1041       fEnergy = TMath::Sqrt(fMUONMASS * fMUONMASS + fPxRec * fPxRec + fPyRec * fPyRec + fPzRec * fPzRec);
1042       fEPvector.SetPxPyPzE(fPxRec, fPyRec, fPzRec, fEnergy);
1043
1044       // total number of muons inside a vertex cut
1045       if((TMath::Abs(fZ)<fZVertexCut) && (TMath::Sqrt(fY*fY+fX*fX)<fRhoVertexCut)) {
1046         nMuons++;
1047         if(fEPvector.Pt() > fLowPtCut) {
1048           nMu1GeV++;
1049           if(fEPvector.Pt() > fHighPtCut) {
1050             nMu3GeV++;
1051             if (fEPvector.Pt() > fVeryHighPtCut) {
1052               nMu10GeV++;
1053             }
1054           }
1055         }
1056       }
1057     }//muon track loop
1058
1059     // Fill the event tags
1060     if(ntrack != 0) meanPt = meanPt/ntrack;
1061
1062     evTag->SetEventId(iEventNumber+1);
1063     evTag->SetPath(filepath);
1064
1065     evTag->SetVertexX(vertexIn->GetXv());
1066     evTag->SetVertexY(vertexIn->GetYv());
1067     evTag->SetVertexZ(vertexIn->GetZv());
1068     evTag->SetVertexZError(vertexIn->GetZRes());
1069     evTag->SetVertexFlag(fVertexflag);
1070
1071     evTag->SetT0VertexZ(esd->GetT0zVertex());
1072
1073     evTag->SetTriggerMask(esd->GetTriggerMask());
1074     evTag->SetTriggerCluster(esd->GetTriggerCluster());
1075
1076     evTag->SetZDCNeutron1Energy(esd->GetZDCN1Energy());
1077     evTag->SetZDCProton1Energy(esd->GetZDCP1Energy());
1078     evTag->SetZDCEMEnergy(esd->GetZDCEMEnergy(0),esd->GetZDCEMEnergy(1));
1079     evTag->SetZDCNeutron1Energy(esd->GetZDCN2Energy());
1080     evTag->SetZDCProton1Energy(esd->GetZDCP2Energy());
1081     evTag->SetNumOfParticipants(esd->GetZDCParticipants());
1082
1083
1084     evTag->SetNumOfTracks(esd->GetNumberOfTracks());
1085     evTag->SetNumOfPosTracks(nPos);
1086     evTag->SetNumOfNegTracks(nNeg);
1087     evTag->SetNumOfNeutrTracks(nNeutr);
1088
1089     evTag->SetNumOfV0s(esd->GetNumberOfV0s());
1090     evTag->SetNumOfCascades(esd->GetNumberOfCascades());
1091     evTag->SetNumOfKinks(esd->GetNumberOfKinks());
1092     evTag->SetNumOfPMDTracks(esd->GetNumberOfPmdTracks());
1093
1094     evTag->SetNumOfProtons(nProtons);
1095     evTag->SetNumOfKaons(nKaons);
1096     evTag->SetNumOfPions(nPions);
1097     evTag->SetNumOfMuons(nMuons);
1098     evTag->SetNumOfElectrons(nElectrons);
1099     evTag->SetNumOfPhotons(nGamas);
1100     evTag->SetNumOfPi0s(nPi0s);
1101     evTag->SetNumOfNeutrons(nNeutrons);
1102     evTag->SetNumOfKaon0s(nK0s);
1103
1104     evTag->SetNumOfChargedAbove1GeV(nCh1GeV);
1105     evTag->SetNumOfChargedAbove3GeV(nCh3GeV);
1106     evTag->SetNumOfChargedAbove10GeV(nCh10GeV);
1107     evTag->SetNumOfMuonsAbove1GeV(nMu1GeV);
1108     evTag->SetNumOfMuonsAbove3GeV(nMu3GeV);
1109     evTag->SetNumOfMuonsAbove10GeV(nMu10GeV);
1110     evTag->SetNumOfElectronsAbove1GeV(nEl1GeV);
1111     evTag->SetNumOfElectronsAbove3GeV(nEl3GeV);
1112     evTag->SetNumOfElectronsAbove10GeV(nEl10GeV);
1113
1114     evTag->SetNumOfPHOSClusters(esd->GetNumberOfPHOSClusters());
1115     evTag->SetNumOfEMCALClusters(esd->GetNumberOfEMCALClusters());
1116
1117     evTag->SetTotalMomentum(totalP);
1118     evTag->SetMeanPt(meanPt);
1119     evTag->SetMaxPt(maxPt);
1120
1121     tag->SetRunId(iInitRunNumber);
1122     if(fIsSim) tag->SetDataType(0);
1123     else tag->SetDataType(1);
1124     tag->AddEventTag(*evTag);
1125   }//event loop
1126   lastEvent = iNumberOfEvents;
1127
1128   t->Delete("");
1129
1130   TString localFileName = "Run"; localFileName += tag->GetRunId();
1131   localFileName += ".Event"; localFileName += firstEvent; localFileName += "_"; localFileName += lastEvent; localFileName += "."; localFileName += Counter;
1132   localFileName += ".ESD.tag.root";
1133
1134   TString fileName;
1135
1136   if(fStorage == 0) {
1137     fileName = localFileName.Data();
1138     AliInfo(Form("Writing tags to local file: %s",fileName.Data()));
1139   }
1140   else if(fStorage == 1) {
1141     TString alienLocation = "/alien";
1142     alienLocation += gGrid->Pwd();
1143     alienLocation += fgridpath.Data();
1144     alienLocation += "/";
1145     alienLocation +=  localFileName;
1146     alienLocation += "?se=";
1147     alienLocation += fSE.Data();
1148     fileName = alienLocation.Data();
1149     AliInfo(Form("Writing tags to grid file: %s",fileName.Data()));
1150   }
1151
1152   TFile* ftag = TFile::Open(fileName, "recreate");
1153   ftag->cd();
1154   ttag.Fill();
1155   tag->Clear();
1156   ttag.Write();
1157   ftag->Close();
1158
1159   delete ftag;
1160   delete esd;
1161
1162   delete tag;
1163 }
1164
1165 //_____________________________________________________________________________
1166 void AliESDTagCreator::CreateESDTags(Int_t fFirstEvent, Int_t fLastEvent, TMap *grpData) {
1167   //GRP
1168   Float_t lhcLuminosity = 0.0;
1169   TString lhcState = "test";
1170   UInt_t detectorMask = 0;
1171
1172   TObjString *s = new TObjString;
1173   s = (TObjString *)grpData->GetValue("fDetectorMask");
1174   detectorMask = atoi(s->GetString().Data());
1175
1176   s = (TObjString *)grpData->GetValue("fAliceStartTime");
1177   Float_t startTime = atof(s->GetString().Data());
1178   TTimeStamp *t1 = new TTimeStamp(startTime);
1179
1180   s = (TObjString *)grpData->GetValue("fAliceStopTime");
1181   Float_t stopTime = atof(s->GetString().Data());
1182   TTimeStamp *t2 = new TTimeStamp(stopTime);
1183
1184   s = (TObjString *)grpData->GetValue("fAliceBeamType");
1185   const char* beamtype = s->GetString().Data();
1186
1187   s = (TObjString *)grpData->GetValue("fAliceBeamEnergy");
1188   Float_t beamenergy = atof(s->GetString().Data());
1189
1190
1191   /////////////
1192   //muon code//
1193   ////////////
1194   Double_t fMUONMASS = 0.105658369;
1195   //Variables
1196   Double_t fX,fY,fZ ;
1197   Double_t fThetaX, fThetaY, fPyz, fChisquare;
1198   Double_t fPxRec,fPyRec, fPzRec, fEnergy;
1199   Int_t fCharge;
1200   TLorentzVector fEPvector;
1201
1202   Float_t fZVertexCut = 10.0;
1203   Float_t fRhoVertexCut = 2.0;
1204
1205   Float_t fLowPtCut = 1.0;
1206   Float_t fHighPtCut = 3.0;
1207   Float_t fVeryHighPtCut = 10.0;
1208   ////////////
1209
1210   Double_t partFrac[5] = {0.01, 0.01, 0.85, 0.10, 0.05};
1211
1212   // Creates the tags for all the events in a given ESD file
1213   Int_t ntrack;
1214   Int_t nProtons, nKaons, nPions, nMuons, nElectrons;
1215   Int_t nPos, nNeg, nNeutr;
1216   Int_t nK0s, nNeutrons, nPi0s, nGamas;
1217   Int_t nCh1GeV, nCh3GeV, nCh10GeV;
1218   Int_t nMu1GeV, nMu3GeV, nMu10GeV;
1219   Int_t nEl1GeV, nEl3GeV, nEl10GeV;
1220   Float_t maxPt = .0, meanPt = .0, totalP = .0;
1221   Int_t fVertexflag;
1222   Int_t iRunNumber = 0;
1223   TString fVertexName("default");
1224
1225   AliInfo(Form("Creating the ESD tags......."));
1226
1227   TFile *file = TFile::Open("AliESDs.root");
1228   if (!file || !file->IsOpen()) {
1229     AliError(Form("opening failed"));
1230     delete file;
1231     return ;
1232   }
1233   Int_t lastEvent = 0;
1234   TTree *b = (TTree*) file->Get("esdTree");
1235   AliESDEvent *esd = new AliESDEvent();
1236   esd->ReadFromTree(b);
1237
1238   b->GetEntry(fFirstEvent);
1239   Int_t iInitRunNumber = esd->GetRunNumber();
1240
1241   Int_t iNumberOfEvents = (Int_t)b->GetEntries();
1242   if(fLastEvent == -1) lastEvent = (Int_t)b->GetEntries();
1243   else lastEvent = fLastEvent;
1244
1245   char fileName[256];
1246   sprintf(fileName, "Run%d.Event%d_%d.ESD.tag.root",
1247           iInitRunNumber,fFirstEvent,lastEvent);
1248   AliInfo(Form("writing tags to file %s", fileName));
1249   AliDebug(1, Form("writing tags to file %s", fileName));
1250
1251   TFile* ftag = TFile::Open(fileName, "recreate");
1252
1253   AliRunTag *tag = new AliRunTag();
1254   AliEventTag *evTag = new AliEventTag();
1255   TTree ttag("T","A Tree with event tags");
1256   TBranch * btag = ttag.Branch("AliTAG", &tag);
1257   btag->SetCompressionLevel(9);
1258
1259   if(fLastEvent != -1) iNumberOfEvents = fLastEvent + 1;
1260   for (Int_t iEventNumber = fFirstEvent; iEventNumber < iNumberOfEvents; iEventNumber++) {
1261     ntrack = 0;
1262     nPos = 0;
1263     nNeg = 0;
1264     nNeutr =0;
1265     nK0s = 0;
1266     nNeutrons = 0;
1267     nPi0s = 0;
1268     nGamas = 0;
1269     nProtons = 0;
1270     nKaons = 0;
1271     nPions = 0;
1272     nMuons = 0;
1273     nElectrons = 0;
1274     nCh1GeV = 0;
1275     nCh3GeV = 0;
1276     nCh10GeV = 0;
1277     nMu1GeV = 0;
1278     nMu3GeV = 0;
1279     nMu10GeV = 0;
1280     nEl1GeV = 0;
1281     nEl3GeV = 0;
1282     nEl10GeV = 0;
1283     maxPt = .0;
1284     meanPt = .0;
1285     totalP = .0;
1286     fVertexflag = 0;
1287
1288     b->GetEntry(iEventNumber);
1289     iRunNumber = esd->GetRunNumber();
1290     if(iRunNumber != iInitRunNumber) AliFatal("Inconsistency of run numbers in the AliESD!!!");
1291     const AliESDVertex * vertexIn = esd->GetVertex();
1292     if (!vertexIn) AliError("ESD has not defined vertex.");
1293     if (vertexIn) fVertexName = vertexIn->GetName();
1294     if(fVertexName != "default") fVertexflag = 1;
1295     for (Int_t iTrackNumber = 0; iTrackNumber < esd->GetNumberOfTracks(); iTrackNumber++) {
1296       AliESDtrack * esdTrack = esd->GetTrack(iTrackNumber);
1297       UInt_t status = esdTrack->GetStatus();
1298
1299       //select only tracks with ITS refit
1300       if ((status&AliESDtrack::kITSrefit)==0) continue;
1301       //select only tracks with TPC refit
1302       if ((status&AliESDtrack::kTPCrefit)==0) continue;
1303
1304       //select only tracks with the "combined PID"
1305       if ((status&AliESDtrack::kESDpid)==0) continue;
1306       Double_t p[3];
1307       esdTrack->GetPxPyPz(p);
1308       Double_t momentum = sqrt(pow(p[0],2) + pow(p[1],2) + pow(p[2],2));
1309       Double_t fPt = sqrt(pow(p[0],2) + pow(p[1],2));
1310       totalP += momentum;
1311       meanPt += fPt;
1312       if(fPt > maxPt) maxPt = fPt;
1313
1314       if(esdTrack->GetSign() > 0) {
1315         nPos++;
1316         if(fPt > fLowPtCut) nCh1GeV++;
1317         if(fPt > fHighPtCut) nCh3GeV++;
1318         if(fPt > fVeryHighPtCut) nCh10GeV++;
1319       }
1320       if(esdTrack->GetSign() < 0) {
1321         nNeg++;
1322         if(fPt > fLowPtCut) nCh1GeV++;
1323         if(fPt > fHighPtCut) nCh3GeV++;
1324         if(fPt > fVeryHighPtCut) nCh10GeV++;
1325       }
1326       if(esdTrack->GetSign() == 0) nNeutr++;
1327
1328       //PID
1329       Double_t prob[5];
1330       esdTrack->GetESDpid(prob);
1331
1332       Double_t rcc = 0.0;
1333       for(Int_t i = 0; i < AliPID::kSPECIES; i++) rcc += prob[i]*partFrac[i];
1334       if(rcc == 0.0) continue;
1335       //Bayes' formula
1336       Double_t w[5];
1337       for(Int_t i = 0; i < AliPID::kSPECIES; i++) w[i] = prob[i]*partFrac[i]/rcc;
1338
1339       //protons
1340       if ((w[4]>w[3])&&(w[4]>w[2])&&(w[4]>w[1])&&(w[4]>w[0])) nProtons++;
1341       //kaons
1342       if ((w[3]>w[4])&&(w[3]>w[2])&&(w[3]>w[1])&&(w[3]>w[0])) nKaons++;
1343       //pions
1344       if ((w[2]>w[4])&&(w[2]>w[3])&&(w[2]>w[1])&&(w[2]>w[0])) nPions++;
1345       //electrons
1346       if ((w[0]>w[4])&&(w[0]>w[3])&&(w[0]>w[2])&&(w[0]>w[1])) {
1347         nElectrons++;
1348         if(fPt > fLowPtCut) nEl1GeV++;
1349         if(fPt > fHighPtCut) nEl3GeV++;
1350         if(fPt > fVeryHighPtCut) nEl10GeV++;
1351       }
1352       ntrack++;
1353     }//track loop
1354
1355     /////////////
1356     //muon code//
1357     ////////////
1358     Int_t nMuonTracks = esd->GetNumberOfMuonTracks();
1359     // loop over all reconstructed tracks (also first track of combination)
1360     for (Int_t iTrack = 0; iTrack <  nMuonTracks;  iTrack++) {
1361       AliESDMuonTrack* muonTrack = esd->GetMuonTrack(iTrack);
1362       if (muonTrack == 0x0) continue;
1363
1364       // Coordinates at vertex
1365       fZ = muonTrack->GetZ();
1366       fY = muonTrack->GetBendingCoor();
1367       fX = muonTrack->GetNonBendingCoor();
1368
1369       fThetaX = muonTrack->GetThetaX();
1370       fThetaY = muonTrack->GetThetaY();
1371
1372       fPyz = 1./TMath::Abs(muonTrack->GetInverseBendingMomentum());
1373       fPzRec = - fPyz / TMath::Sqrt(1.0 + TMath::Tan(fThetaY)*TMath::Tan(fThetaY));
1374       fPxRec = fPzRec * TMath::Tan(fThetaX);
1375       fPyRec = fPzRec * TMath::Tan(fThetaY);
1376       fCharge = Int_t(TMath::Sign(1.,muonTrack->GetInverseBendingMomentum()));
1377
1378       //ChiSquare of the track if needed
1379       fChisquare = muonTrack->GetChi2()/(2.0 * muonTrack->GetNHit() - 5);
1380       fEnergy = TMath::Sqrt(fMUONMASS * fMUONMASS + fPxRec * fPxRec + fPyRec * fPyRec + fPzRec * fPzRec);
1381       fEPvector.SetPxPyPzE(fPxRec, fPyRec, fPzRec, fEnergy);
1382
1383       // total number of muons inside a vertex cut
1384       if((TMath::Abs(fZ)<fZVertexCut) && (TMath::Sqrt(fY*fY+fX*fX)<fRhoVertexCut)) {
1385         nMuons++;
1386         if(fEPvector.Pt() > fLowPtCut) {
1387           nMu1GeV++;
1388           if(fEPvector.Pt() > fHighPtCut) {
1389             nMu3GeV++;
1390             if (fEPvector.Pt() > fVeryHighPtCut) {
1391               nMu10GeV++;
1392             }
1393           }
1394         }
1395       }
1396     }//muon track loop
1397
1398     // Fill the event tags
1399     if(ntrack != 0)
1400       meanPt = meanPt/ntrack;
1401
1402     evTag->SetEventId(iEventNumber+1);
1403     if (vertexIn) {
1404       evTag->SetVertexX(vertexIn->GetXv());
1405       evTag->SetVertexY(vertexIn->GetYv());
1406       evTag->SetVertexZ(vertexIn->GetZv());
1407       evTag->SetVertexZError(vertexIn->GetZRes());
1408     }
1409     evTag->SetVertexFlag(fVertexflag);
1410
1411     evTag->SetT0VertexZ(esd->GetT0zVertex());
1412
1413     evTag->SetTriggerMask(esd->GetTriggerMask());
1414     evTag->SetTriggerCluster(esd->GetTriggerCluster());
1415
1416     evTag->SetZDCNeutron1Energy(esd->GetZDCN1Energy());
1417     evTag->SetZDCProton1Energy(esd->GetZDCP1Energy());
1418     evTag->SetZDCNeutron2Energy(esd->GetZDCN2Energy());
1419     evTag->SetZDCProton2Energy(esd->GetZDCP2Energy());
1420     evTag->SetZDCEMEnergy(esd->GetZDCEMEnergy(0),esd->GetZDCEMEnergy(1));
1421     evTag->SetNumOfParticipants(esd->GetZDCParticipants());
1422
1423
1424     evTag->SetNumOfTracks(esd->GetNumberOfTracks());
1425     evTag->SetNumOfPosTracks(nPos);
1426     evTag->SetNumOfNegTracks(nNeg);
1427     evTag->SetNumOfNeutrTracks(nNeutr);
1428
1429     evTag->SetNumOfV0s(esd->GetNumberOfV0s());
1430     evTag->SetNumOfCascades(esd->GetNumberOfCascades());
1431     evTag->SetNumOfKinks(esd->GetNumberOfKinks());
1432     evTag->SetNumOfPMDTracks(esd->GetNumberOfPmdTracks());
1433
1434     evTag->SetNumOfProtons(nProtons);
1435     evTag->SetNumOfKaons(nKaons);
1436     evTag->SetNumOfPions(nPions);
1437     evTag->SetNumOfMuons(nMuons);
1438     evTag->SetNumOfElectrons(nElectrons);
1439     evTag->SetNumOfPhotons(nGamas);
1440     evTag->SetNumOfPi0s(nPi0s);
1441     evTag->SetNumOfNeutrons(nNeutrons);
1442     evTag->SetNumOfKaon0s(nK0s);
1443
1444     evTag->SetNumOfChargedAbove1GeV(nCh1GeV);
1445     evTag->SetNumOfChargedAbove3GeV(nCh3GeV);
1446     evTag->SetNumOfChargedAbove10GeV(nCh10GeV);
1447     evTag->SetNumOfMuonsAbove1GeV(nMu1GeV);
1448     evTag->SetNumOfMuonsAbove3GeV(nMu3GeV);
1449     evTag->SetNumOfMuonsAbove10GeV(nMu10GeV);
1450     evTag->SetNumOfElectronsAbove1GeV(nEl1GeV);
1451     evTag->SetNumOfElectronsAbove3GeV(nEl3GeV);
1452     evTag->SetNumOfElectronsAbove10GeV(nEl10GeV);
1453
1454     evTag->SetNumOfPHOSClusters(esd->GetNumberOfPHOSClusters());
1455     evTag->SetNumOfEMCALClusters(esd->GetNumberOfEMCALClusters());
1456
1457     evTag->SetTotalMomentum(totalP);
1458     evTag->SetMeanPt(meanPt);
1459     evTag->SetMaxPt(maxPt);
1460
1461     tag->SetLHCTag(lhcLuminosity,lhcState);
1462     tag->SetDetectorTag(detectorMask);
1463
1464     tag->SetRunId(iInitRunNumber);
1465     tag->SetRunStartTime(t1->GetDate());
1466     tag->SetRunStopTime(t2->GetDate());
1467     tag->SetBeamEnergy(beamenergy);
1468     tag->SetBeamType(beamtype);
1469
1470     tag->AddEventTag(*evTag);
1471   }
1472
1473   ftag->cd();
1474   ttag.Fill();
1475   tag->Clear();
1476   ttag.Write();
1477   ftag->Close();
1478   file->cd();
1479   delete tag;
1480   delete evTag;
1481 }
1482
1483 //_____________________________________________________________________________
1484 void AliESDTagCreator::SwitchOffBranches() const {
1485   //
1486   // Switch of branches on user request
1487   TObjArray * tokens = fBranches.Tokenize(" ");
1488   Int_t ntok = tokens->GetEntries();
1489   for (Int_t i = 0; i < ntok; i++)  {
1490     TString str = ((TObjString*) tokens->At(i))->GetString();
1491     fChain->SetBranchStatus(Form("%s%s%s","*", str.Data(), "*"), 0);
1492     AliInfo(Form("Branch %s switched off \n", str.Data()));
1493   }
1494 }