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;
 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;
 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         if(fEPvector.Pt() > fLowPtCut) {
 411           nMu1GeV++;
 412           if(fEPvector.Pt() > fHighPtCut) {
 413             nMu3GeV++;
 414             if (fEPvector.Pt() > fVeryHighPtCut) {
 415               nMu10GeV++;
 416             }
 417           }
 418         }
 419       }
 420     }//muon track loop
 421
 422     // Fill the event tags
 423     if(ntrack != 0) meanPt = meanPt/ntrack;
 424
 425     //AliInfo(Form("====================================="));
 426     //AliInfo(Form("URL: %s - GUID: %s",fturl.Data(),fguid.Data()));
 427     //AliInfo(Form("====================================="));
 428
 429     evTag->SetEventId(iEventNumber+1);
 430     evTag->SetGUID(fguid);
 431     if(fSession == "grid") {
 432       evTag->SetMD5(0);
 433       evTag->SetTURL(fturl);
 434       evTag->SetSize(0);
 435     }
 436     else evTag->SetPath(fturl);
 437
 438     evTag->SetVertexX(vertexIn->GetXv());
 439     evTag->SetVertexY(vertexIn->GetYv());
 440     evTag->SetVertexZ(vertexIn->GetZv());
 441     evTag->SetVertexZError(vertexIn->GetZRes());
 442     evTag->SetVertexFlag(fVertexflag);
 443
 444     evTag->SetT0VertexZ(esd->GetT0zVertex());
 445
 446     evTag->SetTriggerMask(esd->GetTriggerMask());
 447     evTag->SetTriggerCluster(esd->GetTriggerCluster());
 448
 449     evTag->SetZDCNeutron1Energy(esd->GetZDCN1Energy());
 450     evTag->SetZDCProton1Energy(esd->GetZDCP1Energy());
 451     evTag->SetZDCEMEnergy(esd->GetZDCEMEnergy(0),esd->GetZDCEMEnergy(1));
 452     evTag->SetZDCNeutron1Energy(esd->GetZDCN2Energy());
 453     evTag->SetZDCProton1Energy(esd->GetZDCP2Energy());
 454     evTag->SetNumOfParticipants(esd->GetZDCParticipants());
 455     evTag->SetNumOfParticipants(esd->GetZDCParticipants());
 456     evTag->SetNumOfParticipants2(esd->GetZDCParticipants2());
 457
 458
 459     evTag->SetNumOfTracks(esd->GetNumberOfTracks());
 460     evTag->SetNumOfPosTracks(nPos);
 461     evTag->SetNumOfNegTracks(nNeg);
 462     evTag->SetNumOfNeutrTracks(nNeutr);
 463
 464     evTag->SetNumOfV0s(esd->GetNumberOfV0s());
 465     evTag->SetNumOfCascades(esd->GetNumberOfCascades());
 466     evTag->SetNumOfKinks(esd->GetNumberOfKinks());
 467     evTag->SetNumOfPMDTracks(esd->GetNumberOfPmdTracks());
 468
 469     evTag->SetNumOfProtons(nProtons);
 470     evTag->SetNumOfKaons(nKaons);
 471     evTag->SetNumOfPions(nPions);
 472     evTag->SetNumOfMuons(nMuons);
 473     evTag->SetNumOfElectrons(nElectrons);
 474     evTag->SetNumOfPhotons(nGamas);
 475     evTag->SetNumOfPi0s(nPi0s);
 476     evTag->SetNumOfNeutrons(nNeutrons);
 477     evTag->SetNumOfKaon0s(nK0s);
 478
 479     evTag->SetNumOfChargedAbove1GeV(nCh1GeV);
 480     evTag->SetNumOfChargedAbove3GeV(nCh3GeV);
 481     evTag->SetNumOfChargedAbove10GeV(nCh10GeV);
 482     evTag->SetNumOfMuonsAbove1GeV(nMu1GeV);
 483     evTag->SetNumOfMuonsAbove3GeV(nMu3GeV);
 484     evTag->SetNumOfMuonsAbove10GeV(nMu10GeV);
 485     evTag->SetNumOfElectronsAbove1GeV(nEl1GeV);
 486     evTag->SetNumOfElectronsAbove3GeV(nEl3GeV);
 487     evTag->SetNumOfElectronsAbove10GeV(nEl10GeV);
 488
 489     evTag->SetNumOfPHOSClusters(esd->GetNumberOfPHOSClusters());
 490     evTag->SetNumOfEMCALClusters(esd->GetNumberOfEMCALClusters());
 491
 492     evTag->SetTotalMomentum(totalP);
 493     evTag->SetMeanPt(meanPt);
 494     evTag->SetMaxPt(maxPt);
 495
 496     tag->SetRunId(iInitRunNumber);
 497     if(fIsSim) tag->SetDataType(0);
 498     else tag->SetDataType(1);
 499
 500     if(fguid != fTempGuid) {
 501       fTempGuid = fguid;
 502       ttag.Fill();
 503       tag->Clear("");
 504     }
 505     tag->AddEventTag(*evTag);
 506     if(iEventNumber+1 == chain->GetEntries()) {
 507       //AliInfo(Form("File: %s",fturl.Data()));
 508       ttag.Fill();
 509       tag->Clear("");
 510     }
 511   }//event loop
 512   lastEvent = chain->GetEntries();
 513
 514   //gSystem->GetMemInfo(meminfo);
 515   //AliInfo(Form("After the event and track loop - Memory used: %d MB - Increase: %d MB",meminfo->fMemUsed,meminfo->fMemUsed - tempmem));
 516   //tempmem = meminfo->fMemUsed;
 517
 518   //chain->Delete("");
 519
 520   //gSystem->GetMemInfo(meminfo);
 521   //AliInfo(Form("After the t->Delete - Memory used: %d MB - Increase: %d MB",meminfo->fMemUsed,meminfo->fMemUsed - tempmem));
 522   //tempmem = meminfo->fMemUsed;
 523
 524   ftag->cd();
 525   tag->Clear();
 526   ttag.Write();
 527   ftag->Close();
 528
 529   //gSystem->GetMemInfo(meminfo);
 530   //AliInfo(Form("After the file closing - Memory used: %d MB - Increase: %d MB",meminfo->fMemUsed,meminfo->fMemUsed - tempmem));
 531   //tempmem = meminfo->fMemUsed;
 532
 533   delete esd;
 534   delete tag;
 535
 536   //gSystem->GetMemInfo(meminfo);
 537   //AliInfo(Form("After the delete objects - Memory used: %d MB - Increase: %d MB",meminfo->fMemUsed,meminfo->fMemUsed - tempmem));
 538 }
 539
 540 //_____________________________________________________________________________
 541 void AliESDTagCreator::CreateTag(TFile* file, const char *guid, const char *md5, const char *turl, Long64_t size, Int_t Counter) {
 542   //private method that creates tag files
 543   TString fguid = guid;
 544   TString fmd5 = md5;
 545   TString fturl = turl;
 546   /////////////
 547   //muon code//
 548   ////////////
 549   Double_t fMUONMASS = 0.105658369;
 550   //Variables
 551   Double_t fX,fY,fZ ;
 552   Double_t fThetaX, fThetaY, fPyz, fChisquare;
 553   Double_t fPxRec, fPyRec, fPzRec, fEnergy;
 554   Int_t fCharge;
 555   TLorentzVector fEPvector;
 556
 557   Float_t fZVertexCut = 10.0;
 558   Float_t fRhoVertexCut = 2.0;
 559
 560   Float_t fLowPtCut = 1.0;
 561   Float_t fHighPtCut = 3.0;
 562   Float_t fVeryHighPtCut = 10.0;
 563   ////////////
 564
 565   Double_t partFrac[5] = {0.01, 0.01, 0.85, 0.10, 0.05};
 566
 567   // Creates the tags for all the events in a given ESD file
 568   Bool_t fIsSim = kTRUE;
 569   Int_t ntrack;
 570   Int_t nProtons, nKaons, nPions, nMuons, nElectrons;
 571   Int_t nPos, nNeg, nNeutr;
 572   Int_t nK0s, nNeutrons, nPi0s, nGamas;
 573   Int_t nCh1GeV, nCh3GeV, nCh10GeV;
 574   Int_t nMu1GeV, nMu3GeV, nMu10GeV;
 575   Int_t nEl1GeV, nEl3GeV, nEl10GeV;
 576   Float_t maxPt = .0, meanPt = .0, totalP = .0;
 577   Int_t fVertexflag;
 578   Int_t iRunNumber = 0;
 579   TString fVertexName;
 580
 581   AliRunTag *tag = new AliRunTag();
 582   AliEventTag *evTag = new AliEventTag();
 583   TTree ttag("T","A Tree with event tags");
 584   TBranch * btag = ttag.Branch("AliTAG", &tag);
 585   btag->SetCompressionLevel(9);
 586   gSystem->GetMemInfo(meminfo);
 587   AliInfo(Form("After the tag initialization - Memory used: %d MB",meminfo->fMemUsed));
 588   Int_t tempmem = meminfo->fMemUsed;
 589
 590   AliInfo(Form("Creating the ESD tags......."));
 591
 592   Int_t firstEvent = 0,lastEvent = 0;
 593   TTree *t = (TTree*) file->Get("esdTree");
 594   AliESDEvent *esd = new AliESDEvent();
 595   esd->ReadFromTree(t);
 596
 597   gSystem->GetMemInfo(meminfo);
 598   AliInfo(Form("After the esd initialization - Memory used: %d MB - Increase: %d MB",meminfo->fMemUsed,meminfo->fMemUsed - tempmem));
 599   tempmem = meminfo->fMemUsed;
 600
 601   t->GetEntry(0);
 602   Int_t iInitRunNumber = esd->GetRunNumber();
 603
 604   Int_t iNumberOfEvents = (Int_t)t->GetEntries();
 605   for (Int_t iEventNumber = 0; iEventNumber < iNumberOfEvents; iEventNumber++) {
 606     ntrack = 0;
 607     nPos = 0;
 608     nNeg = 0;
 609     nNeutr =0;
 610     nK0s = 0;
 611     nNeutrons = 0;
 612     nPi0s = 0;
 613     nGamas = 0;
 614     nProtons = 0;
 615     nKaons = 0;
 616     nPions = 0;
 617     nMuons = 0;
 618     nElectrons = 0;
 619     nCh1GeV = 0;
 620     nCh3GeV = 0;
 621     nCh10GeV = 0;
 622     nMu1GeV = 0;
 623     nMu3GeV = 0;
 624     nMu10GeV = 0;
 625     nEl1GeV = 0;
 626     nEl3GeV = 0;
 627     nEl10GeV = 0;
 628     maxPt = .0;
 629     meanPt = .0;
 630     totalP = .0;
 631     fVertexflag = 1;
 632
 633     t->GetEntry(iEventNumber);
 634     iRunNumber = esd->GetRunNumber();
 635     if(iRunNumber != iInitRunNumber) AliFatal("Inconsistency of run numbers in the AliESD!!!");
 636     const AliESDVertex * vertexIn = esd->GetVertex();
 637     fVertexName = vertexIn->GetName();
 638     if(fVertexName == "default") fVertexflag = 0;
 639
 640     for (Int_t iTrackNumber = 0; iTrackNumber < esd->GetNumberOfTracks(); iTrackNumber++) {
 641       AliESDtrack * esdTrack = esd->GetTrack(iTrackNumber);
 642       if(esdTrack->GetLabel() != 0) fIsSim = kTRUE;
 643       else if(esdTrack->GetLabel() == 0) fIsSim = kFALSE;
 644       UInt_t status = esdTrack->GetStatus();
 645
 646       //select only tracks with ITS refit
 647       if ((status&AliESDtrack::kITSrefit)==0) continue;
 648       //select only tracks with TPC refit
 649       if ((status&AliESDtrack::kTPCrefit)==0) continue;
 650
 651       //select only tracks with the "combined PID"
 652       if ((status&AliESDtrack::kESDpid)==0) continue;
 653       Double_t p[3];
 654       esdTrack->GetPxPyPz(p);
 655       Double_t momentum = sqrt(pow(p[0],2) + pow(p[1],2) + pow(p[2],2));
 656       Double_t fPt = sqrt(pow(p[0],2) + pow(p[1],2));
 657       totalP += momentum;
 658       meanPt += fPt;
 659       if(fPt > maxPt) maxPt = fPt;
 660
 661       if(esdTrack->GetSign() > 0) {
 662         nPos++;
 663         if(fPt > fLowPtCut) nCh1GeV++;
 664         if(fPt > fHighPtCut) nCh3GeV++;
 665         if(fPt > fVeryHighPtCut) nCh10GeV++;
 666       }
 667       if(esdTrack->GetSign() < 0) {
 668         nNeg++;
 669         if(fPt > fLowPtCut) nCh1GeV++;
 670         if(fPt > fHighPtCut) nCh3GeV++;
 671         if(fPt > fVeryHighPtCut) nCh10GeV++;
 672       }
 673       if(esdTrack->GetSign() == 0) nNeutr++;
 674
 675       //PID
 676       Double_t prob[5];
 677       esdTrack->GetESDpid(prob);
 678
 679       Double_t rcc = 0.0;
 680       for(Int_t i = 0; i < AliPID::kSPECIES; i++) rcc += prob[i]*partFrac[i];
 681       if(rcc == 0.0) continue;
 682       //Bayes' formula
 683       Double_t w[5];
 684       for(Int_t i = 0; i < AliPID::kSPECIES; i++) w[i] = prob[i]*partFrac[i]/rcc;
 685
 686       //protons
 687       if ((w[4]>w[3])&&(w[4]>w[2])&&(w[4]>w[1])&&(w[4]>w[0])) nProtons++;
 688       //kaons
 689       if ((w[3]>w[4])&&(w[3]>w[2])&&(w[3]>w[1])&&(w[3]>w[0])) nKaons++;
 690       //pions
 691       if ((w[2]>w[4])&&(w[2]>w[3])&&(w[2]>w[1])&&(w[2]>w[0])) nPions++;
 692       //electrons
 693       if ((w[0]>w[4])&&(w[0]>w[3])&&(w[0]>w[2])&&(w[0]>w[1])) {
 694         nElectrons++;
 695         if(fPt > fLowPtCut) nEl1GeV++;
 696         if(fPt > fHighPtCut) nEl3GeV++;
 697         if(fPt > fVeryHighPtCut) nEl10GeV++;
 698       }
 699       ntrack++;
 700     }//esd track loop
 701
 702     /////////////
 703     //muon code//
 704     ////////////
 705     Int_t nMuonTracks = esd->GetNumberOfMuonTracks();
 706     // loop over all reconstructed tracks (also first track of combination)
 707     for (Int_t iTrack = 0; iTrack <  nMuonTracks;  iTrack++) {
 708       AliESDMuonTrack* muonTrack = esd->GetMuonTrack(iTrack);
 709       if (muonTrack == 0x0) continue;
 710
 711       // Coordinates at vertex
 712       fZ = muonTrack->GetZ();
 713       fY = muonTrack->GetBendingCoor();
 714       fX = muonTrack->GetNonBendingCoor();
 715
 716       fThetaX = muonTrack->GetThetaX();
 717       fThetaY = muonTrack->GetThetaY();
 718
 719       fPyz = 1./TMath::Abs(muonTrack->GetInverseBendingMomentum());
 720       fPzRec = - fPyz / TMath::Sqrt(1.0 + TMath::Tan(fThetaY)*TMath::Tan(fThetaY));
 721       fPxRec = fPzRec * TMath::Tan(fThetaX);
 722       fPyRec = fPzRec * TMath::Tan(fThetaY);
 723       fCharge = Int_t(TMath::Sign(1.,muonTrack->GetInverseBendingMomentum()));
 724
 725       //ChiSquare of the track if needed
 726       fChisquare = muonTrack->GetChi2()/(2.0 * muonTrack->GetNHit() - 5);
 727       fEnergy = TMath::Sqrt(fMUONMASS * fMUONMASS + fPxRec * fPxRec + fPyRec * fPyRec + fPzRec * fPzRec);
 728       fEPvector.SetPxPyPzE(fPxRec, fPyRec, fPzRec, fEnergy);
 729
 730       // total number of muons inside a vertex cut
 731       if((TMath::Abs(fZ)<fZVertexCut) && (TMath::Sqrt(fY*fY+fX*fX)<fRhoVertexCut)) {
 732         nMuons++;
 733         if(fEPvector.Pt() > fLowPtCut) {
 734           nMu1GeV++;
 735           if(fEPvector.Pt() > fHighPtCut) {
 736             nMu3GeV++;
 737             if (fEPvector.Pt() > fVeryHighPtCut) {
 738               nMu10GeV++;
 739             }
 740           }
 741         }
 742       }
 743     }//muon track loop
 744
 745     // Fill the event tags
 746     if(ntrack != 0) meanPt = meanPt/ntrack;
 747
 748     evTag->SetEventId(iEventNumber+1);
 749     evTag->SetGUID(fguid);
 750     evTag->SetMD5(fmd5);
 751     evTag->SetTURL(fturl);
 752     evTag->SetSize(size);
 753     evTag->SetVertexX(vertexIn->GetXv());
 754     evTag->SetVertexY(vertexIn->GetYv());
 755     evTag->SetVertexZ(vertexIn->GetZv());
 756     evTag->SetVertexZError(vertexIn->GetZRes());
 757     evTag->SetVertexFlag(fVertexflag);
 758
 759     evTag->SetT0VertexZ(esd->GetT0zVertex());
 760
 761     evTag->SetTriggerMask(esd->GetTriggerMask());
 762     evTag->SetTriggerCluster(esd->GetTriggerCluster());
 763
 764     evTag->SetZDCNeutron1Energy(esd->GetZDCN1Energy());
 765     evTag->SetZDCProton1Energy(esd->GetZDCP1Energy());
 766     evTag->SetZDCEMEnergy(esd->GetZDCEMEnergy(0),esd->GetZDCEMEnergy(1));
 767     evTag->SetZDCNeutron1Energy(esd->GetZDCN2Energy());
 768     evTag->SetZDCProton1Energy(esd->GetZDCP2Energy());
 769     evTag->SetNumOfParticipants(esd->GetZDCParticipants());
 770
 771
 772     evTag->SetNumOfTracks(esd->GetNumberOfTracks());
 773     evTag->SetNumOfPosTracks(nPos);
 774     evTag->SetNumOfNegTracks(nNeg);
 775     evTag->SetNumOfNeutrTracks(nNeutr);
 776
 777     evTag->SetNumOfV0s(esd->GetNumberOfV0s());
 778     evTag->SetNumOfCascades(esd->GetNumberOfCascades());
 779     evTag->SetNumOfKinks(esd->GetNumberOfKinks());
 780     evTag->SetNumOfPMDTracks(esd->GetNumberOfPmdTracks());
 781
 782     evTag->SetNumOfProtons(nProtons);
 783     evTag->SetNumOfKaons(nKaons);
 784     evTag->SetNumOfPions(nPions);
 785     evTag->SetNumOfMuons(nMuons);
 786     evTag->SetNumOfElectrons(nElectrons);
 787     evTag->SetNumOfPhotons(nGamas);
 788     evTag->SetNumOfPi0s(nPi0s);
 789     evTag->SetNumOfNeutrons(nNeutrons);
 790     evTag->SetNumOfKaon0s(nK0s);
 791
 792     evTag->SetNumOfChargedAbove1GeV(nCh1GeV);
 793     evTag->SetNumOfChargedAbove3GeV(nCh3GeV);
 794     evTag->SetNumOfChargedAbove10GeV(nCh10GeV);
 795     evTag->SetNumOfMuonsAbove1GeV(nMu1GeV);
 796     evTag->SetNumOfMuonsAbove3GeV(nMu3GeV);
 797     evTag->SetNumOfMuonsAbove10GeV(nMu10GeV);
 798     evTag->SetNumOfElectronsAbove1GeV(nEl1GeV);
 799     evTag->SetNumOfElectronsAbove3GeV(nEl3GeV);
 800     evTag->SetNumOfElectronsAbove10GeV(nEl10GeV);
 801
 802     evTag->SetNumOfPHOSClusters(esd->GetNumberOfPHOSClusters());
 803     evTag->SetNumOfEMCALClusters(esd->GetNumberOfEMCALClusters());
 804
 805     evTag->SetTotalMomentum(totalP);
 806     evTag->SetMeanPt(meanPt);
 807     evTag->SetMaxPt(maxPt);
 808
 809     tag->SetRunId(iInitRunNumber);
 810     if(fIsSim) tag->SetDataType(0);
 811     else tag->SetDataType(1);
 812     tag->AddEventTag(*evTag);
 813   }//event loop
 814   lastEvent = iNumberOfEvents;
 815
 816   gSystem->GetMemInfo(meminfo);
 817   AliInfo(Form("After the event and track loop - Memory used: %d MB - Increase: %d MB",meminfo->fMemUsed,meminfo->fMemUsed - tempmem));
 818   tempmem = meminfo->fMemUsed;
 819   t->Delete("");
 820
 821   gSystem->GetMemInfo(meminfo);
 822   AliInfo(Form("After the t->Delete - Memory used: %d MB - Increase: %d MB",meminfo->fMemUsed,meminfo->fMemUsed - tempmem));
 823   tempmem = meminfo->fMemUsed;
 824
 825   TString localFileName = "Run"; localFileName += tag->GetRunId();
 826   localFileName += ".Event"; localFileName += firstEvent; localFileName += "_"; localFileName += lastEvent; localFileName += "."; localFileName += Counter;
 827   localFileName += ".ESD.tag.root";
 828
 829   TString fileName;
 830
 831   if(fStorage == 0) {
 832     fileName = localFileName.Data();
 833     AliInfo(Form("Writing tags to local file: %s",fileName.Data()));
 834   }
 835   else if(fStorage == 1) {
 836     TString alienLocation = "/alien";
 837     alienLocation += gGrid->Pwd();
 838     alienLocation += fgridpath.Data();
 839     alienLocation += "/";
 840     alienLocation +=  localFileName;
 841     alienLocation += "?se=";
 842     alienLocation += fSE.Data();
 843     fileName = alienLocation.Data();
 844     AliInfo(Form("Writing tags to grid file: %s",fileName.Data()));
 845   }
 846
 847   TFile* ftag = TFile::Open(fileName, "recreate");
 848   ftag->cd();
 849   ttag.Fill();
 850   tag->Clear();
 851   ttag.Write();
 852   ftag->Close();
 853
 854   gSystem->GetMemInfo(meminfo);
 855   AliInfo(Form("After the file closing - Memory used: %d MB - Increase: %d MB",meminfo->fMemUsed,meminfo->fMemUsed - tempmem));
 856   tempmem = meminfo->fMemUsed;
 857
 858   delete ftag;
 859   delete esd;
 860
 861   delete tag;
 862   gSystem->GetMemInfo(meminfo);
 863   AliInfo(Form("After the delete objects - Memory used: %d MB - Increase: %d MB",meminfo->fMemUsed,meminfo->fMemUsed - tempmem));
 864 }
 865
 866 //_____________________________________________________________________________
 867 void AliESDTagCreator::CreateTag(TFile* file, const char *filepath, Int_t Counter) {
 868   //private method that creates tag files
 869   /////////////
 870   //muon code//
 871   ////////////
 872   Double_t fMUONMASS = 0.105658369;
 873   //Variables
 874   Double_t fX,fY,fZ ;
 875   Double_t fThetaX, fThetaY, fPyz, fChisquare;
 876   Double_t fPxRec, fPyRec, fPzRec, fEnergy;
 877   Int_t fCharge;
 878   TLorentzVector fEPvector;
 879
 880   Float_t fZVertexCut = 10.0;
 881   Float_t fRhoVertexCut = 2.0;
 882
 883   Float_t fLowPtCut = 1.0;
 884   Float_t fHighPtCut = 3.0;
 885   Float_t fVeryHighPtCut = 10.0;
 886   ////////////
 887
 888   Double_t partFrac[5] = {0.01, 0.01, 0.85, 0.10, 0.05};
 889
 890   // Creates the tags for all the events in a given ESD file
 891   Bool_t fIsSim = kTRUE;
 892   Int_t ntrack;
 893   Int_t nProtons, nKaons, nPions, nMuons, nElectrons;
 894   Int_t nPos, nNeg, nNeutr;
 895   Int_t nK0s, nNeutrons, nPi0s, nGamas;
 896   Int_t nCh1GeV, nCh3GeV, nCh10GeV;
 897   Int_t nMu1GeV, nMu3GeV, nMu10GeV;
 898   Int_t nEl1GeV, nEl3GeV, nEl10GeV;
 899   Float_t maxPt = .0, meanPt = .0, totalP = .0;
 900   Int_t fVertexflag;
 901   Int_t iRunNumber = 0;
 902   TString fVertexName;
 903
 904   AliRunTag *tag = new AliRunTag();
 905   AliEventTag *evTag = new AliEventTag();
 906   TTree ttag("T","A Tree with event tags");
 907   TBranch * btag = ttag.Branch("AliTAG", &tag);
 908   btag->SetCompressionLevel(9);
 909
 910   AliInfo(Form("Creating the ESD tags......."));
 911
 912   Int_t firstEvent = 0,lastEvent = 0;
 913
 914   TTree *t = (TTree*) file->Get("esdTree");
 915   AliESDEvent *esd = new AliESDEvent();
 916   esd->ReadFromTree(t);
 917
 918   t->GetEntry(0);
 919   Int_t iInitRunNumber = esd->GetRunNumber();
 920
 921   Int_t iNumberOfEvents = (Int_t)t->GetEntries();
 922   for (Int_t iEventNumber = 0; iEventNumber < iNumberOfEvents; iEventNumber++) {
 923     ntrack = 0;
 924     nPos = 0;
 925     nNeg = 0;
 926     nNeutr =0;
 927     nK0s = 0;
 928     nNeutrons = 0;
 929     nPi0s = 0;
 930     nGamas = 0;
 931     nProtons = 0;
 932     nKaons = 0;
 933     nPions = 0;
 934     nMuons = 0;
 935     nElectrons = 0;
 936     nCh1GeV = 0;
 937     nCh3GeV = 0;
 938     nCh10GeV = 0;
 939     nMu1GeV = 0;
 940     nMu3GeV = 0;
 941     nMu10GeV = 0;
 942     nEl1GeV = 0;
 943     nEl3GeV = 0;
 944     nEl10GeV = 0;
 945     maxPt = .0;
 946     meanPt = .0;
 947     totalP = .0;
 948     fVertexflag = 1;
 949
 950     t->GetEntry(iEventNumber);
 951     iRunNumber = esd->GetRunNumber();
 952     if(iRunNumber != iInitRunNumber) AliFatal("Inconsistency of run numbers in the AliESD!!!");
 953     const AliESDVertex * vertexIn = esd->GetVertex();
 954     fVertexName = vertexIn->GetName();
 955     if(fVertexName == "default") fVertexflag = 0;
 956
 957     for (Int_t iTrackNumber = 0; iTrackNumber < esd->GetNumberOfTracks(); iTrackNumber++) {
 958       AliESDtrack * esdTrack = esd->GetTrack(iTrackNumber);
 959       if(esdTrack->GetLabel() != 0) fIsSim = kTRUE;
 960       else if(esdTrack->GetLabel() == 0) fIsSim = kFALSE;
 961       UInt_t status = esdTrack->GetStatus();
 962
 963       //select only tracks with ITS refit
 964       if ((status&AliESDtrack::kITSrefit)==0) continue;
 965       //select only tracks with TPC refit
 966       if ((status&AliESDtrack::kTPCrefit)==0) continue;
 967
 968       //select only tracks with the "combined PID"
 969       if ((status&AliESDtrack::kESDpid)==0) continue;
 970       Double_t p[3];
 971       esdTrack->GetPxPyPz(p);
 972       Double_t momentum = sqrt(pow(p[0],2) + pow(p[1],2) + pow(p[2],2));
 973       Double_t fPt = sqrt(pow(p[0],2) + pow(p[1],2));
 974       totalP += momentum;
 975       meanPt += fPt;
 976       if(fPt > maxPt) maxPt = fPt;
 977
 978       if(esdTrack->GetSign() > 0) {
 979         nPos++;
 980         if(fPt > fLowPtCut) nCh1GeV++;
 981         if(fPt > fHighPtCut) nCh3GeV++;
 982         if(fPt > fVeryHighPtCut) nCh10GeV++;
 983       }
 984       if(esdTrack->GetSign() < 0) {
 985         nNeg++;
 986         if(fPt > fLowPtCut) nCh1GeV++;
 987         if(fPt > fHighPtCut) nCh3GeV++;
 988         if(fPt > fVeryHighPtCut) nCh10GeV++;
 989       }
 990       if(esdTrack->GetSign() == 0) nNeutr++;
 991
 992       //PID
 993       Double_t prob[5];
 994       esdTrack->GetESDpid(prob);
 995
 996       Double_t rcc = 0.0;
 997       for(Int_t i = 0; i < AliPID::kSPECIES; i++) rcc += prob[i]*partFrac[i];
 998       if(rcc == 0.0) continue;
 999       //Bayes' formula
1000       Double_t w[5];
1001       for(Int_t i = 0; i < AliPID::kSPECIES; i++) w[i] = prob[i]*partFrac[i]/rcc;
1002
1003       //protons
1004       if ((w[4]>w[3])&&(w[4]>w[2])&&(w[4]>w[1])&&(w[4]>w[0])) nProtons++;
1005       //kaons
1006       if ((w[3]>w[4])&&(w[3]>w[2])&&(w[3]>w[1])&&(w[3]>w[0])) nKaons++;
1007       //pions
1008       if ((w[2]>w[4])&&(w[2]>w[3])&&(w[2]>w[1])&&(w[2]>w[0])) nPions++;
1009       //electrons
1010       if ((w[0]>w[4])&&(w[0]>w[3])&&(w[0]>w[2])&&(w[0]>w[1])) {
1011         nElectrons++;
1012         if(fPt > fLowPtCut) nEl1GeV++;
1013         if(fPt > fHighPtCut) nEl3GeV++;
1014         if(fPt > fVeryHighPtCut) nEl10GeV++;
1015       }
1016       ntrack++;
1017     }//esd track loop
1018
1019     /////////////
1020     //muon code//
1021     ////////////
1022     Int_t nMuonTracks = esd->GetNumberOfMuonTracks();
1023     // loop over all reconstructed tracks (also first track of combination)
1024     for (Int_t iTrack = 0; iTrack <  nMuonTracks;  iTrack++) {
1025       AliESDMuonTrack* muonTrack = esd->GetMuonTrack(iTrack);
1026       if (muonTrack == 0x0) continue;
1027
1028       // Coordinates at vertex
1029       fZ = muonTrack->GetZ();
1030       fY = muonTrack->GetBendingCoor();
1031       fX = muonTrack->GetNonBendingCoor();
1032
1033       fThetaX = muonTrack->GetThetaX();
1034       fThetaY = muonTrack->GetThetaY();
1035
1036       fPyz = 1./TMath::Abs(muonTrack->GetInverseBendingMomentum());
1037       fPzRec = - fPyz / TMath::Sqrt(1.0 + TMath::Tan(fThetaY)*TMath::Tan(fThetaY));
1038       fPxRec = fPzRec * TMath::Tan(fThetaX);
1039       fPyRec = fPzRec * TMath::Tan(fThetaY);
1040       fCharge = Int_t(TMath::Sign(1.,muonTrack->GetInverseBendingMomentum()));
1041
1042       //ChiSquare of the track if needed
1043       fChisquare = muonTrack->GetChi2()/(2.0 * muonTrack->GetNHit() - 5);
1044       fEnergy = TMath::Sqrt(fMUONMASS * fMUONMASS + fPxRec * fPxRec + fPyRec * fPyRec + fPzRec * fPzRec);
1045       fEPvector.SetPxPyPzE(fPxRec, fPyRec, fPzRec, fEnergy);
1046
1047       // total number of muons inside a vertex cut
1048       if((TMath::Abs(fZ)<fZVertexCut) && (TMath::Sqrt(fY*fY+fX*fX)<fRhoVertexCut)) {
1049         nMuons++;
1050         if(fEPvector.Pt() > fLowPtCut) {
1051           nMu1GeV++;
1052           if(fEPvector.Pt() > fHighPtCut) {
1053             nMu3GeV++;
1054             if (fEPvector.Pt() > fVeryHighPtCut) {
1055               nMu10GeV++;
1056             }
1057           }
1058         }
1059       }
1060     }//muon track loop
1061
1062     // Fill the event tags
1063     if(ntrack != 0) meanPt = meanPt/ntrack;
1064
1065     evTag->SetEventId(iEventNumber+1);
1066     evTag->SetPath(filepath);
1067
1068     evTag->SetVertexX(vertexIn->GetXv());
1069     evTag->SetVertexY(vertexIn->GetYv());
1070     evTag->SetVertexZ(vertexIn->GetZv());
1071     evTag->SetVertexZError(vertexIn->GetZRes());
1072     evTag->SetVertexFlag(fVertexflag);
1073
1074     evTag->SetT0VertexZ(esd->GetT0zVertex());
1075
1076     evTag->SetTriggerMask(esd->GetTriggerMask());
1077     evTag->SetTriggerCluster(esd->GetTriggerCluster());
1078
1079     evTag->SetZDCNeutron1Energy(esd->GetZDCN1Energy());
1080     evTag->SetZDCProton1Energy(esd->GetZDCP1Energy());
1081     evTag->SetZDCEMEnergy(esd->GetZDCEMEnergy(0),esd->GetZDCEMEnergy(1));
1082     evTag->SetZDCNeutron1Energy(esd->GetZDCN2Energy());
1083     evTag->SetZDCProton1Energy(esd->GetZDCP2Energy());
1084     evTag->SetNumOfParticipants(esd->GetZDCParticipants());
1085
1086
1087     evTag->SetNumOfTracks(esd->GetNumberOfTracks());
1088     evTag->SetNumOfPosTracks(nPos);
1089     evTag->SetNumOfNegTracks(nNeg);
1090     evTag->SetNumOfNeutrTracks(nNeutr);
1091
1092     evTag->SetNumOfV0s(esd->GetNumberOfV0s());
1093     evTag->SetNumOfCascades(esd->GetNumberOfCascades());
1094     evTag->SetNumOfKinks(esd->GetNumberOfKinks());
1095     evTag->SetNumOfPMDTracks(esd->GetNumberOfPmdTracks());
1096
1097     evTag->SetNumOfProtons(nProtons);
1098     evTag->SetNumOfKaons(nKaons);
1099     evTag->SetNumOfPions(nPions);
1100     evTag->SetNumOfMuons(nMuons);
1101     evTag->SetNumOfElectrons(nElectrons);
1102     evTag->SetNumOfPhotons(nGamas);
1103     evTag->SetNumOfPi0s(nPi0s);
1104     evTag->SetNumOfNeutrons(nNeutrons);
1105     evTag->SetNumOfKaon0s(nK0s);
1106
1107     evTag->SetNumOfChargedAbove1GeV(nCh1GeV);
1108     evTag->SetNumOfChargedAbove3GeV(nCh3GeV);
1109     evTag->SetNumOfChargedAbove10GeV(nCh10GeV);
1110     evTag->SetNumOfMuonsAbove1GeV(nMu1GeV);
1111     evTag->SetNumOfMuonsAbove3GeV(nMu3GeV);
1112     evTag->SetNumOfMuonsAbove10GeV(nMu10GeV);
1113     evTag->SetNumOfElectronsAbove1GeV(nEl1GeV);
1114     evTag->SetNumOfElectronsAbove3GeV(nEl3GeV);
1115     evTag->SetNumOfElectronsAbove10GeV(nEl10GeV);
1116
1117     evTag->SetNumOfPHOSClusters(esd->GetNumberOfPHOSClusters());
1118     evTag->SetNumOfEMCALClusters(esd->GetNumberOfEMCALClusters());
1119
1120     evTag->SetTotalMomentum(totalP);
1121     evTag->SetMeanPt(meanPt);
1122     evTag->SetMaxPt(maxPt);
1123
1124     tag->SetRunId(iInitRunNumber);
1125     if(fIsSim) tag->SetDataType(0);
1126     else tag->SetDataType(1);
1127     tag->AddEventTag(*evTag);
1128   }//event loop
1129   lastEvent = iNumberOfEvents;
1130
1131   t->Delete("");
1132
1133   TString localFileName = "Run"; localFileName += tag->GetRunId();
1134   localFileName += ".Event"; localFileName += firstEvent; localFileName += "_"; localFileName += lastEvent; localFileName += "."; localFileName += Counter;
1135   localFileName += ".ESD.tag.root";
1136
1137   TString fileName;
1138
1139   if(fStorage == 0) {
1140     fileName = localFileName.Data();
1141     AliInfo(Form("Writing tags to local file: %s",fileName.Data()));
1142   }
1143   else if(fStorage == 1) {
1144     TString alienLocation = "/alien";
1145     alienLocation += gGrid->Pwd();
1146     alienLocation += fgridpath.Data();
1147     alienLocation += "/";
1148     alienLocation +=  localFileName;
1149     alienLocation += "?se=";
1150     alienLocation += fSE.Data();
1151     fileName = alienLocation.Data();
1152     AliInfo(Form("Writing tags to grid file: %s",fileName.Data()));
1153   }
1154
1155   TFile* ftag = TFile::Open(fileName, "recreate");
1156   ftag->cd();
1157   ttag.Fill();
1158   tag->Clear();
1159   ttag.Write();
1160   ftag->Close();
1161
1162   delete ftag;
1163   delete esd;
1164
1165   delete tag;
1166 }
1167
1168 //_____________________________________________________________________________
1169 void AliESDTagCreator::CreateESDTags(Int_t fFirstEvent, Int_t fLastEvent, AliGRPObject *grpData) {
1170   //GRP
1171   Float_t lhcLuminosity = 0.0;
1172   TString lhcState = "test";
1173   //UInt_t detectorMask = 0;
1174   Int_t detectorMask = 0;
1175
1176   detectorMask = grpData->GetDetectorMask();
1177   time_t startTime = grpData->GetTimeStart();
1178   TTimeStamp *t1 = new TTimeStamp(startTime);
1179   time_t endTime = grpData->GetTimeEnd();
1180   TTimeStamp *t2 = new TTimeStamp(endTime);
1181   const char* beamtype = grpData->GetBeamType();
1182   Float_t beamenergy = grpData->GetBeamEnergy();
1183
1184
1185   /////////////
1186   //muon code//
1187   ////////////
1188   Double_t fMUONMASS = 0.105658369;
1189   //Variables
1190   Double_t fX,fY,fZ ;
1191   Double_t fThetaX, fThetaY, fPyz, fChisquare;
1192   Double_t fPxRec,fPyRec, fPzRec, fEnergy;
1193   Int_t fCharge;
1194   TLorentzVector fEPvector;
1195
1196   Float_t fZVertexCut = 10.0;
1197   Float_t fRhoVertexCut = 2.0;
1198
1199   Float_t fLowPtCut = 1.0;
1200   Float_t fHighPtCut = 3.0;
1201   Float_t fVeryHighPtCut = 10.0;
1202   ////////////
1203
1204   Double_t partFrac[5] = {0.01, 0.01, 0.85, 0.10, 0.05};
1205
1206   // Creates the tags for all the events in a given ESD file
1207   Int_t ntrack;
1208   Int_t nProtons, nKaons, nPions, nMuons, nElectrons;
1209   Int_t nPos, nNeg, nNeutr;
1210   Int_t nK0s, nNeutrons, nPi0s, nGamas;
1211   Int_t nCh1GeV, nCh3GeV, nCh10GeV;
1212   Int_t nMu1GeV, nMu3GeV, nMu10GeV;
1213   Int_t nEl1GeV, nEl3GeV, nEl10GeV;
1214   Float_t maxPt = .0, meanPt = .0, totalP = .0;
1215   Int_t fVertexflag;
1216   Int_t iRunNumber = 0;
1217   TString fVertexName("default");
1218
1219   AliInfo(Form("Creating the ESD tags......."));
1220
1221   TFile *file = TFile::Open("AliESDs.root");
1222   if (!file || !file->IsOpen()) {
1223     AliError(Form("opening failed"));
1224     delete file;
1225     return ;
1226   }
1227   Int_t lastEvent = 0;
1228   TTree *b = (TTree*) file->Get("esdTree");
1229   AliESDEvent *esd = new AliESDEvent();
1230   esd->ReadFromTree(b);
1231
1232   b->GetEntry(fFirstEvent);
1233   Int_t iInitRunNumber = esd->GetRunNumber();
1234
1235   Int_t iNumberOfEvents = (Int_t)b->GetEntries();
1236   if(fLastEvent == -1) lastEvent = (Int_t)b->GetEntries();
1237   else lastEvent = fLastEvent;
1238
1239   char fileName[256];
1240   sprintf(fileName, "Run%d.Event%d_%d.ESD.tag.root",
1241           iInitRunNumber,fFirstEvent,lastEvent);
1242   AliInfo(Form("writing tags to file %s", fileName));
1243   AliDebug(1, Form("writing tags to file %s", fileName));
1244
1245   TFile* ftag = TFile::Open(fileName, "recreate");
1246
1247   AliRunTag *tag = new AliRunTag();
1248   AliEventTag *evTag = new AliEventTag();
1249   TTree ttag("T","A Tree with event tags");
1250   TBranch * btag = ttag.Branch("AliTAG", &tag);
1251   btag->SetCompressionLevel(9);
1252
1253   if(fLastEvent != -1) iNumberOfEvents = fLastEvent + 1;
1254   for (Int_t iEventNumber = fFirstEvent; iEventNumber < iNumberOfEvents; iEventNumber++) {
1255     ntrack = 0;
1256     nPos = 0;
1257     nNeg = 0;
1258     nNeutr =0;
1259     nK0s = 0;
1260     nNeutrons = 0;
1261     nPi0s = 0;
1262     nGamas = 0;
1263     nProtons = 0;
1264     nKaons = 0;
1265     nPions = 0;
1266     nMuons = 0;
1267     nElectrons = 0;
1268     nCh1GeV = 0;
1269     nCh3GeV = 0;
1270     nCh10GeV = 0;
1271     nMu1GeV = 0;
1272     nMu3GeV = 0;
1273     nMu10GeV = 0;
1274     nEl1GeV = 0;
1275     nEl3GeV = 0;
1276     nEl10GeV = 0;
1277     maxPt = .0;
1278     meanPt = .0;
1279     totalP = .0;
1280     fVertexflag = 0;
1281
1282     b->GetEntry(iEventNumber);
1283     iRunNumber = esd->GetRunNumber();
1284     if(iRunNumber != iInitRunNumber) AliFatal("Inconsistency of run numbers in the AliESD!!!");
1285     const AliESDVertex * vertexIn = esd->GetVertex();
1286     if (!vertexIn) AliError("ESD has not defined vertex.");
1287     if (vertexIn) fVertexName = vertexIn->GetName();
1288     if(fVertexName != "default") fVertexflag = 1;
1289     for (Int_t iTrackNumber = 0; iTrackNumber < esd->GetNumberOfTracks(); iTrackNumber++) {
1290       AliESDtrack * esdTrack = esd->GetTrack(iTrackNumber);
1291       UInt_t status = esdTrack->GetStatus();
1292
1293       //select only tracks with ITS refit
1294       if ((status&AliESDtrack::kITSrefit)==0) continue;
1295       //select only tracks with TPC refit
1296       if ((status&AliESDtrack::kTPCrefit)==0) continue;
1297
1298       //select only tracks with the "combined PID"
1299       if ((status&AliESDtrack::kESDpid)==0) continue;
1300       Double_t p[3];
1301       esdTrack->GetPxPyPz(p);
1302       Double_t momentum = sqrt(pow(p[0],2) + pow(p[1],2) + pow(p[2],2));
1303       Double_t fPt = sqrt(pow(p[0],2) + pow(p[1],2));
1304       totalP += momentum;
1305       meanPt += fPt;
1306       if(fPt > maxPt) maxPt = fPt;
1307
1308       if(esdTrack->GetSign() > 0) {
1309         nPos++;
1310         if(fPt > fLowPtCut) nCh1GeV++;
1311         if(fPt > fHighPtCut) nCh3GeV++;
1312         if(fPt > fVeryHighPtCut) nCh10GeV++;
1313       }
1314       if(esdTrack->GetSign() < 0) {
1315         nNeg++;
1316         if(fPt > fLowPtCut) nCh1GeV++;
1317         if(fPt > fHighPtCut) nCh3GeV++;
1318         if(fPt > fVeryHighPtCut) nCh10GeV++;
1319       }
1320       if(esdTrack->GetSign() == 0) nNeutr++;
1321
1322       //PID
1323       Double_t prob[5];
1324       esdTrack->GetESDpid(prob);
1325
1326       Double_t rcc = 0.0;
1327       for(Int_t i = 0; i < AliPID::kSPECIES; i++) rcc += prob[i]*partFrac[i];
1328       if(rcc == 0.0) continue;
1329       //Bayes' formula
1330       Double_t w[5];
1331       for(Int_t i = 0; i < AliPID::kSPECIES; i++) w[i] = prob[i]*partFrac[i]/rcc;
1332
1333       //protons
1334       if ((w[4]>w[3])&&(w[4]>w[2])&&(w[4]>w[1])&&(w[4]>w[0])) nProtons++;
1335       //kaons
1336       if ((w[3]>w[4])&&(w[3]>w[2])&&(w[3]>w[1])&&(w[3]>w[0])) nKaons++;
1337       //pions
1338       if ((w[2]>w[4])&&(w[2]>w[3])&&(w[2]>w[1])&&(w[2]>w[0])) nPions++;
1339       //electrons
1340       if ((w[0]>w[4])&&(w[0]>w[3])&&(w[0]>w[2])&&(w[0]>w[1])) {
1341         nElectrons++;
1342         if(fPt > fLowPtCut) nEl1GeV++;
1343         if(fPt > fHighPtCut) nEl3GeV++;
1344         if(fPt > fVeryHighPtCut) nEl10GeV++;
1345       }
1346       ntrack++;
1347     }//track loop
1348
1349     /////////////
1350     //muon code//
1351     ////////////
1352     Int_t nMuonTracks = esd->GetNumberOfMuonTracks();
1353     // loop over all reconstructed tracks (also first track of combination)
1354     for (Int_t iTrack = 0; iTrack <  nMuonTracks;  iTrack++) {
1355       AliESDMuonTrack* muonTrack = esd->GetMuonTrack(iTrack);
1356       if (muonTrack == 0x0) continue;
1357
1358       // Coordinates at vertex
1359       fZ = muonTrack->GetZ();
1360       fY = muonTrack->GetBendingCoor();
1361       fX = muonTrack->GetNonBendingCoor();
1362
1363       fThetaX = muonTrack->GetThetaX();
1364       fThetaY = muonTrack->GetThetaY();
1365
1366       fPyz = 1./TMath::Abs(muonTrack->GetInverseBendingMomentum());
1367       fPzRec = - fPyz / TMath::Sqrt(1.0 + TMath::Tan(fThetaY)*TMath::Tan(fThetaY));
1368       fPxRec = fPzRec * TMath::Tan(fThetaX);
1369       fPyRec = fPzRec * TMath::Tan(fThetaY);
1370       fCharge = Int_t(TMath::Sign(1.,muonTrack->GetInverseBendingMomentum()));
1371
1372       //ChiSquare of the track if needed
1373       fChisquare = muonTrack->GetChi2()/(2.0 * muonTrack->GetNHit() - 5);
1374       fEnergy = TMath::Sqrt(fMUONMASS * fMUONMASS + fPxRec * fPxRec + fPyRec * fPyRec + fPzRec * fPzRec);
1375       fEPvector.SetPxPyPzE(fPxRec, fPyRec, fPzRec, fEnergy);
1376
1377       // total number of muons inside a vertex cut
1378       if((TMath::Abs(fZ)<fZVertexCut) && (TMath::Sqrt(fY*fY+fX*fX)<fRhoVertexCut)) {
1379         nMuons++;
1380         if(fEPvector.Pt() > fLowPtCut) {
1381           nMu1GeV++;
1382           if(fEPvector.Pt() > fHighPtCut) {
1383             nMu3GeV++;
1384             if (fEPvector.Pt() > fVeryHighPtCut) {
1385               nMu10GeV++;
1386             }
1387           }
1388         }
1389       }
1390     }//muon track loop
1391
1392     // Fill the event tags
1393     if(ntrack != 0)
1394       meanPt = meanPt/ntrack;
1395
1396     evTag->SetEventId(iEventNumber+1);
1397     if (vertexIn) {
1398       evTag->SetVertexX(vertexIn->GetXv());
1399       evTag->SetVertexY(vertexIn->GetYv());
1400       evTag->SetVertexZ(vertexIn->GetZv());
1401       evTag->SetVertexZError(vertexIn->GetZRes());
1402     }
1403     evTag->SetVertexFlag(fVertexflag);
1404
1405     evTag->SetT0VertexZ(esd->GetT0zVertex());
1406
1407     evTag->SetTriggerMask(esd->GetTriggerMask());
1408     evTag->SetTriggerCluster(esd->GetTriggerCluster());
1409
1410     evTag->SetZDCNeutron1Energy(esd->GetZDCN1Energy());
1411     evTag->SetZDCProton1Energy(esd->GetZDCP1Energy());
1412     evTag->SetZDCNeutron2Energy(esd->GetZDCN2Energy());
1413     evTag->SetZDCProton2Energy(esd->GetZDCP2Energy());
1414     evTag->SetZDCEMEnergy(esd->GetZDCEMEnergy(0),esd->GetZDCEMEnergy(1));
1415     evTag->SetNumOfParticipants(esd->GetZDCParticipants());
1416
1417
1418     evTag->SetNumOfTracks(esd->GetNumberOfTracks());
1419     evTag->SetNumOfPosTracks(nPos);
1420     evTag->SetNumOfNegTracks(nNeg);
1421     evTag->SetNumOfNeutrTracks(nNeutr);
1422
1423     evTag->SetNumOfV0s(esd->GetNumberOfV0s());
1424     evTag->SetNumOfCascades(esd->GetNumberOfCascades());
1425     evTag->SetNumOfKinks(esd->GetNumberOfKinks());
1426     evTag->SetNumOfPMDTracks(esd->GetNumberOfPmdTracks());
1427
1428     evTag->SetNumOfProtons(nProtons);
1429     evTag->SetNumOfKaons(nKaons);
1430     evTag->SetNumOfPions(nPions);
1431     evTag->SetNumOfMuons(nMuons);
1432     evTag->SetNumOfElectrons(nElectrons);
1433     evTag->SetNumOfPhotons(nGamas);
1434     evTag->SetNumOfPi0s(nPi0s);
1435     evTag->SetNumOfNeutrons(nNeutrons);
1436     evTag->SetNumOfKaon0s(nK0s);
1437
1438     evTag->SetNumOfChargedAbove1GeV(nCh1GeV);
1439     evTag->SetNumOfChargedAbove3GeV(nCh3GeV);
1440     evTag->SetNumOfChargedAbove10GeV(nCh10GeV);
1441     evTag->SetNumOfMuonsAbove1GeV(nMu1GeV);
1442     evTag->SetNumOfMuonsAbove3GeV(nMu3GeV);
1443     evTag->SetNumOfMuonsAbove10GeV(nMu10GeV);
1444     evTag->SetNumOfElectronsAbove1GeV(nEl1GeV);
1445     evTag->SetNumOfElectronsAbove3GeV(nEl3GeV);
1446     evTag->SetNumOfElectronsAbove10GeV(nEl10GeV);
1447
1448     evTag->SetNumOfPHOSClusters(esd->GetNumberOfPHOSClusters());
1449     evTag->SetNumOfEMCALClusters(esd->GetNumberOfEMCALClusters());
1450
1451     evTag->SetTotalMomentum(totalP);
1452     evTag->SetMeanPt(meanPt);
1453     evTag->SetMaxPt(maxPt);
1454
1455     tag->SetLHCTag(lhcLuminosity,lhcState);
1456     tag->SetDetectorTag(detectorMask);
1457
1458     tag->SetRunId(iInitRunNumber);
1459     tag->SetRunStartTime(t1->GetDate());
1460     tag->SetRunStopTime(t2->GetDate());
1461     tag->SetBeamEnergy(beamenergy);
1462     tag->SetBeamType(beamtype);
1463
1464     tag->AddEventTag(*evTag);
1465   }
1466
1467   ftag->cd();
1468   ttag.Fill();
1469   tag->Clear();
1470   ttag.Write();
1471   ftag->Close();
1472   file->cd();
1473   delete tag;
1474   delete evTag;
1475 }
1476
1477 //_____________________________________________________________________________
1478 void AliESDTagCreator::SwitchOffBranches() const {
1479   //
1480   // Switch of branches on user request
1481   TObjArray * tokens = fBranches.Tokenize(" ");
1482   Int_t ntok = tokens->GetEntries();
1483   for (Int_t i = 0; i < ntok; i++)  {
1484     TString str = ((TObjString*) tokens->At(i))->GetString();
1485     fChain->SetBranchStatus(Form("%s%s%s","*", str.Data(), "*"), 0);
1486     AliInfo(Form("Branch %s switched off \n", str.Data()));
1487   }
1488 }