66b28d533e663d6139f3d75159bb6e592d80ae64
[u/mrichter/AliRoot.git] / STEER / AliESDTagCreator.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 //           AliESDTagCreator 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 "AliESDEvent.h"
40 #include "AliESDVertex.h"
41 #include "AliLog.h"
42
43 #include "AliESDTagCreator.h"
44
45
46 ClassImp(AliESDTagCreator)
47
48
49 //______________________________________________________________________________
50   AliESDTagCreator::AliESDTagCreator() :
51     AliTagCreator() {
52   //==============Default constructor for a AliESDTagCreator================
53 }
54
55 //______________________________________________________________________________
56 AliESDTagCreator::~AliESDTagCreator() {
57 //================Default destructor for a AliESDTagCreator===================
58 }
59
60 //______________________________________________________________________________
61 Bool_t AliESDTagCreator::ReadGridCollection(TGridResult *fresult) {
62   // Reads the entry of the TGridResult and creates the tags
63   Int_t nEntries = fresult->GetEntries();
64
65   TString alienUrl;
66   const char *guid;
67   const char *md5;
68   const char *turl;
69   Long64_t size = -1;
70
71   Int_t counter = 0;
72   for(Int_t i = 0; i < nEntries; i++) {
73     alienUrl = fresult->GetKey(i,"turl");
74     guid = fresult->GetKey(i,"guid");
75     if(fresult->GetKey(i,"size")) size = atol (fresult->GetKey(i,"size"));
76     md5 = fresult->GetKey(i,"md5");
77     turl = fresult->GetKey(i,"turl");
78     if(md5 && !strlen(guid)) md5 = 0;
79     if(guid && !strlen(guid)) guid = 0;
80
81     TFile *f = TFile::Open(alienUrl,"READ");
82     CreateTag(f,guid,md5,turl,size,counter);
83     f->Close();
84     delete f;    
85     counter += 1;
86   }//grid result loop
87
88   return kTRUE;
89 }
90
91 //______________________________________________________________________________
92 Bool_t AliESDTagCreator::ReadLocalCollection(const char *localpath) {
93   // Checks the different subdirs of the given local path and in the
94   // case where it finds an AliESDs.root file it creates the tags
95   
96   void *dira =  gSystem->OpenDirectory(localpath);
97   Char_t fPath[256];
98   const char * dirname = 0x0;
99   const char * filename = 0x0;
100   const char * pattern = "AliESDs.root"; 
101
102   Int_t counter = 0;
103   while((dirname = gSystem->GetDirEntry(dira))) {
104     sprintf(fPath,"%s/%s",localpath,dirname);
105     void *dirb =  gSystem->OpenDirectory(fPath);
106     while((filename = gSystem->GetDirEntry(dirb))) {
107       if(strstr(filename,pattern)) {
108         TString fESDFileName;
109         fESDFileName = fPath;
110         fESDFileName += "/";
111         fESDFileName += pattern;
112         TFile *f = TFile::Open(fESDFileName,"READ");
113         CreateTag(f,fESDFileName,counter);
114         f->Close();
115         delete f;        
116         
117         counter += 1;
118       }//pattern check
119     }//child directory's entry loop
120   }//parent directory's entry loop
121
122   return kTRUE;
123 }
124
125 //______________________________________________________________________________
126 Bool_t AliESDTagCreator::ReadCAFCollection(const char *filename) {
127   // Temporary solution for CAF: Takes as an input the ascii file that
128   // lists the ESDs stored in the SE of the CAF and creates the tags.
129
130   // Open the input stream
131   ifstream in;
132   in.open(filename);
133
134   Int_t counter = 0;
135   TString esdfile;
136   // Read the input list of files and add them to the chain
137   while(in.good()) {
138     in >> esdfile;
139     if (!esdfile.Contains("root")) continue; // protection
140     TFile *f = TFile::Open(esdfile,"READ");
141     CreateTag(f,esdfile,counter);
142     f->Close();
143     delete f;    
144     
145     counter += 1;
146   }
147
148   return kTRUE;
149 }
150
151 //_____________________________________________________________________________
152 void AliESDTagCreator::CreateTag(TFile* file, const char *guid, const char *md5, const char *turl, Long64_t size, Int_t Counter) {
153   //private method that creates tag files
154   TString fguid = guid;
155   TString fmd5 = md5;
156   TString fturl = turl;
157   /////////////
158   //muon code//
159   ////////////
160   Double_t fMUONMASS = 0.105658369;
161   //Variables
162   Double_t fX,fY,fZ ;
163   Double_t fThetaX, fThetaY, fPyz, fChisquare;
164   Double_t fPxRec, fPyRec, fPzRec, fEnergy;
165   Int_t fCharge;
166   TLorentzVector fEPvector;
167
168   Float_t fZVertexCut = 10.0; 
169   Float_t fRhoVertexCut = 2.0; 
170
171   Float_t fLowPtCut = 1.0;
172   Float_t fHighPtCut = 3.0;
173   Float_t fVeryHighPtCut = 10.0;
174   ////////////
175
176   Double_t partFrac[5] = {0.01, 0.01, 0.85, 0.10, 0.05};
177
178   // Creates the tags for all the events in a given ESD file
179   Bool_t fIsSim = kTRUE;
180   Int_t ntrack;
181   Int_t nProtons, nKaons, nPions, nMuons, nElectrons;
182   Int_t nPos, nNeg, nNeutr;
183   Int_t nK0s, nNeutrons, nPi0s, nGamas;
184   Int_t nCh1GeV, nCh3GeV, nCh10GeV;
185   Int_t nMu1GeV, nMu3GeV, nMu10GeV;
186   Int_t nEl1GeV, nEl3GeV, nEl10GeV;
187   Float_t maxPt = .0, meanPt = .0, totalP = .0;
188   Int_t fVertexflag;
189   Int_t iRunNumber = 0;
190   TString fVertexName;
191
192   AliRunTag *tag = new AliRunTag();
193   AliEventTag *evTag = new AliEventTag();
194   TTree ttag("T","A Tree with event tags");
195   TBranch * btag = ttag.Branch("AliTAG", &tag);
196   btag->SetCompressionLevel(9);
197   
198   AliInfo(Form("Creating the tags......."));    
199   
200   Int_t firstEvent = 0,lastEvent = 0;
201   TTree *t = (TTree*) file->Get("esdTree");
202   AliESDEvent *esd = new AliESDEvent();
203   esd->ReadFromTree(t);
204   
205   t->GetEntry(0);
206   Int_t iInitRunNumber = esd->GetRunNumber();
207
208   Int_t iNumberOfEvents = (Int_t)t->GetEntries();
209   for (Int_t iEventNumber = 0; iEventNumber < iNumberOfEvents; iEventNumber++) {
210     ntrack = 0;
211     nPos = 0;
212     nNeg = 0;
213     nNeutr =0;
214     nK0s = 0;
215     nNeutrons = 0;
216     nPi0s = 0;
217     nGamas = 0;
218     nProtons = 0;
219     nKaons = 0;
220     nPions = 0;
221     nMuons = 0;
222     nElectrons = 0;       
223     nCh1GeV = 0;
224     nCh3GeV = 0;
225     nCh10GeV = 0;
226     nMu1GeV = 0;
227     nMu3GeV = 0;
228     nMu10GeV = 0;
229     nEl1GeV = 0;
230     nEl3GeV = 0;
231     nEl10GeV = 0;
232     maxPt = .0;
233     meanPt = .0;
234     totalP = .0;
235     fVertexflag = 1;
236     
237     t->GetEntry(iEventNumber);
238     iRunNumber = esd->GetRunNumber();
239     if(iRunNumber != iInitRunNumber) AliFatal("Inconsistency of run numbers in the AliESD!!!");
240     const AliESDVertex * vertexIn = esd->GetVertex();
241     fVertexName = vertexIn->GetName();
242     if(fVertexName == "default") fVertexflag = 0;
243
244     for (Int_t iTrackNumber = 0; iTrackNumber < esd->GetNumberOfTracks(); iTrackNumber++) {
245       AliESDtrack * esdTrack = esd->GetTrack(iTrackNumber);
246       if(esdTrack->GetLabel() != 0) fIsSim = kTRUE;
247       else if(esdTrack->GetLabel() == 0) fIsSim = kFALSE;
248       UInt_t status = esdTrack->GetStatus();
249       
250       //select only tracks with ITS refit
251       if ((status&AliESDtrack::kITSrefit)==0) continue;
252       //select only tracks with TPC refit
253       if ((status&AliESDtrack::kTPCrefit)==0) continue;
254       
255       //select only tracks with the "combined PID"
256       if ((status&AliESDtrack::kESDpid)==0) continue;
257       Double_t p[3];
258       esdTrack->GetPxPyPz(p);
259       Double_t momentum = sqrt(pow(p[0],2) + pow(p[1],2) + pow(p[2],2));
260       Double_t fPt = sqrt(pow(p[0],2) + pow(p[1],2));
261       totalP += momentum;
262       meanPt += fPt;
263       if(fPt > maxPt) maxPt = fPt;
264       
265       if(esdTrack->GetSign() > 0) {
266         nPos++;
267         if(fPt > fLowPtCut) nCh1GeV++;
268         if(fPt > fHighPtCut) nCh3GeV++;
269         if(fPt > fVeryHighPtCut) nCh10GeV++;
270       }
271       if(esdTrack->GetSign() < 0) {
272         nNeg++;
273         if(fPt > fLowPtCut) nCh1GeV++;
274         if(fPt > fHighPtCut) nCh3GeV++;
275         if(fPt > fVeryHighPtCut) nCh10GeV++;
276       }
277       if(esdTrack->GetSign() == 0) nNeutr++;
278       
279       //PID
280       Double_t prob[5];
281       esdTrack->GetESDpid(prob);
282       
283       Double_t rcc = 0.0;
284       for(Int_t i = 0; i < AliPID::kSPECIES; i++) rcc += prob[i]*partFrac[i];
285       if(rcc == 0.0) continue;
286       //Bayes' formula
287       Double_t w[5];
288       for(Int_t i = 0; i < AliPID::kSPECIES; i++) w[i] = prob[i]*partFrac[i]/rcc;
289       
290       //protons
291       if ((w[4]>w[3])&&(w[4]>w[2])&&(w[4]>w[1])&&(w[4]>w[0])) nProtons++;
292       //kaons
293       if ((w[3]>w[4])&&(w[3]>w[2])&&(w[3]>w[1])&&(w[3]>w[0])) nKaons++;
294       //pions
295       if ((w[2]>w[4])&&(w[2]>w[3])&&(w[2]>w[1])&&(w[2]>w[0])) nPions++; 
296       //electrons
297       if ((w[0]>w[4])&&(w[0]>w[3])&&(w[0]>w[2])&&(w[0]>w[1])) {
298         nElectrons++;
299         if(fPt > fLowPtCut) nEl1GeV++;
300         if(fPt > fHighPtCut) nEl3GeV++;
301         if(fPt > fVeryHighPtCut) nEl10GeV++;
302       }   
303       ntrack++;
304     }//esd track loop
305     
306     /////////////
307     //muon code//
308     ////////////
309     Int_t nMuonTracks = esd->GetNumberOfMuonTracks();
310     // loop over all reconstructed tracks (also first track of combination)
311     for (Int_t iTrack = 0; iTrack <  nMuonTracks;  iTrack++) {
312       AliESDMuonTrack* muonTrack = esd->GetMuonTrack(iTrack);
313       if (muonTrack == 0x0) continue;
314       
315       // Coordinates at vertex
316       fZ = muonTrack->GetZ(); 
317       fY = muonTrack->GetBendingCoor();
318       fX = muonTrack->GetNonBendingCoor(); 
319       
320       fThetaX = muonTrack->GetThetaX();
321       fThetaY = muonTrack->GetThetaY();
322       
323       fPyz = 1./TMath::Abs(muonTrack->GetInverseBendingMomentum());
324       fPzRec = - fPyz / TMath::Sqrt(1.0 + TMath::Tan(fThetaY)*TMath::Tan(fThetaY));
325       fPxRec = fPzRec * TMath::Tan(fThetaX);
326       fPyRec = fPzRec * TMath::Tan(fThetaY);
327       fCharge = Int_t(TMath::Sign(1.,muonTrack->GetInverseBendingMomentum()));
328       
329       //ChiSquare of the track if needed
330       fChisquare = muonTrack->GetChi2()/(2.0 * muonTrack->GetNHit() - 5);
331       fEnergy = TMath::Sqrt(fMUONMASS * fMUONMASS + fPxRec * fPxRec + fPyRec * fPyRec + fPzRec * fPzRec);
332       fEPvector.SetPxPyPzE(fPxRec, fPyRec, fPzRec, fEnergy);
333       
334       // total number of muons inside a vertex cut 
335       if((TMath::Abs(fZ)<fZVertexCut) && (TMath::Sqrt(fY*fY+fX*fX)<fRhoVertexCut)) {
336         nMuons++;
337         if(fEPvector.Pt() > fLowPtCut) {
338           nMu1GeV++; 
339           if(fEPvector.Pt() > fHighPtCut) {
340             nMu3GeV++; 
341             if (fEPvector.Pt() > fVeryHighPtCut) {
342               nMu10GeV++;
343             }
344           }
345         }
346       }
347     }//muon track loop
348     
349     // Fill the event tags 
350     if(ntrack != 0) meanPt = meanPt/ntrack;
351     
352     evTag->SetEventId(iEventNumber+1);
353     evTag->SetGUID(fguid);
354     evTag->SetMD5(fmd5);
355     evTag->SetTURL(fturl);
356     evTag->SetSize(size);
357     evTag->SetVertexX(vertexIn->GetXv());
358     evTag->SetVertexY(vertexIn->GetYv());
359     evTag->SetVertexZ(vertexIn->GetZv());
360     evTag->SetVertexZError(vertexIn->GetZRes());
361     evTag->SetVertexFlag(fVertexflag);
362     
363     evTag->SetT0VertexZ(esd->GetT0zVertex());
364     
365     evTag->SetTriggerMask(esd->GetTriggerMask());
366     evTag->SetTriggerCluster(esd->GetTriggerCluster());
367     
368     evTag->SetZDCNeutron1Energy(esd->GetZDCN1Energy());
369     evTag->SetZDCProton1Energy(esd->GetZDCP1Energy());
370     evTag->SetZDCEMEnergy(esd->GetZDCEMEnergy());
371     evTag->SetZDCNeutron1Energy(esd->GetZDCN2Energy());
372     evTag->SetZDCProton1Energy(esd->GetZDCP2Energy());
373     evTag->SetNumOfParticipants(esd->GetZDCParticipants());
374     
375     
376     evTag->SetNumOfTracks(esd->GetNumberOfTracks());
377     evTag->SetNumOfPosTracks(nPos);
378     evTag->SetNumOfNegTracks(nNeg);
379     evTag->SetNumOfNeutrTracks(nNeutr);
380     
381     evTag->SetNumOfV0s(esd->GetNumberOfV0s());
382     evTag->SetNumOfCascades(esd->GetNumberOfCascades());
383     evTag->SetNumOfKinks(esd->GetNumberOfKinks());
384     evTag->SetNumOfPMDTracks(esd->GetNumberOfPmdTracks());
385     
386     evTag->SetNumOfProtons(nProtons);
387     evTag->SetNumOfKaons(nKaons);
388     evTag->SetNumOfPions(nPions);
389     evTag->SetNumOfMuons(nMuons);
390     evTag->SetNumOfElectrons(nElectrons);
391     evTag->SetNumOfPhotons(nGamas);
392     evTag->SetNumOfPi0s(nPi0s);
393     evTag->SetNumOfNeutrons(nNeutrons);
394     evTag->SetNumOfKaon0s(nK0s);
395     
396     evTag->SetNumOfChargedAbove1GeV(nCh1GeV);
397     evTag->SetNumOfChargedAbove3GeV(nCh3GeV);
398     evTag->SetNumOfChargedAbove10GeV(nCh10GeV);
399     evTag->SetNumOfMuonsAbove1GeV(nMu1GeV);
400     evTag->SetNumOfMuonsAbove3GeV(nMu3GeV);
401     evTag->SetNumOfMuonsAbove10GeV(nMu10GeV);
402     evTag->SetNumOfElectronsAbove1GeV(nEl1GeV);
403     evTag->SetNumOfElectronsAbove3GeV(nEl3GeV);
404     evTag->SetNumOfElectronsAbove10GeV(nEl10GeV);
405     
406     evTag->SetNumOfPHOSClusters(esd->GetNumberOfPHOSClusters());
407     evTag->SetNumOfEMCALClusters(esd->GetNumberOfEMCALClusters());
408     
409     evTag->SetTotalMomentum(totalP);
410     evTag->SetMeanPt(meanPt);
411     evTag->SetMaxPt(maxPt);
412     
413     tag->SetRunId(iInitRunNumber);
414     if(fIsSim) tag->SetDataType(0);
415     else tag->SetDataType(1);
416     tag->AddEventTag(*evTag);
417   }//event loop
418   lastEvent = iNumberOfEvents;
419   
420   t->Delete("");
421   
422   ttag.Fill();
423   tag->Clear();
424
425   TString localFileName = "Run"; localFileName += tag->GetRunId(); 
426   localFileName += ".Event"; localFileName += firstEvent; localFileName += "_"; localFileName += lastEvent; localFileName += "."; localFileName += Counter;
427   localFileName += ".ESD.tag.root";
428
429   TString fileName;
430   
431   if(fStorage == 0) {
432     fileName = localFileName.Data();      
433     AliInfo(Form("Writing tags to local file: %s",fileName.Data()));
434   }
435   else if(fStorage == 1) {
436     TString alienLocation = "/alien";
437     alienLocation += gGrid->Pwd();
438     alienLocation += fgridpath.Data();
439     alienLocation += "/";
440     alienLocation +=  localFileName;
441     alienLocation += "?se=";
442     alienLocation += fSE.Data();
443     fileName = alienLocation.Data();
444     AliInfo(Form("Writing tags to grid file: %s",fileName.Data()));
445   }
446
447   TFile* ftag = TFile::Open(fileName, "recreate");
448   ftag->cd();
449   ttag.Write();
450   ftag->Close();
451
452   delete ftag;
453   delete esd;
454
455   delete tag;
456 }
457
458 //_____________________________________________________________________________
459 void AliESDTagCreator::CreateTag(TFile* file, const char *filepath, Int_t Counter) {
460   //private method that creates tag files
461   /////////////
462   //muon code//
463   ////////////
464   Double_t fMUONMASS = 0.105658369;
465   //Variables
466   Double_t fX,fY,fZ ;
467   Double_t fThetaX, fThetaY, fPyz, fChisquare;
468   Double_t fPxRec, fPyRec, fPzRec, fEnergy;
469   Int_t fCharge;
470   TLorentzVector fEPvector;
471
472   Float_t fZVertexCut = 10.0; 
473   Float_t fRhoVertexCut = 2.0; 
474
475   Float_t fLowPtCut = 1.0;
476   Float_t fHighPtCut = 3.0;
477   Float_t fVeryHighPtCut = 10.0;
478   ////////////
479
480   Double_t partFrac[5] = {0.01, 0.01, 0.85, 0.10, 0.05};
481
482   // Creates the tags for all the events in a given ESD file
483   Bool_t fIsSim = kTRUE;
484   Int_t ntrack;
485   Int_t nProtons, nKaons, nPions, nMuons, nElectrons;
486   Int_t nPos, nNeg, nNeutr;
487   Int_t nK0s, nNeutrons, nPi0s, nGamas;
488   Int_t nCh1GeV, nCh3GeV, nCh10GeV;
489   Int_t nMu1GeV, nMu3GeV, nMu10GeV;
490   Int_t nEl1GeV, nEl3GeV, nEl10GeV;
491   Float_t maxPt = .0, meanPt = .0, totalP = .0;
492   Int_t fVertexflag;
493   Int_t iRunNumber = 0;
494   TString fVertexName;
495
496   AliRunTag *tag = new AliRunTag();
497   AliEventTag *evTag = new AliEventTag();
498   TTree ttag("T","A Tree with event tags");
499   TBranch * btag = ttag.Branch("AliTAG", &tag);
500   btag->SetCompressionLevel(9);
501   
502   AliInfo(Form("Creating the tags......."));    
503   
504   Int_t firstEvent = 0,lastEvent = 0;
505   
506   TTree *t = (TTree*) file->Get("esdTree");
507   AliESDEvent *esd = new AliESDEvent();
508   esd->ReadFromTree(t);
509   
510   t->GetEntry(0);
511   Int_t iInitRunNumber = esd->GetRunNumber();
512
513   Int_t iNumberOfEvents = (Int_t)t->GetEntries();
514   for (Int_t iEventNumber = 0; iEventNumber < iNumberOfEvents; iEventNumber++) {
515     ntrack = 0;
516     nPos = 0;
517     nNeg = 0;
518     nNeutr =0;
519     nK0s = 0;
520     nNeutrons = 0;
521     nPi0s = 0;
522     nGamas = 0;
523     nProtons = 0;
524     nKaons = 0;
525     nPions = 0;
526     nMuons = 0;
527     nElectrons = 0;       
528     nCh1GeV = 0;
529     nCh3GeV = 0;
530     nCh10GeV = 0;
531     nMu1GeV = 0;
532     nMu3GeV = 0;
533     nMu10GeV = 0;
534     nEl1GeV = 0;
535     nEl3GeV = 0;
536     nEl10GeV = 0;
537     maxPt = .0;
538     meanPt = .0;
539     totalP = .0;
540     fVertexflag = 1;
541     
542     t->GetEntry(iEventNumber);
543     iRunNumber = esd->GetRunNumber();
544     if(iRunNumber != iInitRunNumber) AliFatal("Inconsistency of run numbers in the AliESD!!!");
545     const AliESDVertex * vertexIn = esd->GetVertex();
546     fVertexName = vertexIn->GetName();
547     if(fVertexName == "default") fVertexflag = 0;
548
549     for (Int_t iTrackNumber = 0; iTrackNumber < esd->GetNumberOfTracks(); iTrackNumber++) {
550       AliESDtrack * esdTrack = esd->GetTrack(iTrackNumber);
551       if(esdTrack->GetLabel() != 0) fIsSim = kTRUE;
552       else if(esdTrack->GetLabel() == 0) fIsSim = kFALSE;
553       UInt_t status = esdTrack->GetStatus();
554       
555       //select only tracks with ITS refit
556       if ((status&AliESDtrack::kITSrefit)==0) continue;
557       //select only tracks with TPC refit
558       if ((status&AliESDtrack::kTPCrefit)==0) continue;
559       
560       //select only tracks with the "combined PID"
561       if ((status&AliESDtrack::kESDpid)==0) continue;
562       Double_t p[3];
563       esdTrack->GetPxPyPz(p);
564       Double_t momentum = sqrt(pow(p[0],2) + pow(p[1],2) + pow(p[2],2));
565       Double_t fPt = sqrt(pow(p[0],2) + pow(p[1],2));
566       totalP += momentum;
567       meanPt += fPt;
568       if(fPt > maxPt) maxPt = fPt;
569       
570       if(esdTrack->GetSign() > 0) {
571         nPos++;
572         if(fPt > fLowPtCut) nCh1GeV++;
573         if(fPt > fHighPtCut) nCh3GeV++;
574         if(fPt > fVeryHighPtCut) nCh10GeV++;
575       }
576       if(esdTrack->GetSign() < 0) {
577         nNeg++;
578         if(fPt > fLowPtCut) nCh1GeV++;
579         if(fPt > fHighPtCut) nCh3GeV++;
580         if(fPt > fVeryHighPtCut) nCh10GeV++;
581       }
582       if(esdTrack->GetSign() == 0) nNeutr++;
583       
584       //PID
585       Double_t prob[5];
586       esdTrack->GetESDpid(prob);
587       
588       Double_t rcc = 0.0;
589       for(Int_t i = 0; i < AliPID::kSPECIES; i++) rcc += prob[i]*partFrac[i];
590       if(rcc == 0.0) continue;
591       //Bayes' formula
592       Double_t w[5];
593       for(Int_t i = 0; i < AliPID::kSPECIES; i++) w[i] = prob[i]*partFrac[i]/rcc;
594       
595       //protons
596       if ((w[4]>w[3])&&(w[4]>w[2])&&(w[4]>w[1])&&(w[4]>w[0])) nProtons++;
597       //kaons
598       if ((w[3]>w[4])&&(w[3]>w[2])&&(w[3]>w[1])&&(w[3]>w[0])) nKaons++;
599       //pions
600       if ((w[2]>w[4])&&(w[2]>w[3])&&(w[2]>w[1])&&(w[2]>w[0])) nPions++; 
601       //electrons
602       if ((w[0]>w[4])&&(w[0]>w[3])&&(w[0]>w[2])&&(w[0]>w[1])) {
603         nElectrons++;
604         if(fPt > fLowPtCut) nEl1GeV++;
605         if(fPt > fHighPtCut) nEl3GeV++;
606         if(fPt > fVeryHighPtCut) nEl10GeV++;
607       }   
608       ntrack++;
609     }//esd track loop
610     
611     /////////////
612     //muon code//
613     ////////////
614     Int_t nMuonTracks = esd->GetNumberOfMuonTracks();
615     // loop over all reconstructed tracks (also first track of combination)
616     for (Int_t iTrack = 0; iTrack <  nMuonTracks;  iTrack++) {
617       AliESDMuonTrack* muonTrack = esd->GetMuonTrack(iTrack);
618       if (muonTrack == 0x0) continue;
619       
620       // Coordinates at vertex
621       fZ = muonTrack->GetZ(); 
622       fY = muonTrack->GetBendingCoor();
623       fX = muonTrack->GetNonBendingCoor(); 
624       
625       fThetaX = muonTrack->GetThetaX();
626       fThetaY = muonTrack->GetThetaY();
627       
628       fPyz = 1./TMath::Abs(muonTrack->GetInverseBendingMomentum());
629       fPzRec = - fPyz / TMath::Sqrt(1.0 + TMath::Tan(fThetaY)*TMath::Tan(fThetaY));
630       fPxRec = fPzRec * TMath::Tan(fThetaX);
631       fPyRec = fPzRec * TMath::Tan(fThetaY);
632       fCharge = Int_t(TMath::Sign(1.,muonTrack->GetInverseBendingMomentum()));
633       
634       //ChiSquare of the track if needed
635       fChisquare = muonTrack->GetChi2()/(2.0 * muonTrack->GetNHit() - 5);
636       fEnergy = TMath::Sqrt(fMUONMASS * fMUONMASS + fPxRec * fPxRec + fPyRec * fPyRec + fPzRec * fPzRec);
637       fEPvector.SetPxPyPzE(fPxRec, fPyRec, fPzRec, fEnergy);
638       
639       // total number of muons inside a vertex cut 
640       if((TMath::Abs(fZ)<fZVertexCut) && (TMath::Sqrt(fY*fY+fX*fX)<fRhoVertexCut)) {
641         nMuons++;
642         if(fEPvector.Pt() > fLowPtCut) {
643           nMu1GeV++; 
644           if(fEPvector.Pt() > fHighPtCut) {
645             nMu3GeV++; 
646             if (fEPvector.Pt() > fVeryHighPtCut) {
647               nMu10GeV++;
648             }
649           }
650         }
651       }
652     }//muon track loop
653     
654     // Fill the event tags 
655     if(ntrack != 0) meanPt = meanPt/ntrack;
656     
657     evTag->SetEventId(iEventNumber+1);
658     evTag->SetPath(filepath);
659  
660     evTag->SetVertexX(vertexIn->GetXv());
661     evTag->SetVertexY(vertexIn->GetYv());
662     evTag->SetVertexZ(vertexIn->GetZv());
663     evTag->SetVertexZError(vertexIn->GetZRes());
664     evTag->SetVertexFlag(fVertexflag);
665     
666     evTag->SetT0VertexZ(esd->GetT0zVertex());
667     
668     evTag->SetTriggerMask(esd->GetTriggerMask());
669     evTag->SetTriggerCluster(esd->GetTriggerCluster());
670     
671     evTag->SetZDCNeutron1Energy(esd->GetZDCN1Energy());
672     evTag->SetZDCProton1Energy(esd->GetZDCP1Energy());
673     evTag->SetZDCEMEnergy(esd->GetZDCEMEnergy());
674     evTag->SetZDCNeutron1Energy(esd->GetZDCN2Energy());
675     evTag->SetZDCProton1Energy(esd->GetZDCP2Energy());
676     evTag->SetNumOfParticipants(esd->GetZDCParticipants());
677     
678     
679     evTag->SetNumOfTracks(esd->GetNumberOfTracks());
680     evTag->SetNumOfPosTracks(nPos);
681     evTag->SetNumOfNegTracks(nNeg);
682     evTag->SetNumOfNeutrTracks(nNeutr);
683     
684     evTag->SetNumOfV0s(esd->GetNumberOfV0s());
685     evTag->SetNumOfCascades(esd->GetNumberOfCascades());
686     evTag->SetNumOfKinks(esd->GetNumberOfKinks());
687     evTag->SetNumOfPMDTracks(esd->GetNumberOfPmdTracks());
688     
689     evTag->SetNumOfProtons(nProtons);
690     evTag->SetNumOfKaons(nKaons);
691     evTag->SetNumOfPions(nPions);
692     evTag->SetNumOfMuons(nMuons);
693     evTag->SetNumOfElectrons(nElectrons);
694     evTag->SetNumOfPhotons(nGamas);
695     evTag->SetNumOfPi0s(nPi0s);
696     evTag->SetNumOfNeutrons(nNeutrons);
697     evTag->SetNumOfKaon0s(nK0s);
698     
699     evTag->SetNumOfChargedAbove1GeV(nCh1GeV);
700     evTag->SetNumOfChargedAbove3GeV(nCh3GeV);
701     evTag->SetNumOfChargedAbove10GeV(nCh10GeV);
702     evTag->SetNumOfMuonsAbove1GeV(nMu1GeV);
703     evTag->SetNumOfMuonsAbove3GeV(nMu3GeV);
704     evTag->SetNumOfMuonsAbove10GeV(nMu10GeV);
705     evTag->SetNumOfElectronsAbove1GeV(nEl1GeV);
706     evTag->SetNumOfElectronsAbove3GeV(nEl3GeV);
707     evTag->SetNumOfElectronsAbove10GeV(nEl10GeV);
708     
709     evTag->SetNumOfPHOSClusters(esd->GetNumberOfPHOSClusters());
710     evTag->SetNumOfEMCALClusters(esd->GetNumberOfEMCALClusters());
711     
712     evTag->SetTotalMomentum(totalP);
713     evTag->SetMeanPt(meanPt);
714     evTag->SetMaxPt(maxPt);
715     
716     tag->SetRunId(iInitRunNumber);
717     if(fIsSim) tag->SetDataType(0);
718     else tag->SetDataType(1);
719     tag->AddEventTag(*evTag);
720   }//event loop
721   lastEvent = iNumberOfEvents;
722   
723   t->Delete("");
724   
725   ttag.Fill();
726   tag->Clear();
727
728   TString localFileName = "Run"; localFileName += tag->GetRunId(); 
729   localFileName += ".Event"; localFileName += firstEvent; localFileName += "_"; localFileName += lastEvent; localFileName += "."; localFileName += Counter;
730   localFileName += ".ESD.tag.root";
731
732   TString fileName;
733   
734   if(fStorage == 0) {
735     fileName = localFileName.Data();      
736     AliInfo(Form("Writing tags to local file: %s",fileName.Data()));
737   }
738   else if(fStorage == 1) {
739     TString alienLocation = "/alien";
740     alienLocation += gGrid->Pwd();
741     alienLocation += fgridpath.Data();
742     alienLocation += "/";
743     alienLocation +=  localFileName;
744     alienLocation += "?se=";
745     alienLocation += fSE.Data();
746     fileName = alienLocation.Data();
747     AliInfo(Form("Writing tags to grid file: %s",fileName.Data()));
748   }
749
750   TFile* ftag = TFile::Open(fileName, "recreate");
751   ftag->cd();
752   ttag.Write();
753   ftag->Close();
754
755   delete ftag;
756   delete esd;
757
758   delete tag;
759 }
760
761 //_____________________________________________________________________________
762 void AliESDTagCreator::CreateESDTags(Int_t fFirstEvent, Int_t fLastEvent) {
763   //GRP
764   Float_t lhcLuminosity = 0.0;
765   TString lhcState = "test";
766   UInt_t detectorMask = 0;
767
768   /////////////
769   //muon code//
770   ////////////
771   Double_t fMUONMASS = 0.105658369;
772   //Variables
773   Double_t fX,fY,fZ ;
774   Double_t fThetaX, fThetaY, fPyz, fChisquare;
775   Double_t fPxRec,fPyRec, fPzRec, fEnergy;
776   Int_t fCharge;
777   TLorentzVector fEPvector;
778
779   Float_t fZVertexCut = 10.0; 
780   Float_t fRhoVertexCut = 2.0; 
781
782   Float_t fLowPtCut = 1.0;
783   Float_t fHighPtCut = 3.0;
784   Float_t fVeryHighPtCut = 10.0;
785   ////////////
786
787   Double_t partFrac[5] = {0.01, 0.01, 0.85, 0.10, 0.05};
788
789   // Creates the tags for all the events in a given ESD file
790   Int_t ntrack;
791   Int_t nProtons, nKaons, nPions, nMuons, nElectrons;
792   Int_t nPos, nNeg, nNeutr;
793   Int_t nK0s, nNeutrons, nPi0s, nGamas;
794   Int_t nCh1GeV, nCh3GeV, nCh10GeV;
795   Int_t nMu1GeV, nMu3GeV, nMu10GeV;
796   Int_t nEl1GeV, nEl3GeV, nEl10GeV;
797   Float_t maxPt = .0, meanPt = .0, totalP = .0;
798   Int_t fVertexflag;
799   Int_t iRunNumber = 0;
800   TString fVertexName("default");
801
802   AliRunTag *tag = new AliRunTag();
803   AliEventTag *evTag = new AliEventTag();
804   TTree ttag("T","A Tree with event tags");
805   TBranch * btag = ttag.Branch("AliTAG", &tag);
806   btag->SetCompressionLevel(9);
807   
808   AliInfo(Form("Creating the tags......."));    
809
810   TFile *file = TFile::Open("AliESDs.root");
811   if (!file || !file->IsOpen()) {
812     AliError(Form("opening failed"));
813     delete file;
814     return ;
815   }  
816   Int_t lastEvent = 0;
817   TTree *b = (TTree*) file->Get("esdTree");
818   AliESDEvent *esd = new AliESDEvent();
819   esd->ReadFromTree(b);
820
821   b->GetEntry(fFirstEvent);
822   Int_t iInitRunNumber = esd->GetRunNumber();
823   
824   Int_t iNumberOfEvents = (Int_t)b->GetEntries();
825   if(fLastEvent != -1) iNumberOfEvents = fLastEvent + 1;
826   for (Int_t iEventNumber = fFirstEvent; iEventNumber < iNumberOfEvents; iEventNumber++) {
827     ntrack = 0;
828     nPos = 0;
829     nNeg = 0;
830     nNeutr =0;
831     nK0s = 0;
832     nNeutrons = 0;
833     nPi0s = 0;
834     nGamas = 0;
835     nProtons = 0;
836     nKaons = 0;
837     nPions = 0;
838     nMuons = 0;
839     nElectrons = 0;       
840     nCh1GeV = 0;
841     nCh3GeV = 0;
842     nCh10GeV = 0;
843     nMu1GeV = 0;
844     nMu3GeV = 0;
845     nMu10GeV = 0;
846     nEl1GeV = 0;
847     nEl3GeV = 0;
848     nEl10GeV = 0;
849     maxPt = .0;
850     meanPt = .0;
851     totalP = .0;
852     fVertexflag = 0;
853
854     b->GetEntry(iEventNumber);
855     iRunNumber = esd->GetRunNumber();
856     if(iRunNumber != iInitRunNumber) AliFatal("Inconsistency of run numbers in the AliESD!!!");
857     const AliESDVertex * vertexIn = esd->GetVertex();
858     if (!vertexIn) AliError("ESD has not defined vertex.");
859     if (vertexIn) fVertexName = vertexIn->GetName();
860     if(fVertexName != "default") fVertexflag = 1;
861     for (Int_t iTrackNumber = 0; iTrackNumber < esd->GetNumberOfTracks(); iTrackNumber++) {
862       AliESDtrack * esdTrack = esd->GetTrack(iTrackNumber);
863       UInt_t status = esdTrack->GetStatus();
864       
865       //select only tracks with ITS refit
866       if ((status&AliESDtrack::kITSrefit)==0) continue;
867       //select only tracks with TPC refit
868       if ((status&AliESDtrack::kTPCrefit)==0) continue;
869       
870       //select only tracks with the "combined PID"
871       if ((status&AliESDtrack::kESDpid)==0) continue;
872       Double_t p[3];
873       esdTrack->GetPxPyPz(p);
874       Double_t momentum = sqrt(pow(p[0],2) + pow(p[1],2) + pow(p[2],2));
875       Double_t fPt = sqrt(pow(p[0],2) + pow(p[1],2));
876       totalP += momentum;
877       meanPt += fPt;
878       if(fPt > maxPt) maxPt = fPt;
879       
880       if(esdTrack->GetSign() > 0) {
881         nPos++;
882         if(fPt > fLowPtCut) nCh1GeV++;
883         if(fPt > fHighPtCut) nCh3GeV++;
884         if(fPt > fVeryHighPtCut) nCh10GeV++;
885       }
886       if(esdTrack->GetSign() < 0) {
887         nNeg++;
888         if(fPt > fLowPtCut) nCh1GeV++;
889         if(fPt > fHighPtCut) nCh3GeV++;
890         if(fPt > fVeryHighPtCut) nCh10GeV++;
891       }
892       if(esdTrack->GetSign() == 0) nNeutr++;
893       
894       //PID
895       Double_t prob[5];
896       esdTrack->GetESDpid(prob);
897       
898       Double_t rcc = 0.0;
899       for(Int_t i = 0; i < AliPID::kSPECIES; i++) rcc += prob[i]*partFrac[i];
900       if(rcc == 0.0) continue;
901       //Bayes' formula
902       Double_t w[5];
903       for(Int_t i = 0; i < AliPID::kSPECIES; i++) w[i] = prob[i]*partFrac[i]/rcc;
904       
905       //protons
906       if ((w[4]>w[3])&&(w[4]>w[2])&&(w[4]>w[1])&&(w[4]>w[0])) nProtons++;
907       //kaons
908       if ((w[3]>w[4])&&(w[3]>w[2])&&(w[3]>w[1])&&(w[3]>w[0])) nKaons++;
909       //pions
910       if ((w[2]>w[4])&&(w[2]>w[3])&&(w[2]>w[1])&&(w[2]>w[0])) nPions++; 
911       //electrons
912       if ((w[0]>w[4])&&(w[0]>w[3])&&(w[0]>w[2])&&(w[0]>w[1])) {
913         nElectrons++;
914         if(fPt > fLowPtCut) nEl1GeV++;
915         if(fPt > fHighPtCut) nEl3GeV++;
916         if(fPt > fVeryHighPtCut) nEl10GeV++;
917       }   
918       ntrack++;
919     }//track loop
920     
921     /////////////
922     //muon code//
923     ////////////
924     Int_t nMuonTracks = esd->GetNumberOfMuonTracks();
925     // loop over all reconstructed tracks (also first track of combination)
926     for (Int_t iTrack = 0; iTrack <  nMuonTracks;  iTrack++) {
927       AliESDMuonTrack* muonTrack = esd->GetMuonTrack(iTrack);
928       if (muonTrack == 0x0) continue;
929       
930       // Coordinates at vertex
931       fZ = muonTrack->GetZ(); 
932       fY = muonTrack->GetBendingCoor();
933       fX = muonTrack->GetNonBendingCoor(); 
934       
935       fThetaX = muonTrack->GetThetaX();
936       fThetaY = muonTrack->GetThetaY();
937       
938       fPyz = 1./TMath::Abs(muonTrack->GetInverseBendingMomentum());
939       fPzRec = - fPyz / TMath::Sqrt(1.0 + TMath::Tan(fThetaY)*TMath::Tan(fThetaY));
940       fPxRec = fPzRec * TMath::Tan(fThetaX);
941       fPyRec = fPzRec * TMath::Tan(fThetaY);
942       fCharge = Int_t(TMath::Sign(1.,muonTrack->GetInverseBendingMomentum()));
943       
944       //ChiSquare of the track if needed
945       fChisquare = muonTrack->GetChi2()/(2.0 * muonTrack->GetNHit() - 5);
946       fEnergy = TMath::Sqrt(fMUONMASS * fMUONMASS + fPxRec * fPxRec + fPyRec * fPyRec + fPzRec * fPzRec);
947       fEPvector.SetPxPyPzE(fPxRec, fPyRec, fPzRec, fEnergy);
948       
949       // total number of muons inside a vertex cut 
950       if((TMath::Abs(fZ)<fZVertexCut) && (TMath::Sqrt(fY*fY+fX*fX)<fRhoVertexCut)) {
951         nMuons++;
952         if(fEPvector.Pt() > fLowPtCut) {
953           nMu1GeV++; 
954           if(fEPvector.Pt() > fHighPtCut) {
955             nMu3GeV++; 
956             if (fEPvector.Pt() > fVeryHighPtCut) {
957               nMu10GeV++;
958             }
959           }
960         }
961       }
962     }//muon track loop
963     
964     // Fill the event tags 
965     if(ntrack != 0)
966       meanPt = meanPt/ntrack;
967     
968     evTag->SetEventId(iEventNumber+1);
969     if (vertexIn) {
970       evTag->SetVertexX(vertexIn->GetXv());
971       evTag->SetVertexY(vertexIn->GetYv());
972       evTag->SetVertexZ(vertexIn->GetZv());
973       evTag->SetVertexZError(vertexIn->GetZRes());
974     }  
975     evTag->SetVertexFlag(fVertexflag);
976
977     evTag->SetT0VertexZ(esd->GetT0zVertex());
978     
979     evTag->SetTriggerMask(esd->GetTriggerMask());
980     evTag->SetTriggerCluster(esd->GetTriggerCluster());
981     
982     evTag->SetZDCNeutron1Energy(esd->GetZDCN1Energy());
983     evTag->SetZDCProton1Energy(esd->GetZDCP1Energy());
984     evTag->SetZDCNeutron2Energy(esd->GetZDCN2Energy());
985     evTag->SetZDCProton2Energy(esd->GetZDCP2Energy());
986     evTag->SetZDCEMEnergy(esd->GetZDCEMEnergy());
987     evTag->SetNumOfParticipants(esd->GetZDCParticipants());
988     
989     
990     evTag->SetNumOfTracks(esd->GetNumberOfTracks());
991     evTag->SetNumOfPosTracks(nPos);
992     evTag->SetNumOfNegTracks(nNeg);
993     evTag->SetNumOfNeutrTracks(nNeutr);
994     
995     evTag->SetNumOfV0s(esd->GetNumberOfV0s());
996     evTag->SetNumOfCascades(esd->GetNumberOfCascades());
997     evTag->SetNumOfKinks(esd->GetNumberOfKinks());
998     evTag->SetNumOfPMDTracks(esd->GetNumberOfPmdTracks());
999     
1000     evTag->SetNumOfProtons(nProtons);
1001     evTag->SetNumOfKaons(nKaons);
1002     evTag->SetNumOfPions(nPions);
1003     evTag->SetNumOfMuons(nMuons);
1004     evTag->SetNumOfElectrons(nElectrons);
1005     evTag->SetNumOfPhotons(nGamas);
1006     evTag->SetNumOfPi0s(nPi0s);
1007     evTag->SetNumOfNeutrons(nNeutrons);
1008     evTag->SetNumOfKaon0s(nK0s);
1009     
1010     evTag->SetNumOfChargedAbove1GeV(nCh1GeV);
1011     evTag->SetNumOfChargedAbove3GeV(nCh3GeV);
1012     evTag->SetNumOfChargedAbove10GeV(nCh10GeV);
1013     evTag->SetNumOfMuonsAbove1GeV(nMu1GeV);
1014     evTag->SetNumOfMuonsAbove3GeV(nMu3GeV);
1015     evTag->SetNumOfMuonsAbove10GeV(nMu10GeV);
1016     evTag->SetNumOfElectronsAbove1GeV(nEl1GeV);
1017     evTag->SetNumOfElectronsAbove3GeV(nEl3GeV);
1018     evTag->SetNumOfElectronsAbove10GeV(nEl10GeV);
1019     
1020     evTag->SetNumOfPHOSClusters(esd->GetNumberOfPHOSClusters());
1021     evTag->SetNumOfEMCALClusters(esd->GetNumberOfEMCALClusters());
1022     
1023     evTag->SetTotalMomentum(totalP);
1024     evTag->SetMeanPt(meanPt);
1025     evTag->SetMaxPt(maxPt);
1026     
1027     tag->SetLHCTag(lhcLuminosity,lhcState);
1028     tag->SetDetectorTag(detectorMask);
1029
1030     tag->SetRunId(iInitRunNumber);
1031     tag->AddEventTag(*evTag);
1032   }
1033   if(fLastEvent == -1) lastEvent = (Int_t)b->GetEntries();
1034   else lastEvent = fLastEvent;
1035         
1036   ttag.Fill();
1037   tag->Clear();
1038
1039   char fileName[256];
1040   sprintf(fileName, "Run%d.Event%d_%d.ESD.tag.root", 
1041           tag->GetRunId(),fFirstEvent,lastEvent );
1042   AliInfo(Form("writing tags to file %s", fileName));
1043   AliDebug(1, Form("writing tags to file %s", fileName));
1044  
1045   TFile* ftag = TFile::Open(fileName, "recreate");
1046   ftag->cd();
1047   ttag.Write();
1048   ftag->Close();
1049   file->cd();
1050   delete tag;
1051   delete evTag;
1052 }
1053
1054 //__________________________________________________________________________
1055 /*void AliTagCreator::CreateAODTags(Int_t fFirstEvent, Int_t fLastEvent) {
1056   //creates tag files for AODs
1057   
1058   Float_t fLowPtCut = 1.0;
1059   Float_t fHighPtCut = 3.0;
1060   Float_t fVeryHighPtCut = 10.0;
1061   ////////////
1062
1063   Double_t partFrac[10] = {0.01, 0.01, 0.85, 0.10, 0.05, 0., 0., 0., 0., 0.};
1064
1065   // Creates the tags for all the events in a given ESD file
1066   Int_t ntrack;
1067   Int_t nProtons, nKaons, nPions, nMuons, nElectrons;
1068   Int_t nPos, nNeg, nNeutr;
1069   Int_t nKinks, nV0s, nCascades;
1070   Int_t nK0s, nNeutrons, nPi0s, nGamas;
1071   Int_t nCh1GeV, nCh3GeV, nCh10GeV;
1072   Int_t nMu1GeV, nMu3GeV, nMu10GeV;
1073   Int_t nEl1GeV, nEl3GeV, nEl10GeV;
1074   Float_t maxPt = .0, meanPt = .0, totalP = .0;
1075
1076   AliRunTag *tag = new AliRunTag();
1077   TTree ttag("T","A Tree with event tags");
1078   TBranch * btag = ttag.Branch("AliTAG", &tag);
1079   btag->SetCompressionLevel(9);
1080   
1081   //reading the esd tag file
1082   TChain *oldTagTree = new TChain("T");
1083   const char * tagPattern = "ESD.tag";
1084   // Open the working directory
1085   void * dirp = gSystem->OpenDirectory(gSystem->pwd());
1086   const char * name = 0x0;
1087   // Add all files matching *pattern* to the chain
1088   while((name = gSystem->GetDirEntry(dirp))) {
1089     if (strstr(name,tagPattern)) oldTagTree->Add(name);  
1090   }//directory loop
1091   AliInfo(Form("Chained tag files: %d",oldTagTree->GetEntries()));
1092
1093   //reading the esd tag file
1094   AliRunTag *oldtag = new AliRunTag();
1095   TString tagFilename;
1096   oldTagTree->SetBranchAddress("AliTAG",&oldtag);
1097   oldTagTree->GetEntry(0);
1098   tag->CopyStandardContent(oldtag);
1099   const TClonesArray *evTagList = oldtag->GetEventTags();
1100
1101   AliInfo(Form("Creating the tags......."));    
1102
1103   TFile *file = TFile::Open("AliAOD.root");
1104   if (!file || !file->IsOpen()) {
1105     AliError(Form("opening failed"));
1106     delete file;
1107     return ;
1108   }
1109   TTree *aodTree = (TTree*)file->Get("AOD");
1110   AliAODEvent *aod = (AliAODEvent*)aodTree->GetUserInfo()->FindObject("AliAODEvent");
1111   TIter next(aod->GetList());
1112   TObject *el;
1113
1114   Int_t lastEvent = 0;  
1115   if(fLastEvent == -1) lastEvent = (Int_t)aodTree->GetEntries();
1116   else lastEvent = fLastEvent;
1117
1118   while((el=(TNamed*)next())) 
1119     aodTree->SetBranchAddress(el->GetName(),aod->GetList()->GetObjectRef(el));
1120   
1121   // loop over events
1122   Int_t nEvents = aodTree->GetEntries();
1123   for (Int_t iEventNumber = 0; iEventNumber < nEvents; iEventNumber++) {
1124     AliEventTag *evTag = (AliEventTag *)evTagList->At(iEventNumber);    
1125     ntrack = 0;
1126     nPos = 0; nNeg = 0; nNeutr =0;
1127     nKinks = 0; nV0s = 0; nCascades = 0;
1128     nK0s = 0; nNeutrons = 0; nPi0s = 0; nGamas = 0;
1129     nProtons = 0;  nKaons = 0; nPions = 0; nMuons = 0; nElectrons = 0;    
1130     nCh1GeV = 0; nCh3GeV = 0; nCh10GeV = 0;
1131     nMu1GeV = 0; nMu3GeV = 0; nMu10GeV = 0;
1132     nEl1GeV = 0; nEl3GeV = 0; nEl10GeV = 0;
1133     maxPt = .0; meanPt = .0; totalP = .0;
1134
1135     // read events
1136     aodTree->GetEvent(iEventNumber);
1137     
1138     // set pointers
1139     aod->GetStdContent();
1140     
1141     Int_t nTracks = aod->GetNTracks();
1142     // loop over vertices
1143     Int_t nVtxs = aod->GetNVertices();
1144     for (Int_t nVtx = 0; nVtx < nVtxs; nVtx++) {
1145       // print track info
1146       AliAODVertex *vertex = aod->GetVertex(nVtx);
1147       if(vertex->GetType() == 1) nKinks += 1;
1148       if(vertex->GetType() == 2) nV0s += 1;
1149       if(vertex->GetType() == 3) nCascades += 1;
1150     }
1151     for (Int_t nTr = 0; nTr < nTracks; nTr++) {
1152       AliAODTrack *track = aod->GetTrack(nTr);
1153       
1154       Double_t fPt = track->Pt();
1155       if(fPt > maxPt) maxPt = fPt;
1156       if(track->Charge() > 0) {
1157         nPos++;
1158         if(fPt > fLowPtCut) nCh1GeV++;
1159         if(fPt > fHighPtCut) nCh3GeV++;
1160         if(fPt > fVeryHighPtCut) nCh10GeV++;
1161       }
1162       if(track->Charge() < 0) {
1163         nNeg++;
1164         if(fPt > fLowPtCut) nCh1GeV++;
1165         if(fPt > fHighPtCut) nCh3GeV++;
1166         if(fPt > fVeryHighPtCut) nCh10GeV++;
1167       }
1168       if(track->Charge() == 0) nNeutr++;
1169
1170       //PID
1171       const Double32_t *prob = track->PID();
1172       Double_t rcc = 0.0;
1173       for(Int_t i = 0; i < AliPID::kSPECIES; i++) rcc += prob[i]*partFrac[i];
1174       if(rcc == 0.0) continue;
1175       //Bayes' formula
1176       Double_t w[10];
1177       for(Int_t i = 0; i < AliPID::kSPECIES; i++) w[i] = prob[i]*partFrac[i]/rcc;
1178       
1179       //protons
1180       if ((w[4]>w[3])&&(w[4]>w[2])&&(w[4]>w[1])&&(w[4]>w[0])) nProtons++;
1181       //kaons
1182       if ((w[3]>w[4])&&(w[3]>w[2])&&(w[3]>w[1])&&(w[3]>w[0])) nKaons++;
1183       //pions
1184       if ((w[2]>w[4])&&(w[2]>w[3])&&(w[2]>w[1])&&(w[2]>w[0])) nPions++; 
1185       //muons
1186       if ((w[1]>w[4])&&(w[1]>w[3])&&(w[1]>w[2])&&(w[1]>w[0])) {
1187         nMuons++;
1188         if(fPt > fLowPtCut) nMu1GeV++;
1189         if(fPt > fHighPtCut) nMu3GeV++;
1190         if(fPt > fVeryHighPtCut) nMu10GeV++;
1191         }
1192       //electrons
1193       if ((w[0]>w[4])&&(w[0]>w[3])&&(w[0]>w[2])&&(w[0]>w[1])) {
1194         nElectrons++;
1195         if(fPt > fLowPtCut) nEl1GeV++;
1196         if(fPt > fHighPtCut) nEl3GeV++;
1197         if(fPt > fVeryHighPtCut) nEl10GeV++;
1198       }
1199
1200       totalP += track->P();
1201       meanPt += fPt;
1202       ntrack++;
1203     }//track loop    
1204     // Fill the event tags 
1205     if(ntrack != 0)
1206       meanPt = meanPt/ntrack;
1207
1208     evTag->SetEventId(iEventNumber+1);
1209         
1210     evTag->SetNumOfTracks(nTracks);
1211     evTag->SetNumOfPosTracks(nPos);
1212     evTag->SetNumOfNegTracks(nNeg);
1213     evTag->SetNumOfNeutrTracks(nNeutr);
1214     
1215     evTag->SetNumOfV0s(nV0s);
1216     evTag->SetNumOfCascades(nCascades);
1217     evTag->SetNumOfKinks(nKinks);
1218     
1219     evTag->SetNumOfProtons(nProtons);
1220     evTag->SetNumOfKaons(nKaons);
1221     evTag->SetNumOfPions(nPions);
1222     evTag->SetNumOfMuons(nMuons);
1223     evTag->SetNumOfElectrons(nElectrons);
1224     evTag->SetNumOfPhotons(nGamas);
1225     evTag->SetNumOfPi0s(nPi0s);
1226     evTag->SetNumOfNeutrons(nNeutrons);
1227     evTag->SetNumOfKaon0s(nK0s);
1228     
1229     evTag->SetNumOfChargedAbove1GeV(nCh1GeV);
1230     evTag->SetNumOfChargedAbove3GeV(nCh3GeV);
1231     evTag->SetNumOfChargedAbove10GeV(nCh10GeV);
1232     evTag->SetNumOfMuonsAbove1GeV(nMu1GeV);
1233     evTag->SetNumOfMuonsAbove3GeV(nMu3GeV);
1234     evTag->SetNumOfMuonsAbove10GeV(nMu10GeV);
1235     evTag->SetNumOfElectronsAbove1GeV(nEl1GeV);
1236     evTag->SetNumOfElectronsAbove3GeV(nEl3GeV);
1237     evTag->SetNumOfElectronsAbove10GeV(nEl10GeV);
1238     
1239     evTag->SetTotalMomentum(totalP);
1240     evTag->SetMeanPt(meanPt);
1241     evTag->SetMaxPt(maxPt);
1242     tag->AddEventTag(*evTag);
1243   }//event loop
1244   if(fLastEvent == -1) lastEvent = (Int_t)aodTree->GetEntries();
1245   else lastEvent = fLastEvent;
1246
1247   ttag.Fill();
1248   tag->Clear();
1249
1250   char fileName[256];
1251   sprintf(fileName, "Run%d.Event%d_%d.AOD.tag.root", 
1252           tag->GetRunId(),fFirstEvent,lastEvent );
1253   AliInfo(Form("writing tags to file %s", fileName));
1254   AliDebug(1, Form("writing tags to file %s", fileName));
1255  
1256   TFile* ftag = TFile::Open(fileName, "recreate");
1257   ftag->cd();
1258   ttag.Write();
1259   ftag->Close();
1260   file->cd();
1261   file->Close();
1262   }*/