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