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::se01"),
  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
 208   //Defining new tag objects
 209   AliRunTag *newTag = new AliRunTag();
 210   TTree ttag("T","A Tree with event tags");
 211   TBranch * btag = ttag.Branch("AliTAG", &newTag);
 212   btag->SetCompressionLevel(9);
 213   for(Int_t iTagFiles = 0; iTagFiles < fgChain->GetEntries(); iTagFiles++) {
 214     fgChain->GetEntry(iTagFiles);
 215     newTag = tag;
 216     ttag.Fill();
 217     newTag->Clear();
 218   }//tag file loop
 219
 220   TString localFileName = "Run"; localFileName += tag->GetRunId();
 221   localFileName += ".Merged"; localFileName += ".ESD.tag.root";
 222
 223   TString filename = 0x0;
 224
 225   if(fStorage == 0) {
 226     filename = localFileName.Data();
 227     AliInfo(Form("Writing merged tags to local file: %s",filename.Data()));
 228   }
 229   else if(fStorage == 1) {
 230     TString alienFileName = "/alien";
 231     alienFileName += gGrid->Pwd();
 232     alienFileName += fgridpath.Data();
 233     alienFileName += "/";
 234     alienFileName +=  localFileName;
 235     alienFileName += "?se=";
 236     alienFileName += fSE.Data();
 237     filename = alienFileName.Data();
 238     AliInfo(Form("Writing merged tags to grid file: %s",filename.Data()));
 239   }
 240
 241   TFile* ftag = TFile::Open(filename, "recreate");
 242   ftag->cd();
 243   ttag.Write();
 244   ftag->Close();
 245
 246   delete newTag;
 247
 248   return kTRUE;
 249 }
 250
 251 //__________________________________________________________________________
 252 Bool_t AliTagCreator::MergeTags(TGridResult *result) {
 253   //Merges the tags that are listed in the TGridResult
 254   AliInfo(Form("Merging tags....."));
 255   TChain *fgChain = new TChain("T");
 256
 257   Int_t nEntries = result->GetEntries();
 258
 259   TString alienUrl;
 260   for(Int_t i = 0; i < nEntries; i++) {
 261     alienUrl = result->GetKey(i,"turl");
 262     fgChain->Add(alienUrl);
 263   }
 264   AliInfo(Form("Chained tag files: %d",fgChain->GetEntries()));
 265   AliRunTag *tag = new AliRunTag;
 266   fgChain->SetBranchAddress("AliTAG",&tag);
 267
 268   //Defining new tag objects
 269   AliRunTag *newTag = new AliRunTag();
 270   TTree ttag("T","A Tree with event tags");
 271   TBranch * btag = ttag.Branch("AliTAG", &newTag);
 272   btag->SetCompressionLevel(9);
 273   for(Int_t iTagFiles = 0; iTagFiles < fgChain->GetEntries(); iTagFiles++) {
 274     fgChain->GetEntry(iTagFiles);
 275     newTag = tag;
 276     ttag.Fill();
 277     newTag->Clear();
 278   }//tag file loop
 279
 280   TString localFileName = "Run"; localFileName += tag->GetRunId();
 281   localFileName += ".Merged"; localFileName += ".ESD.tag.root";
 282
 283   TString filename = 0x0;
 284
 285   if(fStorage == 0) {
 286     filename = localFileName.Data();
 287     AliInfo(Form("Writing merged tags to local file: %s",filename.Data()));
 288   }
 289   else if(fStorage == 1) {
 290     TString alienFileName = "/alien";
 291     alienFileName += gGrid->Pwd();
 292     alienFileName += fgridpath.Data();
 293     alienFileName += "/";
 294     alienFileName +=  localFileName;
 295     alienFileName += "?se=";
 296     alienFileName += fSE.Data();
 297     filename = alienFileName.Data();
 298     AliInfo(Form("Writing merged tags to grid file: %s",filename.Data()));
 299   }
 300
 301   TFile* ftag = TFile::Open(filename, "recreate");
 302   ftag->cd();
 303   ttag.Write();
 304   ftag->Close();
 305
 306   delete newTag;
 307
 308   return kTRUE;
 309 }
 310
 311 //_____________________________________________________________________________
 312 void AliTagCreator::CreateTag(TFile* file, const char *guid, const char *md5, const char *turl, Long64_t size, Int_t Counter) {
 313   //private method that creates tag files
 314   TString fguid = guid;
 315   TString fmd5 = md5;
 316   TString fturl = turl;
 317   /////////////
 318   //muon code//
 319   ////////////
 320   Double_t fMUONMASS = 0.105658369;
 321   //Variables
 322   Double_t fX,fY,fZ ;
 323   Double_t fThetaX, fThetaY, fPyz, fChisquare;
 324   Double_t fPxRec, fPyRec, fPzRec, fEnergy;
 325   Int_t fCharge;
 326   TLorentzVector fEPvector;
 327
 328   Float_t fZVertexCut = 10.0;
 329   Float_t fRhoVertexCut = 2.0;
 330
 331   Float_t fLowPtCut = 1.0;
 332   Float_t fHighPtCut = 3.0;
 333   Float_t fVeryHighPtCut = 10.0;
 334   ////////////
 335
 336   Double_t partFrac[5] = {0.01, 0.01, 0.85, 0.10, 0.05};
 337
 338   // Creates the tags for all the events in a given ESD file
 339   Bool_t fIsSim = kTRUE;
 340   Int_t ntrack;
 341   Int_t nProtons, nKaons, nPions, nMuons, nElectrons;
 342   Int_t nPos, nNeg, nNeutr;
 343   Int_t nK0s, nNeutrons, nPi0s, nGamas;
 344   Int_t nCh1GeV, nCh3GeV, nCh10GeV;
 345   Int_t nMu1GeV, nMu3GeV, nMu10GeV;
 346   Int_t nEl1GeV, nEl3GeV, nEl10GeV;
 347   Float_t maxPt = .0, meanPt = .0, totalP = .0;
 348   Int_t fVertexflag;
 349   Int_t iRunNumber = 0;
 350   TString fVertexName;
 351
 352   AliRunTag *tag = new AliRunTag();
 353   AliEventTag *evTag = new AliEventTag();
 354   TTree ttag("T","A Tree with event tags");
 355   TBranch * btag = ttag.Branch("AliTAG", &tag);
 356   btag->SetCompressionLevel(9);
 357
 358   AliInfo(Form("Creating the tags......."));
 359
 360   Int_t firstEvent = 0,lastEvent = 0;
 361   TTree *t = (TTree*) file->Get("esdTree");
 362   TBranch * b = t->GetBranch("ESD");
 363   AliESD *esd = 0;
 364   b->SetAddress(&esd);
 365
 366   b->GetEntry(0);
 367   Int_t iInitRunNumber = esd->GetRunNumber();
 368
 369   Int_t iNumberOfEvents = b->GetEntries();
 370   for (Int_t iEventNumber = 0; iEventNumber < iNumberOfEvents; iEventNumber++) {
 371     ntrack = 0;
 372     nPos = 0;
 373     nNeg = 0;
 374     nNeutr =0;
 375     nK0s = 0;
 376     nNeutrons = 0;
 377     nPi0s = 0;
 378     nGamas = 0;
 379     nProtons = 0;
 380     nKaons = 0;
 381     nPions = 0;
 382     nMuons = 0;
 383     nElectrons = 0;
 384     nCh1GeV = 0;
 385     nCh3GeV = 0;
 386     nCh10GeV = 0;
 387     nMu1GeV = 0;
 388     nMu3GeV = 0;
 389     nMu10GeV = 0;
 390     nEl1GeV = 0;
 391     nEl3GeV = 0;
 392     nEl10GeV = 0;
 393     maxPt = .0;
 394     meanPt = .0;
 395     totalP = .0;
 396     fVertexflag = 1;
 397
 398     b->GetEntry(iEventNumber);
 399     iRunNumber = esd->GetRunNumber();
 400     if(iRunNumber != iInitRunNumber) AliFatal("Inconsistency of run numbers in the AliESD!!!");
 401     const AliESDVertex * vertexIn = esd->GetVertex();
 402     fVertexName = vertexIn->GetName();
 403     if(fVertexName == "default") fVertexflag = 0;
 404
 405     for (Int_t iTrackNumber = 0; iTrackNumber < esd->GetNumberOfTracks(); iTrackNumber++) {
 406       AliESDtrack * esdTrack = esd->GetTrack(iTrackNumber);
 407       if(esdTrack->GetLabel() != 0) fIsSim = kTRUE;
 408       else if(esdTrack->GetLabel() == 0) fIsSim = kFALSE;
 409       UInt_t status = esdTrack->GetStatus();
 410
 411       //select only tracks with ITS refit
 412       if ((status&AliESDtrack::kITSrefit)==0) continue;
 413       //select only tracks with TPC refit
 414       if ((status&AliESDtrack::kTPCrefit)==0) continue;
 415
 416       //select only tracks with the "combined PID"
 417       if ((status&AliESDtrack::kESDpid)==0) continue;
 418       Double_t p[3];
 419       esdTrack->GetPxPyPz(p);
 420       Double_t momentum = sqrt(pow(p[0],2) + pow(p[1],2) + pow(p[2],2));
 421       Double_t fPt = sqrt(pow(p[0],2) + pow(p[1],2));
 422       totalP += momentum;
 423       meanPt += fPt;
 424       if(fPt > maxPt) maxPt = fPt;
 425
 426       if(esdTrack->GetSign() > 0) {
 427         nPos++;
 428         if(fPt > fLowPtCut) nCh1GeV++;
 429         if(fPt > fHighPtCut) nCh3GeV++;
 430         if(fPt > fVeryHighPtCut) nCh10GeV++;
 431       }
 432       if(esdTrack->GetSign() < 0) {
 433         nNeg++;
 434         if(fPt > fLowPtCut) nCh1GeV++;
 435         if(fPt > fHighPtCut) nCh3GeV++;
 436         if(fPt > fVeryHighPtCut) nCh10GeV++;
 437       }
 438       if(esdTrack->GetSign() == 0) nNeutr++;
 439
 440       //PID
 441       Double_t prob[5];
 442       esdTrack->GetESDpid(prob);
 443
 444       Double_t rcc = 0.0;
 445       for(Int_t i = 0; i < AliPID::kSPECIES; i++) rcc += prob[i]*partFrac[i];
 446       if(rcc == 0.0) continue;
 447       //Bayes' formula
 448       Double_t w[5];
 449       for(Int_t i = 0; i < AliPID::kSPECIES; i++) w[i] = prob[i]*partFrac[i]/rcc;
 450
 451       //protons
 452       if ((w[4]>w[3])&&(w[4]>w[2])&&(w[4]>w[1])&&(w[4]>w[0])) nProtons++;
 453       //kaons
 454       if ((w[3]>w[4])&&(w[3]>w[2])&&(w[3]>w[1])&&(w[3]>w[0])) nKaons++;
 455       //pions
 456       if ((w[2]>w[4])&&(w[2]>w[3])&&(w[2]>w[1])&&(w[2]>w[0])) nPions++;
 457       //electrons
 458       if ((w[0]>w[4])&&(w[0]>w[3])&&(w[0]>w[2])&&(w[0]>w[1])) {
 459         nElectrons++;
 460         if(fPt > fLowPtCut) nEl1GeV++;
 461         if(fPt > fHighPtCut) nEl3GeV++;
 462         if(fPt > fVeryHighPtCut) nEl10GeV++;
 463       }
 464       ntrack++;
 465     }//esd track loop
 466
 467     /////////////
 468     //muon code//
 469     ////////////
 470     Int_t nMuonTracks = esd->GetNumberOfMuonTracks();
 471     // loop over all reconstructed tracks (also first track of combination)
 472     for (Int_t iTrack = 0; iTrack <  nMuonTracks;  iTrack++) {
 473       AliESDMuonTrack* muonTrack = esd->GetMuonTrack(iTrack);
 474       if (muonTrack == 0x0) continue;
 475
 476       // Coordinates at vertex
 477       fZ = muonTrack->GetZ();
 478       fY = muonTrack->GetBendingCoor();
 479       fX = muonTrack->GetNonBendingCoor();
 480
 481       fThetaX = muonTrack->GetThetaX();
 482       fThetaY = muonTrack->GetThetaY();
 483
 484       fPyz = 1./TMath::Abs(muonTrack->GetInverseBendingMomentum());
 485       fPzRec = - fPyz / TMath::Sqrt(1.0 + TMath::Tan(fThetaY)*TMath::Tan(fThetaY));
 486       fPxRec = fPzRec * TMath::Tan(fThetaX);
 487       fPyRec = fPzRec * TMath::Tan(fThetaY);
 488       fCharge = Int_t(TMath::Sign(1.,muonTrack->GetInverseBendingMomentum()));
 489
 490       //ChiSquare of the track if needed
 491       fChisquare = muonTrack->GetChi2()/(2.0 * muonTrack->GetNHit() - 5);
 492       fEnergy = TMath::Sqrt(fMUONMASS * fMUONMASS + fPxRec * fPxRec + fPyRec * fPyRec + fPzRec * fPzRec);
 493       fEPvector.SetPxPyPzE(fPxRec, fPyRec, fPzRec, fEnergy);
 494
 495       // total number of muons inside a vertex cut
 496       if((TMath::Abs(fZ)<fZVertexCut) && (TMath::Sqrt(fY*fY+fX*fX)<fRhoVertexCut)) {
 497         nMuons++;
 498         if(fEPvector.Pt() > fLowPtCut) {
 499           nMu1GeV++;
 500           if(fEPvector.Pt() > fHighPtCut) {
 501             nMu3GeV++;
 502             if (fEPvector.Pt() > fVeryHighPtCut) {
 503               nMu10GeV++;
 504             }
 505           }
 506         }
 507       }
 508     }//muon track loop
 509
 510     // Fill the event tags
 511     if(ntrack != 0) meanPt = meanPt/ntrack;
 512
 513     evTag->SetEventId(iEventNumber+1);
 514     evTag->SetGUID(fguid);
 515     evTag->SetMD5(fmd5);
 516     evTag->SetTURL(fturl);
 517     evTag->SetSize(size);
 518     evTag->SetVertexX(vertexIn->GetXv());
 519     evTag->SetVertexY(vertexIn->GetYv());
 520     evTag->SetVertexZ(vertexIn->GetZv());
 521     evTag->SetVertexZError(vertexIn->GetZRes());
 522     evTag->SetVertexFlag(fVertexflag);
 523
 524     evTag->SetT0VertexZ(esd->GetT0zVertex());
 525
 526     evTag->SetTriggerMask(esd->GetTriggerMask());
 527     evTag->SetTriggerCluster(esd->GetTriggerCluster());
 528
 529     evTag->SetZDCNeutron1Energy(esd->GetZDCN1Energy());
 530     evTag->SetZDCProton1Energy(esd->GetZDCP1Energy());
 531     evTag->SetZDCEMEnergy(esd->GetZDCEMEnergy());
 532     evTag->SetZDCNeutron1Energy(esd->GetZDCN2Energy());
 533     evTag->SetZDCProton1Energy(esd->GetZDCP2Energy());
 534     evTag->SetNumOfParticipants(esd->GetZDCParticipants());
 535
 536
 537     evTag->SetNumOfTracks(esd->GetNumberOfTracks());
 538     evTag->SetNumOfPosTracks(nPos);
 539     evTag->SetNumOfNegTracks(nNeg);
 540     evTag->SetNumOfNeutrTracks(nNeutr);
 541
 542     evTag->SetNumOfV0s(esd->GetNumberOfV0s());
 543     evTag->SetNumOfCascades(esd->GetNumberOfCascades());
 544     evTag->SetNumOfKinks(esd->GetNumberOfKinks());
 545     evTag->SetNumOfPMDTracks(esd->GetNumberOfPmdTracks());
 546
 547     evTag->SetNumOfProtons(nProtons);
 548     evTag->SetNumOfKaons(nKaons);
 549     evTag->SetNumOfPions(nPions);
 550     evTag->SetNumOfMuons(nMuons);
 551     evTag->SetNumOfElectrons(nElectrons);
 552     evTag->SetNumOfPhotons(nGamas);
 553     evTag->SetNumOfPi0s(nPi0s);
 554     evTag->SetNumOfNeutrons(nNeutrons);
 555     evTag->SetNumOfKaon0s(nK0s);
 556
 557     evTag->SetNumOfChargedAbove1GeV(nCh1GeV);
 558     evTag->SetNumOfChargedAbove3GeV(nCh3GeV);
 559     evTag->SetNumOfChargedAbove10GeV(nCh10GeV);
 560     evTag->SetNumOfMuonsAbove1GeV(nMu1GeV);
 561     evTag->SetNumOfMuonsAbove3GeV(nMu3GeV);
 562     evTag->SetNumOfMuonsAbove10GeV(nMu10GeV);
 563     evTag->SetNumOfElectronsAbove1GeV(nEl1GeV);
 564     evTag->SetNumOfElectronsAbove3GeV(nEl3GeV);
 565     evTag->SetNumOfElectronsAbove10GeV(nEl10GeV);
 566
 567     evTag->SetNumOfPHOSClusters(esd->GetNumberOfPHOSClusters());
 568     evTag->SetNumOfEMCALClusters(esd->GetNumberOfEMCALClusters());
 569
 570     evTag->SetTotalMomentum(totalP);
 571     evTag->SetMeanPt(meanPt);
 572     evTag->SetMaxPt(maxPt);
 573
 574     tag->SetRunId(iInitRunNumber);
 575     if(fIsSim) tag->SetDataType(0);
 576     else tag->SetDataType(1);
 577     tag->AddEventTag(*evTag);
 578   }//event loop
 579   lastEvent = iNumberOfEvents;
 580
 581   t->Delete("");
 582
 583   ttag.Fill();
 584   tag->Clear();
 585
 586   TString localFileName = "Run"; localFileName += tag->GetRunId();
 587   localFileName += ".Event"; localFileName += firstEvent; localFileName += "_"; localFileName += lastEvent; localFileName += "."; localFileName += Counter;
 588   localFileName += ".ESD.tag.root";
 589
 590   TString fileName;
 591
 592   if(fStorage == 0) {
 593     fileName = localFileName.Data();
 594     AliInfo(Form("Writing tags to local file: %s",fileName.Data()));
 595   }
 596   else if(fStorage == 1) {
 597     TString alienLocation = "/alien";
 598     alienLocation += gGrid->Pwd();
 599     alienLocation += fgridpath.Data();
 600     alienLocation += "/";
 601     alienLocation +=  localFileName;
 602     alienLocation += "?se=";
 603     alienLocation += fSE.Data();
 604     fileName = alienLocation.Data();
 605     AliInfo(Form("Writing tags to grid file: %s",fileName.Data()));
 606   }
 607
 608   TFile* ftag = TFile::Open(fileName, "recreate");
 609   ftag->cd();
 610   ttag.Write();
 611   ftag->Close();
 612
 613   delete ftag;
 614   delete esd;
 615
 616   delete tag;
 617 }
 618
 619 //_____________________________________________________________________________
 620 void AliTagCreator::CreateTag(TFile* file, const char *filepath, Int_t Counter) {
 621   //private method that creates tag files
 622   /////////////
 623   //muon code//
 624   ////////////
 625   Double_t fMUONMASS = 0.105658369;
 626   //Variables
 627   Double_t fX,fY,fZ ;
 628   Double_t fThetaX, fThetaY, fPyz, fChisquare;
 629   Double_t fPxRec, fPyRec, fPzRec, fEnergy;
 630   Int_t fCharge;
 631   TLorentzVector fEPvector;
 632
 633   Float_t fZVertexCut = 10.0;
 634   Float_t fRhoVertexCut = 2.0;
 635
 636   Float_t fLowPtCut = 1.0;
 637   Float_t fHighPtCut = 3.0;
 638   Float_t fVeryHighPtCut = 10.0;
 639   ////////////
 640
 641   Double_t partFrac[5] = {0.01, 0.01, 0.85, 0.10, 0.05};
 642
 643   // Creates the tags for all the events in a given ESD file
 644   Bool_t fIsSim = kTRUE;
 645   Int_t ntrack;
 646   Int_t nProtons, nKaons, nPions, nMuons, nElectrons;
 647   Int_t nPos, nNeg, nNeutr;
 648   Int_t nK0s, nNeutrons, nPi0s, nGamas;
 649   Int_t nCh1GeV, nCh3GeV, nCh10GeV;
 650   Int_t nMu1GeV, nMu3GeV, nMu10GeV;
 651   Int_t nEl1GeV, nEl3GeV, nEl10GeV;
 652   Float_t maxPt = .0, meanPt = .0, totalP = .0;
 653   Int_t fVertexflag;
 654   Int_t iRunNumber = 0;
 655   TString fVertexName;
 656
 657   AliRunTag *tag = new AliRunTag();
 658   AliEventTag *evTag = new AliEventTag();
 659   TTree ttag("T","A Tree with event tags");
 660   TBranch * btag = ttag.Branch("AliTAG", &tag);
 661   btag->SetCompressionLevel(9);
 662
 663   AliInfo(Form("Creating the tags......."));
 664
 665   Int_t firstEvent = 0,lastEvent = 0;
 666
 667   TTree *t = (TTree*) file->Get("esdTree");
 668   TBranch * b = t->GetBranch("ESD");
 669   AliESD *esd = 0;
 670   b->SetAddress(&esd);
 671
 672   b->GetEntry(0);
 673   Int_t iInitRunNumber = esd->GetRunNumber();
 674
 675   Int_t iNumberOfEvents = b->GetEntries();
 676   for (Int_t iEventNumber = 0; iEventNumber < iNumberOfEvents; iEventNumber++) {
 677     ntrack = 0;
 678     nPos = 0;
 679     nNeg = 0;
 680     nNeutr =0;
 681     nK0s = 0;
 682     nNeutrons = 0;
 683     nPi0s = 0;
 684     nGamas = 0;
 685     nProtons = 0;
 686     nKaons = 0;
 687     nPions = 0;
 688     nMuons = 0;
 689     nElectrons = 0;
 690     nCh1GeV = 0;
 691     nCh3GeV = 0;
 692     nCh10GeV = 0;
 693     nMu1GeV = 0;
 694     nMu3GeV = 0;
 695     nMu10GeV = 0;
 696     nEl1GeV = 0;
 697     nEl3GeV = 0;
 698     nEl10GeV = 0;
 699     maxPt = .0;
 700     meanPt = .0;
 701     totalP = .0;
 702     fVertexflag = 1;
 703
 704     b->GetEntry(iEventNumber);
 705     iRunNumber = esd->GetRunNumber();
 706     if(iRunNumber != iInitRunNumber) AliFatal("Inconsistency of run numbers in the AliESD!!!");
 707     const AliESDVertex * vertexIn = esd->GetVertex();
 708     fVertexName = vertexIn->GetName();
 709     if(fVertexName == "default") fVertexflag = 0;
 710
 711     for (Int_t iTrackNumber = 0; iTrackNumber < esd->GetNumberOfTracks(); iTrackNumber++) {
 712       AliESDtrack * esdTrack = esd->GetTrack(iTrackNumber);
 713       if(esdTrack->GetLabel() != 0) fIsSim = kTRUE;
 714       else if(esdTrack->GetLabel() == 0) fIsSim = kFALSE;
 715       UInt_t status = esdTrack->GetStatus();
 716
 717       //select only tracks with ITS refit
 718       if ((status&AliESDtrack::kITSrefit)==0) continue;
 719       //select only tracks with TPC refit
 720       if ((status&AliESDtrack::kTPCrefit)==0) continue;
 721
 722       //select only tracks with the "combined PID"
 723       if ((status&AliESDtrack::kESDpid)==0) continue;
 724       Double_t p[3];
 725       esdTrack->GetPxPyPz(p);
 726       Double_t momentum = sqrt(pow(p[0],2) + pow(p[1],2) + pow(p[2],2));
 727       Double_t fPt = sqrt(pow(p[0],2) + pow(p[1],2));
 728       totalP += momentum;
 729       meanPt += fPt;
 730       if(fPt > maxPt) maxPt = fPt;
 731
 732       if(esdTrack->GetSign() > 0) {
 733         nPos++;
 734         if(fPt > fLowPtCut) nCh1GeV++;
 735         if(fPt > fHighPtCut) nCh3GeV++;
 736         if(fPt > fVeryHighPtCut) nCh10GeV++;
 737       }
 738       if(esdTrack->GetSign() < 0) {
 739         nNeg++;
 740         if(fPt > fLowPtCut) nCh1GeV++;
 741         if(fPt > fHighPtCut) nCh3GeV++;
 742         if(fPt > fVeryHighPtCut) nCh10GeV++;
 743       }
 744       if(esdTrack->GetSign() == 0) nNeutr++;
 745
 746       //PID
 747       Double_t prob[5];
 748       esdTrack->GetESDpid(prob);
 749
 750       Double_t rcc = 0.0;
 751       for(Int_t i = 0; i < AliPID::kSPECIES; i++) rcc += prob[i]*partFrac[i];
 752       if(rcc == 0.0) continue;
 753       //Bayes' formula
 754       Double_t w[5];
 755       for(Int_t i = 0; i < AliPID::kSPECIES; i++) w[i] = prob[i]*partFrac[i]/rcc;
 756
 757       //protons
 758       if ((w[4]>w[3])&&(w[4]>w[2])&&(w[4]>w[1])&&(w[4]>w[0])) nProtons++;
 759       //kaons
 760       if ((w[3]>w[4])&&(w[3]>w[2])&&(w[3]>w[1])&&(w[3]>w[0])) nKaons++;
 761       //pions
 762       if ((w[2]>w[4])&&(w[2]>w[3])&&(w[2]>w[1])&&(w[2]>w[0])) nPions++;
 763       //electrons
 764       if ((w[0]>w[4])&&(w[0]>w[3])&&(w[0]>w[2])&&(w[0]>w[1])) {
 765         nElectrons++;
 766         if(fPt > fLowPtCut) nEl1GeV++;
 767         if(fPt > fHighPtCut) nEl3GeV++;
 768         if(fPt > fVeryHighPtCut) nEl10GeV++;
 769       }
 770       ntrack++;
 771     }//esd track loop
 772
 773     /////////////
 774     //muon code//
 775     ////////////
 776     Int_t nMuonTracks = esd->GetNumberOfMuonTracks();
 777     // loop over all reconstructed tracks (also first track of combination)
 778     for (Int_t iTrack = 0; iTrack <  nMuonTracks;  iTrack++) {
 779       AliESDMuonTrack* muonTrack = esd->GetMuonTrack(iTrack);
 780       if (muonTrack == 0x0) continue;
 781
 782       // Coordinates at vertex
 783       fZ = muonTrack->GetZ();
 784       fY = muonTrack->GetBendingCoor();
 785       fX = muonTrack->GetNonBendingCoor();
 786
 787       fThetaX = muonTrack->GetThetaX();
 788       fThetaY = muonTrack->GetThetaY();
 789
 790       fPyz = 1./TMath::Abs(muonTrack->GetInverseBendingMomentum());
 791       fPzRec = - fPyz / TMath::Sqrt(1.0 + TMath::Tan(fThetaY)*TMath::Tan(fThetaY));
 792       fPxRec = fPzRec * TMath::Tan(fThetaX);
 793       fPyRec = fPzRec * TMath::Tan(fThetaY);
 794       fCharge = Int_t(TMath::Sign(1.,muonTrack->GetInverseBendingMomentum()));
 795
 796       //ChiSquare of the track if needed
 797       fChisquare = muonTrack->GetChi2()/(2.0 * muonTrack->GetNHit() - 5);
 798       fEnergy = TMath::Sqrt(fMUONMASS * fMUONMASS + fPxRec * fPxRec + fPyRec * fPyRec + fPzRec * fPzRec);
 799       fEPvector.SetPxPyPzE(fPxRec, fPyRec, fPzRec, fEnergy);
 800
 801       // total number of muons inside a vertex cut
 802       if((TMath::Abs(fZ)<fZVertexCut) && (TMath::Sqrt(fY*fY+fX*fX)<fRhoVertexCut)) {
 803         nMuons++;
 804         if(fEPvector.Pt() > fLowPtCut) {
 805           nMu1GeV++;
 806           if(fEPvector.Pt() > fHighPtCut) {
 807             nMu3GeV++;
 808             if (fEPvector.Pt() > fVeryHighPtCut) {
 809               nMu10GeV++;
 810             }
 811           }
 812         }
 813       }
 814     }//muon track loop
 815
 816     // Fill the event tags
 817     if(ntrack != 0) meanPt = meanPt/ntrack;
 818
 819     evTag->SetEventId(iEventNumber+1);
 820     evTag->SetPath(filepath);
 821
 822     evTag->SetVertexX(vertexIn->GetXv());
 823     evTag->SetVertexY(vertexIn->GetYv());
 824     evTag->SetVertexZ(vertexIn->GetZv());
 825     evTag->SetVertexZError(vertexIn->GetZRes());
 826     evTag->SetVertexFlag(fVertexflag);
 827
 828     evTag->SetT0VertexZ(esd->GetT0zVertex());
 829
 830     evTag->SetTriggerMask(esd->GetTriggerMask());
 831     evTag->SetTriggerCluster(esd->GetTriggerCluster());
 832
 833     evTag->SetZDCNeutron1Energy(esd->GetZDCN1Energy());
 834     evTag->SetZDCProton1Energy(esd->GetZDCP1Energy());
 835     evTag->SetZDCEMEnergy(esd->GetZDCEMEnergy());
 836     evTag->SetZDCNeutron1Energy(esd->GetZDCN2Energy());
 837     evTag->SetZDCProton1Energy(esd->GetZDCP2Energy());
 838     evTag->SetNumOfParticipants(esd->GetZDCParticipants());
 839
 840
 841     evTag->SetNumOfTracks(esd->GetNumberOfTracks());
 842     evTag->SetNumOfPosTracks(nPos);
 843     evTag->SetNumOfNegTracks(nNeg);
 844     evTag->SetNumOfNeutrTracks(nNeutr);
 845
 846     evTag->SetNumOfV0s(esd->GetNumberOfV0s());
 847     evTag->SetNumOfCascades(esd->GetNumberOfCascades());
 848     evTag->SetNumOfKinks(esd->GetNumberOfKinks());
 849     evTag->SetNumOfPMDTracks(esd->GetNumberOfPmdTracks());
 850
 851     evTag->SetNumOfProtons(nProtons);
 852     evTag->SetNumOfKaons(nKaons);
 853     evTag->SetNumOfPions(nPions);
 854     evTag->SetNumOfMuons(nMuons);
 855     evTag->SetNumOfElectrons(nElectrons);
 856     evTag->SetNumOfPhotons(nGamas);
 857     evTag->SetNumOfPi0s(nPi0s);
 858     evTag->SetNumOfNeutrons(nNeutrons);
 859     evTag->SetNumOfKaon0s(nK0s);
 860
 861     evTag->SetNumOfChargedAbove1GeV(nCh1GeV);
 862     evTag->SetNumOfChargedAbove3GeV(nCh3GeV);
 863     evTag->SetNumOfChargedAbove10GeV(nCh10GeV);
 864     evTag->SetNumOfMuonsAbove1GeV(nMu1GeV);
 865     evTag->SetNumOfMuonsAbove3GeV(nMu3GeV);
 866     evTag->SetNumOfMuonsAbove10GeV(nMu10GeV);
 867     evTag->SetNumOfElectronsAbove1GeV(nEl1GeV);
 868     evTag->SetNumOfElectronsAbove3GeV(nEl3GeV);
 869     evTag->SetNumOfElectronsAbove10GeV(nEl10GeV);
 870
 871     evTag->SetNumOfPHOSClusters(esd->GetNumberOfPHOSClusters());
 872     evTag->SetNumOfEMCALClusters(esd->GetNumberOfEMCALClusters());
 873
 874     evTag->SetTotalMomentum(totalP);
 875     evTag->SetMeanPt(meanPt);
 876     evTag->SetMaxPt(maxPt);
 877
 878     tag->SetRunId(iInitRunNumber);
 879     if(fIsSim) tag->SetDataType(0);
 880     else tag->SetDataType(1);
 881     tag->AddEventTag(*evTag);
 882   }//event loop
 883   lastEvent = iNumberOfEvents;
 884
 885   t->Delete("");
 886
 887   ttag.Fill();
 888   tag->Clear();
 889
 890   TString localFileName = "Run"; localFileName += tag->GetRunId();
 891   localFileName += ".Event"; localFileName += firstEvent; localFileName += "_"; localFileName += lastEvent; localFileName += "."; localFileName += Counter;
 892   localFileName += ".ESD.tag.root";
 893
 894   TString fileName;
 895
 896   if(fStorage == 0) {
 897     fileName = localFileName.Data();
 898     AliInfo(Form("Writing tags to local file: %s",fileName.Data()));
 899   }
 900   else if(fStorage == 1) {
 901     TString alienLocation = "/alien";
 902     alienLocation += gGrid->Pwd();
 903     alienLocation += fgridpath.Data();
 904     alienLocation += "/";
 905     alienLocation +=  localFileName;
 906     alienLocation += "?se=";
 907     alienLocation += fSE.Data();
 908     fileName = alienLocation.Data();
 909     AliInfo(Form("Writing tags to grid file: %s",fileName.Data()));
 910   }
 911
 912   TFile* ftag = TFile::Open(fileName, "recreate");
 913   ftag->cd();
 914   ttag.Write();
 915   ftag->Close();
 916
 917   delete ftag;
 918   delete esd;
 919
 920   delete tag;
 921 }