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