When calculating a*a-b*b the form (a-b)*(a+b) is usually more numerically stable.
[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     evTag->SetEventId(iEventNumber+1);
432     evTag->SetGUID(fguid);
433     if(fSession == "grid") {
434       evTag->SetMD5(0);
435       evTag->SetTURL(fturl);
436       evTag->SetSize(0);
437     }
438     else evTag->SetPath(fturl);
439     
440     evTag->SetVertexX(vertexIn->GetXv());
441     evTag->SetVertexY(vertexIn->GetYv());
442     evTag->SetVertexZ(vertexIn->GetZv());
443     evTag->SetVertexZError(vertexIn->GetZRes());
444     evTag->SetVertexFlag(fVertexflag);
445     
446     evTag->SetT0VertexZ(esd->GetT0zVertex());
447     
448     evTag->SetTriggerMask(esd->GetTriggerMask());
449     evTag->SetTriggerCluster(esd->GetTriggerCluster());
450     
451     evTag->SetZDCNeutron1Energy(esd->GetZDCN1Energy());
452     evTag->SetZDCProton1Energy(esd->GetZDCP1Energy());
453     evTag->SetZDCEMEnergy(esd->GetZDCEMEnergy(0),esd->GetZDCEMEnergy(1));
454     evTag->SetZDCNeutron1Energy(esd->GetZDCN2Energy());
455     evTag->SetZDCProton1Energy(esd->GetZDCP2Energy());
456     evTag->SetNumOfParticipants(esd->GetZDCParticipants());
457     evTag->SetNumOfParticipants(esd->GetZDCParticipants());
458     evTag->SetNumOfParticipants2(esd->GetZDCParticipants2());
459     
460     
461     evTag->SetNumOfTracks(esd->GetNumberOfTracks());
462     evTag->SetNumOfPosTracks(nPos);
463     evTag->SetNumOfNegTracks(nNeg);
464     evTag->SetNumOfNeutrTracks(nNeutr);
465     
466     evTag->SetNumOfV0s(esd->GetNumberOfV0s());
467     evTag->SetNumOfCascades(esd->GetNumberOfCascades());
468     evTag->SetNumOfKinks(esd->GetNumberOfKinks());
469     evTag->SetNumOfPMDTracks(esd->GetNumberOfPmdTracks());
470     
471     evTag->SetNumOfProtons(nProtons);
472     evTag->SetNumOfKaons(nKaons);
473     evTag->SetNumOfPions(nPions);
474     evTag->SetNumOfMuons(nMuons);
475     evTag->SetNumOfFWMuons(nFWMuons);
476     evTag->SetNumOfElectrons(nElectrons);
477     evTag->SetNumOfPhotons(nGamas);
478     evTag->SetNumOfPi0s(nPi0s);
479     evTag->SetNumOfNeutrons(nNeutrons);
480     evTag->SetNumOfKaon0s(nK0s);
481     
482     evTag->SetNumOfChargedAbove1GeV(nCh1GeV);
483     evTag->SetNumOfChargedAbove3GeV(nCh3GeV);
484     evTag->SetNumOfChargedAbove10GeV(nCh10GeV);
485     evTag->SetNumOfMuonsAbove1GeV(nMu1GeV);
486     evTag->SetNumOfMuonsAbove3GeV(nMu3GeV);
487     evTag->SetNumOfMuonsAbove10GeV(nMu10GeV);
488     evTag->SetNumOfElectronsAbove1GeV(nEl1GeV);
489     evTag->SetNumOfElectronsAbove3GeV(nEl3GeV);
490     evTag->SetNumOfElectronsAbove10GeV(nEl10GeV);
491     
492     evTag->SetNumOfPHOSClusters(esd->GetNumberOfPHOSClusters());
493     evTag->SetNumOfEMCALClusters(esd->GetNumberOfEMCALClusters());
494     
495     evTag->SetTotalMomentum(totalP);
496     evTag->SetMeanPt(meanPt);
497     evTag->SetMaxPt(maxPt);
498     
499     tag->SetRunId(iInitRunNumber);
500     if(fIsSim) tag->SetDataType(0);
501     else tag->SetDataType(1);
502
503     if(fguid != fTempGuid) {
504       fTempGuid = fguid;
505       ttag.Fill();
506       tag->Clear("");
507     }
508     tag->AddEventTag(*evTag);
509     if(iEventNumber+1 == chain->GetEntries()) {
510       //AliInfo(Form("File: %s",fturl.Data()));
511       ttag.Fill();
512       tag->Clear("");
513     }      
514   }//event loop
515   lastEvent = chain->GetEntries();
516   
517   //gSystem->GetMemInfo(meminfo);
518   //AliInfo(Form("After the event and track loop - Memory used: %d MB - Increase: %d MB",meminfo->fMemUsed,meminfo->fMemUsed - tempmem));
519   //tempmem = meminfo->fMemUsed;
520
521   //chain->Delete("");
522   
523   //gSystem->GetMemInfo(meminfo);
524   //AliInfo(Form("After the t->Delete - Memory used: %d MB - Increase: %d MB",meminfo->fMemUsed,meminfo->fMemUsed - tempmem));
525   //tempmem = meminfo->fMemUsed;
526
527   ftag->cd();
528   tag->Clear();
529   ttag.Write();
530   ftag->Close();
531
532   //gSystem->GetMemInfo(meminfo);
533   //AliInfo(Form("After the file closing - Memory used: %d MB - Increase: %d MB",meminfo->fMemUsed,meminfo->fMemUsed - tempmem));
534   //tempmem = meminfo->fMemUsed;
535
536   delete esd;
537   delete tag;
538
539   //gSystem->GetMemInfo(meminfo);
540   //AliInfo(Form("After the delete objects - Memory used: %d MB - Increase: %d MB",meminfo->fMemUsed,meminfo->fMemUsed - tempmem));
541 }
542
543 //_____________________________________________________________________________
544 void AliESDTagCreator::CreateTag(TFile* file, const char *guid, const char *md5, const char *turl, Long64_t size, Int_t Counter) {
545   //private method that creates tag files
546   TString fguid = guid;
547   TString fmd5 = md5;
548   TString fturl = turl;
549   /////////////
550   //muon code//
551   ////////////
552   Double_t fMUONMASS = 0.105658369;
553   //Variables
554   Double_t fX,fY,fZ ;
555   Double_t fThetaX, fThetaY, fPyz, fChisquare;
556   Double_t fPxRec, fPyRec, fPzRec, fEnergy;
557   Int_t fCharge;
558   TLorentzVector fEPvector;
559
560   Float_t fZVertexCut = 10.0; 
561   Float_t fRhoVertexCut = 2.0; 
562
563   Float_t fLowPtCut = 1.0;
564   Float_t fHighPtCut = 3.0;
565   Float_t fVeryHighPtCut = 10.0;
566   ////////////
567
568   Double_t partFrac[5] = {0.01, 0.01, 0.85, 0.10, 0.05};
569
570   // Creates the tags for all the events in a given ESD file
571   Bool_t fIsSim = kTRUE;
572   Int_t ntrack;
573   Int_t nProtons, nKaons, nPions, nMuons, nElectrons, nFWMuons;
574   Int_t nPos, nNeg, nNeutr;
575   Int_t nK0s, nNeutrons, nPi0s, nGamas;
576   Int_t nCh1GeV, nCh3GeV, nCh10GeV;
577   Int_t nMu1GeV, nMu3GeV, nMu10GeV;
578   Int_t nEl1GeV, nEl3GeV, nEl10GeV;
579   Float_t maxPt = .0, meanPt = .0, totalP = .0;
580   Int_t fVertexflag;
581   Int_t iRunNumber = 0;
582   TString fVertexName;
583
584   AliRunTag *tag = new AliRunTag();
585   AliEventTag *evTag = new AliEventTag();
586   TTree ttag("T","A Tree with event tags");
587   TBranch * btag = ttag.Branch("AliTAG", &tag);
588   btag->SetCompressionLevel(9);
589   gSystem->GetMemInfo(meminfo);
590   AliInfo(Form("After the tag initialization - Memory used: %d MB",meminfo->fMemUsed));
591   Int_t tempmem = meminfo->fMemUsed;
592   
593   AliInfo(Form("Creating the ESD tags......."));        
594   
595   Int_t firstEvent = 0,lastEvent = 0;
596   TTree *t = (TTree*) file->Get("esdTree");
597   AliESDEvent *esd = new AliESDEvent();
598   esd->ReadFromTree(t);
599   
600   gSystem->GetMemInfo(meminfo);
601   AliInfo(Form("After the esd initialization - Memory used: %d MB - Increase: %d MB",meminfo->fMemUsed,meminfo->fMemUsed - tempmem));
602   tempmem = meminfo->fMemUsed;
603
604   t->GetEntry(0);
605   Int_t iInitRunNumber = esd->GetRunNumber();
606
607   Int_t iNumberOfEvents = (Int_t)t->GetEntries();
608   for (Int_t iEventNumber = 0; iEventNumber < iNumberOfEvents; iEventNumber++) {
609     ntrack = 0;
610     nPos = 0;
611     nNeg = 0;
612     nNeutr =0;
613     nK0s = 0;
614     nNeutrons = 0;
615     nPi0s = 0;
616     nGamas = 0;
617     nProtons = 0;
618     nKaons = 0;
619     nPions = 0;
620     nMuons = 0;
621     nFWMuons = 0;
622     nElectrons = 0;       
623     nCh1GeV = 0;
624     nCh3GeV = 0;
625     nCh10GeV = 0;
626     nMu1GeV = 0;
627     nMu3GeV = 0;
628     nMu10GeV = 0;
629     nEl1GeV = 0;
630     nEl3GeV = 0;
631     nEl10GeV = 0;
632     maxPt = .0;
633     meanPt = .0;
634     totalP = .0;
635     fVertexflag = 1;
636     
637     t->GetEntry(iEventNumber);
638     iRunNumber = esd->GetRunNumber();
639     if(iRunNumber != iInitRunNumber) AliFatal("Inconsistency of run numbers in the AliESD!!!");
640     const AliESDVertex * vertexIn = esd->GetVertex();
641     fVertexName = vertexIn->GetName();
642     if(fVertexName == "default") fVertexflag = 0;
643
644     for (Int_t iTrackNumber = 0; iTrackNumber < esd->GetNumberOfTracks(); iTrackNumber++) {
645       AliESDtrack * esdTrack = esd->GetTrack(iTrackNumber);
646       if(esdTrack->GetLabel() != 0) fIsSim = kTRUE;
647       else if(esdTrack->GetLabel() == 0) fIsSim = kFALSE;
648       UInt_t status = esdTrack->GetStatus();
649       
650       //select only tracks with ITS refit
651       if ((status&AliESDtrack::kITSrefit)==0) continue;
652       //select only tracks with TPC refit
653       if ((status&AliESDtrack::kTPCrefit)==0) continue;
654       
655       //select only tracks with the "combined PID"
656       if ((status&AliESDtrack::kESDpid)==0) continue;
657       Double_t p[3];
658       esdTrack->GetPxPyPz(p);
659       Double_t pt2 = p[0]*p[0]+p[1]*p[1];
660       Double_t momentum = TMath::Sqrt(pt2+p[2]*p[2]);
661       Double_t fPt = TMath::Sqrt(pt2);
662       totalP += momentum;
663       meanPt += fPt;
664       if(fPt > maxPt) maxPt = fPt;
665       
666       if(esdTrack->GetSign() > 0) {
667         nPos++;
668         if(fPt > fLowPtCut) nCh1GeV++;
669         if(fPt > fHighPtCut) nCh3GeV++;
670         if(fPt > fVeryHighPtCut) nCh10GeV++;
671       }
672       if(esdTrack->GetSign() < 0) {
673         nNeg++;
674         if(fPt > fLowPtCut) nCh1GeV++;
675         if(fPt > fHighPtCut) nCh3GeV++;
676         if(fPt > fVeryHighPtCut) nCh10GeV++;
677       }
678       if(esdTrack->GetSign() == 0) nNeutr++;
679       
680       //PID
681       Double_t prob[5];
682       esdTrack->GetESDpid(prob);
683       
684       Double_t rcc = 0.0;
685       for(Int_t i = 0; i < AliPID::kSPECIES; i++) rcc += prob[i]*partFrac[i];
686       if(rcc == 0.0) continue;
687       //Bayes' formula
688       Double_t w[5];
689       for(Int_t i = 0; i < AliPID::kSPECIES; i++) w[i] = prob[i]*partFrac[i]/rcc;
690       
691       //protons
692       if ((w[4]>w[3])&&(w[4]>w[2])&&(w[4]>w[1])&&(w[4]>w[0])) nProtons++;
693       //kaons
694       if ((w[3]>w[4])&&(w[3]>w[2])&&(w[3]>w[1])&&(w[3]>w[0])) nKaons++;
695       //pions
696       if ((w[2]>w[4])&&(w[2]>w[3])&&(w[2]>w[1])&&(w[2]>w[0])) nPions++; 
697       //electrons
698       if ((w[0]>w[4])&&(w[0]>w[3])&&(w[0]>w[2])&&(w[0]>w[1])) {
699         nElectrons++;
700         if(fPt > fLowPtCut) nEl1GeV++;
701         if(fPt > fHighPtCut) nEl3GeV++;
702         if(fPt > fVeryHighPtCut) nEl10GeV++;
703       }   
704       ntrack++;
705     }//esd track loop
706     
707     /////////////
708     //muon code//
709     ////////////
710     Int_t nMuonTracks = esd->GetNumberOfMuonTracks();
711     // loop over all reconstructed tracks (also first track of combination)
712     for (Int_t iTrack = 0; iTrack <  nMuonTracks;  iTrack++) {
713       AliESDMuonTrack* muonTrack = esd->GetMuonTrack(iTrack);
714       if (muonTrack == 0x0) continue;
715       
716       // Coordinates at vertex
717       fZ = muonTrack->GetZ(); 
718       fY = muonTrack->GetBendingCoor();
719       fX = muonTrack->GetNonBendingCoor(); 
720       
721       fThetaX = muonTrack->GetThetaX();
722       fThetaY = muonTrack->GetThetaY();
723       
724       fPyz = 1./TMath::Abs(muonTrack->GetInverseBendingMomentum());
725       fPzRec = - fPyz / TMath::Sqrt(1.0 + TMath::Tan(fThetaY)*TMath::Tan(fThetaY));
726       fPxRec = fPzRec * TMath::Tan(fThetaX);
727       fPyRec = fPzRec * TMath::Tan(fThetaY);
728       fCharge = Int_t(TMath::Sign(1.,muonTrack->GetInverseBendingMomentum()));
729       
730       //ChiSquare of the track if needed
731       fChisquare = muonTrack->GetChi2()/(2.0 * muonTrack->GetNHit() - 5);
732       fEnergy = TMath::Sqrt(fMUONMASS * fMUONMASS + fPxRec * fPxRec + fPyRec * fPyRec + fPzRec * fPzRec);
733       fEPvector.SetPxPyPzE(fPxRec, fPyRec, fPzRec, fEnergy);
734       
735       // total number of muons inside a vertex cut 
736       if((TMath::Abs(fZ)<fZVertexCut) && (TMath::Sqrt(fY*fY+fX*fX)<fRhoVertexCut)) {
737         nMuons++;
738         nFWMuons++;
739         if(fEPvector.Pt() > fLowPtCut) {
740           nMu1GeV++; 
741           if(fEPvector.Pt() > fHighPtCut) {
742             nMu3GeV++; 
743             if (fEPvector.Pt() > fVeryHighPtCut) {
744               nMu10GeV++;
745             }
746           }
747         }
748       }
749     }//muon track loop
750     
751     // Fill the event tags 
752     if(ntrack != 0) meanPt = meanPt/ntrack;
753     
754     evTag->SetEventId(iEventNumber+1);
755     evTag->SetGUID(fguid);
756     evTag->SetMD5(fmd5);
757     evTag->SetTURL(fturl);
758     evTag->SetSize(size);
759     evTag->SetVertexX(vertexIn->GetXv());
760     evTag->SetVertexY(vertexIn->GetYv());
761     evTag->SetVertexZ(vertexIn->GetZv());
762     evTag->SetVertexZError(vertexIn->GetZRes());
763     evTag->SetVertexFlag(fVertexflag);
764     
765     evTag->SetT0VertexZ(esd->GetT0zVertex());
766     
767     evTag->SetTriggerMask(esd->GetTriggerMask());
768     evTag->SetTriggerCluster(esd->GetTriggerCluster());
769     
770     evTag->SetZDCNeutron1Energy(esd->GetZDCN1Energy());
771     evTag->SetZDCProton1Energy(esd->GetZDCP1Energy());
772     evTag->SetZDCEMEnergy(esd->GetZDCEMEnergy(0),esd->GetZDCEMEnergy(1));
773     evTag->SetZDCNeutron1Energy(esd->GetZDCN2Energy());
774     evTag->SetZDCProton1Energy(esd->GetZDCP2Energy());
775     evTag->SetNumOfParticipants(esd->GetZDCParticipants());
776     
777     
778     evTag->SetNumOfTracks(esd->GetNumberOfTracks());
779     evTag->SetNumOfPosTracks(nPos);
780     evTag->SetNumOfNegTracks(nNeg);
781     evTag->SetNumOfNeutrTracks(nNeutr);
782     
783     evTag->SetNumOfV0s(esd->GetNumberOfV0s());
784     evTag->SetNumOfCascades(esd->GetNumberOfCascades());
785     evTag->SetNumOfKinks(esd->GetNumberOfKinks());
786     evTag->SetNumOfPMDTracks(esd->GetNumberOfPmdTracks());
787     
788     evTag->SetNumOfProtons(nProtons);
789     evTag->SetNumOfKaons(nKaons);
790     evTag->SetNumOfPions(nPions);
791     evTag->SetNumOfMuons(nMuons);
792     evTag->SetNumOfFWMuons(nFWMuons);
793     evTag->SetNumOfElectrons(nElectrons);
794     evTag->SetNumOfPhotons(nGamas);
795     evTag->SetNumOfPi0s(nPi0s);
796     evTag->SetNumOfNeutrons(nNeutrons);
797     evTag->SetNumOfKaon0s(nK0s);
798     
799     evTag->SetNumOfChargedAbove1GeV(nCh1GeV);
800     evTag->SetNumOfChargedAbove3GeV(nCh3GeV);
801     evTag->SetNumOfChargedAbove10GeV(nCh10GeV);
802     evTag->SetNumOfMuonsAbove1GeV(nMu1GeV);
803     evTag->SetNumOfMuonsAbove3GeV(nMu3GeV);
804     evTag->SetNumOfMuonsAbove10GeV(nMu10GeV);
805     evTag->SetNumOfElectronsAbove1GeV(nEl1GeV);
806     evTag->SetNumOfElectronsAbove3GeV(nEl3GeV);
807     evTag->SetNumOfElectronsAbove10GeV(nEl10GeV);
808     
809     evTag->SetNumOfPHOSClusters(esd->GetNumberOfPHOSClusters());
810     evTag->SetNumOfEMCALClusters(esd->GetNumberOfEMCALClusters());
811     
812     evTag->SetTotalMomentum(totalP);
813     evTag->SetMeanPt(meanPt);
814     evTag->SetMaxPt(maxPt);
815     
816     tag->SetRunId(iInitRunNumber);
817     if(fIsSim) tag->SetDataType(0);
818     else tag->SetDataType(1);
819     tag->AddEventTag(*evTag);
820   }//event loop
821   lastEvent = iNumberOfEvents;
822   
823   gSystem->GetMemInfo(meminfo);
824   AliInfo(Form("After the event and track loop - Memory used: %d MB - Increase: %d MB",meminfo->fMemUsed,meminfo->fMemUsed - tempmem));
825   tempmem = meminfo->fMemUsed;
826   t->Delete("");
827   
828   gSystem->GetMemInfo(meminfo);
829   AliInfo(Form("After the t->Delete - Memory used: %d MB - Increase: %d MB",meminfo->fMemUsed,meminfo->fMemUsed - tempmem));
830   tempmem = meminfo->fMemUsed;
831
832   TString localFileName = "Run"; localFileName += tag->GetRunId(); 
833   localFileName += ".Event"; localFileName += firstEvent; localFileName += "_"; localFileName += lastEvent; localFileName += "."; localFileName += Counter;
834   localFileName += ".ESD.tag.root";
835
836   TString fileName;
837   
838   if(fStorage == 0) {
839     fileName = localFileName.Data();      
840     AliInfo(Form("Writing tags to local file: %s",fileName.Data()));
841   }
842   else if(fStorage == 1) {
843     TString alienLocation = "/alien";
844     alienLocation += gGrid->Pwd();
845     alienLocation += fgridpath.Data();
846     alienLocation += "/";
847     alienLocation +=  localFileName;
848     alienLocation += "?se=";
849     alienLocation += fSE.Data();
850     fileName = alienLocation.Data();
851     AliInfo(Form("Writing tags to grid file: %s",fileName.Data()));
852   }
853
854   TFile* ftag = TFile::Open(fileName, "recreate");
855   ftag->cd();
856   ttag.Fill();
857   tag->Clear();
858   ttag.Write();
859   ftag->Close();
860
861   gSystem->GetMemInfo(meminfo);
862   AliInfo(Form("After the file closing - Memory used: %d MB - Increase: %d MB",meminfo->fMemUsed,meminfo->fMemUsed - tempmem));
863   tempmem = meminfo->fMemUsed;
864
865   delete ftag;
866   delete esd;
867
868   delete tag;
869   gSystem->GetMemInfo(meminfo);
870   AliInfo(Form("After the delete objects - Memory used: %d MB - Increase: %d MB",meminfo->fMemUsed,meminfo->fMemUsed - tempmem));
871 }
872
873 //_____________________________________________________________________________
874 void AliESDTagCreator::CreateTag(TFile* file, const char *filepath, Int_t Counter) {
875   //private method that creates tag files
876   /////////////
877   //muon code//
878   ////////////
879   Double_t fMUONMASS = 0.105658369;
880   //Variables
881   Double_t fX,fY,fZ ;
882   Double_t fThetaX, fThetaY, fPyz, fChisquare;
883   Double_t fPxRec, fPyRec, fPzRec, fEnergy;
884   Int_t fCharge;
885   TLorentzVector fEPvector;
886
887   Float_t fZVertexCut = 10.0; 
888   Float_t fRhoVertexCut = 2.0; 
889
890   Float_t fLowPtCut = 1.0;
891   Float_t fHighPtCut = 3.0;
892   Float_t fVeryHighPtCut = 10.0;
893   ////////////
894
895   Double_t partFrac[5] = {0.01, 0.01, 0.85, 0.10, 0.05};
896
897   // Creates the tags for all the events in a given ESD file
898   Bool_t fIsSim = kTRUE;
899   Int_t ntrack;
900   Int_t nProtons, nKaons, nPions, nMuons, nElectrons, nFWMuons;
901   Int_t nPos, nNeg, nNeutr;
902   Int_t nK0s, nNeutrons, nPi0s, nGamas;
903   Int_t nCh1GeV, nCh3GeV, nCh10GeV;
904   Int_t nMu1GeV, nMu3GeV, nMu10GeV;
905   Int_t nEl1GeV, nEl3GeV, nEl10GeV;
906   Float_t maxPt = .0, meanPt = .0, totalP = .0;
907   Int_t fVertexflag;
908   Int_t iRunNumber = 0;
909   TString fVertexName;
910
911   AliRunTag *tag = new AliRunTag();
912   AliEventTag *evTag = new AliEventTag();
913   TTree ttag("T","A Tree with event tags");
914   TBranch * btag = ttag.Branch("AliTAG", &tag);
915   btag->SetCompressionLevel(9);
916   
917   AliInfo(Form("Creating the ESD tags......."));        
918   
919   Int_t firstEvent = 0,lastEvent = 0;
920   
921   TTree *t = (TTree*) file->Get("esdTree");
922   AliESDEvent *esd = new AliESDEvent();
923   esd->ReadFromTree(t);
924   
925   t->GetEntry(0);
926   Int_t iInitRunNumber = esd->GetRunNumber();
927
928   Int_t iNumberOfEvents = (Int_t)t->GetEntries();
929   for (Int_t iEventNumber = 0; iEventNumber < iNumberOfEvents; iEventNumber++) {
930     ntrack = 0;
931     nPos = 0;
932     nNeg = 0;
933     nNeutr =0;
934     nK0s = 0;
935     nNeutrons = 0;
936     nPi0s = 0;
937     nGamas = 0;
938     nProtons = 0;
939     nKaons = 0;
940     nPions = 0;
941     nMuons = 0;
942     nFWMuons = 0;
943     nElectrons = 0;       
944     nCh1GeV = 0;
945     nCh3GeV = 0;
946     nCh10GeV = 0;
947     nMu1GeV = 0;
948     nMu3GeV = 0;
949     nMu10GeV = 0;
950     nEl1GeV = 0;
951     nEl3GeV = 0;
952     nEl10GeV = 0;
953     maxPt = .0;
954     meanPt = .0;
955     totalP = .0;
956     fVertexflag = 1;
957     
958     t->GetEntry(iEventNumber);
959     iRunNumber = esd->GetRunNumber();
960     if(iRunNumber != iInitRunNumber) AliFatal("Inconsistency of run numbers in the AliESD!!!");
961     const AliESDVertex * vertexIn = esd->GetVertex();
962     fVertexName = vertexIn->GetName();
963     if(fVertexName == "default") fVertexflag = 0;
964
965     for (Int_t iTrackNumber = 0; iTrackNumber < esd->GetNumberOfTracks(); iTrackNumber++) {
966       AliESDtrack * esdTrack = esd->GetTrack(iTrackNumber);
967       if(esdTrack->GetLabel() != 0) fIsSim = kTRUE;
968       else if(esdTrack->GetLabel() == 0) fIsSim = kFALSE;
969       UInt_t status = esdTrack->GetStatus();
970       
971       //select only tracks with ITS refit
972       if ((status&AliESDtrack::kITSrefit)==0) continue;
973       //select only tracks with TPC refit
974       if ((status&AliESDtrack::kTPCrefit)==0) continue;
975       
976       //select only tracks with the "combined PID"
977       if ((status&AliESDtrack::kESDpid)==0) continue;
978       Double_t p[3];
979       esdTrack->GetPxPyPz(p);
980       Double_t pt2 = p[0]*p[0]+p[1]*p[1];
981       Double_t momentum = TMath::Sqrt(pt2+p[2]*p[2]);
982       Double_t fPt = TMath::Sqrt(pt2);
983       totalP += momentum;
984       meanPt += fPt;
985       if(fPt > maxPt) maxPt = fPt;
986       
987       if(esdTrack->GetSign() > 0) {
988         nPos++;
989         if(fPt > fLowPtCut) nCh1GeV++;
990         if(fPt > fHighPtCut) nCh3GeV++;
991         if(fPt > fVeryHighPtCut) nCh10GeV++;
992       }
993       if(esdTrack->GetSign() < 0) {
994         nNeg++;
995         if(fPt > fLowPtCut) nCh1GeV++;
996         if(fPt > fHighPtCut) nCh3GeV++;
997         if(fPt > fVeryHighPtCut) nCh10GeV++;
998       }
999       if(esdTrack->GetSign() == 0) nNeutr++;
1000       
1001       //PID
1002       Double_t prob[5];
1003       esdTrack->GetESDpid(prob);
1004       
1005       Double_t rcc = 0.0;
1006       for(Int_t i = 0; i < AliPID::kSPECIES; i++) rcc += prob[i]*partFrac[i];
1007       if(rcc == 0.0) continue;
1008       //Bayes' formula
1009       Double_t w[5];
1010       for(Int_t i = 0; i < AliPID::kSPECIES; i++) w[i] = prob[i]*partFrac[i]/rcc;
1011       
1012       //protons
1013       if ((w[4]>w[3])&&(w[4]>w[2])&&(w[4]>w[1])&&(w[4]>w[0])) nProtons++;
1014       //kaons
1015       if ((w[3]>w[4])&&(w[3]>w[2])&&(w[3]>w[1])&&(w[3]>w[0])) nKaons++;
1016       //pions
1017       if ((w[2]>w[4])&&(w[2]>w[3])&&(w[2]>w[1])&&(w[2]>w[0])) nPions++; 
1018       //electrons
1019       if ((w[0]>w[4])&&(w[0]>w[3])&&(w[0]>w[2])&&(w[0]>w[1])) {
1020         nElectrons++;
1021         if(fPt > fLowPtCut) nEl1GeV++;
1022         if(fPt > fHighPtCut) nEl3GeV++;
1023         if(fPt > fVeryHighPtCut) nEl10GeV++;
1024       }   
1025       ntrack++;
1026     }//esd track loop
1027     
1028     /////////////
1029     //muon code//
1030     ////////////
1031     Int_t nMuonTracks = esd->GetNumberOfMuonTracks();
1032     // loop over all reconstructed tracks (also first track of combination)
1033     for (Int_t iTrack = 0; iTrack <  nMuonTracks;  iTrack++) {
1034       AliESDMuonTrack* muonTrack = esd->GetMuonTrack(iTrack);
1035       if (muonTrack == 0x0) continue;
1036       
1037       // Coordinates at vertex
1038       fZ = muonTrack->GetZ(); 
1039       fY = muonTrack->GetBendingCoor();
1040       fX = muonTrack->GetNonBendingCoor(); 
1041       
1042       fThetaX = muonTrack->GetThetaX();
1043       fThetaY = muonTrack->GetThetaY();
1044       
1045       fPyz = 1./TMath::Abs(muonTrack->GetInverseBendingMomentum());
1046       fPzRec = - fPyz / TMath::Sqrt(1.0 + TMath::Tan(fThetaY)*TMath::Tan(fThetaY));
1047       fPxRec = fPzRec * TMath::Tan(fThetaX);
1048       fPyRec = fPzRec * TMath::Tan(fThetaY);
1049       fCharge = Int_t(TMath::Sign(1.,muonTrack->GetInverseBendingMomentum()));
1050       
1051       //ChiSquare of the track if needed
1052       fChisquare = muonTrack->GetChi2()/(2.0 * muonTrack->GetNHit() - 5);
1053       fEnergy = TMath::Sqrt(fMUONMASS * fMUONMASS + fPxRec * fPxRec + fPyRec * fPyRec + fPzRec * fPzRec);
1054       fEPvector.SetPxPyPzE(fPxRec, fPyRec, fPzRec, fEnergy);
1055       
1056       // total number of muons inside a vertex cut 
1057       if((TMath::Abs(fZ)<fZVertexCut) && (TMath::Sqrt(fY*fY+fX*fX)<fRhoVertexCut)) {
1058         nMuons++;
1059         nFWMuons++;
1060         if(fEPvector.Pt() > fLowPtCut) {
1061           nMu1GeV++; 
1062           if(fEPvector.Pt() > fHighPtCut) {
1063             nMu3GeV++; 
1064             if (fEPvector.Pt() > fVeryHighPtCut) {
1065               nMu10GeV++;
1066             }
1067           }
1068         }
1069       }
1070     }//muon track loop
1071     
1072     // Fill the event tags 
1073     if(ntrack != 0) meanPt = meanPt/ntrack;
1074     
1075     evTag->SetEventId(iEventNumber+1);
1076     evTag->SetPath(filepath);
1077  
1078     evTag->SetVertexX(vertexIn->GetXv());
1079     evTag->SetVertexY(vertexIn->GetYv());
1080     evTag->SetVertexZ(vertexIn->GetZv());
1081     evTag->SetVertexZError(vertexIn->GetZRes());
1082     evTag->SetVertexFlag(fVertexflag);
1083     
1084     evTag->SetT0VertexZ(esd->GetT0zVertex());
1085     
1086     evTag->SetTriggerMask(esd->GetTriggerMask());
1087     evTag->SetTriggerCluster(esd->GetTriggerCluster());
1088     
1089     evTag->SetZDCNeutron1Energy(esd->GetZDCN1Energy());
1090     evTag->SetZDCProton1Energy(esd->GetZDCP1Energy());
1091     evTag->SetZDCEMEnergy(esd->GetZDCEMEnergy(0),esd->GetZDCEMEnergy(1));
1092     evTag->SetZDCNeutron1Energy(esd->GetZDCN2Energy());
1093     evTag->SetZDCProton1Energy(esd->GetZDCP2Energy());
1094     evTag->SetNumOfParticipants(esd->GetZDCParticipants());
1095     
1096     
1097     evTag->SetNumOfTracks(esd->GetNumberOfTracks());
1098     evTag->SetNumOfPosTracks(nPos);
1099     evTag->SetNumOfNegTracks(nNeg);
1100     evTag->SetNumOfNeutrTracks(nNeutr);
1101     
1102     evTag->SetNumOfV0s(esd->GetNumberOfV0s());
1103     evTag->SetNumOfCascades(esd->GetNumberOfCascades());
1104     evTag->SetNumOfKinks(esd->GetNumberOfKinks());
1105     evTag->SetNumOfPMDTracks(esd->GetNumberOfPmdTracks());
1106     
1107     evTag->SetNumOfProtons(nProtons);
1108     evTag->SetNumOfKaons(nKaons);
1109     evTag->SetNumOfPions(nPions);
1110     evTag->SetNumOfMuons(nMuons);
1111     evTag->SetNumOfFWMuons(nFWMuons);
1112     evTag->SetNumOfElectrons(nElectrons);
1113     evTag->SetNumOfPhotons(nGamas);
1114     evTag->SetNumOfPi0s(nPi0s);
1115     evTag->SetNumOfNeutrons(nNeutrons);
1116     evTag->SetNumOfKaon0s(nK0s);
1117     
1118     evTag->SetNumOfChargedAbove1GeV(nCh1GeV);
1119     evTag->SetNumOfChargedAbove3GeV(nCh3GeV);
1120     evTag->SetNumOfChargedAbove10GeV(nCh10GeV);
1121     evTag->SetNumOfMuonsAbove1GeV(nMu1GeV);
1122     evTag->SetNumOfMuonsAbove3GeV(nMu3GeV);
1123     evTag->SetNumOfMuonsAbove10GeV(nMu10GeV);
1124     evTag->SetNumOfElectronsAbove1GeV(nEl1GeV);
1125     evTag->SetNumOfElectronsAbove3GeV(nEl3GeV);
1126     evTag->SetNumOfElectronsAbove10GeV(nEl10GeV);
1127     
1128     evTag->SetNumOfPHOSClusters(esd->GetNumberOfPHOSClusters());
1129     evTag->SetNumOfEMCALClusters(esd->GetNumberOfEMCALClusters());
1130     
1131     evTag->SetTotalMomentum(totalP);
1132     evTag->SetMeanPt(meanPt);
1133     evTag->SetMaxPt(maxPt);
1134     
1135     tag->SetRunId(iInitRunNumber);
1136     if(fIsSim) tag->SetDataType(0);
1137     else tag->SetDataType(1);
1138     tag->AddEventTag(*evTag);
1139   }//event loop
1140   lastEvent = iNumberOfEvents;
1141   
1142   t->Delete("");
1143   
1144   TString localFileName = "Run"; localFileName += tag->GetRunId(); 
1145   localFileName += ".Event"; localFileName += firstEvent; localFileName += "_"; localFileName += lastEvent; localFileName += "."; localFileName += Counter;
1146   localFileName += ".ESD.tag.root";
1147
1148   TString fileName;
1149   
1150   if(fStorage == 0) {
1151     fileName = localFileName.Data();      
1152     AliInfo(Form("Writing tags to local file: %s",fileName.Data()));
1153   }
1154   else if(fStorage == 1) {
1155     TString alienLocation = "/alien";
1156     alienLocation += gGrid->Pwd();
1157     alienLocation += fgridpath.Data();
1158     alienLocation += "/";
1159     alienLocation +=  localFileName;
1160     alienLocation += "?se=";
1161     alienLocation += fSE.Data();
1162     fileName = alienLocation.Data();
1163     AliInfo(Form("Writing tags to grid file: %s",fileName.Data()));
1164   }
1165
1166   TFile* ftag = TFile::Open(fileName, "recreate");
1167   ftag->cd();
1168   ttag.Fill();
1169   tag->Clear();
1170   ttag.Write();
1171   ftag->Close();
1172
1173   delete ftag;
1174   delete esd;
1175
1176   delete tag;
1177 }
1178
1179 //_____________________________________________________________________________
1180 void AliESDTagCreator::CreateESDTags(Int_t fFirstEvent, Int_t fLastEvent, AliGRPObject *grpData, ULong_t * qa, Bool_t * es, Int_t qalength, Int_t eslength) {
1181   //GRP
1182   Float_t lhcLuminosity = 0.0;
1183   TString lhcState = "test";
1184   //UInt_t detectorMask = 0;
1185   Int_t detectorMask = 0;
1186
1187   detectorMask = grpData->GetDetectorMask();
1188   time_t startTime = grpData->GetTimeStart();
1189   TTimeStamp *t1 = new TTimeStamp(startTime);
1190   time_t endTime = grpData->GetTimeEnd();
1191   TTimeStamp *t2 = new TTimeStamp(endTime);
1192   const char* beamtype = grpData->GetBeamType();
1193   Float_t beamenergy = grpData->GetBeamEnergy();
1194
1195
1196   /////////////
1197   //muon code//
1198   ////////////
1199   Double_t fMUONMASS = 0.105658369;
1200   //Variables
1201   Double_t fX,fY,fZ ;
1202   Double_t fThetaX, fThetaY, fPyz, fChisquare;
1203   Double_t fPxRec,fPyRec, fPzRec, fEnergy;
1204   Int_t fCharge;
1205   TLorentzVector fEPvector;
1206
1207   Float_t fZVertexCut = 10.0; 
1208   Float_t fRhoVertexCut = 2.0; 
1209
1210   Float_t fLowPtCut = 1.0;
1211   Float_t fHighPtCut = 3.0;
1212   Float_t fVeryHighPtCut = 10.0;
1213   ////////////
1214
1215   Double_t partFrac[5] = {0.01, 0.01, 0.85, 0.10, 0.05};
1216
1217   // Creates the tags for all the events in a given ESD file
1218   Int_t ntrack;
1219   Int_t nProtons, nKaons, nPions, nMuons, nElectrons, nFWMuons;
1220   Int_t nPos, nNeg, nNeutr;
1221   Int_t nK0s, nNeutrons, nPi0s, nGamas;
1222   Int_t nCh1GeV, nCh3GeV, nCh10GeV;
1223   Int_t nMu1GeV, nMu3GeV, nMu10GeV;
1224   Int_t nEl1GeV, nEl3GeV, nEl10GeV;
1225   Float_t maxPt = .0, meanPt = .0, totalP = .0;
1226   Int_t fVertexflag;
1227   Int_t iRunNumber = 0;
1228   TString fVertexName("default");
1229   
1230   AliInfo(Form("Creating the ESD tags......."));        
1231
1232   TFile *file = TFile::Open("AliESDs.root");
1233   if (!file || !file->IsOpen()) {
1234     AliError(Form("opening failed"));
1235     delete file;
1236     return ;
1237   }  
1238   Int_t lastEvent = 0;
1239   TTree *b = (TTree*) file->Get("esdTree");
1240   AliESDEvent *esd = new AliESDEvent();
1241   esd->ReadFromTree(b);
1242
1243   b->GetEntry(fFirstEvent);
1244   Int_t iInitRunNumber = esd->GetRunNumber();
1245   
1246   Int_t iNumberOfEvents = (Int_t)b->GetEntries();
1247   if(fLastEvent == -1) lastEvent = (Int_t)b->GetEntries();
1248   else lastEvent = fLastEvent;
1249
1250   char fileName[256];
1251   sprintf(fileName, "Run%d.Event%d_%d.ESD.tag.root", 
1252           iInitRunNumber,fFirstEvent,lastEvent);
1253   AliInfo(Form("writing tags to file %s", fileName));
1254   AliDebug(1, Form("writing tags to file %s", fileName));
1255  
1256   TFile* ftag = TFile::Open(fileName, "recreate");
1257  
1258   AliRunTag *tag = new AliRunTag();
1259   AliEventTag *evTag = new AliEventTag();
1260   TTree ttag("T","A Tree with event tags");
1261   TBranch * btag = ttag.Branch("AliTAG", &tag);
1262   btag->SetCompressionLevel(9);
1263
1264   if(fLastEvent != -1) iNumberOfEvents = fLastEvent + 1;
1265   for (Int_t iEventNumber = fFirstEvent; iEventNumber < iNumberOfEvents; iEventNumber++) {
1266     ntrack = 0;
1267     nPos = 0;
1268     nNeg = 0;
1269     nNeutr =0;
1270     nK0s = 0;
1271     nNeutrons = 0;
1272     nPi0s = 0;
1273     nGamas = 0;
1274     nProtons = 0;
1275     nKaons = 0;
1276     nPions = 0;
1277     nMuons = 0;
1278     nFWMuons = 0;
1279     nElectrons = 0;       
1280     nCh1GeV = 0;
1281     nCh3GeV = 0;
1282     nCh10GeV = 0;
1283     nMu1GeV = 0;
1284     nMu3GeV = 0;
1285     nMu10GeV = 0;
1286     nEl1GeV = 0;
1287     nEl3GeV = 0;
1288     nEl10GeV = 0;
1289     maxPt = .0;
1290     meanPt = .0;
1291     totalP = .0;
1292     fVertexflag = 0;
1293
1294     b->GetEntry(iEventNumber);
1295     iRunNumber = esd->GetRunNumber();
1296     if(iRunNumber != iInitRunNumber) AliFatal("Inconsistency of run numbers in the AliESD!!!");
1297     const AliESDVertex * vertexIn = esd->GetVertex();
1298     if (!vertexIn) AliError("ESD has not defined vertex.");
1299     if (vertexIn) fVertexName = vertexIn->GetName();
1300     if(fVertexName != "default") fVertexflag = 1;
1301     for (Int_t iTrackNumber = 0; iTrackNumber < esd->GetNumberOfTracks(); iTrackNumber++) {
1302       AliESDtrack * esdTrack = esd->GetTrack(iTrackNumber);
1303       UInt_t status = esdTrack->GetStatus();
1304       
1305       //select only tracks with ITS refit
1306       if ((status&AliESDtrack::kITSrefit)==0) continue;
1307       //select only tracks with TPC refit
1308       if ((status&AliESDtrack::kTPCrefit)==0) continue;
1309       
1310       //select only tracks with the "combined PID"
1311       if ((status&AliESDtrack::kESDpid)==0) continue;
1312       Double_t p[3];
1313       esdTrack->GetPxPyPz(p);
1314       Double_t pt2 = p[0]*p[0]+p[1]*p[1];
1315       Double_t momentum = TMath::Sqrt(pt2+p[2]*p[2]);
1316       Double_t fPt = TMath::Sqrt(pt2);
1317       totalP += momentum;
1318       meanPt += fPt;
1319       if(fPt > maxPt) maxPt = fPt;
1320       
1321       if(esdTrack->GetSign() > 0) {
1322         nPos++;
1323         if(fPt > fLowPtCut) nCh1GeV++;
1324         if(fPt > fHighPtCut) nCh3GeV++;
1325         if(fPt > fVeryHighPtCut) nCh10GeV++;
1326       }
1327       if(esdTrack->GetSign() < 0) {
1328         nNeg++;
1329         if(fPt > fLowPtCut) nCh1GeV++;
1330         if(fPt > fHighPtCut) nCh3GeV++;
1331         if(fPt > fVeryHighPtCut) nCh10GeV++;
1332       }
1333       if(esdTrack->GetSign() == 0) nNeutr++;
1334       
1335       //PID
1336       Double_t prob[5];
1337       esdTrack->GetESDpid(prob);
1338       
1339       Double_t rcc = 0.0;
1340       for(Int_t i = 0; i < AliPID::kSPECIES; i++) rcc += prob[i]*partFrac[i];
1341       if(rcc == 0.0) continue;
1342       //Bayes' formula
1343       Double_t w[5];
1344       for(Int_t i = 0; i < AliPID::kSPECIES; i++) w[i] = prob[i]*partFrac[i]/rcc;
1345       
1346       //protons
1347       if ((w[4]>w[3])&&(w[4]>w[2])&&(w[4]>w[1])&&(w[4]>w[0])) nProtons++;
1348       //kaons
1349       if ((w[3]>w[4])&&(w[3]>w[2])&&(w[3]>w[1])&&(w[3]>w[0])) nKaons++;
1350       //pions
1351       if ((w[2]>w[4])&&(w[2]>w[3])&&(w[2]>w[1])&&(w[2]>w[0])) nPions++; 
1352       //electrons
1353       if ((w[0]>w[4])&&(w[0]>w[3])&&(w[0]>w[2])&&(w[0]>w[1])) {
1354         nElectrons++;
1355         if(fPt > fLowPtCut) nEl1GeV++;
1356         if(fPt > fHighPtCut) nEl3GeV++;
1357         if(fPt > fVeryHighPtCut) nEl10GeV++;
1358       }   
1359       ntrack++;
1360     }//track loop
1361     
1362     /////////////
1363     //muon code//
1364     ////////////
1365     Int_t nMuonTracks = esd->GetNumberOfMuonTracks();
1366     // loop over all reconstructed tracks (also first track of combination)
1367     for (Int_t iTrack = 0; iTrack <  nMuonTracks;  iTrack++) {
1368       AliESDMuonTrack* muonTrack = esd->GetMuonTrack(iTrack);
1369       if (muonTrack == 0x0) continue;
1370       
1371       // Coordinates at vertex
1372       fZ = muonTrack->GetZ(); 
1373       fY = muonTrack->GetBendingCoor();
1374       fX = muonTrack->GetNonBendingCoor(); 
1375       
1376       fThetaX = muonTrack->GetThetaX();
1377       fThetaY = muonTrack->GetThetaY();
1378       
1379       fPyz = 1./TMath::Abs(muonTrack->GetInverseBendingMomentum());
1380       fPzRec = - fPyz / TMath::Sqrt(1.0 + TMath::Tan(fThetaY)*TMath::Tan(fThetaY));
1381       fPxRec = fPzRec * TMath::Tan(fThetaX);
1382       fPyRec = fPzRec * TMath::Tan(fThetaY);
1383       fCharge = Int_t(TMath::Sign(1.,muonTrack->GetInverseBendingMomentum()));
1384       
1385       //ChiSquare of the track if needed
1386       fChisquare = muonTrack->GetChi2()/(2.0 * muonTrack->GetNHit() - 5);
1387       fEnergy = TMath::Sqrt(fMUONMASS * fMUONMASS + fPxRec * fPxRec + fPyRec * fPyRec + fPzRec * fPzRec);
1388       fEPvector.SetPxPyPzE(fPxRec, fPyRec, fPzRec, fEnergy);
1389       
1390       // total number of muons inside a vertex cut 
1391       if((TMath::Abs(fZ)<fZVertexCut) && (TMath::Sqrt(fY*fY+fX*fX)<fRhoVertexCut)) {
1392         nMuons++;
1393         nFWMuons++;
1394         if(fEPvector.Pt() > fLowPtCut) {
1395           nMu1GeV++; 
1396           if(fEPvector.Pt() > fHighPtCut) {
1397             nMu3GeV++; 
1398             if (fEPvector.Pt() > fVeryHighPtCut) {
1399               nMu10GeV++;
1400             }
1401           }
1402         }
1403       }
1404     }//muon track loop
1405     
1406     // Fill the event tags 
1407     if(ntrack != 0)
1408       meanPt = meanPt/ntrack;
1409     
1410     evTag->SetEventId(iEventNumber+1);
1411     if (vertexIn) {
1412       evTag->SetVertexX(vertexIn->GetXv());
1413       evTag->SetVertexY(vertexIn->GetYv());
1414       evTag->SetVertexZ(vertexIn->GetZv());
1415       evTag->SetVertexZError(vertexIn->GetZRes());
1416     }  
1417     evTag->SetVertexFlag(fVertexflag);
1418
1419     evTag->SetT0VertexZ(esd->GetT0zVertex());
1420     
1421     evTag->SetTriggerMask(esd->GetTriggerMask());
1422     evTag->SetTriggerCluster(esd->GetTriggerCluster());
1423     
1424     evTag->SetZDCNeutron1Energy(esd->GetZDCN1Energy());
1425     evTag->SetZDCProton1Energy(esd->GetZDCP1Energy());
1426     evTag->SetZDCNeutron2Energy(esd->GetZDCN2Energy());
1427     evTag->SetZDCProton2Energy(esd->GetZDCP2Energy());
1428     evTag->SetZDCEMEnergy(esd->GetZDCEMEnergy(0),esd->GetZDCEMEnergy(1));
1429     evTag->SetNumOfParticipants(esd->GetZDCParticipants());
1430     
1431     
1432     evTag->SetNumOfTracks(esd->GetNumberOfTracks());
1433     evTag->SetNumOfPosTracks(nPos);
1434     evTag->SetNumOfNegTracks(nNeg);
1435     evTag->SetNumOfNeutrTracks(nNeutr);
1436     
1437     evTag->SetNumOfV0s(esd->GetNumberOfV0s());
1438     evTag->SetNumOfCascades(esd->GetNumberOfCascades());
1439     evTag->SetNumOfKinks(esd->GetNumberOfKinks());
1440     evTag->SetNumOfPMDTracks(esd->GetNumberOfPmdTracks());
1441     
1442     evTag->SetNumOfProtons(nProtons);
1443     evTag->SetNumOfKaons(nKaons);
1444     evTag->SetNumOfPions(nPions);
1445     evTag->SetNumOfMuons(nMuons);
1446     evTag->SetNumOfFWMuons(nFWMuons);
1447     evTag->SetNumOfElectrons(nElectrons);
1448     evTag->SetNumOfPhotons(nGamas);
1449     evTag->SetNumOfPi0s(nPi0s);
1450     evTag->SetNumOfNeutrons(nNeutrons);
1451     evTag->SetNumOfKaon0s(nK0s);
1452     
1453     evTag->SetNumOfChargedAbove1GeV(nCh1GeV);
1454     evTag->SetNumOfChargedAbove3GeV(nCh3GeV);
1455     evTag->SetNumOfChargedAbove10GeV(nCh10GeV);
1456     evTag->SetNumOfMuonsAbove1GeV(nMu1GeV);
1457     evTag->SetNumOfMuonsAbove3GeV(nMu3GeV);
1458     evTag->SetNumOfMuonsAbove10GeV(nMu10GeV);
1459     evTag->SetNumOfElectronsAbove1GeV(nEl1GeV);
1460     evTag->SetNumOfElectronsAbove3GeV(nEl3GeV);
1461     evTag->SetNumOfElectronsAbove10GeV(nEl10GeV);
1462     
1463     evTag->SetNumOfPHOSClusters(esd->GetNumberOfPHOSClusters());
1464     evTag->SetNumOfEMCALClusters(esd->GetNumberOfEMCALClusters());
1465     
1466     evTag->SetTotalMomentum(totalP);
1467     evTag->SetMeanPt(meanPt);
1468     evTag->SetMaxPt(maxPt);
1469     
1470     tag->SetLHCTag(lhcLuminosity,lhcState);
1471     tag->SetDetectorTag(detectorMask);
1472
1473     tag->SetRunId(iInitRunNumber);
1474     tag->SetRunStartTime(t1->GetDate());
1475     tag->SetRunStopTime(t2->GetDate());
1476     tag->SetBeamEnergy(beamenergy);
1477     tag->SetBeamType(beamtype);
1478     
1479     //QA setting 
1480     tag->SetQA(qa, qalength) ; 
1481     tag->SetEventSpecies(es, eslength) ;
1482
1483     tag->AddEventTag(*evTag);
1484   }
1485         
1486   ftag->cd();
1487   ttag.Fill();
1488   tag->Clear();
1489   ttag.Write();
1490   ftag->Close();
1491   file->cd();
1492   delete tag;
1493   delete evTag;
1494 }
1495
1496 //_____________________________________________________________________________
1497 void AliESDTagCreator::SwitchOffBranches() const {
1498   //
1499   // Switch of branches on user request
1500   TObjArray * tokens = fBranches.Tokenize(" ");
1501   Int_t ntok = tokens->GetEntries();
1502   for (Int_t i = 0; i < ntok; i++)  {
1503     TString str = ((TObjString*) tokens->At(i))->GetString();
1504     fChain->SetBranchStatus(Form("%s%s%s","*", str.Data(), "*"), 0);
1505     AliInfo(Form("Branch %s switched off \n", str.Data()));
1506   }
1507 }