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