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