]>
Commit | Line | Data |
---|---|---|
d8770a0f | 1 | AliRunLoader *gAL=0; |
2 | Int_t gEvt=0; Int_t gMaxEvt=0; | |
3df18479 | 3 | TObjArray *pNmean,*pQthre; |
d8770a0f | 4 | TTree *gEsdTr; |
5 | //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ | |
e4290ede | 6 | void HESDfromKin(const char *name="default") |
d8770a0f | 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(); | |
6e06db1d | 18 | AliESDEvent *pEsd = new AliESDEvent(); |
d8770a0f | 19 | TFile *pEsdFl=TFile::Open("AliESDs.root","recreate"); |
20 | gEsdTr=new TTree("esdTree","Sim ESD from kinematics"); | |
e4290ede | 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 | ||
d8770a0f | 25 | } else return; |
26 | ||
3df18479 | 27 | if(!OpenCalib()) {Printf("Problems in OpenCalib!Bye.");return;} |
d8770a0f | 28 | |
c0d7adf8 | 29 | TString ttl=name; |
30 | Bool_t htaCheck=ttl.Contains("HTA"); | |
9785d5fb | 31 | // if(!htaCheck) SimEsd(pHL,pEsd); else SimEsdHidden(pHL,pEsd); |
32 | SimEsd(pHL,pEsd,htaCheck); | |
e4290ede | 33 | |
34 | pEsdFl->cd(); | |
d8770a0f | 35 | pEsdFl->Write();pEsdFl->Close(); |
36 | } | |
37 | //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ | |
9785d5fb | 38 | void SimEsd(AliLoader *pHL,AliESDEvent *pEsd,Bool_t htaCheck) |
d8770a0f | 39 | { |
9785d5fb | 40 | if(htaCheck) { |
41 | TFile *fout = new TFile("HTA.root","recreate"); | |
42 | TH1F *hdC = new TH1F("dC" ,";delta Cerenkov (rad)",100,-0.2,0.2); | |
43 | TH1F *hCer = new TH1F("Cer" ,"Theta Cerenkov (rad)",250,0.,0.75); | |
44 | TH2F *htvsp = new TH2F("tvsp",";momentum (GeV/c);theta Cerenkov (rad)",100,0.,5.,1000,0.,0.75); | |
45 | TH1F *hdth = new TH1F("dth" ,";Delta theta Trk (mrad)",100,-250,250); | |
46 | TH1F *hdph = new TH1F("dph" ,";Delta phi Trk (mrad)",100,-500,500); | |
47 | Double_t rd=TMath::RadToDeg(); | |
48 | Printf("----------------------------------------------"); | |
49 | Printf("| SimHTA:Utility to embed ESD from kinematics|"); | |
50 | Printf("| with Hidden Track Algorithm (HTA) |"); | |
51 | Printf("----------------------------------------------"); | |
52 | } else { | |
53 | Printf("-----------------------------------------------"); | |
54 | Printf("| SimESD: Utility to embed ESD from kinematics|"); | |
55 | Printf("-----------------------------------------------"); | |
56 | } | |
d8770a0f | 57 | AliHMPIDTracker::SetFieldMap(gAL->GetAliRun()->Field(),kTRUE); |
58 | AliHMPID *pH=(AliHMPID*)gAL->GetAliRun()->GetDetector("HMPID"); | |
d8770a0f | 59 | Int_t iNevt=gAL->GetNumberOfEvents(); |
1b46ef29 | 60 | pEsd->SetMagneticField(AliHMPIDTracker::GetBz()); |
d8770a0f | 61 | for(Int_t iEvt=0;iEvt<iNevt;iEvt++){//events loop |
d8770a0f | 62 | gAL->GetEvent(iEvt); |
63 | pHL->TreeR()->GetEntry(0); | |
64 | AliStack *pStack=gAL->Stack(); | |
9785d5fb | 65 | Int_t nTrkHMPID=0; |
d8770a0f | 66 | for(Int_t i=0;i<pStack->GetNtrack();i++){ |
9785d5fb | 67 | if(!pStack->IsPhysicalPrimary(i)) continue; |
d8770a0f | 68 | TParticle *pTrack=pStack->Particle(i); |
1b46ef29 | 69 | if(pTrack->GetPDG()->Charge()==0) continue; |
e4290ede | 70 | AliESDtrack trk(pTrack); |
39cd22e6 | 71 | Float_t xPc,yPc,xRa,yRa,thRa,phRa; |
72 | Int_t iCh=AliHMPIDTracker::IntTrkCha(&trk,xPc,yPc,xRa,yRa,thRa,phRa); //get chamber intersected by this track | |
73 | if(iCh<0) { | |
74 | trk.SetHMPIDtrk(0,0,0,0); //no intersection found | |
75 | trk.SetHMPIDcluIdx (99,99999); //chamber not found, mip not yet considered | |
76 | trk.SetHMPIDsignal(AliHMPIDRecon::kNotPerformed); //ring reconstruction not yet performed | |
77 | continue; //no intersection at all, go after next track | |
78 | } | |
9785d5fb | 79 | nTrkHMPID++; |
39cd22e6 | 80 | trk.SetHMPIDcluIdx (iCh,99999); //chamber not found, mip not yet considered |
81 | trk.SetHMPIDtrk(xRa,yRa,thRa,phRa); //store initial infos | |
d8770a0f | 82 | pEsd->AddTrack(&trk); |
d8770a0f | 83 | }// track loop |
ac66a50d | 84 | |
9785d5fb | 85 | if(!(iEvt%50)) Printf("Number of events processed: %i with tracks %i in HMPID",iEvt,nTrkHMPID); |
86 | if(!htaCheck) AliHMPIDTracker::Recon(pEsd,pH->CluLst(),pNmean,pQthre); | |
ac66a50d | 87 | // Perform PID |
88 | ||
9785d5fb | 89 | |
ac66a50d | 90 | for(Int_t iTrk=0;iTrk<pEsd->GetNumberOfTracks();iTrk++){ //ESD tracks loop |
9785d5fb | 91 | |
ac66a50d | 92 | AliESDtrack *pTrk = pEsd->GetTrack(iTrk); //get reconstructed track |
9785d5fb | 93 | |
94 | Float_t radX,radY,thetaTrk,phiTrk; | |
95 | pTrk->GetHMPIDtrk(radX,radY,thetaTrk,phiTrk); | |
96 | ||
97 | if(htaCheck) { | |
98 | Int_t iCh = pTrk->GetHMPIDcluIdx(); //chamber not found, mip not yet considered | |
99 | iCh/=1000000; | |
9a573d52 | 100 | Printf("simulated track theta %f phi %f",thetaTrk*rd,phiTrk*rd); |
9785d5fb | 101 | TObjArray *pClus = pH->CluLst(); |
102 | ||
103 | if(AliHMPIDTracker::ReconHiddenTrk(pEsd,pClus,pNmean,pQthre)!=0) continue; | |
104 | ||
105 | Double_t thetaCerSim = TMath::ACos(pTrack->Energy()/(1.292*pTrack->P())); | |
106 | Printf("-----------------------------------------------------------"); | |
107 | Printf(" theta Cerenkov simulated %f",thetaCerSim); | |
108 | Printf(" theta Cerenkov reconstructed %f",pTrk->GetHMPIDsignal()); | |
109 | Float_t thetaFit,phiFit; | |
110 | pTrk->GetHMPIDtrk(radX,radY,thetaFit,phiFit); | |
111 | // Printf("reconstr. track theta %f phi %f",thetaFit*rd,phiFit*rd); | |
112 | ||
113 | // fill useful histos | |
114 | hdC->Fill(pTrk->GetHMPIDsignal()-thetaCerSim); | |
115 | hCer->Fill(pTrk->GetHMPIDsignal()); | |
116 | htvsp->Fill(pTrk->P(),pTrk->GetHMPIDsignal()); | |
117 | hdth->Fill((thetaFit-thetaTrk)*1000); | |
118 | hdph->Fill((phiFit-phiTrk)*1000); | |
119 | } | |
120 | ||
ac66a50d | 121 | AliHMPIDPid pID; |
122 | Double_t prob[5]; | |
123 | pID.FindPid(pTrk,5,prob); | |
d35d51b7 | 124 | pTrk->SetHMPIDpid(prob); |
ac66a50d | 125 | // 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); |
126 | } | |
127 | ||
d8770a0f | 128 | gEsdTr->Fill(); |
129 | pEsd->Reset(); | |
ac66a50d | 130 | |
d8770a0f | 131 | }// event loop |
132 | Printf("Events processed %i",iEvt); | |
133 | gAL->UnloadHeader(); gAL->UnloadKinematics(); | |
3b49956b | 134 | }//Esd() |
d8770a0f | 135 | //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
6e06db1d | 136 | void SimEsdHidden(AliLoader *pHL,AliESDEvent *pEsd) |
3b49956b | 137 | { |
1b46ef29 | 138 | TFile *fout = new TFile("HTA.root","recreate"); |
139 | TH1F *hdC = new TH1F("dC" ,";delta Cerenkov (rad)",100,-0.2,0.2); | |
140 | TH1F *hCer = new TH1F("Cer" ,"Theta Cerenkov (rad)",250,0.,0.75); | |
141 | TH2F *htvsp = new TH2F("tvsp",";momentum (GeV/c);theta Cerenkov (rad)",100,0.,5.,1000,0.,0.75); | |
142 | TH1F *hdth = new TH1F("dth" ,";Delta theta Trk (mrad)",100,-250,250); | |
143 | TH1F *hdph = new TH1F("dph" ,";Delta phi Trk (mrad)",100,-500,500); | |
3b49956b | 144 | Double_t rd=TMath::RadToDeg(); |
c0d7adf8 | 145 | Printf("----------------------------------------------"); |
146 | Printf("| SimHTA:Utility to embed ESD from kinematics|"); | |
147 | Printf("| with Hidden Track Algorithm (HTA) |"); | |
148 | Printf("----------------------------------------------"); | |
d8770a0f | 149 | AliHMPIDTracker::SetFieldMap(gAL->GetAliRun()->Field(),kTRUE); |
150 | AliHMPID *pH=(AliHMPID*)gAL->GetAliRun()->GetDetector("HMPID"); | |
d8770a0f | 151 | Int_t iNevt=gAL->GetNumberOfEvents(); |
9785d5fb | 152 | pEsd->SetMagneticField(AliHMPIDTracker::GetBz()); |
1b46ef29 | 153 | |
d8770a0f | 154 | Printf("Number of events to process: %i",iNevt); |
1b46ef29 | 155 | |
d8770a0f | 156 | for(Int_t iEvt=0;iEvt<iNevt;iEvt++){//events loop |
9785d5fb | 157 | if(!(iEvt%1)) Printf("Events processed %i",iEvt); |
3b49956b | 158 | gAL->GetEvent(iEvt); |
d8770a0f | 159 | pHL->TreeR()->GetEntry(0); |
3b49956b | 160 | AliStack *pStack=gAL->Stack(); |
1b46ef29 | 161 | |
3b49956b | 162 | for(Int_t i=0;i<pStack->GetNtrack();i++){ |
1b46ef29 | 163 | |
9785d5fb | 164 | (!pStack->IsPhysicalPrimary(i)) continue; |
3b49956b | 165 | TParticle *pTrack=pStack->Particle(i); |
1b46ef29 | 166 | |
9785d5fb | 167 | if(pTrack->GetPDG()->Charge()==0) continue; |
3b49956b | 168 | //find the chamber that intersects HMPID |
9785d5fb | 169 | pTrack->Print(); |
3b49956b | 170 | AliESDtrack trk(pTrack); |
39cd22e6 | 171 | Float_t xPc,yPc,xRa,yRa,thRa,phRa; |
172 | Int_t iCh=AliHMPIDTracker::IntTrkCha(&trk,xPc,yPc,xRa,yRa,thRa,phRa); //get chamber intersected by this track | |
3b49956b | 173 | if(iCh<0) continue; //no intersection at all, go after next track |
39cd22e6 | 174 | trk.SetHMPIDtrk(xRa,yRa,thRa,phRa); //store initial infos |
3b49956b | 175 | Float_t radX,radY,thetaTrk,phiTrk; |
176 | trk.GetHMPIDtrk(radX,radY,thetaTrk,phiTrk); | |
3df18479 | 177 | // Printf("simulated track theta %f phi %f",thetaTrk*rd,phiTrk*rd); |
3b49956b | 178 | TObjArray *pClus = pH->CluLst(); |
1b46ef29 | 179 | |
39cd22e6 | 180 | if(AliHMPIDTracker::ReconHiddenTrk(iCh,0,&trk,(TClonesArray *)pClus->At(iCh),pNmean,pQthre)!=0) continue; |
1b46ef29 | 181 | |
3b49956b | 182 | pEsd->AddTrack(&trk); |
183 | Double_t thetaCerSim = TMath::ACos(pTrack->Energy()/(1.292*pTrack->P())); | |
3df18479 | 184 | // Printf(" theta Cerenkov simulated %f",thetaCerSim); |
185 | // Printf(" theta Cerenkov reconstructed %f",trk.GetHMPIDsignal()); | |
1b46ef29 | 186 | Float_t thetaFit,phiFit; |
187 | trk.GetHMPIDtrk(radX,radY,thetaFit,phiFit); | |
188 | // Printf("reconstr. track theta %f phi %f",thetaFit*rd,phiFit*rd); | |
189 | ||
190 | // fill useful histos | |
3df18479 | 191 | hdC->Fill(trk.GetHMPIDsignal()-thetaCerSim); |
192 | hCer->Fill(trk.GetHMPIDsignal()); | |
193 | htvsp->Fill(pTrack->P(),trk.GetHMPIDsignal()); | |
1b46ef29 | 194 | hdth->Fill((thetaFit-thetaTrk)*1000); |
195 | hdph->Fill((phiFit-phiTrk)*1000); | |
3b49956b | 196 | }// track loop |
d8770a0f | 197 | pEsd->SetMagneticField(AliHMPIDTracker::GetBz()); |
198 | gEsdTr->Fill(); | |
199 | pEsd->Reset(); | |
200 | }// event loop | |
201 | Printf("Events processed %i",iEvt); | |
1b46ef29 | 202 | fout->Write(); |
203 | fout->Close(); | |
204 | delete fout; | |
d8770a0f | 205 | gAL->UnloadHeader(); gAL->UnloadKinematics(); |
3b49956b | 206 | }//EsdHidden() |
d8770a0f | 207 | //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
3df18479 | 208 | Bool_t OpenCalib() |
d8770a0f | 209 | { |
210 | AliCDBManager* pCDB = AliCDBManager::Instance(); | |
211 | pCDB->SetDefaultStorage("local://$HOME"); | |
212 | AliCDBEntry *pQthreEnt=pCDB->Get("HMPID/Calib/Qthre",0); | |
213 | AliCDBEntry *pNmeanEnt=pCDB->Get("HMPID/Calib/Nmean",0); | |
214 | ||
3df18479 | 215 | if(!pQthreEnt || !pNmeanEnt) return kFALSE; |
d8770a0f | 216 | |
217 | pNmean=(TObjArray*)pNmeanEnt->GetObject(); | |
3df18479 | 218 | pQthre=(TObjArray*)pQthreEnt->GetObject(); |
219 | ||
220 | if(!pQthre || !pNmean) return kFALSE; | |
221 | return kTRUE; | |
d8770a0f | 222 | } |
223 | //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |