fe8dd235fb284a383836900df5b0b62a97893b82
[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
31 //ROOT-AliEn
32 #include <TGrid.h>
33 #include <TGridResult.h>
34
35 //AliRoot
36 #include "AliRunTag.h"
37 #include "AliEventTag.h"
38 #include "AliPID.h"
39 #include "AliAODEvent.h"
40 #include "AliAODVertex.h"
41 #include "AliLog.h"
42
43 #include "AliAODTagCreator.h"
44
45
46 ClassImp(AliAODTagCreator)
47
48
49 //______________________________________________________________________________
50   AliAODTagCreator::AliAODTagCreator() :
51     AliTagCreator(), 
52     fChain(0), 
53     fAODEvent(0), 
54     fTreeT(0), 
55     fRunTag(0),
56     fTreeTEsd(0), 
57     fRunTagEsd(0)  
58 {
59   //==============Default constructor for a AliAODTagCreator================
60 }
61
62 //______________________________________________________________________________
63 AliAODTagCreator::~AliAODTagCreator() {
64 //================Default destructor for a AliAODTagCreator===================
65   delete fChain;
66   delete fAODEvent;
67 }
68
69 //______________________________________________________________________________
70 Bool_t AliAODTagCreator::ReadGridCollection(TGridResult *fresult) {
71   // Reads the entry of the TGridResult and creates the tags
72   Int_t nEntries = fresult->GetEntries();
73
74   TString alienUrl;
75   const char* guid;
76   const char* md5;
77   const char* turl;
78   Long64_t size = -1;
79
80   fChain = new TChain("aodTree");
81  
82   for(Int_t i = 0; i < nEntries; i++) {
83     alienUrl = fresult->GetKey(i,"turl");
84     guid = fresult->GetKey(i,"guid");
85     if(fresult->GetKey(i,"size")) size = atol (fresult->GetKey(i,"size"));
86     md5 = fresult->GetKey(i,"md5");
87     turl = fresult->GetKey(i,"turl");
88     if(md5 && !strlen(guid)) md5 = 0;
89     if(guid && !strlen(guid)) guid = 0;
90     
91     fChain->Add(alienUrl);
92   }//grid result loop
93   
94   AliInfo(Form("AOD chain created......."));    
95   AliInfo(Form("Chain entries: %d",fChain->GetEntries()));      
96
97   CreateTag(fChain, "grid");
98   
99   return kTRUE;
100 }
101
102 //______________________________________________________________________________
103 Bool_t AliAODTagCreator::ReadLocalCollection(const char *localpath) {
104   // Checks the different subdirs of the given local path and in the
105   // case where it finds an AliAODs.root file it creates the tags
106   
107   void *dira =  gSystem->OpenDirectory(localpath);
108   Char_t fPath[256];
109   const char * dirname = 0x0;
110   const char * filename = 0x0;
111   const char * pattern = "AliAODs.root"; 
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
191     char fileName[256];
192     sprintf(fileName, "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
222     //reading the esd tag file                                                  
223     fTreeTEsd = new TChain("T");
224     const char * tagPattern = "ESD.tag";
225     // Open the working directory
226     void * dirp = gSystem->OpenDirectory(gSystem->pwd());
227     const char * name = 0x0;
228     // Add all files matching *pattern* to the chain
229     while((name = gSystem->GetDirEntry(dirp))) {
230         if (strstr(name,tagPattern)) fTreeTEsd->Add(name);
231     }//directory loop
232     AliInfo(Form("Chained tag files: %d", fTreeTEsd->GetEntries()));
233   
234     
235     fChain = chain;
236     
237     TString fSession = type;
238     TString fguid, fmd5, fturl;
239     fAODEvent = new AliAODEvent();
240     fAODEvent->ReadFromTree(fChain);
241
242     Int_t firstEvent = 0;
243     
244     TString localFileName = "Run"; 
245     localFileName += fAODEvent->GetRunNumber(); 
246     localFileName += ".Event"; 
247     localFileName += firstEvent; 
248     localFileName += "_"; 
249     localFileName += chain->GetEntries(); //localFileName += "."; localFileName += Counter;
250     localFileName += ".AOD.tag.root";
251     
252     TString fileName;
253     
254     if(fStorage == 0) {
255         fileName = localFileName.Data();      
256         AliInfo(Form("Writing tags to local file: %s",fileName.Data()));
257     }
258     else if(fStorage == 1) {
259         TString alienLocation = "/alien";
260         alienLocation += gGrid->Pwd();
261         alienLocation += fgridpath.Data();
262         alienLocation += "/";
263         alienLocation +=  localFileName;
264         alienLocation += "?se=";
265         alienLocation += fSE.Data();
266         fileName = alienLocation.Data();
267         AliInfo(Form("Writing tags to grid file: %s",fileName.Data()));
268     }
269     
270     TFile* ftag = TFile::Open(fileName, "recreate");
271     
272     fRunTag = new AliRunTag();
273     fTreeT  = new TTree("T", "A Tree with event tags");
274     TBranch * btag = fTreeT->Branch("AliTAG", &fRunTag);
275     btag->SetCompressionLevel(9);
276     
277     // Access information from esd tag
278     fRunTagEsd = new AliRunTag();
279     fTreeTEsd->SetBranchAddress("AliTAG",&fRunTagEsd);
280     // Creating new information of aod
281     AliInfo(Form("Creating the AOD tags......."));
282     CreateTags(type);
283     ftag->cd();
284     fRunTag->Clear();
285     fTreeT->Write();
286     ftag->Close();
287 }
288
289 void AliAODTagCreator::CreateTags(const char* type)
290 {
291     // Event loop for tag creation
292     TString fturl;
293     TString fguid;
294     Int_t   oldRun = -1;
295     fChain->GetEntry(0);
296     TFile *f = fChain->GetFile();
297     TString ftempGuid = f->GetUUID().AsString();
298     // Loop over events 
299     Int_t nEvents = fChain->GetEntries();
300     Int_t ntags    = 0;
301     Int_t tagentry = 0;
302     const TClonesArray *evTagList = 0;
303
304     for (Int_t iEventNumber = 0; iEventNumber < nEvents; iEventNumber++) {
305         // Copy old tag information
306         if (iEventNumber >= ntags) {
307             fTreeTEsd->GetEntry(tagentry++);
308             fRunTag->CopyStandardContent(fRunTagEsd);
309             evTagList = fRunTagEsd->GetEventTags();
310             ntags += evTagList->GetEntries();
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         
330         fAODEvent->GetStdContent();
331         
332         // Fill the event tag from the aod informatiom
333         FillEventTag(fAODEvent, evTag);
334         // Set the event and input file references
335         evTag->SetEventId(iEventNumber+1);
336         evTag->SetGUID(fguid);
337         if(!strcmp(type,"grid")) {
338             evTag->SetMD5(0);
339             evTag->SetTURL(fturl);
340             evTag->SetSize(0);
341             }
342         else evTag->SetPath(fturl);
343
344         // Check if a new run has to be created
345         // File has changed
346         if(fguid != ftempGuid) {
347             ftempGuid = fguid;
348             fTreeT->Fill();
349             fRunTag->Clear("");
350         }
351         // Run# has changed
352         if (oldRun != (fAODEvent->GetRunNumber()))
353         {
354             oldRun = fAODEvent->GetRunNumber();
355             fTreeT->Fill();
356             fRunTag->Clear("");
357         }
358         
359         // Add the event tag
360         fRunTag->AddEventTag(*evTag);
361         delete evTag;
362         // Last event
363         if(iEventNumber+1 == fChain->GetEntries()) {
364             fTreeT->Fill();
365             fRunTag->Clear("");
366         }      
367     }//event loop
368 }
369
370
371
372 void AliAODTagCreator::FillEventTag(AliAODEvent* aod, AliEventTag* evTag)
373 {
374 // 
375 // Fill the event tag information       
376     //
377     fAODEvent = aod;
378     
379     //
380     Float_t fLowPtCut      =  1.0;
381     Float_t fHighPtCut     =  3.0;
382     Float_t fVeryHighPtCut = 10.0;
383     ////////////                            
384     Double_t partFrac[10] = {0.01, 0.01, 0.85, 0.10, 0.05, 0., 0., 0., 0., 0.};
385     
386     // Creates the tags for all the events in a given AOD file
387     Int_t   ntrack = 0;
388     Int_t   nPos = 0, nNeg = 0, nNeutr =0;
389     Int_t   nKinks = 0, nV0s = 0, nCascades = 0;
390     Int_t   nK0s = 0, nNeutrons = 0, nPi0s = 0, nGamas = 0;
391     Int_t   nProtons = 0,  nKaons = 0, nPions = 0, nMuons = 0, nElectrons = 0, nFWMuons = 0;
392     Int_t   nCh1GeV = 0, nCh3GeV = 0, nCh10GeV = 0;
393     Int_t   nMu1GeV = 0, nMu3GeV = 0, nMu10GeV = 0;
394     Int_t   nEl1GeV = 0, nEl3GeV = 0, nEl10GeV = 0;
395     Float_t maxPt =  .0, meanPt = .0, totalP =  .0;
396
397     Int_t nTracks = fAODEvent->GetNTracks();
398     // loop over vertices 
399     Int_t nVtxs = fAODEvent->GetNumberOfVertices();
400     for (Int_t nVtx = 0; nVtx < nVtxs; nVtx++) {
401         AliAODVertex *vertex = fAODEvent->GetVertex(nVtx);
402         if(vertex->GetType() == 1) nKinks    += 1;
403         if(vertex->GetType() == 2) nV0s      += 1;
404         if(vertex->GetType() == 3) nCascades += 1;
405     }
406     for (Int_t nTr = 0; nTr < nTracks; nTr++) {
407         AliAODTrack *track = fAODEvent->GetTrack(nTr);
408         
409         Double_t fPt = track->Pt();
410         if(fPt > maxPt) maxPt = fPt;
411         if(track->Charge() > 0) {
412             nPos++;
413             if(fPt > fLowPtCut)      nCh1GeV++;
414             if(fPt > fHighPtCut)     nCh3GeV++;
415             if(fPt > fVeryHighPtCut) nCh10GeV++;
416         }
417         if(track->Charge() < 0) {
418             nNeg++;
419             if(fPt > fLowPtCut)      nCh1GeV++;
420             if(fPt > fHighPtCut)     nCh3GeV++;
421             if(fPt > fVeryHighPtCut) nCh10GeV++;
422             }
423         if(track->Charge() == 0) nNeutr++;
424         //PID
425         const Double32_t *prob = track->PID();
426         Double_t rcc = 0.0;
427         for(Int_t i = 0; i < AliPID::kSPECIES; i++) rcc += prob[i]*partFrac[i];
428         if(rcc == 0.0) continue;
429         //Bayes' formula  
430         Double_t w[10];
431         for(Int_t i = 0; i < AliPID::kSPECIES; i++) w[i] = prob[i]*partFrac[i]/rcc;
432         
433         //protons
434         if ((w[4]>w[3])&&(w[4]>w[2])&&(w[4]>w[1])&&(w[4]>w[0])) nProtons++;
435         //kaons
436         if ((w[3]>w[4])&&(w[3]>w[2])&&(w[3]>w[1])&&(w[3]>w[0])) nKaons++;
437         //pions
438         if ((w[2]>w[4])&&(w[2]>w[3])&&(w[2]>w[1])&&(w[2]>w[0])) nPions++;
439         //muons 
440         if ((w[1]>w[4])&&(w[1]>w[3])&&(w[1]>w[2])&&(w[1]>w[0])) {
441             nMuons++;
442             if(fPt > fLowPtCut)      nMu1GeV++;
443             if(fPt > fHighPtCut)     nMu3GeV++;
444             if(fPt > fVeryHighPtCut) nMu10GeV++;
445         }
446         //electrons  
447         if ((w[0]>w[4])&&(w[0]>w[3])&&(w[0]>w[2])&&(w[0]>w[1])) {
448             nElectrons++;
449             if(fPt > fLowPtCut) nEl1GeV++;
450             if(fPt > fHighPtCut) nEl3GeV++;
451             if(fPt > fVeryHighPtCut) nEl10GeV++;
452         }
453         totalP += track->P();
454         meanPt += fPt;
455         ntrack++;
456         // forward muons (in the dimuon spectrometer)
457         if(track->IsMuonTrack()) nFWMuons++;   
458                                   
459     }//track loop
460     //
461     // Fill the event tags  
462     if(ntrack != 0)
463         meanPt = meanPt/ntrack;
464     
465     
466     evTag->SetNumOfTracks(nTracks);
467     evTag->SetNumOfPosTracks(nPos);
468     evTag->SetNumOfNegTracks(nNeg);
469     evTag->SetNumOfNeutrTracks(nNeutr);
470     
471     evTag->SetNumOfV0s(nV0s);
472     evTag->SetNumOfCascades(nCascades);
473     evTag->SetNumOfKinks(nKinks);
474     
475     evTag->SetNumOfProtons(nProtons);
476     evTag->SetNumOfKaons(nKaons);
477     evTag->SetNumOfPions(nPions);
478     evTag->SetNumOfMuons(nMuons);
479     evTag->SetNumOfFWMuons(nFWMuons);
480     evTag->SetNumOfElectrons(nElectrons);
481     evTag->SetNumOfPhotons(nGamas);
482     evTag->SetNumOfPi0s(nPi0s);
483     evTag->SetNumOfNeutrons(nNeutrons);
484     evTag->SetNumOfKaon0s(nK0s);
485     
486     evTag->SetNumOfChargedAbove1GeV(nCh1GeV);
487     evTag->SetNumOfChargedAbove3GeV(nCh3GeV);
488     evTag->SetNumOfChargedAbove10GeV(nCh10GeV);
489     evTag->SetNumOfMuonsAbove1GeV(nMu1GeV);
490     evTag->SetNumOfMuonsAbove3GeV(nMu3GeV);
491     evTag->SetNumOfMuonsAbove10GeV(nMu10GeV);
492     evTag->SetNumOfElectronsAbove1GeV(nEl1GeV);
493     evTag->SetNumOfElectronsAbove3GeV(nEl3GeV);
494     evTag->SetNumOfElectronsAbove10GeV(nEl10GeV);
495         
496     evTag->SetTotalMomentum(totalP);
497     evTag->SetMeanPt(meanPt);
498     evTag->SetMaxPt(maxPt);
499 }