first checkin of Claus' code with a few modifications
[u/mrichter/AliRoot.git] / PWG0 / esdTrackCuts / testESDtrackCuts.C
1 testESDtrackCuts(Char_t* dataDir=, Int_t nRuns=10) {
2
3   Char_t str[256];
4
5   gSystem->Load("libESDtrackQuality.so");
6
7   // ########################################################
8   // definition of ESD track cuts
9
10   ESDtrackQualityCuts* trackCuts = new ESDtrackQualityCuts();
11   trackCuts->DefineHistograms(4);
12
13   trackCuts->SetMinNClustersTPC(50);
14   trackCuts->SetMaxChi2PerClusterTPC(3.5);
15   trackCuts->SetMaxCovDiagonalElements(2,2,0.5,0.5,2);
16   trackCuts->SetRequireTPCRefit(kTRUE);
17
18   trackCuts->SetMinNsigmaToVertex(3);
19   trackCuts->SetAcceptKingDaughters(kFALSE);
20
21   AliLog::SetClassDebugLevel("ESDtrackQualityCuts",1);
22
23   // ########################################################
24   // definition of used pointers
25   TFile* esdFile;
26   TTree* esdTree;
27   TBranch* esdBranch;
28
29   AliESD* esd = 0;
30
31   // ########################################################
32   // get the data dir  
33   Char_t execDir[256];
34   sprintf(execDir,gSystem->pwd());
35   TSystemDirectory* baseDir = new TSystemDirectory(".",dataDir);
36   TList* dirList            = baseDir->GetListOfFiles();
37   Int_t nDirs               = dirList->GetEntries();
38   // go back to the dir where this script is executed
39   gSystem->cd(execDir);
40   
41   // ########################################################
42   // loop over runs
43   Int_t nRunCounter = 0;
44   for (Int_t r=1; r<=nDirs; r++) {
45
46     TSystemFile* presentDir = (TSystemFile*)dirList->At(r);
47     if (!presentDir->IsDirectory())
48       continue;
49     // first check that the files are there
50     sprintf(str,"%s/%s",dataDir, presentDir->GetName());
51     if ((!gSystem->Which(str,"galice.root")) ||
52         (!gSystem->Which(str,"AliESDs.root"))) 
53       continue;
54     
55     if (nRunCounter++ >= nRuns)
56       break;    
57     
58     cout << "run #" << nRunCounter << endl;
59
60     // #########################################################
61     // setup galice and runloader
62     if (gAlice) {
63       delete gAlice->GetRunLoader();
64       delete gAlice;
65       gAlice=0;
66     }
67
68     sprintf(str,"%s/run%d/galice.root",dataDir,r);
69     AliRunLoader* runLoader = AliRunLoader::Open(str);
70
71     runLoader->LoadgAlice();
72     gAlice = runLoader->GetAliRun();
73     runLoader->LoadHeader();
74
75     // #########################################################
76     // open esd file and get the tree
77
78     sprintf(str,"%s/run%d/AliESDs.root",dataDir,r);
79     // close it first to avoid memory leak
80     if (esdFile)
81       if (esdFile->IsOpen())
82         esdFile->Close();
83
84     esdFile = TFile::Open(str);
85     esdTree = (TTree*)esdFile->Get("esdTree");
86     if (!esdTree)
87       continue;
88     esdBranch = esdTree->GetBranch("ESD");
89     esdBranch->SetAddress(&esd);
90     if (!esdBranch)
91       continue;
92
93     // ########################################################
94     // Magnetic field
95     AliTracker::SetFieldMap(gAlice->Field(),kTRUE); // kTRUE means uniform magnetic field
96
97     // ########################################################
98     // getting number of events
99     Int_t nEvents    = (Int_t)runLoader->GetNumberOfEvents();
100     Int_t nESDEvents = esdBranch->GetEntries();
101     
102     if (nEvents!=nESDEvents) 
103       cout << " Warning: Different number of events from runloader and esdtree!!!" << nEvents << " / " << nESDEvents << endl;
104     
105     // ########################################################
106     // loop over number of events
107     cout << " looping over events..." << endl;
108     for(Int_t i=1; i<nEvents; i++) {
109       
110       esdBranch->GetEntry(i);
111       runLoader->GetEvent(i);
112       
113       // ########################################################
114       // get the EDS vertex
115       AliESDVertex* vtxESD = esd->GetVertex();
116       
117       Double_t vtxSigma[3];
118       vtxESD->GetSigmaXYZ(vtxSigma);
119       
120       // ########################################################
121       // loop over esd tracks      
122       Int_t nTracks = esd->GetNumberOfTracks();      
123
124       for (Int_t t=0; t<nTracks; t++) {
125         AliESDtrack* esdTrack = esd->GetTrack(t);      
126         
127         //trackCuts->AcceptTrack(esdTrack, vtxESD, esd->GetMagneticField());
128         trackCuts->AcceptTrack(esdTrack);
129         
130       } // end of track loop
131     } // end  of event loop
132   } // end of run loop
133
134   TFile* fout = new TFile("out.root","RECREATE");
135
136   trackCuts->SaveHistograms("esd_track_cuts");
137   
138   fout->Write();
139   fout->Close();
140
141 }