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