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