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