91a746a5913b8fd0d7259bca9ce77db9577256f5
[u/mrichter/AliRoot.git] / 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: %l",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: %l",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: %l",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: %l", 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 }