Calculate crosstalk correction in separate function. Do calculation in 2 itterations...
[u/mrichter/AliRoot.git] / test / merge / test.C
1 // Usage in compiled mode
2 // gSystem->SetIncludePath("-I$ROOTSYS/include -I$ALICE_ROOT/include");  
3 // gROOT->LoadMacro("test.C+");
4 // test()
5
6 #if !defined(__CINT__) || defined(__MAKECINT__)
7
8 // Root include files
9 #include <Riostream.h>
10 #include <TFile.h>
11 #include <TTree.h>
12 #include <TBranch.h>
13 #include <TStopwatch.h>
14 #include <TObject.h>
15 #include <TParticle.h>
16
17 // AliRoot include files
18 #include "AliESDEvent.h"
19 #include "AliRunLoader.h"
20 #include "AliRun.h"
21 #include "AliStack.h"
22
23 #endif
24
25 void test(const char * sdir ="signal",
26           const char * bdir ="backgr") {
27
28   TStopwatch timer;
29   timer.Start();
30
31   TString name;
32
33   // Signal file, tree, and branch
34   name = sdir;
35   name += "/AliESDs.root";
36   TFile * fSig = TFile::Open(name.Data());
37   TTree * tSig = (TTree*)fSig->Get("esdTree");
38
39   AliESDEvent * esdSig = new AliESDEvent();// The signal ESD object is put here
40   esdSig->ReadFromTree(tSig);
41
42   // Run loader (signal events)
43   name = sdir;
44   name += "/galice.root";
45   AliRunLoader* rlSig = AliRunLoader::Open(name.Data());
46
47   // Run loader (underlying events)
48   name = bdir;
49   name += "/galice.root";
50   AliRunLoader* rlUnd = AliRunLoader::Open(name.Data(),"Underlying");
51
52   // gAlice
53   rlSig->LoadgAlice();
54   rlUnd->LoadgAlice();
55   gAlice = rlSig->GetAliRun();
56
57   // Now load kinematics and event header
58   rlSig->LoadKinematics();
59   rlSig->LoadHeader();
60   rlUnd->LoadKinematics();
61   rlUnd->LoadHeader();
62
63   // Loop on events: check that MC and data contain the same number of events
64   Long64_t nevSig = rlSig->GetNumberOfEvents();
65   Long64_t nevUnd = rlUnd->GetNumberOfEvents();
66   Long64_t nSigPerUnd = nevSig/nevUnd;
67
68   cout << nevSig << " signal events" << endl;
69   cout << nevUnd << " underlying events" << endl;
70   cout << nSigPerUnd << " signal events per one underlying" << endl;
71
72   for (Int_t iev=0; iev<nevSig; iev++) {
73     cout << "Signal event " << iev << endl;
74     Int_t ievUnd = iev/nSigPerUnd;
75     cout << "Underlying event " << ievUnd << endl;
76
77     // Get signal ESD
78     tSig->GetEntry(iev);
79     // Get signal kinematics
80     rlSig->GetEvent(iev);
81     // Get underlying kinematics
82     rlUnd->GetEvent(ievUnd);
83
84     // Particle stack
85     AliStack * stackSig = rlSig->Stack();
86     Int_t nPartSig = stackSig->GetNtrack();
87     AliStack * stackUnd = rlUnd->Stack();
88     Int_t nPartUnd = stackUnd->GetNtrack();
89
90     Int_t nrec = esdSig->GetNumberOfTracks();
91     cout << nrec << " reconstructed tracks" << endl;
92     for(Int_t irec=0; irec<nrec; irec++) {
93       AliESDtrack * track = esdSig->GetTrack(irec);
94       UInt_t label = TMath::Abs(track->GetTPCLabel());
95       if (label>=10000000) {
96         // Underlying event. 10000000 is the
97         // value of fkMASKSTEP in AliRunDigitizer
98 //      cout << " Track from the underlying event" << endl;
99         label %=10000000;
100         if (label>=nPartUnd) continue;
101         TParticle * part = stackUnd->Particle(label);
102         if(part) part->Print();
103       }
104       else {
105         cout << " Track " << label << " from the signal event" << endl;
106         if (label>=nPartSig) {
107           cout <<"Strange, label outside the range "<< endl;
108           continue;
109         }
110         TParticle * part = stackSig->Particle(label);
111         if(part) part->Print();
112       }
113
114     }
115
116   }
117
118   fSig->Close();
119
120   timer.Stop();
121   timer.Print();
122 }