Possibility to set a default momentum (needed for the reconstruction without the...
[u/mrichter/AliRoot.git] / STEER / AliTagCreator.cxx
CommitLineData
b45e5084 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// AliTagCreator 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
c336835b 23#include <Riostream.h>
b45e5084 24#include <TFile.h>
25#include <TString.h>
cb1645b7 26#include <TTree.h>
4fd84456 27#include <TSystem.h>
28#include <TChain.h>
2bdb9d38 29#include <TLorentzVector.h>
b45e5084 30
31//ROOT-AliEn
32#include <TGrid.h>
b45e5084 33#include <TGridResult.h>
b45e5084 34
35//AliRoot
36#include "AliRunTag.h"
37#include "AliEventTag.h"
38#include "AliESD.h"
39#include "AliESDVertex.h"
b45e5084 40#include "AliLog.h"
41
42
43#include "AliTagCreator.h"
44
45
46ClassImp(AliTagCreator)
47
48
49//______________________________________________________________________________
442646f2 50 AliTagCreator::AliTagCreator() :
51 TObject(),
52 fSE("ALICE::CERN::se01"),
53 fgridpath(""),
54 fStorage(0)
55{
b45e5084 56 //==============Default constructor for a AliTagCreator==================
b45e5084 57}
58
59//______________________________________________________________________________
c336835b 60AliTagCreator::~AliTagCreator() {
b45e5084 61//================Default destructor for a AliTagCreator=======================
62}
63
64//______________________________________________________________________________
c336835b 65void AliTagCreator::SetStorage(Int_t storage) {
cb1645b7 66 // Sets correctly the storage: 0 for local, 1 for GRID
28afeb2e 67 fStorage = storage;
68 if(fStorage == 0)
69 AliInfo(Form("Tags will be stored locally...."));
70 if(fStorage == 1)
71 AliInfo(Form("Tags will be stored in the grid...."));
72 if((fStorage != 0)&&(fStorage != 1))
73 {
74 AliInfo(Form("Storage was not properly set!!!"));
75 abort();
76 }
77}
78
b45e5084 79
80//______________________________________________________________________________
c336835b 81Bool_t AliTagCreator::ReadGridCollection(TGridResult *fresult) {
cb1645b7 82 // Reads the entry of the TGridResult and creates the tags
83 Int_t nEntries = fresult->GetEntries();
b45e5084 84
cb1645b7 85 TString alienUrl;
b45e5084 86 const char *guid;
e16601cf 87 const char *md5;
88 const char *turl;
89 Long64_t size = -1;
90
b45e5084 91 Int_t counter = 0;
c336835b 92 for(Int_t i = 0; i < nEntries; i++) {
93 alienUrl = fresult->GetKey(i,"turl");
94 guid = fresult->GetKey(i,"guid");
95 if(fresult->GetKey(i,"size")) size = atol (fresult->GetKey(i,"size"));
96 md5 = fresult->GetKey(i,"md5");
97 turl = fresult->GetKey(i,"turl");
98 if(md5 && !strlen(guid)) md5 = 0;
99 if(guid && !strlen(guid)) guid = 0;
100
101 TFile *f = TFile::Open(alienUrl,"READ");
102 CreateTag(f,guid,md5,turl,size,counter);
103 f->Close();
104 delete f;
105 counter += 1;
106 }//grid result loop
107
108 return kTRUE;
109}
110
111//______________________________________________________________________________
112Bool_t AliTagCreator::ReadLocalCollection(const char *localpath) {
113 // Checks the different subdirs of the given local path and in the
114 // case where it finds an AliESDs.root file it creates the tags
115
116 void *dira = gSystem->OpenDirectory(localpath);
117 Char_t fPath[256];
118 const char * dirname = 0x0;
119 const char * filename = 0x0;
120 const char * pattern = "AliESDs.root";
121
122 Int_t counter = 0;
2005260b 123 while((dirname = gSystem->GetDirEntry(dira))) {
c336835b 124 sprintf(fPath,"%s/%s",localpath,dirname);
125 void *dirb = gSystem->OpenDirectory(fPath);
2005260b 126 while((filename = gSystem->GetDirEntry(dirb))) {
c336835b 127 if(strstr(filename,pattern)) {
128 TString fESDFileName;
129 fESDFileName = fPath;
130 fESDFileName += "/";
131 fESDFileName += pattern;
132 TFile *f = TFile::Open(fESDFileName,"READ");
133 CreateTag(f,fESDFileName,counter);
134 f->Close();
135 delete f;
136
137 counter += 1;
138 }//pattern check
139 }//child directory's entry loop
140 }//parent directory's entry loop
141
142 return kTRUE;
143}
144
145//______________________________________________________________________________
146Bool_t AliTagCreator::ReadCAFCollection(const char *filename) {
147 // Temporary solution for CAF: Takes as an input the ascii file that
148 // lists the ESDs stored in the SE of the CAF and creates the tags.
149
150 // Open the input stream
151 ifstream in;
152 in.open(filename);
153
154 Int_t counter = 0;
155 TString esdfile;
156 // Read the input list of files and add them to the chain
157 while(in.good()) {
158 in >> esdfile;
159 if (!esdfile.Contains("root")) continue; // protection
160 TFile *f = TFile::Open(esdfile,"READ");
161 CreateTag(f,esdfile,counter);
162 f->Close();
163 delete f;
164
165 counter += 1;
166 }
b45e5084 167
168 return kTRUE;
169}
170
c336835b 171
4fd84456 172//__________________________________________________________________________
c336835b 173Bool_t AliTagCreator::MergeTags() {
4fd84456 174 //Merges the tags and stores the merged tag file
175 //locally if fStorage=0 or in the grid if fStorage=1
176 AliInfo(Form("Merging tags....."));
177 TChain *fgChain = new TChain("T");
178
179 if(fStorage == 0) {
180 const char * tagPattern = "tag";
181 // Open the working directory
182 void * dirp = gSystem->OpenDirectory(gSystem->pwd());
183 const char * name = 0x0;
184 // Add all files matching *pattern* to the chain
185 while((name = gSystem->GetDirEntry(dirp))) {
22c08f55 186 if (strstr(name,tagPattern)) fgChain->Add(name);
4fd84456 187 }//directory loop
188 AliInfo(Form("Chained tag files: %d",fgChain->GetEntries()));
189 }//local mode
190
191 else if(fStorage == 1) {
192 TString alienLocation = gGrid->Pwd();
193 alienLocation += fgridpath.Data();
194 alienLocation += "/";
195
196 TGridResult *tagresult = gGrid->Query(alienLocation,"*tag.root","","");
197 Int_t nEntries = tagresult->GetEntries();
198 for(Int_t i = 0; i < nEntries; i++) {
199 TString alienUrl = tagresult->GetKey(i,"turl");
200 fgChain->Add(alienUrl);
201 }//grid result loop
202 AliInfo(Form("Chained tag files: %d",fgChain->GetEntries()));
203 }//grid mode
204
205 AliRunTag *tag = new AliRunTag;
4fd84456 206 fgChain->SetBranchAddress("AliTAG",&tag);
207
208 //Defining new tag objects
209 AliRunTag *newTag = new AliRunTag();
210 TTree ttag("T","A Tree with event tags");
211 TBranch * btag = ttag.Branch("AliTAG", &newTag);
212 btag->SetCompressionLevel(9);
213 for(Int_t iTagFiles = 0; iTagFiles < fgChain->GetEntries(); iTagFiles++) {
214 fgChain->GetEntry(iTagFiles);
22c08f55 215 newTag = tag;
442646f2 216 ttag.Fill();
217 newTag->Clear();
218 }//tag file loop
219
220 TString localFileName = "Run"; localFileName += tag->GetRunId();
221 localFileName += ".Merged"; localFileName += ".ESD.tag.root";
222
223 TString filename = 0x0;
224
225 if(fStorage == 0) {
226 filename = localFileName.Data();
227 AliInfo(Form("Writing merged tags to local file: %s",filename.Data()));
228 }
229 else if(fStorage == 1) {
230 TString alienFileName = "/alien";
231 alienFileName += gGrid->Pwd();
232 alienFileName += fgridpath.Data();
233 alienFileName += "/";
234 alienFileName += localFileName;
235 alienFileName += "?se=";
236 alienFileName += fSE.Data();
237 filename = alienFileName.Data();
238 AliInfo(Form("Writing merged tags to grid file: %s",filename.Data()));
239 }
240
241 TFile* ftag = TFile::Open(filename, "recreate");
242 ftag->cd();
243 ttag.Write();
244 ftag->Close();
245
442646f2 246 delete newTag;
247
248 return kTRUE;
249}
250
251//__________________________________________________________________________
252Bool_t AliTagCreator::MergeTags(TGridResult *result) {
253 //Merges the tags that are listed in the TGridResult
254 AliInfo(Form("Merging tags....."));
255 TChain *fgChain = new TChain("T");
256
257 Int_t nEntries = result->GetEntries();
258
259 TString alienUrl;
260 for(Int_t i = 0; i < nEntries; i++) {
261 alienUrl = result->GetKey(i,"turl");
262 fgChain->Add(alienUrl);
263 }
264 AliInfo(Form("Chained tag files: %d",fgChain->GetEntries()));
265 AliRunTag *tag = new AliRunTag;
442646f2 266 fgChain->SetBranchAddress("AliTAG",&tag);
267
268 //Defining new tag objects
269 AliRunTag *newTag = new AliRunTag();
270 TTree ttag("T","A Tree with event tags");
271 TBranch * btag = ttag.Branch("AliTAG", &newTag);
272 btag->SetCompressionLevel(9);
273 for(Int_t iTagFiles = 0; iTagFiles < fgChain->GetEntries(); iTagFiles++) {
274 fgChain->GetEntry(iTagFiles);
22c08f55 275 newTag = tag;
4fd84456 276 ttag.Fill();
277 newTag->Clear();
278 }//tag file loop
279
280 TString localFileName = "Run"; localFileName += tag->GetRunId();
281 localFileName += ".Merged"; localFileName += ".ESD.tag.root";
282
4fd84456 283 TString filename = 0x0;
284
285 if(fStorage == 0) {
286 filename = localFileName.Data();
287 AliInfo(Form("Writing merged tags to local file: %s",filename.Data()));
288 }
c336835b 289 else if(fStorage == 1) {
290 TString alienFileName = "/alien";
291 alienFileName += gGrid->Pwd();
292 alienFileName += fgridpath.Data();
293 alienFileName += "/";
294 alienFileName += localFileName;
295 alienFileName += "?se=";
296 alienFileName += fSE.Data();
4fd84456 297 filename = alienFileName.Data();
298 AliInfo(Form("Writing merged tags to grid file: %s",filename.Data()));
299 }
300
301 TFile* ftag = TFile::Open(filename, "recreate");
302 ftag->cd();
303 ttag.Write();
304 ftag->Close();
305
4fd84456 306 delete newTag;
307
308 return kTRUE;
309}
b45e5084 310
311//_____________________________________________________________________________
c336835b 312void AliTagCreator::CreateTag(TFile* file, const char *guid, const char *md5, const char *turl, Long64_t size, Int_t Counter) {
313 //private method that creates tag files
314 TString fguid = guid;
315 TString fmd5 = md5;
316 TString fturl = turl;
2bdb9d38 317 /////////////
318 //muon code//
319 ////////////
56982dd1 320 Double_t fMUONMASS = 0.105658369;
2bdb9d38 321 //Variables
56982dd1 322 Double_t fX,fY,fZ ;
323 Double_t fThetaX, fThetaY, fPyz, fChisquare;
324 Double_t fPxRec, fPyRec, fPzRec, fEnergy;
325 Int_t fCharge;
326 TLorentzVector fEPvector;
327
328 Float_t fZVertexCut = 10.0;
329 Float_t fRhoVertexCut = 2.0;
330
331 Float_t fLowPtCut = 1.0;
332 Float_t fHighPtCut = 3.0;
333 Float_t fVeryHighPtCut = 10.0;
2bdb9d38 334 ////////////
335
336 Double_t partFrac[5] = {0.01, 0.01, 0.85, 0.10, 0.05};
337
cb1645b7 338 // Creates the tags for all the events in a given ESD file
c336835b 339 Bool_t fIsSim = kTRUE;
b45e5084 340 Int_t ntrack;
cb1645b7 341 Int_t nProtons, nKaons, nPions, nMuons, nElectrons;
342 Int_t nPos, nNeg, nNeutr;
343 Int_t nK0s, nNeutrons, nPi0s, nGamas;
344 Int_t nCh1GeV, nCh3GeV, nCh10GeV;
345 Int_t nMu1GeV, nMu3GeV, nMu10GeV;
346 Int_t nEl1GeV, nEl3GeV, nEl10GeV;
347 Float_t maxPt = .0, meanPt = .0, totalP = .0;
d1a50cb5 348 Int_t fVertexflag;
ea865c85 349 Int_t iRunNumber = 0;
d1a50cb5 350 TString fVertexName;
b45e5084 351
352 AliRunTag *tag = new AliRunTag();
b45e5084 353 AliEventTag *evTag = new AliEventTag();
354 TTree ttag("T","A Tree with event tags");
4fd84456 355 TBranch * btag = ttag.Branch("AliTAG", &tag);
b45e5084 356 btag->SetCompressionLevel(9);
357
358 AliInfo(Form("Creating the tags......."));
359
360 Int_t firstEvent = 0,lastEvent = 0;
361 TTree *t = (TTree*) file->Get("esdTree");
362 TBranch * b = t->GetBranch("ESD");
363 AliESD *esd = 0;
364 b->SetAddress(&esd);
365
ea865c85 366 b->GetEntry(0);
367 Int_t iInitRunNumber = esd->GetRunNumber();
368
5904bcea 369 Int_t iNumberOfEvents = (Int_t)b->GetEntries();
2bdb9d38 370 for (Int_t iEventNumber = 0; iEventNumber < iNumberOfEvents; iEventNumber++) {
371 ntrack = 0;
372 nPos = 0;
373 nNeg = 0;
374 nNeutr =0;
375 nK0s = 0;
376 nNeutrons = 0;
377 nPi0s = 0;
378 nGamas = 0;
379 nProtons = 0;
380 nKaons = 0;
381 nPions = 0;
382 nMuons = 0;
383 nElectrons = 0;
384 nCh1GeV = 0;
385 nCh3GeV = 0;
386 nCh10GeV = 0;
387 nMu1GeV = 0;
388 nMu3GeV = 0;
389 nMu10GeV = 0;
390 nEl1GeV = 0;
391 nEl3GeV = 0;
392 nEl10GeV = 0;
393 maxPt = .0;
394 meanPt = .0;
395 totalP = .0;
d1a50cb5 396 fVertexflag = 1;
2bdb9d38 397
398 b->GetEntry(iEventNumber);
ea865c85 399 iRunNumber = esd->GetRunNumber();
400 if(iRunNumber != iInitRunNumber) AliFatal("Inconsistency of run numbers in the AliESD!!!");
2bdb9d38 401 const AliESDVertex * vertexIn = esd->GetVertex();
d1a50cb5 402 fVertexName = vertexIn->GetName();
403 if(fVertexName == "default") fVertexflag = 0;
404
2bdb9d38 405 for (Int_t iTrackNumber = 0; iTrackNumber < esd->GetNumberOfTracks(); iTrackNumber++) {
406 AliESDtrack * esdTrack = esd->GetTrack(iTrackNumber);
c336835b 407 if(esdTrack->GetLabel() != 0) fIsSim = kTRUE;
408 else if(esdTrack->GetLabel() == 0) fIsSim = kFALSE;
2bdb9d38 409 UInt_t status = esdTrack->GetStatus();
b45e5084 410
2bdb9d38 411 //select only tracks with ITS refit
412 if ((status&AliESDtrack::kITSrefit)==0) continue;
413 //select only tracks with TPC refit
414 if ((status&AliESDtrack::kTPCrefit)==0) continue;
b45e5084 415
2bdb9d38 416 //select only tracks with the "combined PID"
417 if ((status&AliESDtrack::kESDpid)==0) continue;
418 Double_t p[3];
419 esdTrack->GetPxPyPz(p);
420 Double_t momentum = sqrt(pow(p[0],2) + pow(p[1],2) + pow(p[2],2));
421 Double_t fPt = sqrt(pow(p[0],2) + pow(p[1],2));
422 totalP += momentum;
423 meanPt += fPt;
424 if(fPt > maxPt) maxPt = fPt;
b45e5084 425
2bdb9d38 426 if(esdTrack->GetSign() > 0) {
427 nPos++;
56982dd1 428 if(fPt > fLowPtCut) nCh1GeV++;
429 if(fPt > fHighPtCut) nCh3GeV++;
430 if(fPt > fVeryHighPtCut) nCh10GeV++;
2bdb9d38 431 }
432 if(esdTrack->GetSign() < 0) {
433 nNeg++;
56982dd1 434 if(fPt > fLowPtCut) nCh1GeV++;
435 if(fPt > fHighPtCut) nCh3GeV++;
436 if(fPt > fVeryHighPtCut) nCh10GeV++;
2bdb9d38 437 }
438 if(esdTrack->GetSign() == 0) nNeutr++;
b45e5084 439
2bdb9d38 440 //PID
441 Double_t prob[5];
442 esdTrack->GetESDpid(prob);
b45e5084 443
2bdb9d38 444 Double_t rcc = 0.0;
445 for(Int_t i = 0; i < AliPID::kSPECIES; i++) rcc += prob[i]*partFrac[i];
446 if(rcc == 0.0) continue;
447 //Bayes' formula
448 Double_t w[5];
449 for(Int_t i = 0; i < AliPID::kSPECIES; i++) w[i] = prob[i]*partFrac[i]/rcc;
b45e5084 450
2bdb9d38 451 //protons
452 if ((w[4]>w[3])&&(w[4]>w[2])&&(w[4]>w[1])&&(w[4]>w[0])) nProtons++;
453 //kaons
454 if ((w[3]>w[4])&&(w[3]>w[2])&&(w[3]>w[1])&&(w[3]>w[0])) nKaons++;
455 //pions
456 if ((w[2]>w[4])&&(w[2]>w[3])&&(w[2]>w[1])&&(w[2]>w[0])) nPions++;
457 //electrons
458 if ((w[0]>w[4])&&(w[0]>w[3])&&(w[0]>w[2])&&(w[0]>w[1])) {
459 nElectrons++;
56982dd1 460 if(fPt > fLowPtCut) nEl1GeV++;
461 if(fPt > fHighPtCut) nEl3GeV++;
462 if(fPt > fVeryHighPtCut) nEl10GeV++;
2bdb9d38 463 }
464 ntrack++;
465 }//esd track loop
466
467 /////////////
468 //muon code//
469 ////////////
470 Int_t nMuonTracks = esd->GetNumberOfMuonTracks();
471 // loop over all reconstructed tracks (also first track of combination)
472 for (Int_t iTrack = 0; iTrack < nMuonTracks; iTrack++) {
473 AliESDMuonTrack* muonTrack = esd->GetMuonTrack(iTrack);
474 if (muonTrack == 0x0) continue;
b45e5084 475
2bdb9d38 476 // Coordinates at vertex
56982dd1 477 fZ = muonTrack->GetZ();
478 fY = muonTrack->GetBendingCoor();
479 fX = muonTrack->GetNonBendingCoor();
b45e5084 480
56982dd1 481 fThetaX = muonTrack->GetThetaX();
482 fThetaY = muonTrack->GetThetaY();
b45e5084 483
56982dd1 484 fPyz = 1./TMath::Abs(muonTrack->GetInverseBendingMomentum());
485 fPzRec = - fPyz / TMath::Sqrt(1.0 + TMath::Tan(fThetaY)*TMath::Tan(fThetaY));
486 fPxRec = fPzRec * TMath::Tan(fThetaX);
487 fPyRec = fPzRec * TMath::Tan(fThetaY);
488 fCharge = Int_t(TMath::Sign(1.,muonTrack->GetInverseBendingMomentum()));
b45e5084 489
2bdb9d38 490 //ChiSquare of the track if needed
56982dd1 491 fChisquare = muonTrack->GetChi2()/(2.0 * muonTrack->GetNHit() - 5);
492 fEnergy = TMath::Sqrt(fMUONMASS * fMUONMASS + fPxRec * fPxRec + fPyRec * fPyRec + fPzRec * fPzRec);
493 fEPvector.SetPxPyPzE(fPxRec, fPyRec, fPzRec, fEnergy);
b45e5084 494
2bdb9d38 495 // total number of muons inside a vertex cut
56982dd1 496 if((TMath::Abs(fZ)<fZVertexCut) && (TMath::Sqrt(fY*fY+fX*fX)<fRhoVertexCut)) {
2bdb9d38 497 nMuons++;
56982dd1 498 if(fEPvector.Pt() > fLowPtCut) {
2bdb9d38 499 nMu1GeV++;
56982dd1 500 if(fEPvector.Pt() > fHighPtCut) {
2bdb9d38 501 nMu3GeV++;
56982dd1 502 if (fEPvector.Pt() > fVeryHighPtCut) {
2bdb9d38 503 nMu10GeV++;
504 }
505 }
506 }
507 }
508 }//muon track loop
509
510 // Fill the event tags
c336835b 511 if(ntrack != 0) meanPt = meanPt/ntrack;
2bdb9d38 512
513 evTag->SetEventId(iEventNumber+1);
c336835b 514 evTag->SetGUID(fguid);
515 evTag->SetMD5(fmd5);
516 evTag->SetTURL(fturl);
2bdb9d38 517 evTag->SetSize(size);
518 evTag->SetVertexX(vertexIn->GetXv());
519 evTag->SetVertexY(vertexIn->GetYv());
520 evTag->SetVertexZ(vertexIn->GetZv());
d1a50cb5 521 evTag->SetVertexZError(vertexIn->GetZRes());
522 evTag->SetVertexFlag(fVertexflag);
2bdb9d38 523
524 evTag->SetT0VertexZ(esd->GetT0zVertex());
525
8bd8ac26 526 evTag->SetTriggerMask(esd->GetTriggerMask());
527 evTag->SetTriggerCluster(esd->GetTriggerCluster());
2bdb9d38 528
32a5cab4 529 evTag->SetZDCNeutron1Energy(esd->GetZDCN1Energy());
530 evTag->SetZDCProton1Energy(esd->GetZDCP1Energy());
2bdb9d38 531 evTag->SetZDCEMEnergy(esd->GetZDCEMEnergy());
32a5cab4 532 evTag->SetZDCNeutron1Energy(esd->GetZDCN2Energy());
533 evTag->SetZDCProton1Energy(esd->GetZDCP2Energy());
2bdb9d38 534 evTag->SetNumOfParticipants(esd->GetZDCParticipants());
535
536
537 evTag->SetNumOfTracks(esd->GetNumberOfTracks());
538 evTag->SetNumOfPosTracks(nPos);
539 evTag->SetNumOfNegTracks(nNeg);
540 evTag->SetNumOfNeutrTracks(nNeutr);
541
542 evTag->SetNumOfV0s(esd->GetNumberOfV0s());
543 evTag->SetNumOfCascades(esd->GetNumberOfCascades());
544 evTag->SetNumOfKinks(esd->GetNumberOfKinks());
545 evTag->SetNumOfPMDTracks(esd->GetNumberOfPmdTracks());
546
547 evTag->SetNumOfProtons(nProtons);
548 evTag->SetNumOfKaons(nKaons);
549 evTag->SetNumOfPions(nPions);
550 evTag->SetNumOfMuons(nMuons);
551 evTag->SetNumOfElectrons(nElectrons);
552 evTag->SetNumOfPhotons(nGamas);
553 evTag->SetNumOfPi0s(nPi0s);
554 evTag->SetNumOfNeutrons(nNeutrons);
555 evTag->SetNumOfKaon0s(nK0s);
556
557 evTag->SetNumOfChargedAbove1GeV(nCh1GeV);
558 evTag->SetNumOfChargedAbove3GeV(nCh3GeV);
559 evTag->SetNumOfChargedAbove10GeV(nCh10GeV);
560 evTag->SetNumOfMuonsAbove1GeV(nMu1GeV);
561 evTag->SetNumOfMuonsAbove3GeV(nMu3GeV);
562 evTag->SetNumOfMuonsAbove10GeV(nMu10GeV);
563 evTag->SetNumOfElectronsAbove1GeV(nEl1GeV);
564 evTag->SetNumOfElectronsAbove3GeV(nEl3GeV);
565 evTag->SetNumOfElectronsAbove10GeV(nEl10GeV);
566
85c60a8e 567 evTag->SetNumOfPHOSClusters(esd->GetNumberOfPHOSClusters());
568 evTag->SetNumOfEMCALClusters(esd->GetNumberOfEMCALClusters());
2bdb9d38 569
570 evTag->SetTotalMomentum(totalP);
571 evTag->SetMeanPt(meanPt);
572 evTag->SetMaxPt(maxPt);
573
ea865c85 574 tag->SetRunId(iInitRunNumber);
c336835b 575 if(fIsSim) tag->SetDataType(0);
576 else tag->SetDataType(1);
2bdb9d38 577 tag->AddEventTag(*evTag);
578 }//event loop
cb1645b7 579 lastEvent = iNumberOfEvents;
b45e5084 580
581 t->Delete("");
2bdb9d38 582
b45e5084 583 ttag.Fill();
584 tag->Clear();
585
cb1645b7 586 TString localFileName = "Run"; localFileName += tag->GetRunId();
587 localFileName += ".Event"; localFileName += firstEvent; localFileName += "_"; localFileName += lastEvent; localFileName += "."; localFileName += Counter;
588 localFileName += ".ESD.tag.root";
b45e5084 589
cb1645b7 590 TString fileName;
28afeb2e 591
c336835b 592 if(fStorage == 0) {
593 fileName = localFileName.Data();
594 AliInfo(Form("Writing tags to local file: %s",fileName.Data()));
595 }
596 else if(fStorage == 1) {
597 TString alienLocation = "/alien";
598 alienLocation += gGrid->Pwd();
599 alienLocation += fgridpath.Data();
600 alienLocation += "/";
601 alienLocation += localFileName;
602 alienLocation += "?se=";
603 alienLocation += fSE.Data();
604 fileName = alienLocation.Data();
605 AliInfo(Form("Writing tags to grid file: %s",fileName.Data()));
606 }
b45e5084 607
cb1645b7 608 TFile* ftag = TFile::Open(fileName, "recreate");
b45e5084 609 ftag->cd();
610 ttag.Write();
611 ftag->Close();
612
613 delete ftag;
614 delete esd;
615
616 delete tag;
b45e5084 617}
618
c336835b 619//_____________________________________________________________________________
620void AliTagCreator::CreateTag(TFile* file, const char *filepath, Int_t Counter) {
621 //private method that creates tag files
622 /////////////
623 //muon code//
624 ////////////
625 Double_t fMUONMASS = 0.105658369;
626 //Variables
627 Double_t fX,fY,fZ ;
628 Double_t fThetaX, fThetaY, fPyz, fChisquare;
629 Double_t fPxRec, fPyRec, fPzRec, fEnergy;
630 Int_t fCharge;
631 TLorentzVector fEPvector;
632
633 Float_t fZVertexCut = 10.0;
634 Float_t fRhoVertexCut = 2.0;
635
636 Float_t fLowPtCut = 1.0;
637 Float_t fHighPtCut = 3.0;
638 Float_t fVeryHighPtCut = 10.0;
639 ////////////
640
641 Double_t partFrac[5] = {0.01, 0.01, 0.85, 0.10, 0.05};
642
643 // Creates the tags for all the events in a given ESD file
644 Bool_t fIsSim = kTRUE;
645 Int_t ntrack;
646 Int_t nProtons, nKaons, nPions, nMuons, nElectrons;
647 Int_t nPos, nNeg, nNeutr;
648 Int_t nK0s, nNeutrons, nPi0s, nGamas;
649 Int_t nCh1GeV, nCh3GeV, nCh10GeV;
650 Int_t nMu1GeV, nMu3GeV, nMu10GeV;
651 Int_t nEl1GeV, nEl3GeV, nEl10GeV;
652 Float_t maxPt = .0, meanPt = .0, totalP = .0;
653 Int_t fVertexflag;
654 Int_t iRunNumber = 0;
655 TString fVertexName;
656
657 AliRunTag *tag = new AliRunTag();
658 AliEventTag *evTag = new AliEventTag();
659 TTree ttag("T","A Tree with event tags");
660 TBranch * btag = ttag.Branch("AliTAG", &tag);
661 btag->SetCompressionLevel(9);
662
663 AliInfo(Form("Creating the tags......."));
664
665 Int_t firstEvent = 0,lastEvent = 0;
666
667 TTree *t = (TTree*) file->Get("esdTree");
668 TBranch * b = t->GetBranch("ESD");
669 AliESD *esd = 0;
670 b->SetAddress(&esd);
671
672 b->GetEntry(0);
673 Int_t iInitRunNumber = esd->GetRunNumber();
674
5904bcea 675 Int_t iNumberOfEvents = (Int_t)b->GetEntries();
c336835b 676 for (Int_t iEventNumber = 0; iEventNumber < iNumberOfEvents; iEventNumber++) {
677 ntrack = 0;
678 nPos = 0;
679 nNeg = 0;
680 nNeutr =0;
681 nK0s = 0;
682 nNeutrons = 0;
683 nPi0s = 0;
684 nGamas = 0;
685 nProtons = 0;
686 nKaons = 0;
687 nPions = 0;
688 nMuons = 0;
689 nElectrons = 0;
690 nCh1GeV = 0;
691 nCh3GeV = 0;
692 nCh10GeV = 0;
693 nMu1GeV = 0;
694 nMu3GeV = 0;
695 nMu10GeV = 0;
696 nEl1GeV = 0;
697 nEl3GeV = 0;
698 nEl10GeV = 0;
699 maxPt = .0;
700 meanPt = .0;
701 totalP = .0;
702 fVertexflag = 1;
703
704 b->GetEntry(iEventNumber);
705 iRunNumber = esd->GetRunNumber();
706 if(iRunNumber != iInitRunNumber) AliFatal("Inconsistency of run numbers in the AliESD!!!");
707 const AliESDVertex * vertexIn = esd->GetVertex();
708 fVertexName = vertexIn->GetName();
709 if(fVertexName == "default") fVertexflag = 0;
710
711 for (Int_t iTrackNumber = 0; iTrackNumber < esd->GetNumberOfTracks(); iTrackNumber++) {
712 AliESDtrack * esdTrack = esd->GetTrack(iTrackNumber);
713 if(esdTrack->GetLabel() != 0) fIsSim = kTRUE;
714 else if(esdTrack->GetLabel() == 0) fIsSim = kFALSE;
715 UInt_t status = esdTrack->GetStatus();
716
717 //select only tracks with ITS refit
718 if ((status&AliESDtrack::kITSrefit)==0) continue;
719 //select only tracks with TPC refit
720 if ((status&AliESDtrack::kTPCrefit)==0) continue;
721
722 //select only tracks with the "combined PID"
723 if ((status&AliESDtrack::kESDpid)==0) continue;
724 Double_t p[3];
725 esdTrack->GetPxPyPz(p);
726 Double_t momentum = sqrt(pow(p[0],2) + pow(p[1],2) + pow(p[2],2));
727 Double_t fPt = sqrt(pow(p[0],2) + pow(p[1],2));
728 totalP += momentum;
729 meanPt += fPt;
730 if(fPt > maxPt) maxPt = fPt;
731
732 if(esdTrack->GetSign() > 0) {
733 nPos++;
734 if(fPt > fLowPtCut) nCh1GeV++;
735 if(fPt > fHighPtCut) nCh3GeV++;
736 if(fPt > fVeryHighPtCut) nCh10GeV++;
737 }
738 if(esdTrack->GetSign() < 0) {
739 nNeg++;
740 if(fPt > fLowPtCut) nCh1GeV++;
741 if(fPt > fHighPtCut) nCh3GeV++;
742 if(fPt > fVeryHighPtCut) nCh10GeV++;
743 }
744 if(esdTrack->GetSign() == 0) nNeutr++;
745
746 //PID
747 Double_t prob[5];
748 esdTrack->GetESDpid(prob);
749
750 Double_t rcc = 0.0;
751 for(Int_t i = 0; i < AliPID::kSPECIES; i++) rcc += prob[i]*partFrac[i];
752 if(rcc == 0.0) continue;
753 //Bayes' formula
754 Double_t w[5];
755 for(Int_t i = 0; i < AliPID::kSPECIES; i++) w[i] = prob[i]*partFrac[i]/rcc;
756
757 //protons
758 if ((w[4]>w[3])&&(w[4]>w[2])&&(w[4]>w[1])&&(w[4]>w[0])) nProtons++;
759 //kaons
760 if ((w[3]>w[4])&&(w[3]>w[2])&&(w[3]>w[1])&&(w[3]>w[0])) nKaons++;
761 //pions
762 if ((w[2]>w[4])&&(w[2]>w[3])&&(w[2]>w[1])&&(w[2]>w[0])) nPions++;
763 //electrons
764 if ((w[0]>w[4])&&(w[0]>w[3])&&(w[0]>w[2])&&(w[0]>w[1])) {
765 nElectrons++;
766 if(fPt > fLowPtCut) nEl1GeV++;
767 if(fPt > fHighPtCut) nEl3GeV++;
768 if(fPt > fVeryHighPtCut) nEl10GeV++;
769 }
770 ntrack++;
771 }//esd track loop
772
773 /////////////
774 //muon code//
775 ////////////
776 Int_t nMuonTracks = esd->GetNumberOfMuonTracks();
777 // loop over all reconstructed tracks (also first track of combination)
778 for (Int_t iTrack = 0; iTrack < nMuonTracks; iTrack++) {
779 AliESDMuonTrack* muonTrack = esd->GetMuonTrack(iTrack);
780 if (muonTrack == 0x0) continue;
781
782 // Coordinates at vertex
783 fZ = muonTrack->GetZ();
784 fY = muonTrack->GetBendingCoor();
785 fX = muonTrack->GetNonBendingCoor();
786
787 fThetaX = muonTrack->GetThetaX();
788 fThetaY = muonTrack->GetThetaY();
789
790 fPyz = 1./TMath::Abs(muonTrack->GetInverseBendingMomentum());
791 fPzRec = - fPyz / TMath::Sqrt(1.0 + TMath::Tan(fThetaY)*TMath::Tan(fThetaY));
792 fPxRec = fPzRec * TMath::Tan(fThetaX);
793 fPyRec = fPzRec * TMath::Tan(fThetaY);
794 fCharge = Int_t(TMath::Sign(1.,muonTrack->GetInverseBendingMomentum()));
795
796 //ChiSquare of the track if needed
797 fChisquare = muonTrack->GetChi2()/(2.0 * muonTrack->GetNHit() - 5);
798 fEnergy = TMath::Sqrt(fMUONMASS * fMUONMASS + fPxRec * fPxRec + fPyRec * fPyRec + fPzRec * fPzRec);
799 fEPvector.SetPxPyPzE(fPxRec, fPyRec, fPzRec, fEnergy);
800
801 // total number of muons inside a vertex cut
802 if((TMath::Abs(fZ)<fZVertexCut) && (TMath::Sqrt(fY*fY+fX*fX)<fRhoVertexCut)) {
803 nMuons++;
804 if(fEPvector.Pt() > fLowPtCut) {
805 nMu1GeV++;
806 if(fEPvector.Pt() > fHighPtCut) {
807 nMu3GeV++;
808 if (fEPvector.Pt() > fVeryHighPtCut) {
809 nMu10GeV++;
810 }
811 }
812 }
813 }
814 }//muon track loop
815
816 // Fill the event tags
817 if(ntrack != 0) meanPt = meanPt/ntrack;
818
819 evTag->SetEventId(iEventNumber+1);
820 evTag->SetPath(filepath);
821
822 evTag->SetVertexX(vertexIn->GetXv());
823 evTag->SetVertexY(vertexIn->GetYv());
824 evTag->SetVertexZ(vertexIn->GetZv());
825 evTag->SetVertexZError(vertexIn->GetZRes());
826 evTag->SetVertexFlag(fVertexflag);
827
828 evTag->SetT0VertexZ(esd->GetT0zVertex());
829
830 evTag->SetTriggerMask(esd->GetTriggerMask());
831 evTag->SetTriggerCluster(esd->GetTriggerCluster());
832
833 evTag->SetZDCNeutron1Energy(esd->GetZDCN1Energy());
834 evTag->SetZDCProton1Energy(esd->GetZDCP1Energy());
835 evTag->SetZDCEMEnergy(esd->GetZDCEMEnergy());
836 evTag->SetZDCNeutron1Energy(esd->GetZDCN2Energy());
837 evTag->SetZDCProton1Energy(esd->GetZDCP2Energy());
838 evTag->SetNumOfParticipants(esd->GetZDCParticipants());
839
840
841 evTag->SetNumOfTracks(esd->GetNumberOfTracks());
842 evTag->SetNumOfPosTracks(nPos);
843 evTag->SetNumOfNegTracks(nNeg);
844 evTag->SetNumOfNeutrTracks(nNeutr);
845
846 evTag->SetNumOfV0s(esd->GetNumberOfV0s());
847 evTag->SetNumOfCascades(esd->GetNumberOfCascades());
848 evTag->SetNumOfKinks(esd->GetNumberOfKinks());
849 evTag->SetNumOfPMDTracks(esd->GetNumberOfPmdTracks());
850
851 evTag->SetNumOfProtons(nProtons);
852 evTag->SetNumOfKaons(nKaons);
853 evTag->SetNumOfPions(nPions);
854 evTag->SetNumOfMuons(nMuons);
855 evTag->SetNumOfElectrons(nElectrons);
856 evTag->SetNumOfPhotons(nGamas);
857 evTag->SetNumOfPi0s(nPi0s);
858 evTag->SetNumOfNeutrons(nNeutrons);
859 evTag->SetNumOfKaon0s(nK0s);
860
861 evTag->SetNumOfChargedAbove1GeV(nCh1GeV);
862 evTag->SetNumOfChargedAbove3GeV(nCh3GeV);
863 evTag->SetNumOfChargedAbove10GeV(nCh10GeV);
864 evTag->SetNumOfMuonsAbove1GeV(nMu1GeV);
865 evTag->SetNumOfMuonsAbove3GeV(nMu3GeV);
866 evTag->SetNumOfMuonsAbove10GeV(nMu10GeV);
867 evTag->SetNumOfElectronsAbove1GeV(nEl1GeV);
868 evTag->SetNumOfElectronsAbove3GeV(nEl3GeV);
869 evTag->SetNumOfElectronsAbove10GeV(nEl10GeV);
870
871 evTag->SetNumOfPHOSClusters(esd->GetNumberOfPHOSClusters());
872 evTag->SetNumOfEMCALClusters(esd->GetNumberOfEMCALClusters());
873
874 evTag->SetTotalMomentum(totalP);
875 evTag->SetMeanPt(meanPt);
876 evTag->SetMaxPt(maxPt);
877
878 tag->SetRunId(iInitRunNumber);
879 if(fIsSim) tag->SetDataType(0);
880 else tag->SetDataType(1);
881 tag->AddEventTag(*evTag);
882 }//event loop
883 lastEvent = iNumberOfEvents;
884
885 t->Delete("");
886
887 ttag.Fill();
888 tag->Clear();
889
890 TString localFileName = "Run"; localFileName += tag->GetRunId();
891 localFileName += ".Event"; localFileName += firstEvent; localFileName += "_"; localFileName += lastEvent; localFileName += "."; localFileName += Counter;
892 localFileName += ".ESD.tag.root";
893
894 TString fileName;
895
896 if(fStorage == 0) {
897 fileName = localFileName.Data();
898 AliInfo(Form("Writing tags to local file: %s",fileName.Data()));
899 }
900 else if(fStorage == 1) {
901 TString alienLocation = "/alien";
902 alienLocation += gGrid->Pwd();
903 alienLocation += fgridpath.Data();
904 alienLocation += "/";
905 alienLocation += localFileName;
906 alienLocation += "?se=";
907 alienLocation += fSE.Data();
908 fileName = alienLocation.Data();
909 AliInfo(Form("Writing tags to grid file: %s",fileName.Data()));
910 }
911
912 TFile* ftag = TFile::Open(fileName, "recreate");
913 ftag->cd();
914 ttag.Write();
915 ftag->Close();
916
917 delete ftag;
918 delete esd;
919
920 delete tag;
921}