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