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