]>
Commit | Line | Data |
---|---|---|
1 | AliRunLoader *gAL=0; | |
2 | Int_t gEvt=0; Int_t gMaxEvt=0; | |
3 | TObjArray *pNmean,*pQthre; | |
4 | TTree *gEsdTr; | |
5 | //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ | |
6 | void HESDfromKin(const char *name="default") | |
7 | {//simulate ESD from kinematics | |
8 | ||
9 | if(gSystem->IsFileInIncludePath("galice.root")){// tries to open session | |
10 | if(gAlice) delete gAlice; //in case we execute this in aliroot delete default AliRun object | |
11 | gAL=AliRunLoader::Open(); //try to open galice.root from current dir | |
12 | gAL->LoadgAlice(); //take new AliRun object from galice.root | |
13 | if(gAL->LoadHeader()) return; | |
14 | if(gAL->LoadKinematics()) return; | |
15 | ||
16 | AliLoader *pHL=gAL->GetDetectorLoader("HMPID"); | |
17 | pHL->LoadRecPoints(); | |
18 | AliESDEvent *pEsd = new AliESDEvent(); | |
19 | TFile *pEsdFl=TFile::Open("AliESDs.root","recreate"); | |
20 | gEsdTr=new TTree("esdTree","Sim ESD from kinematics"); | |
21 | pEsd->CreateStdContent(); pEsd->WriteToTree(gEsdTr); //clm: new ESD write schema: see Task Force meeting 20th June, 2007 | |
22 | gEsdTr->GetUserInfo()->Add(pEsd); //clm: TList has to be created for ReadFromTree method -- this was not needed by the old ESD | |
23 | ||
24 | ||
25 | } else return; | |
26 | ||
27 | if(!OpenCalib()) {Printf("Problems in OpenCalib!Bye.");return;} | |
28 | ||
29 | TString ttl=name; | |
30 | Bool_t htaCheck=ttl.Contains("HTA"); | |
31 | if(!htaCheck) SimEsd(pHL,pEsd); else SimEsdHidden(pHL,pEsd); | |
32 | ||
33 | pEsdFl->cd(); | |
34 | pEsdFl->Write();pEsdFl->Close(); | |
35 | } | |
36 | //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ | |
37 | void SimEsd(AliLoader *pHL,AliESDEvent *pEsd) | |
38 | { | |
39 | Printf("-----------------------------------------------"); | |
40 | Printf("| SimESD: Utility to embed ESD from kinematics|"); | |
41 | Printf("-----------------------------------------------"); | |
42 | AliHMPIDTracker::SetFieldMap(gAL->GetAliRun()->Field(),kTRUE); | |
43 | AliHMPID *pH=(AliHMPID*)gAL->GetAliRun()->GetDetector("HMPID"); | |
44 | Int_t mtid=-1; | |
45 | Int_t iNevt=gAL->GetNumberOfEvents(); | |
46 | pEsd->SetMagneticField(AliHMPIDTracker::GetBz()); | |
47 | Printf("Number of events to process: %i",iNevt); | |
48 | for(Int_t iEvt=0;iEvt<iNevt;iEvt++){//events loop | |
49 | if(!(iEvt%50)) Printf("Events processed %i",iEvt); | |
50 | gAL->GetEvent(iEvt); | |
51 | pHL->TreeR()->GetEntry(0); | |
52 | AliStack *pStack=gAL->Stack(); | |
53 | for(Int_t i=0;i<pStack->GetNtrack();i++){ | |
54 | TParticle *pTrack=pStack->Particle(i); | |
55 | mtid=pTrack->GetFirstMother(); | |
56 | if(mtid>=0) continue; // only primaries | |
57 | if(pTrack->GetPDG()->Charge()==0) continue; | |
58 | AliESDtrack trk(pTrack); | |
59 | Float_t xPc,yPc,xRa,yRa,thRa,phRa; | |
60 | Int_t iCh=AliHMPIDTracker::IntTrkCha(&trk,xPc,yPc,xRa,yRa,thRa,phRa); //get chamber intersected by this track | |
61 | if(iCh<0) { | |
62 | trk.SetHMPIDtrk(0,0,0,0); //no intersection found | |
63 | trk.SetHMPIDcluIdx (99,99999); //chamber not found, mip not yet considered | |
64 | trk.SetHMPIDsignal(AliHMPIDRecon::kNotPerformed); //ring reconstruction not yet performed | |
65 | continue; //no intersection at all, go after next track | |
66 | } | |
67 | // Printf(" particle analyzed %d with mtid %d is %s hitting chamber %d",i,mtid,pTrack->GetName(),iCh); | |
68 | trk.SetHMPIDcluIdx (iCh,99999); //chamber not found, mip not yet considered | |
69 | trk.SetHMPIDtrk(xRa,yRa,thRa,phRa); //store initial infos | |
70 | pEsd->AddTrack(&trk); | |
71 | }// track loop | |
72 | ||
73 | AliHMPIDTracker::Recon(pEsd,pH->CluLst(),pNmean,pQthre); | |
74 | // Perform PID | |
75 | ||
76 | for(Int_t iTrk=0;iTrk<pEsd->GetNumberOfTracks();iTrk++){ //ESD tracks loop | |
77 | AliESDtrack *pTrk = pEsd->GetTrack(iTrk); //get reconstructed track | |
78 | AliHMPIDPid pID; | |
79 | Double_t prob[5]; | |
80 | pID.FindPid(pTrk,5,prob); | |
81 | pTrk->SetHMPIDpid(prob); | |
82 | // Printf(" Prob e- %6.2f mu %6.2f pi %6.2f k %6.2f p %6.2f",prob[0]*100,prob[1]*100,prob[2]*100,prob[3]*100,prob[4]*100); | |
83 | } | |
84 | ||
85 | gEsdTr->Fill(); | |
86 | pEsd->Reset(); | |
87 | ||
88 | }// event loop | |
89 | Printf("Events processed %i",iEvt); | |
90 | gAL->UnloadHeader(); gAL->UnloadKinematics(); | |
91 | }//Esd() | |
92 | //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ | |
93 | void SimEsdHidden(AliLoader *pHL,AliESDEvent *pEsd) | |
94 | { | |
95 | TFile *fout = new TFile("HTA.root","recreate"); | |
96 | TH1F *hdC = new TH1F("dC" ,";delta Cerenkov (rad)",100,-0.2,0.2); | |
97 | TH1F *hCer = new TH1F("Cer" ,"Theta Cerenkov (rad)",250,0.,0.75); | |
98 | TH2F *htvsp = new TH2F("tvsp",";momentum (GeV/c);theta Cerenkov (rad)",100,0.,5.,1000,0.,0.75); | |
99 | TH1F *hdth = new TH1F("dth" ,";Delta theta Trk (mrad)",100,-250,250); | |
100 | TH1F *hdph = new TH1F("dph" ,";Delta phi Trk (mrad)",100,-500,500); | |
101 | Double_t rd=TMath::RadToDeg(); | |
102 | Printf("----------------------------------------------"); | |
103 | Printf("| SimHTA:Utility to embed ESD from kinematics|"); | |
104 | Printf("| with Hidden Track Algorithm (HTA) |"); | |
105 | Printf("----------------------------------------------"); | |
106 | AliHMPIDTracker::SetFieldMap(gAL->GetAliRun()->Field(),kTRUE); | |
107 | AliHMPID *pH=(AliHMPID*)gAL->GetAliRun()->GetDetector("HMPID"); | |
108 | Int_t mtid=-1; | |
109 | Int_t iNevt=gAL->GetNumberOfEvents(); | |
110 | ||
111 | Printf("Number of events to process: %i",iNevt); | |
112 | ||
113 | for(Int_t iEvt=0;iEvt<iNevt;iEvt++){//events loop | |
114 | if(!(iEvt%50)) Printf("Events processed %i",iEvt); | |
115 | gAL->GetEvent(iEvt); | |
116 | pHL->TreeR()->GetEntry(0); | |
117 | AliStack *pStack=gAL->Stack(); | |
118 | ||
119 | for(Int_t i=0;i<pStack->GetNtrack();i++){ | |
120 | ||
121 | TParticle *pTrack=pStack->Particle(i); | |
122 | mtid=pTrack->GetFirstMother(); | |
123 | ||
124 | if(mtid>=0) continue; // only primaries | |
125 | ||
126 | //find the chamber that intersects HMPID | |
127 | AliESDtrack trk(pTrack); | |
128 | Float_t xPc,yPc,xRa,yRa,thRa,phRa; | |
129 | Int_t iCh=AliHMPIDTracker::IntTrkCha(&trk,xPc,yPc,xRa,yRa,thRa,phRa); //get chamber intersected by this track | |
130 | if(iCh<0) continue; //no intersection at all, go after next track | |
131 | trk.SetHMPIDtrk(xRa,yRa,thRa,phRa); //store initial infos | |
132 | Float_t radX,radY,thetaTrk,phiTrk; | |
133 | trk.GetHMPIDtrk(radX,radY,thetaTrk,phiTrk); | |
134 | // Printf("simulated track theta %f phi %f",thetaTrk*rd,phiTrk*rd); | |
135 | TObjArray *pClus = pH->CluLst(); | |
136 | ||
137 | if(AliHMPIDTracker::ReconHiddenTrk(iCh,0,&trk,(TClonesArray *)pClus->At(iCh),pNmean,pQthre)!=0) continue; | |
138 | ||
139 | pEsd->AddTrack(&trk); | |
140 | Double_t thetaCerSim = TMath::ACos(pTrack->Energy()/(1.292*pTrack->P())); | |
141 | // Printf(" theta Cerenkov simulated %f",thetaCerSim); | |
142 | // Printf(" theta Cerenkov reconstructed %f",trk.GetHMPIDsignal()); | |
143 | Float_t thetaFit,phiFit; | |
144 | trk.GetHMPIDtrk(radX,radY,thetaFit,phiFit); | |
145 | // Printf("reconstr. track theta %f phi %f",thetaFit*rd,phiFit*rd); | |
146 | ||
147 | // fill useful histos | |
148 | hdC->Fill(trk.GetHMPIDsignal()-thetaCerSim); | |
149 | hCer->Fill(trk.GetHMPIDsignal()); | |
150 | htvsp->Fill(pTrack->P(),trk.GetHMPIDsignal()); | |
151 | hdth->Fill((thetaFit-thetaTrk)*1000); | |
152 | hdph->Fill((phiFit-phiTrk)*1000); | |
153 | }// track loop | |
154 | pEsd->SetMagneticField(AliHMPIDTracker::GetBz()); | |
155 | gEsdTr->Fill(); | |
156 | pEsd->Reset(); | |
157 | }// event loop | |
158 | Printf("Events processed %i",iEvt); | |
159 | fout->Write(); | |
160 | fout->Close(); | |
161 | delete fout; | |
162 | gAL->UnloadHeader(); gAL->UnloadKinematics(); | |
163 | }//EsdHidden() | |
164 | //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ | |
165 | Bool_t OpenCalib() | |
166 | { | |
167 | AliCDBManager* pCDB = AliCDBManager::Instance(); | |
168 | pCDB->SetDefaultStorage("local://$HOME"); | |
169 | AliCDBEntry *pQthreEnt=pCDB->Get("HMPID/Calib/Qthre",0); | |
170 | AliCDBEntry *pNmeanEnt=pCDB->Get("HMPID/Calib/Nmean",0); | |
171 | ||
172 | if(!pQthreEnt || !pNmeanEnt) return kFALSE; | |
173 | ||
174 | pNmean=(TObjArray*)pNmeanEnt->GetObject(); | |
175 | pQthre=(TObjArray*)pQthreEnt->GetObject(); | |
176 | ||
177 | if(!pQthre || !pNmean) return kFALSE; | |
178 | return kTRUE; | |
179 | } | |
180 | //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |