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