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