STEER/AliAODTagCreator.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 //           AliAODTagCreator 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 #include <TRefArray.h>
  31
  32 //ROOT-AliEn
  33 #include <TGrid.h>
  34 #include <TGridResult.h>
  35
  36 //AliRoot
  37 #include "AliRunTag.h"
  38 #include "AliEventTag.h"
  39 #include "AliPID.h"
  40 #include "AliAODEvent.h"
  41 #include "AliAODVertex.h"
  42 #include "AliLog.h"
  43
  44 #include "AliAODTagCreator.h"
  45
  46
  47 ClassImp(AliAODTagCreator)
  48
  49
  50 //______________________________________________________________________________
  51   AliAODTagCreator::AliAODTagCreator() :
  52     AliTagCreator(),
  53     fChain(0),
  54     fAODEvent(0),
  55     fTreeT(0),
  56     fRunTag(0),
  57     fTreeTEsd(0),
  58     fRunTagEsd(0)
  59 {
  60   //==============Default constructor for a AliAODTagCreator================
  61 }
  62
  63 //______________________________________________________________________________
  64 AliAODTagCreator::~AliAODTagCreator() {
  65 //================Default destructor for a AliAODTagCreator===================
  66   delete fChain;
  67   delete fAODEvent;
  68 }
  69
  70 //______________________________________________________________________________
  71 Bool_t AliAODTagCreator::ReadGridCollection(TGridResult *fresult) {
  72   // Reads the entry of the TGridResult and creates the tags
  73   Int_t nEntries = fresult->GetEntries();
  74
  75   TString alienUrl;
  76   const char* guid;
  77   const char* md5;
  78   const char* turl;
  79   Long64_t size = -1;
  80
  81   fChain = new TChain("aodTree");
  82
  83   for(Int_t i = 0; i < nEntries; i++) {
  84     alienUrl = fresult->GetKey(i,"turl");
  85     guid = fresult->GetKey(i,"guid");
  86     if(fresult->GetKey(i,"size")) size = atol (fresult->GetKey(i,"size"));
  87     md5 = fresult->GetKey(i,"md5");
  88     turl = fresult->GetKey(i,"turl");
  89     if(md5 && !strlen(guid)) md5 = 0;
  90     if(guid && !strlen(guid)) guid = 0;
  91
  92     fChain->Add(alienUrl);
  93   }//grid result loop
  94
  95   AliInfo(Form("AOD chain created......."));
  96   AliInfo(Form("Chain entries: %d",fChain->GetEntries()));
  97
  98   CreateTag(fChain, "grid");
  99
 100   return kTRUE;
 101 }
 102
 103 //______________________________________________________________________________
 104 Bool_t AliAODTagCreator::ReadLocalCollection(const char *localpath, const char* pattern) {
 105   // Checks the different subdirs of the given local path and in the
 106   // case where it finds an AliAOD.root file it creates the tags
 107
 108   void *dira =  gSystem->OpenDirectory(localpath);
 109   Char_t fPath[256];
 110   const char * dirname  = 0x0;
 111   const char * filename = 0x0;
 112
 113   fChain = new TChain("aodTree");
 114
 115   while((dirname = gSystem->GetDirEntry(dira))) {
 116     sprintf(fPath,"%s/%s",localpath,dirname);
 117     void *dirb =  gSystem->OpenDirectory(fPath);
 118     while((filename = gSystem->GetDirEntry(dirb))) {
 119         TString bstr = dirname;
 120         if(bstr.Contains("..")) continue;
 121         if(strstr(filename,pattern)) {
 122             TString aodFileName;
 123             aodFileName = fPath;
 124             aodFileName += "/";
 125             aodFileName += pattern;
 126             fChain->Add(aodFileName);
 127         } //pattern check
 128     } //child directory's entry loop
 129   } //parent directory's entry loop
 130
 131   AliInfo(Form("AOD chain created......."));
 132   AliInfo(Form("Chain entries: %d",fChain->GetEntries()));
 133
 134   CreateTag(fChain, "local");
 135
 136   return kTRUE;
 137 }
 138
 139 //______________________________________________________________________________
 140 Bool_t AliAODTagCreator::ReadCAFCollection(const char *filename) {
 141   // Temporary solution for CAF: Takes as an input the ascii file that
 142   // lists the AODs stored in the SE of the CAF and creates the tags.
 143
 144   // Open the input stream
 145   ifstream in;
 146   in.open(filename);
 147
 148   TString aodfile;
 149
 150   fChain = new TChain("aodTree");
 151
 152   // Read the input list of files and add them to the chain
 153   while(in.good()) {
 154     in >> aodfile;
 155     if (!aodfile.Contains("root")) continue; // protection
 156     fChain->Add(aodfile);
 157   }
 158
 159   AliInfo(Form("AOD chain created......."));
 160   AliInfo(Form("Chain entries: %d",fChain->GetEntries()));
 161
 162   CreateTag(fChain, "proof");
 163
 164   return kTRUE;
 165 }
 166
 167 //__________________________________________________________________________
 168 void AliAODTagCreator::CreateAODTags(Int_t fFirstEvent, Int_t fLastEvent, TList */*grpList*/) {
 169   // Creates tag files for AODs
 170
 171     AliInfo(Form("Creating the AOD tags......."));
 172
 173     TFile *file = TFile::Open("AliAOD.root");
 174     if (!file || !file->IsOpen()) {
 175         AliError(Form("opening failed"));
 176         delete file;
 177         return ;
 178     }
 179
 180     fChain = new TChain("aodTree");
 181     fChain->Add("AliAOD.root");
 182
 183     fAODEvent = new AliAODEvent();
 184     fAODEvent->ReadFromTree(fChain);
 185
 186     Int_t lastEvent = 0;
 187     if(fLastEvent == -1) lastEvent = (Int_t)fChain->GetEntries();
 188     else lastEvent = fLastEvent;
 189
 190     char fileName[256];
 191     sprintf(fileName, "Run%d.Event%d_%d.AOD.tag.root",
 192             fAODEvent->GetRunNumber(), fFirstEvent, lastEvent );
 193     AliInfo(Form("writing tags to file %s", fileName));
 194     AliDebug(1, Form("writing tags to file %s", fileName));
 195
 196     TFile* ftag = TFile::Open(fileName, "recreate");
 197
 198     fRunTag = new AliRunTag();
 199     fTreeT = new TTree("T","A Tree with event tags");
 200     TBranch * btag = fTreeT->Branch("AliTAG", &fRunTag);
 201     btag->SetCompressionLevel(9);
 202
 203     CreateTags();
 204
 205     ftag->cd();
 206     fTreeT->Fill();
 207     fRunTag->Clear();
 208     fTreeT->Write();
 209     ftag->Close();
 210     file->cd();
 211     file->Close();
 212 }
 213
 214 //_____________________________________________________________________________
 215 void AliAODTagCreator::CreateTag(TChain* chain, const char *type) {
 216
 217     // Private method that creates tag files
 218     //
 219
 220     //reading the esd tag file
 221     fTreeTEsd = new TChain("T");
 222     const char * tagPattern = "ESD.tag";
 223     // Open the working directory
 224     void * dirp = gSystem->OpenDirectory(gSystem->pwd());
 225     const char * name = 0x0;
 226     // Add all files matching *pattern* to the chain
 227     while((name = gSystem->GetDirEntry(dirp))) {
 228         if (strstr(name,tagPattern)) fTreeTEsd->Add(name);
 229     }//directory loop
 230     AliInfo(Form("Chained tag files: %d", fTreeTEsd->GetEntries()));
 231
 232     fChain = chain;
 233
 234     TString fSession = type;
 235     TString fguid, fmd5, fturl;
 236     fAODEvent = new AliAODEvent();
 237     fAODEvent->ReadFromTree(fChain);
 238
 239     Int_t firstEvent = 0;
 240
 241     TString localFileName = "Run";
 242     localFileName += fAODEvent->GetRunNumber();
 243     localFileName += ".Event";
 244     localFileName += firstEvent;
 245     localFileName += "_";
 246     localFileName += chain->GetEntries(); //localFileName += "."; localFileName += Counter;
 247     localFileName += ".AOD.tag.root";
 248
 249     TString fileName;
 250
 251     if(fStorage == 0) {
 252         fileName = localFileName.Data();
 253         AliInfo(Form("Writing tags to local file: %s",fileName.Data()));
 254     }
 255     else if(fStorage == 1) {
 256         TString alienLocation = "/alien";
 257         alienLocation += gGrid->Pwd();
 258         alienLocation += fgridpath.Data();
 259         alienLocation += "/";
 260         alienLocation +=  localFileName;
 261         alienLocation += "?se=";
 262         alienLocation += fSE.Data();
 263         fileName = alienLocation.Data();
 264         AliInfo(Form("Writing tags to grid file: %s",fileName.Data()));
 265     }
 266
 267     TFile* ftag = TFile::Open(fileName, "recreate");
 268
 269     fRunTag = new AliRunTag();
 270     fTreeT  = new TTree("T", "A Tree with event tags");
 271     TBranch * btag = fTreeT->Branch("AliTAG", &fRunTag);
 272     btag->SetCompressionLevel(9);
 273
 274     // Access information from esd tag
 275     fRunTagEsd = new AliRunTag();
 276     fTreeTEsd->SetBranchAddress("AliTAG",&fRunTagEsd);
 277     // Creating new information of aod
 278     AliInfo(Form("Creating the AOD tags......."));
 279     CreateTags(type);
 280     ftag->cd();
 281     fRunTag->Clear();
 282     fTreeT->Write();
 283     ftag->Close();
 284 }
 285
 286 void AliAODTagCreator::CreateTags(const char* type)
 287 {
 288     // Event loop for tag creation
 289     TString fturl;
 290     TString fguid;
 291     Int_t   oldRun = -1;
 292     fChain->GetEntry(0);
 293     TFile *f = fChain->GetFile();
 294     TString ftempGuid = f->GetUUID().AsString();
 295     // Loop over events
 296     Int_t nEvents = fChain->GetEntries();
 297     Int_t ntags    = 0;
 298     Int_t tagentry = 0;
 299     const TClonesArray *evTagList = 0;
 300
 301     for (Int_t iEventNumber = 0; iEventNumber < nEvents; iEventNumber++) {
 302         // Copy old tag information
 303         if (iEventNumber >= ntags) {
 304             fTreeTEsd->GetEntry(tagentry++);
 305             fRunTag->CopyStandardContent(fRunTagEsd);
 306             evTagList = fRunTagEsd->GetEventTags();
 307             ntags += evTagList->GetEntries();
 308         }
 309
 310         // Create a new Tag
 311         AliEventTag* evTag = new AliEventTag();
 312         // Read event
 313         fChain->GetEntry(iEventNumber);
 314         if (iEventNumber == 0) oldRun = fAODEvent->GetRunNumber();
 315         // Reference to the input file
 316         TFile *file = fChain->GetFile();
 317         const TUrl *url = file->GetEndpointUrl();
 318         fguid = file->GetUUID().AsString();
 319
 320         if (!strcmp(type,"grid")) {
 321             TString fturltemp = "alien://"; fturltemp += url->GetFile();
 322             fturl = fturltemp(0,fturltemp.Index(".root",5,0,TString::kExact)+5);
 323         } else {
 324             fturl = url->GetFile();
 325         }
 326
 327         fAODEvent->GetStdContent();
 328
 329         // Fill the event tag from the aod informatiom
 330         FillEventTag(fAODEvent, evTag);
 331         // Set the event and input file references
 332         //evTag->SetEventId(iEventNumber+1);
 333         evTag->SetGUID(fguid);
 334         if(!strcmp(type,"grid")) {
 335             evTag->SetMD5("");
 336             evTag->SetTURL(fturl);
 337             evTag->SetSize(0);
 338             }
 339         else evTag->SetPath(fturl);
 340
 341         // Check if a new run has to be created
 342         // File has changed
 343         if(fguid != ftempGuid) {
 344             ftempGuid = fguid;
 345             fTreeT->Fill();
 346             fRunTag->Clear("");
 347         }
 348         // Run# has changed
 349         if (oldRun != (fAODEvent->GetRunNumber()))
 350         {
 351             oldRun = fAODEvent->GetRunNumber();
 352             fTreeT->Fill();
 353             fRunTag->Clear("");
 354         }
 355
 356         // Add the event tag
 357         fRunTag->AddEventTag(*evTag);
 358         delete evTag;
 359         // Last event
 360         if(iEventNumber+1 == fChain->GetEntries()) {
 361             fTreeT->Fill();
 362             fRunTag->Clear("");
 363         }
 364     }//event loop
 365 }
 366
 367
 368
 369 void AliAODTagCreator::FillEventTag(AliAODEvent* aod, AliEventTag* evTag)
 370 {
 371 //
 372 // Fill the event tag information
 373     //
 374     fAODEvent = aod;
 375
 376     //
 377     Float_t fLowPtCut      =  1.0;
 378     Float_t fHighPtCut     =  3.0;
 379     Float_t fVeryHighPtCut = 10.0;
 380     ////////////
 381     Double_t partFrac[10] = {0.01, 0.01, 0.85, 0.10, 0.05, 0., 0., 0., 0., 0.};
 382
 383     // Creates the tags for all the events in a given AOD file
 384     Int_t   ntrack = 0;
 385     Int_t   nPos = 0, nNeg = 0, nNeutr =0;
 386     Int_t   nKinks = 0, nV0s = 0, nCascades = 0;
 387     Int_t   nK0s = 0, nNeutrons = 0, nPi0s = 0, nGamas = 0;
 388     Int_t   nProtons = 0,  nKaons = 0, nPions = 0, nMuons = 0, nElectrons = 0, nFWMuons = 0;
 389     Int_t   nCh1GeV = 0, nCh3GeV = 0, nCh10GeV = 0;
 390     Int_t   nMu1GeV = 0, nMu3GeV = 0, nMu10GeV = 0;
 391     Int_t   nEl1GeV = 0, nEl3GeV = 0, nEl10GeV = 0;
 392     Float_t maxPt =  .0, etamaxPt = -999., phimaxPt = -999., meanPt = .0, totalP =  .0;
 393
 394     TRefArray tmp;
 395
 396
 397     // Primary Vertex
 398     AliAODVertex *pVertex = fAODEvent->GetPrimaryVertex();
 399     if (pVertex) {
 400         evTag->SetVertexX(pVertex->GetX());
 401         evTag->SetVertexY(pVertex->GetY());
 402         evTag->SetVertexZ(pVertex->GetZ());
 403         Double_t covmatrix[6];
 404         pVertex->GetCovarianceMatrix(covmatrix);
 405         evTag->SetVertexZError(sqrt(covmatrix[5]));
 406     }
 407     // loop over vertices
 408     Int_t nVtxs = fAODEvent->GetNumberOfVertices();
 409     for (Int_t nVtx = 0; nVtx < nVtxs; nVtx++) {
 410         AliAODVertex *vertex = fAODEvent->GetVertex(nVtx);
 411         if(vertex->GetType() == 1) nKinks    += 1;
 412         if(vertex->GetType() == 2) nV0s      += 1;
 413         if(vertex->GetType() == 3) nCascades += 1;
 414     }
 415     Int_t nTracks = fAODEvent->GetNTracks();
 416     for (Int_t nTr = 0; nTr < nTracks; nTr++) {
 417         AliAODTrack *track = fAODEvent->GetTrack(nTr);
 418
 419         Double_t fPt = track->Pt();
 420         if(fPt > maxPt) {
 421             maxPt = fPt;
 422             etamaxPt = track->Eta();
 423             phimaxPt = track->Phi();
 424         }
 425
 426         if(track->Charge() > 0) {
 427             nPos++;
 428             if(fPt > fLowPtCut)      nCh1GeV++;
 429             if(fPt > fHighPtCut)     nCh3GeV++;
 430             if(fPt > fVeryHighPtCut) nCh10GeV++;
 431         }
 432         if(track->Charge() < 0) {
 433             nNeg++;
 434             if(fPt > fLowPtCut)      nCh1GeV++;
 435             if(fPt > fHighPtCut)     nCh3GeV++;
 436             if(fPt > fVeryHighPtCut) nCh10GeV++;
 437             }
 438         if(track->Charge() == 0) nNeutr++;
 439         //PID
 440         const Double32_t *prob = track->PID();
 441         Double_t rcc = 0.0;
 442         for(Int_t i = 0; i < AliPID::kSPECIES; i++) rcc += prob[i]*partFrac[i];
 443         if(rcc == 0.0) continue;
 444         //Bayes' formula
 445         Double_t w[10];
 446         for(Int_t i = 0; i < AliPID::kSPECIES; i++) w[i] = prob[i]*partFrac[i]/rcc;
 447
 448         //protons
 449         if ((w[4]>w[3])&&(w[4]>w[2])&&(w[4]>w[1])&&(w[4]>w[0])) nProtons++;
 450         //kaons
 451         if ((w[3]>w[4])&&(w[3]>w[2])&&(w[3]>w[1])&&(w[3]>w[0])) nKaons++;
 452         //pions
 453         if ((w[2]>w[4])&&(w[2]>w[3])&&(w[2]>w[1])&&(w[2]>w[0])) nPions++;
 454         //muons
 455         if ((w[1]>w[4])&&(w[1]>w[3])&&(w[1]>w[2])&&(w[1]>w[0])) {
 456             nMuons++;
 457             if(fPt > fLowPtCut)      nMu1GeV++;
 458             if(fPt > fHighPtCut)     nMu3GeV++;
 459             if(fPt > fVeryHighPtCut) nMu10GeV++;
 460         }
 461         //electrons
 462         if ((w[0]>w[4])&&(w[0]>w[3])&&(w[0]>w[2])&&(w[0]>w[1])) {
 463             nElectrons++;
 464             if(fPt > fLowPtCut) nEl1GeV++;
 465             if(fPt > fHighPtCut) nEl3GeV++;
 466             if(fPt > fVeryHighPtCut) nEl10GeV++;
 467         }
 468         totalP += track->P();
 469         meanPt += fPt;
 470         ntrack++;
 471         // forward muons (in the dimuon spectrometer)
 472         if(track->IsMuonTrack()) nFWMuons++;
 473
 474     }//track loop
 475     //
 476     // Fill the event tags
 477     if(ntrack != 0)
 478         meanPt = meanPt/ntrack;
 479
 480
 481     evTag->SetNumOfTracks(nTracks);
 482     evTag->SetNumOfPosTracks(nPos);
 483     evTag->SetNumOfNegTracks(nNeg);
 484     evTag->SetNumOfNeutrTracks(nNeutr);
 485
 486     evTag->SetNumOfV0s(nV0s);
 487     evTag->SetNumOfCascades(nCascades);
 488     evTag->SetNumOfKinks(nKinks);
 489
 490     evTag->SetNumOfProtons(nProtons);
 491     evTag->SetNumOfKaons(nKaons);
 492     evTag->SetNumOfPions(nPions);
 493     evTag->SetNumOfMuons(nMuons);
 494     evTag->SetNumOfFWMuons(nFWMuons);
 495     evTag->SetNumOfElectrons(nElectrons);
 496     evTag->SetNumOfPhotons(nGamas);
 497     evTag->SetNumOfPi0s(nPi0s);
 498     evTag->SetNumOfNeutrons(nNeutrons);
 499     evTag->SetNumOfKaon0s(nK0s);
 500
 501     evTag->SetNumOfChargedAbove1GeV(nCh1GeV);
 502     evTag->SetNumOfChargedAbove3GeV(nCh3GeV);
 503     evTag->SetNumOfChargedAbove10GeV(nCh10GeV);
 504     evTag->SetNumOfMuonsAbove1GeV(nMu1GeV);
 505     evTag->SetNumOfMuonsAbove3GeV(nMu3GeV);
 506     evTag->SetNumOfMuonsAbove10GeV(nMu10GeV);
 507     evTag->SetNumOfElectronsAbove1GeV(nEl1GeV);
 508     evTag->SetNumOfElectronsAbove3GeV(nEl3GeV);
 509     evTag->SetNumOfElectronsAbove10GeV(nEl10GeV);
 510
 511     tmp.Clear();
 512     evTag->SetNumOfPHOSClusters(fAODEvent->GetPHOSClusters(&tmp));
 513     tmp.Clear();
 514     evTag->SetNumOfEMCALClusters(fAODEvent->GetEMCALClusters(&tmp));
 515
 516     evTag->SetTotalMomentum(totalP);
 517     evTag->SetMeanPt(meanPt);
 518     evTag->SetMaxPt(maxPt);
 519     evTag->SetEtaMaxPt(etamaxPt);
 520     evTag->SetPhiMaxPt(phimaxPt);
 521 }