STEER/AliTagCreator.cxx

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