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