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