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